!  dblas_prb.f90  15 February 2001
!
program dblas_prb
!
!*******************************************************************************
!
!! DBLAS_PRB tests the DBLAS routines.
!
  character ( len = 8 ) date
  character ( len = 10 ) time
!
  call date_and_time ( date, time )

  write ( *, * ) ' '
  write ( *, * ) 'DBLAS_PRB'
  write ( *, * ) '  Tests for DBLAS,'
  write ( *, * ) '  the Double Precision basic linear algebra subprograms.'
  write ( *, * ) ' '
  write ( *, * ) '  Today''s date: ', date
  write ( *, * ) '  Today''s time: ', time
 
  call test05
  call test99
  call test01
  call test02
  call test03
  call test04
 
  write ( *, * ) ' '
  write ( *, * ) 'DBLAS_PRB'
  write ( *, * ) '  Normal end of DBLAS tests.'

  stop
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 demonstrates IDAMAX.
!
  integer, parameter :: n = 11
!
  integer i
  integer i1
  integer incx
  integer idamax
  double precision x(n)
!
  write ( *, * ) ' '
  write ( *, * ) 'TEST05'
  write ( *, * ) '  IDAMAX returns the index of maximum magnitude;'
  write ( *, * ) ' '
 
  do i = 1, n
    x(i) = dble ( mod ( 7 * i, 11 ) ) - dble ( n / 2 )
  end do

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

  incx = 1

  i1 = idamax ( n, x, incx )

  write ( *, * ) ' '
  write ( *, * ) '  The index of maximum magnitude = ', i1

  return
end
subroutine test99
!
!*******************************************************************************
!
!! TEST99 tests SASUM.
!
  integer, parameter :: lda = 6
  integer, parameter :: ma = 5
  integer, parameter :: na = 4
  integer, parameter :: nx = 10
!
  double precision a(lda,na)
  double precision dasum
  integer i
  integer j
  double precision x(nx)
!
  do i = 1, nx
    x(i) = (-1.0D+00)**i * real ( 2 * i )
  end do

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

  write ( *, * ) ' '
  write ( *, * ) '  DASUM ( NX, X, 1 ) = ', dasum ( nx, x, 1 )
  write ( *, * ) '  DASUM ( NX/2, X, 2 ) = ', dasum ( nx/2, x, 2 )
  write ( *, * ) '  DASUM ( 2, X, NX/2 ) = ', dasum ( 2, x, nx/2 )

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

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

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

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

  return
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 demonstrates DDOT.
!
  integer, parameter :: n = 5
  integer, parameter :: lda = 10
  integer, parameter :: ldb = 7
  integer, parameter :: ldc = 6
!
  double precision a(lda,lda)
  double precision b(ldb,ldb)
  double precision c(ldc,ldc)
  integer i
  integer j
  double precision ddot
  double precision sum
  double precision x(n)
  double precision y(n)
!
  write ( *, * ) ' '
  write ( *, * ) 'TEST01'
  write ( *, * ) '  DDOT computes the dot product of vectors.'
  write ( *, * ) ' '
 
  do i = 1, n
    x(i) = dble ( i )
    y(i) = dble ( i )
    do j = 1, n
      a(i,j) = dble ( i + j )
      b(i,j) = dble ( 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) = dble ( i )
    y(i) = - dble ( i )
  end do

  sum = ddot ( 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) = dble ( i )
  end do

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

  sum = ddot ( 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 = ddot ( 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 ddot:
!
  do i = 1, n
    do j = 1, n
      c(i,j) = ddot ( n, a(i,1), lda, b(1,j), 1 )
    end do
  end do
 
  write ( *, * ) ' '
  write ( *, * ) '  Matrix product computed with ddot:'
  write ( *, * ) ' '
  do i = 1, n
    write ( *, '(5g14.6)' ) c(i,1:n)
  end do
 
  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 demonstrates DNRM2.
!
  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.
!
  double precision a(lda,lda)
  integer i
  integer incx
  integer j
  double precision dnrm2
  double precision sum
  double precision x(n)
!
  write ( *, * ) ' '
  write ( *, * ) 'TEST02'
  write ( *, * ) '  DNRM2 computes the Euclidean norm of a vector.'
  write ( *, * ) ' '
!
!  Compute the euclidean norm of a vector:
!
  do i = 1, n
    x(i) = dble ( i )
  end do

  write ( *, * ) ' '
  write ( *, * ) '  The vector X:'
  write ( *, * ) ' '
  do i = 1, n
    write ( *, '( i6, f8.4 )' ) i, x(i)
  end do
  incx = 1
  write ( *, * ) '  The 2-norm of X is ', dnrm2 ( 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) = dble ( i + j )
    end do
  end do

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

  incx = 1
  sum = dnrm2 ( n, a(1,2), incx )
  write ( *, * ) '  The 2-norm of column 2 of A is ', sum
 
  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests IDAMAX.
!! TEST03 tests DAXPY.
!! TEST03 tests DDOT.
!! TEST03 tests DSCAL.
!
  integer, parameter :: n = 10
  integer, parameter :: lda = n
!
  double precision a(lda,n)
  double precision b(n)
  integer i
  integer info
  integer ipvt(n)
  integer idamax
  integer j
  integer k
  integer l
  double precision t
!
  write ( *, * ) ' '
  write ( *, * ) 'TEST03'
  write ( *, * ) '  Use IDAMAX, DAXPY, DDOT and DSCAL'
  write ( *, * ) '  in a Gauss elimination routine.'
  write ( *, * ) ' '
!
!  Set the matrix.
!
  do i = 1, n
    do j = 1, n

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

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

  info = 0

  do k = 1, n - 1

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

    if ( a(l,k) == 0.0D+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.0D+00 / a(k,k)
      call dscal ( 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 daxpy ( 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.0D+00 ) then
    info = n
  end if

  if ( info /= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) '  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 daxpy ( 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 daxpy ( k-1, t, a(1,k), 1, b(1), 1 )
  end do

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

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests a Cholesky solver.
!
!
!  Solve a positive definite symmetric linear system A*X = B.
!
  integer, parameter :: n = 20
!
  double precision a(n,n)
  integer i
  integer info
  integer j
  double precision x(n)
!
  write ( *, * ) ' '
  write ( *, * ) 'TEST04'
  write ( *, * ) '  Demonstrate Cholesky factor/solve routines'
  write ( *, * ) '  built out of level 1, 2 and 3 blas'
  write ( *, * ) ' '
!
!  Set the entries of the matrix.
!
  do i = 1, n
 
    a(i,1:n) = 0.0D+00

    a(i,i) = 2.0D+00

    if ( i > 1 ) then
      a(i,i-1) = -1.0D+00
    end if

    if ( i < n ) then
      a(i,i+1) = -1.0D+00
    end if

  end do
!
!  Set the right hand side vector.
!
  x(1:n-1) = 0.0D+00
  x(n) = dble ( n + 1 )
!
!  Factor the matrix.
!
  call dlltb ( n, a, n, info )
 
  if ( info /= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) '  The matrix is singular!'
    return
  end if
!
!  Solve the system.
!
  call dllts ( n, a, n, x )

  write ( *, * ) ' '
  write ( *, * ) '  First 10 entries of solution:'
  write ( *, * ) ' '
  write ( *, '(5g14.6)' ) x(1:10)
 
  return
end
subroutine dlltb ( n, a, lda, info )
!
!*******************************************************************************
!
!! DLLTB 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.
!
!  DLLT is an equivalent routine, but it solves the problem
!  at a lower level.  
!
!  DLLTB solves the problem in chunks,
!  which is hoped to allow for greater optimization.
!
  integer lda
  integer n

  integer, parameter :: nb = 64
!
  double precision 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 dsyrk ( 'lower', 'no transpose', jb, j-1, -1.0D+00, a(j,1), lda, &
      1.0D+00, a(j,j), lda )
!
!  Factorize diagonal block and test for non-positive-definiteness.
!
    call dllt ( 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 dgemm ( 'no transpose', 'transpose', n-j-jb+1, jb, j-1, &
        -1.0D+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 dtrsm ( 'right', 'lower', 'transpose', 'non-unit', n-j-jb+1, &
        jb, 1.0D+00, a(j,j), lda, a(j+jb,j), lda )

    end if
 
  end do
 
  return
end
subroutine dllt ( n, a, lda, info )
!
!*******************************************************************************
!
!! DLLT computes an l*transpose(l) factorization of a symmetric positive-
!  definite matrix a.
!
!  dllt uses level 2 and level 1 blas routines.
!
  integer lda
  integer n
!
  double precision a(lda,n)
  integer info
  integer j
  double precision ddot
!
  info = 0
 
  do j = 1, n
!
!  Update A(J,J).
!
    a(j,j) = a(j,j) - ddot ( 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.0D+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 dgemv ( 'no transpose', n-j, j-1, -1.0D+00, a(j+1,1), lda, &
        a(j,1), lda, 1.0D+00, a(j+1,j), 1 )
!
!  Compute elements J+1 to N of J-th column of L.
!
      call dscal ( n-j, 1.0D+00/a(j,j), a(j+1,j), 1 )

    end if
 
  end do
 
  return
end
subroutine dllts ( n, a, lda, b )
!
!*******************************************************************************
!
!! DLLTS solves A*X=B, for a positive definite symmetric matrix.
!
!  The matrix A should have been factored by DLLT or DLLTB.
!
  integer lda
  integer n
!
  double precision a(lda,n)
  double precision b(n)
  integer k
  double precision ddot
  double precision t
!
  do k = 1, n
    t = ddot ( k-1, a(k,1), lda, b(1), 1 )
    b(k) = ( b(k) - t ) / a(k,k)
  end do
 
  do k = n, 1, -1
    b(k) = b(k) / a(k,k)
    t = - b(k)
    call daxpy ( k-1, t, a(k,1), lda, b(1), 1 )
  end do
 
  return
end