program lapacks_prb
!
!***********************************************************************
!
!! LAPACKS_PRB calls the single precision real LAPACK test routines.
!
  implicit none
!
  call timestamp ( )

  write ( *, * ) ' '
  write ( *, * ) 'LAPACKS_PRB'
  write ( *, * ) '  Tests for the LAPACKS linear algebra library.'
  write ( *, * ) '  Single precision version.'

  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'LAPACKS_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!***********************************************************************
!
!! TEST01 tests SGECON.
!! TEST01 tests SGETRF.
!! TEST01 tests SGETRS.
!! TEST01 tests SGETRI.
!
  implicit none
!
  integer, parameter :: n = 3

  integer, parameter :: lda = n
  integer, parameter :: lwork = 4 * n
!
  real a(lda,n)
  real anorm
  real b(n)
  integer i
  integer info
  integer ipiv(n)
  integer iwork(n)
  integer j
  real rcond
  real work(4*n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  SGETRF factors a general matrix;'
  write ( *, '(a)' ) '  SGETRS solves a linear system;'
  write ( *, '(a)' ) '  SGECON computes the condition number;'
  write ( *, '(a)' ) '  SGETRI computes the inverse.'
!
!  Set the matrix.
!
  a(1,1) = 1.0
  a(1,2) = 2.0
  a(1,3) = 3.0

  a(2,1) = 4.0
  a(2,2) = 5.0
  a(2,3) = 6.0

  a(3,1) = 7.0
  a(3,2) = 8.0
  a(3,3) = 0.0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  System matrix A:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:n)
  end do
!
!  Compute the L-Infinity norm of the matrix.
!
  call rmat_norm_li ( n, n, anorm )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Matrix L-infinity norm is ', anorm
!
!  Set the right hand side.
!
  b(1:3) = (/ 14.0, 32.0, 23.0 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Right hand side b:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  Factor the matrix.
!
  call sgetrf ( n, n, a, lda, ipiv, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  SGETRF returned INFO = ', info
    write ( *, '(a)' ) '  The matrix is numerically singular.'
    return
  end if
!
!  Get the condition number.
!
  call sgecon ( 'I', n, a, lda, anorm, rcond, work, iwork, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Condition number calculation failed!'
    write ( *, '(a,i6)' ) '  INFO = ', info
    return
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Matrix condition number = ', rcond
!
!  Solve the linear system.
!
  call sgetrs ( 'n', n, 1, a, lda, ipiv, b, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Solution procedure failed!'
    write ( *, '(a,i6)' ) '  INFO = ', info
    return
  end if
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Solution X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  Compute the inverse matrix.
!
  call sgetri ( n, a, lda, ipiv, work, lwork, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Inversion procedure failed!'
    write ( *, '(a,i6)' ) '  INFO = ', info
    return
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Inverse matrix:'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    write ( *, '(5g14.6)' ) a(i,1:n)
  end do
 
  return
end
subroutine test02
!
!***********************************************************************
!
!! TEST02 tests SGETRF.
!! TEST02 tests SGETRS.
!
!  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)
!
!  For this problem, no pivoting is required.
!
  implicit none
!
  integer, parameter :: n = 25
  integer, parameter :: lda = n
!
  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer ipiv(n)
  integer j
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  SGETRF factors a general matrix;'
  write ( *, '(a)' ) '  SGETRS solves a linear system;'
  write ( *, '(a)' ) ' '
!
!  Assign values to matrix A and right hand side b.
!
  do i = 1, n
    do j = 1, n
      if ( i == j ) then
        a(i,j) = n
      else
        a(i,j) = -1.0
      end if
    end do
  end do

  b(1:n) = 1.0
!
!  Factor the matrix.
!
  call sgetrf ( n, n, a, lda, ipiv, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Matrix is singular, INFO = ', info
    return
  end if
!
!  Solve the linear system.
!
  call sgetrs ( 'n', n, 1, a, lda, ipiv, b, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Solution procedure failed, INFO = ', info
    return
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First five entries of solution:'
  write ( *, '(a)' ) '  (all should be 1)'
  write ( *, '(a)' ) ' '
  write ( *, '(5g14.6)' ) b(1:5)
 
  return
end
subroutine test03
!
!***********************************************************************
!
!! TEST03 tests SPOTRF.
!! TEST03 tests SPOTRI.
!
  implicit none
!
  integer, parameter :: n = 25

  integer, parameter :: lda = n
!
  real a(lda,n)
  integer i
  integer info
  integer j
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  For a positive definite symmetric matrix,'
  write ( *, '(a)' ) '  SPOTRF factors;'
  write ( *, '(a)' ) '  SPOTRI computes the inverse.'
  write ( *, '(a)' ) ' '
!
!  Zero out the matrix.
!
  a(1:n,1:n) = 0.0
!
!  Subdiagonal.
!
  do i = 2, n
    a(i,i-1) = -1.0
  end do
!
!  Diagonal.
!
  do i = 1, n
    a(i,i) = 2.0
  end do
!
!  Superdiagonal.
!
  do i = 1, n - 1
    a(i,i+1) = -1.0
  end do
!
!  Factor the matrix.
!
  call spotrf ( 'u', n, a, lda, info )
 
  if ( info /= 0 ) then
    write ( *, '(a,i6)' ) '  SPOTRF returns INFO = ', info
    return
  end if
!
!  Compute the inverse matrix.
!
  call spotri ( 'u', n, a, lda, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Inversion procedure failed, INFO = ', info
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First row of inverse:'
  write ( *, '(a)' ) ' '
  write ( *, '(5g14.6)' ) a(1,1:n)
 
  return
end
subroutine test04
!
!***********************************************************************
!
!! TEST04 tests SPBTRF.
!! TEST04 tests SPBTRS.
!
!  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)
!
  implicit none
!
  integer, parameter :: n = 25
  integer, parameter :: nband = 1

  integer, parameter :: lda = nband + 1
!
  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer j
  integer nrhs
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  For a positive definite symmetric band matrix:'
  write ( *, '(a)' ) '  SPBTRF factors;'
  write ( *, '(a)' ) '  SPBTRS solves linear systems.'
  write ( *, '(a)' ) ' '
!
!  Zero out the matrix.
!
  a(1:lda,1:n) = 0.0
!
!  Super (and sub) diagonal.
!
  do i = 2, n
    a(1,i) = -1.0
  end do
!
!  Diagonal.
!
  do i = 1, n
    a(2,i) = 2.0
  end do
!
!  Set the right hand side.
!
  b(1) = 1.0
  do i = 2, n-1
    b(i) = 0.0
  end do
  b(n) = 1.0
!
!  Factor the matrix.
!
  call spbtrf ( 'u', n, nband, a, lda, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Factorization failed, INFO = ', info
    return
  end if
!
!  Solve the linear system.
!
  nrhs = 1
  call spbtrs ( 'u', n, nband, nrhs, a, lda, b, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Solution failed, INFO = ', info
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First 5 entries of solution (All should be 1):'
  write ( *, '(5g14.6)' ) b(1:5)
 
  return
end
subroutine test05
!
!***********************************************************************
!
!! TEST05 tests SGBTRF.
!! TEST05 tests SGBTRS.
!
!  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)
!
  implicit none
!
  integer, parameter :: n = 25
  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 ipiv(n)
  integer j
  integer m
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  SGBTRF factors a general band matrix.'
  write ( *, '(a)' ) '  SGBTRS solves a factored system.'
  write ( *, '(a)' ) ' '
!
!  Assign values to matrix A and right hand side b.
!
  b(1) = 1.0
  b(2:n-1) = 0.0
  b(n) = 1.0
!
!  Zero out the matrix.
!
  a(1:lda,1:n) = 0.0

  m = ml + mu + 1
!
!  Superdiagonal
!
  do j = 2, n
    a(m-1,j) = -1.0
  end do
!
!  Diagonal
!
  do j = 1, n
    a(m,j) = 2.0
  end do
!
!  Subdiagonal
!
  do j = 1, n - 1
    a(m+1,j) = -1.0
  end do
!
!  Factor the matrix.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Bandwidth is ', m
  write ( *, '(a)' ) ' '

  call sgbtrf ( n, n, ml, mu, a, lda, ipiv, info )

  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Factorization failed, INFO = ', info
    return
  end if
!
!  Solve the linear system.
!
  call sgbtrs ( 'n', n, ml, mu, 1, a, lda, ipiv, b, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Solution failed, INFO = ', info
    return
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First 5 entries of solution (All should be 1):'
  write ( *, '(5g14.6)' ) b(1:5)

  return
end
subroutine test06
!
!***********************************************************************
!
!! TEST06 tests SGBTRF.
!! TEST06 tests SGBTRS.
!
  implicit none
!
  integer, parameter :: n = 25
  integer, parameter :: ml = 3 
  integer, parameter :: mu = 3

  integer, parameter :: lda = 2 * ml + mu + 1
!
  real a(lda,n)
  real b(n)
  integer i
  integer ihi
  integer ilo
  integer info
  integer ipiv(n)
  integer j
  integer m
  real temp
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  SGBTRF factors a general band matrix.'
  write ( *, '(a)' ) '  SGBTRS solves a factored system.'
  write ( *, '(a)' ) ' '
!
!  Assign values to matrix A and right hand side b.
!
!  We want to try a problem with a significant bandwidth.
!
  m = ml + mu + 1

  do j = 1, n
    ilo = max ( 1, j-mu )
    ihi = min ( n, j+ml )
 
    temp = 0.0
    do i = ilo, ihi
      a(i-j+m,j) = -1.0
      temp = temp - 1.0
    end do
    temp = temp + 1.0
 
    a(m,j) = 4.0 - temp
  end do
!
!  Right hand side.
!
  b(1:n) = 4.0
!
!  Factor the matrix.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Bandwidth is ', m
 
  call sgbtrf ( n, n, ml, mu, a, lda, ipiv, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Factorization failed, INFO = ', info
    return
  end if
!
!  Solve the linear system.
!
  call sgbtrs ( 'n', n, ml, mu, 1, a, lda, ipiv, b, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Solution failed, INFO = ', info
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First 5 entries of solution (All should be 1):'
  write ( *, '(5g14.6)' ) b(1:5)
 
  return
end
subroutine test07
!
!***********************************************************************
!
!! TEST07 tests SGTSV.
!
  implicit none
!
  integer, parameter :: n = 100
!
  real b(n)
  real c(n-1)
  real d(n)
  real e(n-1)
  integer i
  integer info
  integer ldb
  integer nrhs
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  SGTSV factors and solves a linear system'
  write ( *, '(a)' ) '  with a general tridiagonal matrix.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The system is of order N = ', n
  write ( *, '(a)' ) ' '
!
!  Right hand side.
! 
  b(1:n-1) = 0.0
  b(n) = n + 1
!
!  Subdiagonal.
! 
  c(1:n) = -1.0
!
!  Diagonal.
!
  d(1:n) = 2.0
!
!  Superdiagonal.
!
  e(1:n) = -1.0
 
  nrhs = 1
  ldb = n
!
!  Factor and solve the linear system.
!
  call sgtsv ( n, nrhs, c, d, e, b, ldb, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Solution procedure failed.'
    write ( *, '(a,i6)' ) '  INFO = ', info
    return
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First 5 entries of solution:'
  write ( *, '(a)' ) '  (Should be 1, 2, 3, ...)'
  write ( *, '(5g14.6)' ) b(1:5)
 
  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  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
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
function rmat_norm_li ( m, n, a )
!
!*******************************************************************************
!
!! RMAT_NORM_LI returns the matrix L-infinity norm of an M by N matrix.
!
!
!  Definition:
!
!    The matrix L-infinity norm is defined as:
!
!      RMAT_NORM_LI =  Max ( 1 <= I <= M ) Sum ( 1 <= J <= N ) abs ( A(I,J) ).
!
!  Note:
!
!    The matrix L-infinity norm is derived from the vector L-infinity norm,
!    and satisifies:
!
!      rvec_norm_li ( A*x ) <= rmat_norm_li ( A ) * rvec_norm_li ( x ).
!
!  Modified:
!
!    24 March 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the number of rows in A.
!
!    Input, integer N, the number of columns in A.
!
!    Input, real A(M,N), the matrix whose L-infinity norm is desired.
!
!    Output, real RMAT_NORM_LI, the L-infinity norm of A.
!
  implicit none
!
  integer m
  integer n
!
  real a(m,n)
  integer i
  real rmat_norm_li
!
  rmat_norm_li = 0.0E+00

  do i = 1, m
    rmat_norm_li = max ( rmat_norm_li, sum ( abs ( a(i,1:n) ) ) )
  end do

  return
end