program sblas_prb
!
!*******************************************************************************
!
!! SBLAS_PRB tests the SBLAS routines.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SBLAS_PRB'
  write ( *, '(a)' ) '  Tests for SBLAS,'
  write ( *, '(a)' ) '  the Single Precision basic linear algebra subprograms.'
 
  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SBLAS_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 demonstrates ISAMAX.
!! TEST01 demonstrates ISAMIN.
!! TEST01 demonstrates ISMAX.
!! TEST01 demonstrates ISMIN.
!! TEST01 demonstrates SAMAX.
!! TEST01 demonstrates SAMIN.
!! TEST01 demonstrates SMAX.
!! TEST01 demonstrates SMIN.
!
  implicit none
!
  integer, parameter :: n = 11
!
  integer i
  integer i1
  integer i2
  integer i3
  integer i4
  integer incx
  integer isamax
  integer isamin
  integer ismin
  integer ismax
  real s1
  real s2
  real s3
  real s4
  real samax
  real samin
  real smax
  real smin
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  ISAMAX returns the index of maximum magnitude;'
  write ( *, '(a)' ) '  ISAMIN returns the index of minimum magnitude;'
  write ( *, '(a)' ) '  ISMAX returns the index of maximum value;'
  write ( *, '(a)' ) '  ISMIN returns the index of minimum value;'
  write ( *, '(a)' ) '  SAMAX returns the maximum magnitude;'
  write ( *, '(a)' ) '  SAMIN returns the minimum magnitude;'
  write ( *, '(a)' ) '  SMAX returns the maximum value;'
  write ( *, '(a)' ) '  SMIN returns the minimum value;'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    x(i) = real ( mod ( 7 * i, 11 ) ) - real ( n / 2 )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The vector X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '( i6, f8.4 )' ) i, x(i)
  end do

  incx = 1

  i1 = isamax ( n, x, incx )
  i2 = isamin ( n, x, incx )
  i3 = ismax ( n, x, incx )
  i4 = ismin ( n, x, incx )
  s1 = samax ( n, x, incx )
  s2 = samin ( n, x, incx )
  s3 = smax ( n, x, incx )
  s4 = smin ( n, x, incx )

  write ( *, '(a)' ) ' '
  write ( *, * ) '  The index of maximum magnitude = ', i1
  write ( *, * ) '  The maximum magnitude = ', s1
  write ( *, * ) '  The index of minimum magnitude = ', i2
  write ( *, * ) '  The minimum magnitude = ', s2
  write ( *, * ) '  The index of maximum value = ', i3
  write ( *, * ) '  The maximum value = ', s3
  write ( *, * ) '  The index of minimum value = ', i4
  write ( *, * ) '  The minimum value = ', s4

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests SASUM.
!
  implicit none
!
  integer, parameter :: lda = 6
  integer, parameter :: ma = 5
  integer, parameter :: na = 4
  integer, parameter :: nx = 10
!
  real a(lda,na)
  integer i
  integer j
  real sasum
  real x(nx)
!
  do i = 1, nx
    x(i) = (-1.0E+00)**i * real ( 2 * i )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Demonstrate the use of SASUM'
  write ( *, '(a)' ) '  which adds the absolute values of elements of a vector.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  X = '
  do i = 1, nx
    write ( *, * ) i, x(i)
  end do

  write ( *, * ) ' '
  write ( *, * ) '  SASUM ( NX, X, 1 )   = ', sasum ( nx,   x, 1 )
  write ( *, * ) '  SASUM ( NX/2, X, 2 ) = ', sasum ( nx/2, x, 2 )
  write ( *, * ) '  SASUM ( 2, X, NX/2 ) = ', sasum ( 2,    x, nx/2 )

  a(1:lda,1:na) = 0.0E+00

  do i = 1, ma
    do j = 1, na
      a(i,j) = (-1.0E+00)**(i+j) * real ( 10 * i + j )
    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Demonstrate with a matrix A:'
  write ( *, '(a)' ) ' '
  do i = 1, ma
    write ( *, * ) a(i,1:na)
  end do

  write ( *, '(a)' ) ' '
  write ( *, * ) '  SASUM ( MA, A(1,2), 1 )   = ', sasum ( ma, a(1,2), 1 )
  write ( *, * ) '  SASUM ( NA, A(2,1), LDA ) = ', sasum ( na, a(2,1), lda )

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 demonstrates SDOT.
!
  implicit none
!
  integer, parameter :: n = 5
  integer, parameter :: lda = 10
  integer, parameter :: ldb = 7
  integer, parameter :: ldc = 6
!
  real a(lda,lda)
  real b(ldb,ldb)
  real c(ldc,ldc)
  integer i
  integer j
  real sdot
  real sum
  real x(n)
  real y(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  SDOT computes the dot product of vectors.'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    do j = 1, n
      a(i,j) = real ( i + j )
      b(i,j) = real ( i - j )
    end do
  end do
!
!  To compute a simple dot product of two vectors, use a
!  call like this:
!
  do i = 1, n
    x(i) = real ( i )
    y(i) = - real ( i )
  end do

  sum = sdot ( n, x, 1, y, 1 )
  write ( *, * ) '  Dot product of X and Y is ', sum
!
!  To multiply a ROW of a matrix A times a vector X, we need to
!  specify the increment between successive entries of the row of A:
!
  do i = 1, n
    x(i) = real ( i )
  end do

  do i = 1, n
    do j = 1, n
      a(i,j) = real ( i + j )
    end do
  end do

  sum = sdot ( n, a(2,1), lda, x, 1 )
  write ( *, * ) '  Product of row 2 of A and X is ', sum
!
!  Product of a column of A and a vector is simpler:
!
  sum = sdot ( n, a(1,2), 1, x, 1 )
  write ( *, * ) '  Product of column 2 of A and X is ', sum
!
!  Here's how matrix multiplication, c = a*b, could be done
!  with SDOT:
!
  do i = 1, n
    do j = 1, n
      c(i,j) = sdot ( n, a(i,1), lda, b(1,j), 1 )
    end do
  end do
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Matrix product computed with SDOT:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(5g14.6)' ) c(i,1:n)
  end do
 
  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 demonstrates SNRM2.
!
  implicit none
!
  integer, parameter :: n = 5
  integer, parameter :: lda = n + 5
!
!  these parameters illustrate the fact that matrices are typically
!  dimensioned with more space than the user requires.
!
  real a(lda,lda)
  integer i
  integer incx
  integer j
  real snrm2
  real sum
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  SNRM2 computes the Euclidean norm of a vector.'
  write ( *, '(a)' ) ' '
!
!  Compute the euclidean norm of a vector:
!
  do i = 1, n
    x(i) = real ( i )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The vector X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '( i6, f8.4 )' ) i, x(i)
  end do
  incx = 1
  write ( *, * ) '  The 2-norm of X is ', snrm2 ( n, x, incx )
!
!  Compute the euclidean norm of a row or column of a matrix:
!
  do i = 1, n
    do j = 1, n
      a(i,j) = real ( i + j )
    end do
  end do

  incx = lda
  sum = snrm2 ( n, a(2,1), incx )
  write ( *, * ) '  The 2-norm of row 2 of A is ', sum

  incx = 1
  sum = snrm2 ( n, a(1,2), incx )
  write ( *, * ) '  The 2-norm of column 2 of A is ', sum
 
  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests ISAMAX.
!! TEST05 tests SAXPY.
!! TEST05 tests SDOT.
!! TEST05 tests SSCAL.
!
  implicit none
!
  integer, parameter :: n = 10
  integer, parameter :: lda = n
!
  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer ipvt(n)
  integer isamax
  integer j
  integer k
  integer l
  real t
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  Use ISAMAX, SAXPY, SDOT and SSCAL'
  write ( *, '(a)' ) '  in a Gauss elimination routine.'
  write ( *, '(a)' ) ' '
!
!  Set the matrix.
!
  do i = 1, n
    do j = 1, n

      if ( i == j ) then
        a(i,j) = 2.0E+00
      else if ( i == j + 1 ) then
        a(i,j) = - 1.0E+00
      else if ( i == j - 1 ) then
        a(i,j) = - 1.0E+00
      else
        a(i,j) = 0.0E+00
      end if

    end do
  end do
!
!  Set the right hand side.
!
  b(1:n-1) = 0.0E+00
  b(n) = real ( n ) + 1.0E+00

  info = 0

  do k = 1, n - 1

    l = isamax ( n-k+1, a(k,k), 1 ) + k - 1
    ipvt(k) = l

    if ( a(l,k) == 0.0E+00 ) then

      info = k

    else

      if ( l /= k ) then
        t = a(l,k)
        a(l,k) = a(k,k)
        a(k,k) = t
      end if

      t = -1.0E+00 / a(k,k)
      call sscal ( n-k, t, a(k+1,k), 1 )

      do j = k+1, n
        t = a(l,j)
        if ( l /= k ) then
          a(l,j) = a(k,j)
          a(k,j) = t
        end if
        call saxpy ( n-k, t, a(k+1,k), 1, a(k+1,j), 1 )
      end do

    end if

  end do

  ipvt(n) = n
  if (a(n,n) == 0.0E+00 ) then
    info = n
  end if

  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST05 - Fatal error!'
    write ( *, '(a)' ) '  The matrix is singular!'
    return
  end if

  do k = 1, n-1
    l = ipvt(k)
    t = b(l)
    if ( l /= k ) then
      b(l) = b(k)
      b(k) = t
    end if
    call saxpy ( n-k, t, a(k+1,k), 1, b(k+1), 1 )
  end do

  do k = n, 1, -1
    b(k) = b(k) / a(k,k)
    t = - b(k)
    call saxpy ( k-1, t, a(1,k), 1, b(1), 1 )
  end do

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

  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests a Cholesky solver.
!
!
!  Solve a positive definite symmetric linear system A*X = B.
!
  implicit none
!
  integer, parameter :: n = 20
!
  real a(n,n)
  integer i
  integer info
  integer j
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  SLLTB is a Cholesky factor routine,'
  write ( *, '(a)' ) '  SLLTS is a Cholesky solve routine,'
  write ( *, '(a)' ) '  built out of level 1, 2 and 3 BLAS.'
  write ( *, '(a)' ) ' '
!
!  Set the entries of the matrix.
!
  do i = 1, n
 
    a(i,1:n) = 0.0E+00

    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
!
!  Set the right hand side vector.
!
  x(1:n-1) = 0.0E+00
  x(n) = real ( n + 1 )
!
!  Factor the matrix.
!
  call slltb ( n, a, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST06 - Fatal error!'
    write ( *, '(a)' ) '  The matrix is singular!'
    return
  end if
!
!  Solve the system.
!
  call sllts ( n, a, n, x )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First 10 entries of solution:'
  write ( *, '(a)' ) ' '
  write ( *, '(5g14.6)' ) x(1:10)
 
  return
end
subroutine slltb ( n, a, lda, info )
!
!*******************************************************************************
!
!! SLLTB Cholesky factors a positive definite symmetric matrix.
!
!  The factorization has the form A = L*transpose(L).
!
!  The parameter NB determines the 'blocking factor' used
!  by the level 3 BLAS routines.
!
!  SLLT is an equivalent routine, but it solves the problem
!  at a lower level.  
!
!  SLLTB solves the problem in chunks,
!  which is hoped to allow for greater optimization.
!
  implicit none
!
  integer lda
  integer n

  integer, parameter :: nb = 64
!
  real a(lda,*)
  integer info
  integer j
  integer jb
!
  info = 0
 
  do j = 1, n, nb
 
    jb = min ( nb, n-j+1 )
!
!  Update diagonal block
!
    call ssyrk ( 'lower', 'no transpose', jb, j-1, -1.0E+00, a(j,1), lda, &
      1.0E+00, a(j,j), lda )
!
!  Factorize diagonal block and test for non-positive-definiteness.
!
    call sllt ( jb, a(j,j), lda, info )
 
    if ( info /= 0 )then
      info = info + j - 1
      return
    end if
 
    if ( j+jb <= n ) then
!
!  Update subdiagonal block.
!
      call sgemm ( 'no transpose', 'transpose', n-j-jb+1, jb, j-1, &
        -1.0E+00, a(j+jb,1), lda, a(j,1), lda, 1.0E+00, a(j+jb,j), lda )
!
!  Compute the subdiagonal block of L.
!
      call strsm ( 'right', 'lower', 'transpose', 'non-unit', n-j-jb+1, &
        jb, 1.0E+00, a(j,j), lda, a(j+jb,j), lda )

    end if
 
  end do
 
  return
end
subroutine sllt ( n, a, lda, info )
!
!*******************************************************************************
!
!! SLLT computes an l*transpose(l) factorization of a symmetric positive-
!  definite matrix a.
!
!  sllt uses level 2 and level 1 blas routines.
!
  implicit none
!
  integer lda
  integer n
!
  real a(lda,n)
  integer info
  integer j
  real sdot
!
  info = 0
 
  do j = 1, n
!
!  Update A(J,J).
!
    a(j,j) = a(j,j) - sdot ( j-1, a(j,1), lda, a(j,1), lda )
!
!  Compute L(J,J) and test for non-positive-definiteness.
!
    if ( a(j,j) <= 0.0E+00 ) then
      info = j
      return
    end if
 
    a(j,j) = sqrt ( a(j,j) )
!
!  Update elements J+1 to N of J-th column.
!
    if ( j  <  n ) then

      call sgemv ( 'no transpose', n-j, j-1, -1.0E+00, a(j+1,1), lda, &
        a(j,1), lda, 1.0E+00, a(j+1,j), 1 )
!
!  Compute elements J+1 to N of J-th column of L.
!
      call sscal ( n-j, 1.0E+00/a(j,j), a(j+1,j), 1 )

    end if
 
  end do
 
  return
end
subroutine sllts ( n, a, lda, b )
!
!*******************************************************************************
!
!! SLLTS solves A*X=B, for a positive definite symmetric matrix.
!
!  The matrix A should have been factored by SLLT or SLLTB.
!
  implicit none
!
  integer lda
  integer n
!
  real a(lda,n)
  real b(n)
  integer k
  real sdot
!
  do k = 1, n
    b(k) = ( b(k) - sdot ( k-1, a(k,1), lda, b(1), 1 ) ) / a(k,k)
  end do
 
  do k = n, 1, -1
    b(k) = b(k) / a(k,k)
    call saxpy ( k-1, -b(k), a(k,1), lda, b(1), 1 )
  end do
 
  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