program linpacks_prb
!
!*******************************************************************************
!
!! LINPACKS_PRB calls each of the LINPACKS test routines.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'LINPACKS_PRB'
  write ( *, '(a)' ) '  Tests for LINPACKS.'
  write ( *, '(a)' ) ' '

  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09
  call test10
  call test11

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'LINPACKS_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'
 
  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests SGEFA.
!! TEST01 tests SGESL.
!
!  This test shows the simplest use of LINPACKS:
!
!    Solve A*x = b
!
!  where A is a given matrix, and B a right hand side.
!
!  We will also assume that A is stored in the simplest
!  possible way.
!
  implicit none
!
  integer, parameter :: n = 3
  integer, parameter :: lda = n

  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer ipvt(n)
  integer j
  integer job
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  SGEFA/SGESL'
  write ( *, '(a)' ) '  General storage, factor/solve'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Set the values of the matrix A.
!
  a(1,1) = 1.0E+00
  a(1,2) = 2.0E+00
  a(1,3) = 3.0E+00
 
  a(2,1) = 4.0E+00
  a(2,2) = 5.0E+00
  a(2,3) = 6.0E+00
 
  a(3,1) = 7.0E+00
  a(3,2) = 8.0E+00
  a(3,3) = 0.0E+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The matrix A is '
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:n)
  end do
!
!  Set the values of the right hand side vector B.
!
  b(1:n) = (/ 6.0E+00, 15.0E+00, 15.0E+00 /)
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The right hand side B is '
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  Call SGEFA to factor the matrix A.
!
  call sgefa ( a, lda, n, ipvt, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  SGEFA returned an error flag INFO = ', info
    return
  end if
!
!  If no error occurred, now use SGESL to solve the system.
!
  job = 0
  call sgesl ( a, lda, n, ipvt, b, job )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)') '  SGESL returns the solution:'
  write ( *, '(a)' ) '  (Should be (1,1,1))'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
 
  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 solves a large linear system using SGEFA and SGESL.
!
  implicit none
!
  integer, parameter :: n = 100
  integer, parameter :: lda = n
!
  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer ipvt(n)
  integer j
  integer job
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  SGEFA/SGESL'
  write ( *, '(a)' ) '  General storage factor/solve.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Assign values to the matrix A and the right hand side B.
!
!  The problem is just an enlarged version of the
!  problem for N = 5, which is:
!
!  Matrix A is ( n -1 -1 -1 -1)    Right hand side B is  (1)
!              (-1  n -1 -1 -1)                          (1)
!              (-1 -1  n -1 -1)                          (1)
!              (-1 -1 -1  n -1)                          (1)
!              (-1 -1 -1 -1  n)                          (1)
!
!  Solution is   (1)
!                (1)
!                (1)
!                (1)
!                (1)
!
  b(1:n) = 1.0E+00

  a(1:n,1:n) = -1.0E+00
  do i = 1, n
    a(i,i) = n
  end do
!
!  Call SGEFA to factor the matrix A.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGEFA to factor the matrix.'
 
  call sgefa ( a, lda, n, ipvt, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  SGEFA returned an error flag INFO = ', info
    return
  end if
!
!  Call SGESL to solve the linear system.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGESL to solve the factored system.'
 
  job = 0
  call sgesl ( a, lda, n, ipvt, b, job )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first five entries of the solution:'
  write ( *, '(a)' ) '  (All of them should be 1.)'
  write ( *, '(a)' ) ' '
  do i = 1, 5
    write ( *, '(g14.6)' ) b(i)
  end do
 
  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests SPOFA and SPODI, for a positive definite symmetric matrix.
!
!  SPOFA factors a positive definite symmetric matrix,
!  and the SPODI can compute the determinant, and/or the inverse.
!
  implicit none
!
  integer, parameter :: n = 5
  integer, parameter :: lda = n
!
  real a(lda,n)
  real det(2)
  integer i
  integer info
  integer j
  integer job
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  SPOFA/SPODI,'
  write ( *, '(a)' ) '  factor/inverse for positive definite matrices.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Set the matrix A.
!
  a(1:n,1:n) = 0.0E+00

  do i = 1, n
    a(i,i) = 2.0E+00
    if ( i > 1 ) then
      a(i,i-1) = -1.0E+00
    end if
    if ( i < n ) then
      a(i,i+1) = -1.0E+00
    end if
  end do
!
!  Factor the matrix using SPOFA.
!
  call spofa ( a, lda, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Error, SPOFA returns INFO = ', info
    return
  end if
!
!  Invert the matrix using SPODI.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Calling SPODI for determinant and inverse.'
 
  job = 11
  call spodi ( a, lda, n, det, job )
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a,g14.6)' ) &
    '  Determinant  = ', det(1), ' * 10 ** ', det(2)
!
!  SPODI produces only the 'upper half triangle' of the inverse,
!  which is actually symmetric.  Thus, the lower half could be
!  produced by copying from the upper half.  However, the first row
!  of A, as returned, is exactly the first row of the inverse.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First row of inverse:'
  write ( *, '(a)' ) ' '
  write ( *, '(5g14.6)' ) a(1,1:n)
 
  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests SPBFA and SPBSL, for a positive definite band matrix.
!
  implicit none
!
  integer, parameter :: n = 10
  integer, parameter :: lda = 2
!
  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer m
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  SPBFA and SPBSL,'
  write ( *, '(a)' ) '  Factor/solve for positive definite band matrices.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Assign values to matrix A and right hand side B.
!
!  The problem is just an enlarged version of the
!  problem for N = 5, which is:
!
!  Matrix A is ( 2 -1  0  0  0)    right hand side B is  (1)
!              (-1  2 -1  0  0)                          (0)
!              ( 0 -1  2 -1  0)                          (0)
!              ( 0  0 -1  2 -1)                          (0)
!              ( 0  0  0 -1  2)                          (1)
!
!
!  solution is   (1)
!                (1)
!                (1)
!                (1)
!                (1)
!
!  Set the right hand side.
!
  b(1) = 1.0E+00
  b(2:n-1) = 0.0E+00
  b(n) = 1.0E+00
!
!  Set the number of nonzero diagonals.
!
  m = 1
!
!  Set the value of the subdiagonal and diagonal.
!
  a(1,1:n) = -1.0E+00
  a(2,1:n) = 2.0E+00
!
!  Call SPBFA to factor the matrix A.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SPBFA to factor the matrix.'
 
  call spbfa ( a, lda, n, m, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Error!  SPBFA returns INFO = ',info
    return
  end if
!
!  Call SPBSL to solve the linear system A*X = B.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SPBSL to solve the system.'
 
  call spbsl ( a, lda, n, m, b )
!
!  Print the result.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first 5 entries of the solution:'
  write ( *, '(a)' ) '  (All should be 1):'
  write ( *, '(a)' ) ' '
 
  do i = 1, 5
    write ( *, '(g14.6)' ) b(i)
  end do
 
  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests SGBFA and SGBSL, for a general banded matrix.
!
!
!  The problem solved here is a larger version of this one:
!
!  Matrix A is ( 2 -1  0  0  0)    right hand side B is  (1)
!              (-1  2 -1  0  0)                          (0)
!              ( 0 -1  2 -1  0)                          (0)
!              ( 0  0 -1  2 -1)                          (0)
!              ( 0  0  0 -1  2)                          (1)
!
!
!  Solution is   (1)
!                (1)
!                (1)
!                (1)
!                (1)
!
!
!  N is the number of equations.
!
!  ML is the number of subdiagonals,
!  MU the number of superdiagonals.
!
!  LDA is the leading dimension of the array used to store the
!  matrix, which must be at least 2*ML+MU+1.
!
  implicit none
!
  integer, parameter :: n = 10
  integer, parameter :: ml = 1
  integer, parameter :: mu = 1
  integer, parameter :: lda = 2*ml+mu+1
!
  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer ipivot(n)
  integer j
  integer job
  integer m
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  SGBFA factors a general band matrix.'
  write ( *, '(a)' ) '  SGBSL solves a linear system associated with'
  write ( *, '(a)' ) '    a factored general band matrix.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Set the right hand side B.
!
  b(1) = 1
  b(2:n-1) = 0.0E+00
  b(n) = 1.0E+00
!
!  Set the matrix A.
!
  m = ml + mu + 1
  write ( *, '(a,i6)' ) '  The bandwidth of the matrix is ', m
 
  a(m-1,2:n  ) = -1.0E+00
  a(m,  1:n  ) =  2.0E+00
  a(m+1,1:n-1) = -1.0E+00
!
!  Call SGBFA to factor the matrix A.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Calling SGBFA to factor the matrix.'
 
  call sgbfa ( a, lda, n, ml, mu, ipivot, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Error!  SGBFA returns INFO = ', info
    return
  end if
!
!  Call SGBSL to solve the linear system.  The solution
!  is returned in B, that is, it overwrites the right hand side.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Calling SGBSL to solve the linear system.'
 
  job = 0
  call sgbsl ( a, lda, n, ml, mu, ipivot, b, job )
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first 5 entries of the solution:'
  write ( *, '(a)' ) '  (All should be 1):'
  write ( *, '(a)' ) ' '
 
  do i = 1, 5
    write ( *, '(g14.6)' ) b(i)
  end do
 
  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests SGBFA and SGBSL, for general banded matrices.
!
  implicit none
!
  integer, parameter :: n = 100
  integer, parameter :: ml = 25
  integer, parameter :: mu = 25
  integer, parameter :: lda = 2*ml+mu+1
!
  real a(lda,n)
  real b(n)
  integer i
  integer ihi
  integer ilo
  integer info
  integer ipivot(n)
  integer j
  integer job
  integer m
  real temp
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  SGBFA and SGBSL,'
  write ( *, '(a)' ) '  general band matrix factor/linear solve.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Assign values to matrix A and right hand side B.
!
!  We want to try a problem with a significant bandwidth.
!
  m = ml+mu+1
  write ( *, '(a,i6)' ) '  The bandwidth of the matrix is ', m
 
  do j = 1, n
 
    ilo = max(1,j-mu)
    ihi = min(n,j+ml)
 
    temp = 0.0E+00
    do i = ilo,ihi
      a(i-j+m,j) = -1.0E+00
      temp = temp - 1.0E+00
    end do
 
    temp = temp + 1.0E+00
    a(m,j) = 4.0E+00 - temp
    b(j) = 4.0E+00
 
  end do
!
!  Call SGBFA to factor the matrix A.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGBFA to factor the matrix.'
 
  call sgbfa ( a, lda, n, ml, mu, ipivot, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Error!  SGBFA returns INFO = ', info
    return
  end if
!
!  Call SGBSL to solve the linear system.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGBSL to solve the linear system.'
 
  job = 0
  call sgbsl ( a, lda, n, ml, mu, ipivot, b, job )
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first 5 entries of the solution:'
  write ( *, '(a)' ) '  (All should be 1):'
  write ( *, '(a)' ) ' '
  do i = 1, 5
    write ( *, '(g14.6)' ) b(i)
  end do
 
  return
end
subroutine test07
!
!*******************************************************************************
!
!! TEST07 tests SSIFA.
!! TEST07 tests SSISL.
!! TEST07 tests SSIDI.
!
  implicit none
!
  integer, parameter :: n = 100
  integer, parameter :: lda = n
!
  real a(lda,n)
  real b(n)
  real det(2)
  integer i
  integer inert(3)
  integer info
  integer ipvt(n)
  integer j
  integer job
  real work(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  SSIFA factors a symmetric indefinite matrix,'
  write ( *, '(a)' ) '  SSISL solves a linear system associated with a'
  write ( *, '(a)' ) '    factored symmetric indefinite matrix.'
  write ( *, '(a)' ) '  SSIDI computes the determinant of a factored'
  write ( *, '(a)' ) '    symmetric indefinite matrix.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Assign values to the matrix A and the right hand side B.
!
  b(1:n-1) = 0.0E+00
  b(n)= n+1

  a(1:n,1:n) = 0.0E+00

  do i = 1, n
    a(i,i) = 2.0E+00
    if ( i < n ) then
      a(i,i+1) = -1.0E+00
    end if
  end do
!
!  Call SSIFA to factor the matrix A.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SSIFA to factor the matrix.'
 
  call ssifa ( a, lda, n, ipvt, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Error.  SSIFA returns INFO = ', info
    return
  end if
!
!  Call SSISL to solve the linear system.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SSISL to solve the linear system.'
 
  call ssisl ( a, lda, n, ipvt, b )
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The first 5 entries of the solution:'
  write ( *, '(a)' ) '  (Should be (1,2,3,4,5))'
  write ( *, '(a)' ) ' '
 
  do i = 1, 5
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  Call SSIDI for the determinant.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SSISL for the determinant.'
 
  job = 010

  call ssidi ( a, lda, n, ipvt, det, inert, work, job )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a,g14.6)' ) &
    '  The determinant is ', det(1), ' * 10** ', det(2)

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08 tests the QR factor routine SQRDC.
!! TEST08 tests the QR solve routine SQRSL.
!
  implicit none
!
  integer, parameter :: n = 3
  integer, parameter :: p = 3
  integer, parameter :: lda = n
!
  real a(lda,p)
  real b(lda,p)
  integer i
  integer info
  integer ipvt(p)
  integer j
  integer job
  integer k
  real q(n,n)
  real qraux(p)
  real qty(n)
  real qy(n)
  real r(n,p)
  real rsd(n)
  real work(p)
  real xb(n)
  real y(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  SQRDC computes the QR decomposition of a '
  write ( *, '(a)' ) '  matrix, but does not return Q and R explicitly.'
  write ( *, '(a)' ) '  Show how Q and R can be recovered using SQRSL.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix dimensions are '
  write ( *, '(a,i6)' ) '  N = ', n
  write ( *, '(a,i6)' ) '  P = ', p
!
!  Set the matrix A.
!
  a(1,1) = 1.0E+00
  a(2,1) = 1.0E+00
  a(3,1) = 0.0E+00
 
  a(1,2) = 1.0E+00
  a(2,2) = 0.0E+00
  a(3,2) = 1.0E+00
 
  a(1,3) = 0.0E+00
  a(2,3) = 1.0E+00
  a(3,3) = 1.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The original matrix A:'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:p)
  end do
!
!  Call SQRDC to compute the QR decomposition of A.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SQRDC to compute the QR factors.'
 
  job = 0
  call sqrdc ( a, lda, n, p, qraux, ipvt, work, job )
!
!  Print out what SQRDC has stored in A...
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The packed matrix A which describes Q and R:'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:p)
  end do
!
!  ...and in QRAUX.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The QRAUX vector, containing some additional'
  write ( *, '(a)' ) '  information defining Q:'
  write ( *, '(a)' ) ' '
 
  write ( *, '(5g14.6)' ) qraux(1:n)
!
!  Print out the resulting R factor.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The R factor:'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    do j = 1, p
 
      if ( j < i ) then
        r(i,j) = 0.0E+00
      else
        r(i,j) = a(i,j)
      end if
 
    end do
 
    write ( *, '(5g14.6)' ) r(i,1:p)
 
  end do
!
!  Call SQRSL to extract the information about the Q matrix.
!  We do this, essentially, by asking SQRSL to tell us the
!  value of Q*Y, where Y is a column of the identity matrix.
!
  job = 10000
 
  do i = 1, n
!
!  Set the vector Y.
!
    y(1:n) = 0.0E+00
 
    y(i) = 1.0E+00
!
!  Ask SQRSL to tell us what Q*Y is.
!
    call sqrsl ( a, lda, n, p, qraux, y, qy, qty, b, rsd, xb, job, info )
 
    if ( info /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) '  Error!  SQRSL returns INFO = ',info
      return
    end if
!
!  Copy QY into the appropriate column of Q.
!
    q(1:n,i) = qy(1:n)
 
  end do
!
!  Now print out the Q matrix we have extracted.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The Q factor:'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(5g14.6)' ) q(i,1:n)
  end do
!
!  Compute Q*R to verify that it equals A.
!
  b(1:n,1:p) = matmul ( q(1:n,1:n), r(1:n,1:p) )
!
!  Print the result.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The product, Q*R, = '
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(5g14.6)' ) b(i,1:p)
  end do
 
  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09 demonstrates SGTSL.
!
  implicit none
!
  integer, parameter :: n = 100
!
  real b(n)
  real c(n)
  real d(n)
  real e(n)
  integer i
  integer info
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  SGTSL,'
  write ( *, '(a)' ) '    general tridiagonal factor/solve routine.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Set up the linear system, by storing the values of the
!  subdiagonal, diagonal, and superdiagonal in C, D, and E,
!  and the right hand side in B.
!
  c(1) = 0.0E+00
  c(2:n) = -1.0E+00

  d(1:n) = 2.0E+00

  e(1:n-1) = -1.0E+00
  e(n) = 0.0E+00

  b(1:n-1) = 0.0E+00
  b(n) = real ( n + 1 )
!
!  Call SGTSL to factor and solve the system in one step.
!
  call sgtsl ( n, c, d, e, b, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST09 - Error!'
    write ( *, '(a,i6)' ) '  SGTSL returns nonzero INFO = ', info
    return
  end if
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First 5 entries of solution:'
  write ( *, '(a)' ) '  (Should be (1,2,3,4,5))'
  write ( *, '(a)' ) ' '
 
  do i = 1, 5
    write ( *, '(g14.6)' ) b(i)
  end do
 
  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 tests SGECO.
!! TEST10 tests SGEFA.
!! TEST10 tests SGESL.
!! TEST10 tests SGEDI.
!
  implicit none
!
  integer, parameter :: lda = 10
!
  real a(lda,lda)
  real b(lda)
  real det(2)
  integer i
  integer ipvt(lda)
  integer j
  integer job
  integer n
  real rcond
  real work(lda)
  real z(lda)
!
  n = 3

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  Demonstrate the use of LINPACKS routines'
  write ( *, '(a)' ) '  SGECO, SGEFA, SGESL and SGEDI.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Set the values of the matrix A.
!
  a(1,1) = 1.0E+00
  a(1,2) = 2.0E+00
  a(1,3) = 3.0E+00
 
  a(2,1) = 4.0E+00
  a(2,2) = 5.0E+00
  a(2,3) = 6.0E+00
 
  a(3,1) = 7.0E+00
  a(3,2) = 8.0E+00
  a(3,3) = 0.0E+00
!
!  Call SGECO to factor the matrix and compute its condition number.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGECO to factor the matrix, and compute'
  write ( *, '(a)' ) '  its condition number.'
 
  call sgeco ( a, lda, n, ipvt, rcond, z )
 
  write ( *, '(a,g14.6)' ) 'The reciprocal matrix condition number = ',rcond
 
  if ( rcond + 1.0E+00 == 1.0E+00 ) then
    write ( *, '(a)') ' '
    write ( *, '(a)' ) '  Error!  The matrix is nearly singular!'
    return
  end if
!
!  Set a right hand side.
!
  b(1:3) = (/ 6.0E+00, 15.0E+00, 15.0E+00 /)
!
!  Solve the linear system.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGESL to solve a linear system.'
 
  job = 0
  call sgesl ( a, lda, n, ipvt, b, job )
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Solution returned by SGESL'
  write ( *, '(a)' ) '  (Should be (1,1,1))'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  A second right hand side can be solved without refactoring a.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGESL for a new right hand '
  write ( *, '(a)' ) '  side for the same, factored matrix.'
!
!  Set the right hand side.
!
  b(1:3) = (/ 1.0E+00, 4.0E+00, 7.0E+00 /)
!
!  Solve the system.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGESL to solve a linear system.'
 
  job = 0
  call sgesl ( a, lda, n, ipvt, b, job )
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Solution returned by SGESL'
  write ( *, '(a)' ) '  (should be (1,0,0))'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  The transposed problem transpose(A)*X = B can be solved by SGESL
!  as well, without any refactoring.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGESL for transposed problem.'
!
!  Set the right hand side.
!
  b(1:3) = (/ 6.0E+00, 6.0E+00, -3.0E+00 /)
!
!  Solve the transposed system.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGESL to solve a transposed linear system.'
 
  job = 1
  call sgesl ( a, lda, n, ipvt, b, job )
!
!  Print the results.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Solution returned by SGESL'
  write ( *, '(a)' ) '  (should be (-1,0,1))'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  Now call SGEDI for the inverse and determinant.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SGEDI to get inverse and determinant'
 
  job = 11
  call sgedi ( a, lda, n, ipvt, det, work, job )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a,g14.6)' ) &
    '  The determinant = ',det(1),' * 10 ** ',det(2)
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The inverse matrix:'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(5g14.6)') a(i,1:n)
  end do
 
  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11 tests SCHDC.
!
  implicit none
!
  integer, parameter :: n = 4
  integer, parameter :: lda = n

  real a(lda,n)
  integer i
  integer info
  integer ipvt(n)
  integer j
  integer job
  integer work(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  SCHDC computes the Cholesky decomposition of'
  write ( *, '(a)' ) '  a general storage positive definite symmetric matrix:'
  write ( *, '(a)' ) '    A = U'' * U'
  write ( *, '(a)' ) '  where U is an upper triangular matrix.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The matrix order is N = ', n
!
!  Choose a desired factor U:
!
  call random_number ( harvest = a(1:n,1:n) )

  do i = 1, n
    a(i,1:i-1) = 0.0E+00
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The desired Cholesky factor U:'
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:n)
  end do
!
!  Construct the matrix A.
!
  a(1:n,1:n) = matmul ( transpose ( a(1:n,1:n) ), a(1:n,1:n) )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The matrix A:'
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:n)
  end do
!
!  Call SCHDC to factor the matrix A.
!
  job = 0

  call schdc ( a, lda, n, work, ipvt, job, info )

  if ( info /= n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST11'
    write ( *, '(a,i6)' ) '  SCHDC returned INFO = ', info
    write ( *, '(a)' ) '  This means the matrix is not positive definite.'
  end if
!
!  Only for clarity in printing, we zero out the strict lower triangle.
!
  do i = 2, n
    a(i,1:i-1) = 0.0E+00
  end do
!
!  Now we print the factorization matrix U.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The computed Cholesky factor U:'
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:n)
  end do
!
!  Compute the Cholesky product.
!
  a(1:n,1:n) = matmul ( transpose ( a(1:n,1:n) ), a(1:n,1:n) )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The product U''*U: '
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:n)
  end do

  return
end