program main

!*****************************************************************************80
!
!! qr_solve_test() tests qr_solve().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 October 2012
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'qr_solve_test():'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Test qr_solve().'
  write ( *, '(a)' ) '  QR_SOLVE needs the R8LIB library.'
  write ( *, '(a)' ) '  This test also needs the TEST_LS library.'

  call normal_solve_test ( )
  call qr_solve_test ( )
  call svd_solve_test ( )
  call dqrls_test ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'qr_solve_test():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop 0
end
subroutine normal_solve_test ( )

!*****************************************************************************80
!
!! NORMAL_SOLVE_TEST tests NORMAL_SOLVE.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 April 2012
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), allocatable :: a(:,:)
  real ( kind = rk8 ), allocatable :: b(:)
  real ( kind = rk8 ) b_norm
  integer flag
  integer m
  integer n
  integer prob
  integer prob_num
  real ( kind = rk8 ), allocatable :: r1(:)
  real ( kind = rk8 ) r1_norm
  real ( kind = rk8 ), allocatable :: r2(:)
  real ( kind = rk8 ) r2_norm
  real ( kind = rk8 ) r8vec_norm
  real ( kind = rk8 ) r8vec_norm_affine
  real ( kind = rk8 ) x_diff_norm
  real ( kind = rk8 ), allocatable :: x1(:)
  real ( kind = rk8 ) x1_norm
  real ( kind = rk8 ), allocatable :: x2(:)
  real ( kind = rk8 ) x2_norm

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'normal_solve_test():'
  write ( *, '(a)' ) '  normal_solve() is a function with a simple interface which'
  write ( *, '(a)' ) '  solves a linear system A*x = b in the least squares sense.'
  write ( *, '(a)' ) '  Compare a tabulated solution X1 to the NORMAL_SOLVE result X2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  NORMAL_SOLVE cannot be applied when N < M,'
  write ( *, '(a)' ) '  or if the matrix does not have full column rank.'

  call p00_prob_num ( prob_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  Number of problems = ', prob_num
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Index     M     N     ||B||         ||X1 - X2||   ||X1||       ||X2||        ||R1||        ||R2||'
  write ( *, '(a)' ) ' '

  do prob = 1, prob_num
!
!  Get problem size.
!
    call p00_m ( prob, m )
    call p00_n ( prob, n )
!
!  Allocate space.
!
    allocate ( a(1:m,1:n) )
    allocate ( b(1:m) )
    allocate ( r1(1:m) )
    allocate ( r2(1:m) )
    allocate ( x1(1:n) )
    allocate ( x2(1:n) )
!
!  Retrieve problem data.
!
    call p00_a ( prob, m, n, a )
    call p00_b ( prob, m, b )
    call p00_x ( prob, n, x1 )

    b_norm = r8vec_norm ( m, b )
    x1_norm = r8vec_norm ( n, x1 )
    r1(1:m) = matmul ( a(1:m,1:n), x1(1:n) ) - b(1:m)
    r1_norm = r8vec_norm ( m, r1 )
!
!  Use NORMAL_SOLVE on the problem.
!
    call normal_solve ( m, n, a, b, x2, flag )

    if ( flag /= 0 ) then

      write ( *, &
        '(2x,i5,2x,i4,2x,i4,2x,g12.4,2x,a,2x,g12.4,2x,a,2x,g12.4,2x,a)' ) &
        prob, m, n, b_norm, '------------', x1_norm, '------------', &
        r1_norm, '------------'

    else

      x2_norm = r8vec_norm ( n, x2 )
      r2(1:m) = matmul ( a(1:m,1:n), x2(1:n) ) - b(1:m)
      r2_norm = r8vec_norm ( m, r2 )
!
!  Compare tabulated and computed solutions.
!
      x_diff_norm = r8vec_norm_affine ( n, x1, x2 )
!
!  Report results for this problem.
!
      write ( *, '(2x,i5,2x,i4,2x,i4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4)' ) &
        prob, m, n, b_norm, x_diff_norm, x1_norm, x2_norm, r1_norm, r2_norm

    end if
!
!  Deallocate memory.
!
    deallocate ( a )
    deallocate ( b )
    deallocate ( r1 )
    deallocate ( r2 )
    deallocate ( x1 )
    deallocate ( x2 )

  end do

  return
end
subroutine qr_solve_test ( )

!*****************************************************************************80
!
!! QR_SOLVE_TEST tests QR_SOLVE.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 April 2012
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), allocatable :: a(:,:)
  real ( kind = rk8 ), allocatable :: b(:)
  real ( kind = rk8 ) b_norm
  integer m
  integer n
  integer prob
  integer prob_num
  real ( kind = rk8 ), allocatable :: r1(:)
  real ( kind = rk8 ) r1_norm
  real ( kind = rk8 ), allocatable :: r2(:)
  real ( kind = rk8 ) r2_norm
  real ( kind = rk8 ) r8vec_norm
  real ( kind = rk8 ) r8vec_norm_affine
  real ( kind = rk8 ) x_diff_norm
  real ( kind = rk8 ), allocatable :: x1(:)
  real ( kind = rk8 ) x1_norm
  real ( kind = rk8 ), allocatable :: x2(:)
  real ( kind = rk8 ) x2_norm

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QR_SOLVE_TEST'
  write ( *, '(a)' ) '  QR_SOLVE is a function with a simple interface which'
  write ( *, '(a)' ) '  solves a linear system A*x = b in the least squares sense.'
  write ( *, '(a)' ) '  Compare a tabulated solution X1 to the QR_SOLVE result X2.'

  call p00_prob_num ( prob_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  Number of problems = ', prob_num
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Index     M     N     ||B||         ||X1 - X2||' // &
    '   ||X1||       ||X2||        ||R1||        ||R2||'
  write ( *, '(a)' ) ' '

  do prob = 1, prob_num
!
!  Get problem size.
!
    call p00_m ( prob, m )
    call p00_n ( prob, n )
!
!  Allocate space.
!
    allocate ( a(1:m,1:n) )
    allocate ( b(1:m) )
    allocate ( r1(1:m) )
    allocate ( r2(1:m) )
    allocate ( x1(1:n) )
    allocate ( x2(1:n) )
!
!  Retrieve problem data.
!
    call p00_a ( prob, m, n, a )
    call p00_b ( prob, m, b )
    call p00_x ( prob, n, x1 )

    b_norm = r8vec_norm ( m, b )
    x1_norm = r8vec_norm ( n, x1 )
    r1(1:m) = matmul ( a(1:m,1:n), x1(1:n) ) - b(1:m)
    r1_norm = r8vec_norm ( m, r1 )
!
!  Use QR_SOLVE on the problem.
!
    call qr_solve ( m, n, a, b, x2 )

    x2_norm = r8vec_norm ( n, x2 )
    r2(1:m) = matmul ( a(1:m,1:n), x2(1:n) ) - b(1:m)
    r2_norm = r8vec_norm ( m, r2 )
!
!  Compare tabulated and computed solutions.
!
    x_diff_norm = r8vec_norm_affine ( n, x1, x2 )
!
!  Report results for this problem.
!
    write ( *, '(2x,i5,2x,i4,2x,i4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4)' ) &
      prob, m, n, b_norm, x_diff_norm, x1_norm, x2_norm, r1_norm, r2_norm
!
!  Deallocate memory.
!
    deallocate ( a )
    deallocate ( b )
    deallocate ( r1 )
    deallocate ( r2 )
    deallocate ( x1 )
    deallocate ( x2 )

  end do

  return
end
subroutine r8mat_cholesky_factor ( n, a, c, flag )

!*****************************************************************************80
!
!! r8mat_cholesky_factor() computes the Cholesky factor of a symmetric matrix.
!
!  Discussion:
!
!    The matrix must be symmetric and positive semidefinite.
!
!    For a positive semidefinite symmetric matrix A, the Cholesky factorization
!    is a lower triangular matrix L such that:
!
!      A = L * L'
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    08 April 2009
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) A(N,N), the N by N matrix.
!
!  Output:
!
!    real ( kind = rk8 ) C(N,N), the lower triangular Cholesky factor.
!
!    integer FLAG:
!    0, no error occurred.
!    1, the matrix is not positive definite.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n,n)
  real ( kind = rk8 ) c(n,n)
  integer flag
  integer i
  integer j
  real ( kind = rk8 ) sum2

  flag = 0

  c(1:n,1:n) = a(1:n,1:n)

  do j = 1, n

    c(1:j-1,j) = 0.0D+00

    do i = j, n

      sum2 = c(j,i) - dot_product ( c(j,1:j-1), c(i,1:j-1) )

      if ( i == j ) then
        if ( sum2 <= 0.0D+00 ) then
          flag = 1
          return
        else
          c(i,j) = sqrt ( sum2 )
        end if
      else
        if ( c(j,j) /= 0.0D+00 ) then
          c(i,j) = sum2 / c(j,j)
        else
          c(i,j) = 0.0D+00
        end if
      end if

    end do

  end do

  return
end
subroutine r8mat_cholesky_solve ( n, l, b, x )

!*****************************************************************************80
!
!! r8mat_cholesky_solve() solves a Cholesky factored linear system A * x = b.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!    This routine works with the lower triangular Cholesky factor A = L * L'.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) L(N,N), the N by N lower Cholesky factor of the
!    system matrix A.
!
!    real ( kind = rk8 ) B(N), the right hand side of the linear system.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), the solution of the linear system.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) b(n)
  real ( kind = rk8 ) l(n,n)
  real ( kind = rk8 ) x(n)
!
!  Solve L * y = b.
!
  call r8mat_l_solve ( n, l, b, x )
!
!  Solve L' * x = y.
!
  call r8mat_lt_solve ( n, l, x, x )

  return
end
subroutine r8mat_l_solve ( n, a, b, x )

!*****************************************************************************80
!
!! r8mat_l_solve() solves a lower triangular linear system.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) A(N,N), the N by N lower triangular matrix.
!
!    real ( kind = rk8 ) B(N), the right hand side of the linear system.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), the solution of the linear system.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n,n)
  real ( kind = rk8 ) b(n)
  integer i
  real ( kind = rk8 ) x(n)
!
!  Solve L * x = b.
!
  do i = 1, n
    x(i) = ( b(i) - dot_product ( a(i,1:i-1), x(1:i-1) ) ) / a(i,i)
  end do

  return
end
subroutine r8mat_lt_solve ( n, a, b, x )

!*****************************************************************************80
!
!! r8mat_lt_solve() solves a transposed lower triangular linear system.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!    Given the lower triangular matrix A, the linear system to be solved is:
!
!      A' * x = b
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of rows and columns.
!
!    real ( kind = rk8 ) A(N,N), the N by N lower triangular matrix.
!
!    real ( kind = rk8 ) B(N), the right hand side of the linear system.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), the solution of the linear system.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n,n)
  real ( kind = rk8 ) b(n)
  integer i
  real ( kind = rk8 ) x(n)
!
!  Solve L'*x = b.
!
  do i = n, 1, -1
    x(i) = ( b(i) - dot_product ( x(i+1:n), a(i+1:n,i) ) ) / a(i,i)
  end do

  return
end
subroutine r8mat_print ( m, n, a, title )

!*****************************************************************************80
!
!! r8mat_print() prints a real matrix.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    12 September 2004
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer M, the number of rows in A.
!
!    integer N, the number of columns in A.
!
!    real ( kind = rk8 ) A(M,N), the matrix.
!
!    character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer m
  integer n

  real ( kind = rk8 ) a(m,n)
  character ( len = * ) title

  call r8mat_print_some ( m, n, a, 1, 1, m, n, title )

  return
end
subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

!*****************************************************************************80
!
!! r8mat_print_some prints some of a real matrix.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 September 2009
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer M, N, the number of rows and columns.
!
!    real ( kind = rk8 ) A(M,N), an M by N matrix to be printed.
!
!    integer ILO, JLO, the first row and column to print.
!
!    integer IHI, JHI, the last row and column to print.
!
!    character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: incx = 5
  integer m
  integer n

  real ( kind = rk8 ) a(m,n)
  character ( len = 14 ) ctemp(incx)
  integer i
  integer i2hi
  integer i2lo
  integer ihi
  integer ilo
  integer inc
  integer j
  integer j2
  integer j2hi
  integer j2lo
  integer jhi
  integer jlo
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )

  if ( m <= 0 .or. n <= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  (None)'
    return
  end if

  do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx

    j2hi = j2lo + incx - 1
    j2hi = min ( j2hi, n )
    j2hi = min ( j2hi, jhi )

    inc = j2hi + 1 - j2lo

    write ( *, '(a)' ) ' '

    do j = j2lo, j2hi
      j2 = j + 1 - j2lo
      write ( ctemp(j2), '(i8,6x)' ) j
    end do

    write ( *, '(''  Col   '',5a14)' ) ctemp(1:inc)
    write ( *, '(a)' ) '  Row'
    write ( *, '(a)' ) ' '

    i2lo = max ( ilo, 1 )
    i2hi = min ( ihi, m )

    do i = i2lo, i2hi

      do j2 = 1, inc

        j = j2lo - 1 + j2

        write ( ctemp(j2), '(g14.6)' ) a(i,j)

      end do

      write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc )

    end do

  end do

  return
end
function r8vec_norm_affine ( n, v0, v1 )

!*****************************************************************************80
!
!! r8vec_norm_affine() returns the affine norm of an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    The affine vector L2 norm is defined as:
!
!      R8VEC_NORM_AFFINE(V0,V1)
!        = sqrt ( sum ( 1 <= I <= N ) ( V1(I) - V0(I) )^2
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 October 2010
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the order of the vectors.
!
!    real ( kind = rk8 ) V0(N), the base vector.
!
!    real ( kind = rk8 ) V1(N), the vector whose affine norm is desired.
!
!  Output:
!
!    real ( kind = rk8 ) R8VEC_NORM_AFFINE, the L2 norm of V1-V0.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) r8vec_norm_affine
  real ( kind = rk8 ) v0(n)
  real ( kind = rk8 ) v1(n)

  r8vec_norm_affine = sqrt ( sum ( ( v0(1:n) - v1(1:n) ) ** 2 ) )

  return
end
subroutine r8vec_print ( n, a, title )

!*****************************************************************************80
!
!! r8vec_print() prints an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 September 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of components of the vector.
!
!    real ( kind = rk8 ) A(N), the vector to be printed.
!
!    character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n)
  integer i
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i)
  end do

  return
end
subroutine svd_solve_test ( )

!*****************************************************************************80
!
!! SVD_SOLVE_TEST tests SVD_SOLVE.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 April 2012
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), allocatable :: a(:,:)
  real ( kind = rk8 ), allocatable :: b(:)
  real ( kind = rk8 ) b_norm
  integer m
  integer n
  integer prob
  integer prob_num
  real ( kind = rk8 ), allocatable :: r1(:)
  real ( kind = rk8 ) r1_norm
  real ( kind = rk8 ), allocatable :: r2(:)
  real ( kind = rk8 ) r2_norm
  real ( kind = rk8 ) r8vec_norm
  real ( kind = rk8 ) r8vec_norm_affine
  real ( kind = rk8 ) x_diff_norm
  real ( kind = rk8 ), allocatable :: x1(:)
  real ( kind = rk8 ) x1_norm
  real ( kind = rk8 ), allocatable :: x2(:)
  real ( kind = rk8 ) x2_norm

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SVD_SOLVE_TEST'
  write ( *, '(a)' ) '  SVD_SOLVE is a function with a simple interface which'
  write ( *, '(a)' ) '  solves a linear system A*x = b in the least squares sense.'
  write ( *, '(a)' ) '  Compare a tabulated solution X1 to the SVD_SOLVE result X2.'

  call p00_prob_num ( prob_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  Number of problems = ', prob_num
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Index     M     N     ||B||         ||X1 - X2||   ||X1||       ||X2||        ||R1||        ||R2||'
  write ( *, '(a)' ) ' '

  do prob = 1, prob_num
!
!  Get problem size.
!
    call p00_m ( prob, m )
    call p00_n ( prob, n )
!
!  Allocate space.
!
    allocate ( a(1:m,1:n) )
    allocate ( b(1:m) )
    allocate ( r1(1:m) )
    allocate ( r2(1:m) )
    allocate ( x1(1:n) )
    allocate ( x2(1:n) )
!
!  Retrieve problem data.
!
    call p00_a ( prob, m, n, a )
    call p00_b ( prob, m, b )
    call p00_x ( prob, n, x1 )

    b_norm = r8vec_norm ( m, b )
    x1_norm = r8vec_norm ( n, x1 )
    r1(1:m) = matmul ( a(1:m,1:n), x1(1:n) ) - b(1:m)
    r1_norm = r8vec_norm ( m, r1 )
!
!  Use SVD_SOLVE on the problem.
!
    call svd_solve ( m, n, a, b, x2 )

    x2_norm = r8vec_norm ( n, x2 )
    r2(1:m) = matmul ( a(1:m,1:n), x2(1:n) ) - b(1:m)
    r2_norm = r8vec_norm ( m, r2 )
!
!  Compare tabulated and computed solutions.
!
    x_diff_norm = r8vec_norm_affine ( n, x1, x2 )
!
!  Report results for this problem.
!
    write ( *, '(2x,i5,2x,i4,2x,i4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4,2x,g12.4)' ) &
      prob, m, n, b_norm, x_diff_norm, x1_norm, x2_norm, r1_norm, r2_norm
!
!  Deallocate memory.
!
    deallocate ( a )
    deallocate ( b )
    deallocate ( r1 )
    deallocate ( r2 )
    deallocate ( x1 )
    deallocate ( x2 )

  end do

  return
end
subroutine dqrls_test ( )

!*****************************************************************************80
!
!! DQRLS_TEST tests DQRLS.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 February 2012
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: m = 5
  integer, parameter :: n = 3

  real ( kind = rk8 ) a(m,n)
  real ( kind = rk8 ) b(m)
  integer i
  integer ind
  integer itask
  integer j
  integer jpvt(n)
  integer kr
  real ( kind = rk8 ) qraux(n)
  real ( kind = rk8 ) tol
  real ( kind = rk8 ) work(n)
  real ( kind = rk8 ) x(n)
!
!  Set up least-squares problem
!  quadratic model, equally-spaced points
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DQRLS_TEST'
  write ( *, '(a)' ) '  DQRLS solves a linear system A*x = b '
  write ( *, '(a)' ) '  in the least squares sense.'

  do i = 1, m
    a(i,1) = 1.0D+00
    do j = 2, n
      a(i,j) = a(i,j-1) * real ( i, kind = rk8 )
    end do
  end do

  b(1:5) = (/ 1.0D+00, 2.3D+00, 4.6D+00, 3.1D+00, 1.2D+00 /)

  tol = 1.0D-06

  call r8mat_print ( m, n, a, '  Coefficient matrix A:' )

  call r8vec_print ( m, b, '  Right hand side b:' )
!
!  Solve least-squares problem.
!
  itask = 1
  call dqrls ( a, m, m, n, tol, kr, b, x, b, jpvt, qraux, work, itask, ind )
!
!  Print results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  Error code =', ind
  write ( *, '(a,i4)' ) '  Estimated matrix rank =', kr

  call r8vec_print ( n, x, '  Least squares solution x:' )

  call r8vec_print ( m, b, '  Residuals A*x-b' )

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end