有个家伙不会在HPCC上连接编译LAPACK, 稍微记一下.
module load LAPACK
module load BLAS
gfortran lapack_prb.f90 -o test.exe -llapack -lblas
用的是下面的代码进行测试(网上给的可能会有错). BLAS对下面代码也是需要的. 否则报错如dpbtrf.f:(.text+0x78e): undefined reference to `dtrsm_`.
program main
!*****************************************************************************80
!
!! MAIN is the main program for LAPACK_PRB.
!
! Discussion:
!
! LAPACK_PRB tests the LAPACK library.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 19 January 2011
!
! Author:
!
! John Burkardt
!
implicit none
call timestamp ( )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'LAPACK_PRB'
write ( *, '(a)' ) ' FORTRAN90 version'
write ( *, '(a)' ) ' Test the LAPACK library.'
call test01 ( )
call test02 ( )
call test03 ( )
call test04 ( )
call test045 ( )
call test05 ( )
call test06 ( )
call test07 ( )
call test08 ( )
call test09 ( )
call test10 ( )
call test11 ( )
call test12 ( )
call test13 ( )
call test14 ( )
call test15 ( )
call test16 ( )
call test17 ( )
!
! Terminate.
!
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'LAPACK_PRB'
write ( *, '(a)' ) ' Normal end of execution.'
write ( *, '(a)' ) ' '
call timestamp ( )
stop
end
subroutine test01 ( )
!*****************************************************************************80
!
!! TEST01 tests DGBTRF and DGBTRS.
!
! Discussion:
!
! 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)
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 25
integer ( kind = 4 ), parameter :: ml = 1
integer ( kind = 4 ), parameter :: mu = 1
integer ( kind = 4 ), parameter :: lda = 2 * ml + mu + 1
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) b(n)
integer ( kind = 4 ) info
integer ( kind = 4 ) ipiv(n)
integer ( kind = 4 ) m
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST01'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general band storage mode (GB):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGBTRF factors a general band matrix.'
write ( *, '(a)' ) ' DGBTRS solves a factored system.'
write ( *, '(a)' ) ' '
!
! Assign values to matrix A and right hand side b.
!
b(1) = 1.0D+00
b(2:n-1) = 0.0D+00
b(n) = 1.0D+00
!
! Zero out the matrix.
!
a(1:lda,1:n) = 0.0D+00
m = ml + mu + 1
!
! Superdiagonal,
! Diagonal,
! Subdiagonal.
!
a(m-1,2:n) = -1.0D+00
a(m,1:n) = 2.0D+00
a(m+1,1:n-1) = -1.0D+00
!
! Factor the matrix.
!
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Bandwidth is ', m
write ( *, '(a)' ) ' '
call dgbtrf ( n, n, ml, mu, a, lda, ipiv, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST01'
write ( *, '(a,i8)' ) ' Factorization failed, INFO = ', info
return
end if
!
! Solve the linear system.
!
call dgbtrs ( 'n', n, ml, mu, 1, a, lda, ipiv, b, n, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST01'
write ( *, '(a,i8)' ) ' Solution failed, INFO = ', info
return
end if
call r8vec_print_some ( n, b, 1, 5, ' Partial solution (all should be 1)' )
return
end
subroutine test02 ( )
!*****************************************************************************80
!
!! TEST02 tests DGBTRF and DGBTRS.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 25
integer ( kind = 4 ), parameter :: ml = 3
integer ( kind = 4 ), parameter :: mu = 3
integer ( kind = 4 ), parameter :: lda = 2 * ml + mu + 1
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) b(n)
integer ( kind = 4 ) i
integer ( kind = 4 ) ihi
integer ( kind = 4 ) ilo
integer ( kind = 4 ) info
integer ( kind = 4 ) ipiv(n)
integer ( kind = 4 ) j
integer ( kind = 4 ) m
real ( kind = 8 ) temp
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST02'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general band storage mode (GB):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGBTRF factors a general band matrix.'
write ( *, '(a)' ) ' DGBTRS 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.0D+00
do i = ilo, ihi
a(i-j+m,j) = -1.0D+00
temp = temp - 1.0D+00
end do
temp = temp + 1.0D+00
a(m,j) = 4.0D+00 - temp
end do
!
! Right hand side.
!
b(1:n) = 4.0D+00
!
! Factor the matrix.
!
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Bandwidth is ', m
call dgbtrf ( n, n, ml, mu, a, lda, ipiv, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Factorization failed, INFO = ', info
return
end if
!
! Solve the linear system.
!
call dgbtrs ( 'n', n, ml, mu, 1, a, lda, ipiv, b, n, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Solution failed, INFO = ', info
end if
call r8vec_print_some ( n, b, 1, 5, ' Partial solution (all should be 1)' )
return
end
subroutine test03 ( )
!*****************************************************************************80
!
!! TEST03 tests DGECON and DGETRF.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 3
integer ( kind = 4 ), parameter :: lda = n
integer ( kind = 4 ), parameter :: lwork = 4 * n
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) anorm
integer ( kind = 4 ) info
integer ( kind = 4 ) ipiv(n)
integer ( kind = 4 ) iwork(n)
real ( kind = 8 ) rcond
real ( kind = 8 ) r8mat_norm_li
real ( kind = 8 ) work(lwork)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST03'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGETRF computes the LU factorization;'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGECON computes the condition number '
write ( *, '(a)' ) ' of a factored matrix'
!
! Set the matrix.
!
a(1,1) = 1.0D+00
a(1,2) = 2.0D+00
a(1,3) = 3.0D+00
a(2,1) = 4.0D+00
a(2,2) = 5.0D+00
a(2,3) = 6.0D+00
a(3,1) = 7.0D+00
a(3,2) = 8.0D+00
a(3,3) = 0.0D+00
call r8mat_print ( n, n, a, ' The matrix A:' )
!
! Compute the L-Infinity norm of the matrix.
!
anorm = r8mat_norm_li ( n, n, a )
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' Matrix L-infinity norm is ', anorm
!
! Factor the matrix.
!
call dgetrf ( n, n, a, lda, ipiv, info )
!
! Get the condition number.
!
call dgecon ( 'I', n, a, lda, anorm, rcond, work, iwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Condition number calculation failed!'
write ( *, '(a,i8)' ) ' INFO = ', info
return
end if
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' Matrix reciprocal condition number = ', rcond
return
end
subroutine test04 ( )
!*****************************************************************************80
!
!! TEST04 tests DGEQRF and DORGQR.
!
! Discussion:
!
! DGEQRF computes the QR factorization of an M by N matrix A:
!
! A(MxN) = Q(MxK) * R(KxN)
!
! where K = min ( M, N ).
!
! DORGQR computes the explicit form of the Q factor.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 27 June 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: m = 8
integer ( kind = 4 ), parameter :: n = 6
integer ( kind = 4 ), parameter :: k = min ( m, n )
integer ( kind = 4 ), parameter :: lwork = n
real ( kind = 8 ) a(m,n)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) lda
real ( kind = 8 ) q(m,k)
real ( kind = 8 ) r(k,n)
integer ( kind = 4 ) seed
real ( kind = 8 ) tau(k)
real ( kind = 8 ) work(lwork)
seed = 123456789
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST04'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRF computes the QR factorization:'
write ( *, '(a)' ) ' A = Q * R'
write ( *, '(a)' ) ' DORGQR computes the explicit form of the Q factor.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' In this case, our M x N matrix A has more rows'
write ( *, '(a)' ) ' than columns:'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' M = ', m
write ( *, '(a,i8)' ) ' N = ', n
!
! Set A.
!
call r8mat_uniform_01 ( m, n, seed, a )
call r8mat_print ( m, n, a, ' The matrix A:' )
!
! Compute the QR factorization.
!
lda = m
call dgeqrf ( m, n, a, lda, tau, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DGEQRF returned nonzero INFO = ', info
return
end if
r(1:k,1:n) = 0.0D+00
do i = 1, k
r(i,i:n) = a(i,i:n)
end do
!
! Construct Q explicitly.
!
q(1:m,1:k) = a(1:m,1:k)
call dorgqr ( m, k, k, q, lda, tau, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DORGQR returned nonzero INFO = ', info
return
end if
call r8mat_print ( m, k, q, ' The Q factor:' )
call r8mat_print ( k, n, r, ' The R factor:' )
a(1:m,1:n) = matmul ( q(1:m,1:k), r(1:k,1:n) )
call r8mat_print ( m, n, a, ' The product Q * R:' )
return
end
subroutine test045 ( )
!*****************************************************************************80
!
!! TEST045 tests DGEQRF and DORMGQR.
!
! Discussion:
!
! We want to solve the MxN linear system A*x=b using the QR approach:
!
! Factor A:
!
! A = Q * R (step 1)
!
! Transform:
!
! A * x = b
! ==> Q * R * x = b
! ==> R * x = Q' * b (step 2)
! ==> x = inv(R) * Q' * b. (step 3)
!
! Step 1) DGEQRF computes the QR factorization of an M by N matrix A:
! A(MxN) = Q(MxK) * R(KxN) where K = min ( M, N ).
!
! Step 2) DORMQR can multiply Q' * b, putting the result back into b.
!
! Step 3) We could call a LAPACK routine to solve the upper
! triangular
! system R * x = Q' * b. Instead, we will try this part ourselves.
!
!
! LAPACK makes this process tricky because of two things it does
! for efficiency:
!
! *) LAPACK computes the Q and R factors in a
! compressed and encoded form, overwriting the matrix A and
! storing some extra information in a vector called TAU.
!
! *) LAPACK defines K = min ( M, N ), and
! does NOT compute the QR factorization as an MxM Q
! times an MxN R. Instead, it computes an MxK Q times
! a KxN R. This saves it worrying about zeroes, but it
! means the programmer has to worry about proper storage
! and correct dimensioning.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 19 January 2011
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: m = 8
integer ( kind = 4 ), parameter :: n = 6
integer ( kind = 4 ), parameter :: k = min ( m, n )
integer ( kind = 4 ), parameter :: lwork = n
real ( kind = 8 ) a(m,n)
real ( kind = 8 ) b(m)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) j
integer ( kind = 4 ) lda
real ( kind = 8 ) q(m,k)
real ( kind = 8 ) r(k,n)
integer ( kind = 4 ) seed
real ( kind = 8 ) tau(k)
real ( kind = 8 ) work(lwork)
real ( kind = 8 ) x(n)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST045'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRF computes the QR factorization:'
write ( *, '(a)' ) ' A = Q * R'
write ( *, '(a)' ) ' DORMQR can compute Q'' * b.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' We use these routines to carry out a QR'
write ( *, '(a)' ) ' solve of an M by N linear system A * x = b.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' In this case, our M x N matrix A has more rows'
write ( *, '(a)' ) ' than columns:'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' M = ', m
write ( *, '(a,i8)' ) ' N = ', n
!
! STEP 0: Set random A, simple x, and compute b.
!
seed = 123456789
call r8mat_uniform_01 ( m, n, seed, a )
do i = 1, n
x(i) = real ( i, kind = 8 )
end do
b(1:m) = matmul ( a(1:m,1:n), x(1:n) )
!
! Wipe out X so we believe it when we get it back...
!
x(1:n) = 0.0D+00
call r8mat_print ( m, n, a, ' The matrix A:' )
!
! STEP 1: Compute the QR factorization.
!
lda = m
call dgeqrf ( m, n, a, lda, tau, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DGEQRF returned nonzero INFO = ', info
return
end if
!
! Extract and save the K by N matrix R.
!
r(1:k,1:n) = 0.0D+00
do i = 1, k
r(i,i:n) = a(i,i:n)
end do
!
! STEP 2: Multiply Q' * b to get the N by 1 result in b.
!
call dormqr ( 'L', 'T', m, 1, k, a, m, tau, b, m, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DORMQR returned nonzero INFO = ', info
return
end if
!
! STEP 3: Compute inv(R) * Q' * b, or, equivalently,
! solve R * x = Q' * b.
!
do j = n, 1, -1
x(j) = b(j) / r(j,j)
do i = 1, n - 1
b(i) = b(i) - r(i,j) * x(j)
end do
end do
call r8vec_print ( n, x, ' The solution X:' )
return
end
subroutine test05 ( )
!*****************************************************************************80
!
!! TEST05 tests DGEQRF and DORGQR.
!
! Discussion:
!
! DGEQRF computes the QR factorization of an M by N matrix A:
!
! A(MxN) = Q(MxK) * R(KxN)
!
! where K = min ( M, N ).
!
! In order to get an M x M matrix Q, we are going to pad our
! original matrix A with extra columns of zeroes!
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: m = 8
integer ( kind = 4 ), parameter :: n = 6
integer ( kind = 4 ), parameter :: n2 = 8
integer ( kind = 4 ), parameter :: k = min ( m, n2 )
integer ( kind = 4 ), parameter :: lwork = n2
real ( kind = 8 ) a(m,n)
real ( kind = 8 ) a2(m,n2)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) lda
real ( kind = 8 ) q(m,k)
real ( kind = 8 ) r(k,n2)
integer ( kind = 4 ) seed
real ( kind = 8 ) tau(k)
real ( kind = 8 ) work(lwork)
seed = 123456789
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST05'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRF computes the QR factorization:'
write ( *, '(a)' ) ' A = Q * R'
write ( *, '(a)' ) ' DORGQR computes the explicit form of the Q factor.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' In this case, our M x N matrix A has more rows'
write ( *, '(a)' ) ' than columns:'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' M = ', m
write ( *, '(a,i8)' ) ' N = ', n
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Normally, LAPACK will only return an M x min(M,N)'
write ( *, '(a)' ) ' portion of Q. When N < M, we lose information.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Here, we force the computation of the full Q by making'
write ( *, '(a)' ) ' a copy of A with N-M extra zero columns.'
write ( *, '(a)' ) ' '
!
! Set A.
!
call r8mat_uniform_01 ( m, n, seed, a )
call r8mat_print ( m, n, a, ' The matrix A:' )
!
! Copy A, and pad with extra columns.
!
a2(1:m,1:n) = a(1:m,1:n)
a2(1:m,n+1:m) = 0.0D+00
!
! Compute the QR factorization.
!
lda = m
call dgeqrf ( m, n2, a2, lda, tau, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DGEQRF returned nonzero INFO = ', info
return
end if
r(1:k,1:n2) = 0.0D+00
do i = 1, k
r(i,i:n2) = a2(i,i:n2)
end do
!
! Construct Q explicitly.
!
q(1:m,1:k) = a2(1:m,1:k)
call dorgqr ( m, k, k, q, lda, tau, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DORGQR returned nonzero INFO = ', info
return
end if
call r8mat_print ( m, k, q, ' The Q factor:' )
call r8mat_print ( k, n, r, ' The R factor:' )
a2(1:m,1:n2) = matmul ( q(1:m,1:k), r(1:k,1:n2) )
call r8mat_print ( m, n2, a2, ' The product Q * R:' )
return
end
subroutine test06 ( )
!*****************************************************************************80
!
!! TEST06 tests DGEQRF and DORGQR.
!
! Discussion:
!
! DGEQRF computes the QR factorization of an M by N matrix A:
!
! A(MxN) = Q(MxK) * R(KxN)
!
! where K = min ( M, N ).
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: m = 6
integer ( kind = 4 ), parameter :: n = 8
integer ( kind = 4 ), parameter :: k = min ( m, n )
integer ( kind = 4 ), parameter :: lwork = n
real ( kind = 8 ) a(m,n)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) lda
real ( kind = 8 ) q(m,k)
real ( kind = 8 ) r(k,n)
integer ( kind = 4 ) seed
real ( kind = 8 ) tau(k)
real ( kind = 8 ) work(lwork)
seed = 123456789
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST06'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRF computes the QR factorization:'
write ( *, '(a)' ) ' A = Q * R'
write ( *, '(a)' ) ' DORGQR computes the explicit form of the Q factor.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' In this case, our M x N matrix A has more columns'
write ( *, '(a)' ) ' than rows:'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' M = ', m
write ( *, '(a,i8)' ) ' N = ', n
!
! Set A.
!
call r8mat_uniform_01 ( m, n, seed, a )
call r8mat_print ( m, n, a, ' The matrix A:' )
!
! Compute the QR factorization.
!
lda = m
call dgeqrf ( m, n, a, lda, tau, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DGEQRF returned nonzero INFO = ', info
return
end if
r(1:k,1:n) = 0.0D+00
do i = 1, k
r(i,i:n) = a(i,i:n)
end do
!
! Construct Q explicitly.
!
q(1:m,1:k) = a(1:m,1:k)
call dorgqr ( m, k, k, q, lda, tau, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DORGQR returned nonzero INFO = ', info
return
end if
call r8mat_print ( m, k, q, ' The Q factor:' )
call r8mat_print ( k, n, r, ' The R factor:' )
a(1:m,1:n) = matmul ( q(1:m,1:k), r(1:k,1:n) )
call r8mat_print ( m, n, a, ' The product Q*R' )
return
end
subroutine test07 ( )
!*****************************************************************************80
!
!! TEST07 tests DGESVD.
!
! Discussion:
!
! DGESVD computes the singular value decomposition:
!
! A = U * S * V'
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 08 February 2007
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: m = 6
integer ( kind = 4 ), parameter :: n = 4
integer ( kind = 4 ), parameter :: lwork = 3*min(m,n) + max (max(m,n), 2*min(m,n) )
real ( kind = 8 ) a(m,n)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) lda
integer ( kind = 4 ) ldu
integer ( kind = 4 ) ldvt
character jobu
character jobvt
real ( kind = 8 ) s(min(m,n))
integer ( kind = 4 ) seed
real ( kind = 8 ) sigma(m,n)
real ( kind = 8 ) u(m,m)
real ( kind = 8 ) vt(n,n)
real ( kind = 8 ) work(lwork)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST07'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGESVD computes the singular value decomposition:'
write ( *, '(a)' ) ' A = U * S * V'''
!
! Set A.
!
seed = 123456789
call r8mat_uniform_01 ( m, n, seed, a )
call r8mat_print ( m, n, a, ' The matrix A:' )
!
! Compute the singular values and singular vectors.
!
jobu = 'A'
jobvt = 'A'
lda = m
ldu = m
ldvt = n
call dgesvd ( jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, &
lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DGESVD returned nonzero INFO = ', info
return
end if
call r8vec_print ( min ( m, n ), s, ' Singular values' )
call r8mat_print ( m, m, u, ' Left singular vectors U:' )
call r8mat_print ( n, n, vt, ' Right singular vectors V'':' )
sigma(1:m,1:n) = 0.0D+00
do i = 1, min ( m, n )
sigma(i,i) = s(i)
end do
a(1:m,1:n) = matmul ( u(1:m,1:m), matmul ( sigma(1:m,1:n), vt(1:n,1:n)) )
call r8mat_print ( m, n, a, ' The product U * S * V'':' )
return
end
subroutine test08 ( )
!*****************************************************************************80
!
!! TEST08 tests DGETRF and DGETRI.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 3
integer ( kind = 4 ), parameter :: lda = n
integer ( kind = 4 ), parameter :: lwork = n
real ( kind = 8 ) a(lda,n)
integer ( kind = 4 ) info
integer ( kind = 4 ) ipiv(n)
real ( kind = 8 ) work(lwork)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST08'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGETRF factors a general matrix;'
write ( *, '(a)' ) ' DGETRI computes the inverse.'
!
! Set the matrix.
!
a(1,1) = 1.0D+00
a(1,2) = 2.0D+00
a(1,3) = 3.0D+00
a(2,1) = 4.0D+00
a(2,2) = 5.0D+00
a(2,3) = 6.0D+00
a(3,1) = 7.0D+00
a(3,2) = 8.0D+00
a(3,3) = 0.0D+00
call r8mat_print ( n, n, a, ' The matrix A:' )
!
! Factor the matrix.
!
call dgetrf ( n, n, a, lda, ipiv, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DGETRF returned INFO = ', info
write ( *, '(a)' ) ' The matrix is numerically singular.'
return
end if
!
! Compute the inverse matrix.
!
call dgetri ( n, a, lda, ipiv, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' The inversion procedure failed!'
write ( *, '(a,i8)' ) ' INFO = ', info
return
end if
call r8mat_print ( n, n, a, ' The inverse matrix:' )
return
end
subroutine test09 ( )
!*****************************************************************************80
!
!! TEST09 tests DGETRF and DGETRS.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 3
integer ( kind = 4 ), parameter :: lda = n
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) b(n)
integer ( kind = 4 ) info
integer ( kind = 4 ) ipiv(n)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST09'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGETRF computes the LU factorization;'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGETRS solves linear systems using the LU factors;'
!
! Set the matrix.
!
a(1,1) = 1.0D+00
a(1,2) = 2.0D+00
a(1,3) = 3.0D+00
a(2,1) = 4.0D+00
a(2,2) = 5.0D+00
a(2,3) = 6.0D+00
a(3,1) = 7.0D+00
a(3,2) = 8.0D+00
a(3,3) = 0.0D+00
call r8mat_print ( n, n, a, ' The matrix A:' )
!
! Factor the matrix.
!
call dgetrf ( n, n, a, lda, ipiv, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DGETRF returned INFO = ', info
write ( *, '(a)' ) ' The matrix is numerically singular.'
return
end if
!
! Set the right hand side.
!
b(1:3) = (/ 14.0D+00, 32.0D+00, 23.0D+00 /)
call r8vec_print ( n, b, ' Right hand side B' )
!
! Solve the linear system.
!
call dgetrs ( 'n', n, 1, a, lda, ipiv, b, n, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Solution procedure failed!'
write ( *, '(a,i8)' ) ' INFO = ', info
return
end if
call r8vec_print ( n, b, ' The solution X' )
return
end
subroutine test10 ( )
!*****************************************************************************80
!
!! TEST10 tests DGETRF and DGETRS.
!
! Discussion:
!
! 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.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 25
integer ( kind = 4 ), parameter :: lda = n
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) b(n)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) ipiv(n)
integer ( kind = 4 ) j
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST10'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general storage mode (GE):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGETRF factors a general matrix;'
write ( *, '(a)' ) ' DGETRS 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) = dble ( n )
else
a(i,j) = -1.0D+00
end if
end do
end do
b(1:n) = 1.0D+00
!
! Factor the matrix.
!
call dgetrf ( n, n, a, lda, ipiv, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Matrix is singular, INFO = ', info
return
end if
!
! Solve the linear system.
!
call dgetrs ( 'n', n, 1, a, lda, ipiv, b, n, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Solution procedure failed, INFO = ', info
return
end if
call r8vec_print_some ( n, b, 1, 5, ' Partial solution (all should be 1)' )
return
end
subroutine test11 ( )
!*****************************************************************************80
!
!! TEST11 tests DGTSV.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 100
real ( kind = 8 ) b(n)
real ( kind = 8 ) c(n-1)
real ( kind = 8 ) d(n)
real ( kind = 8 ) e(n-1)
integer ( kind = 4 ) info
integer ( kind = 4 ) ldb
integer ( kind = 4 ) nrhs
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST11'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in general tridiagonal storage mode (GT):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGTSV factors and solves a linear system'
write ( *, '(a)' ) ' with a general tridiagonal matrix.'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' The system is of order N = ', n
write ( *, '(a)' ) ' '
!
! Right hand side.
!
b(1:n-1) = 0.0D+00
b(n) = n + 1
!
! Subdiagonal.
! Diagonal.
! Superdiagonal.
!
c(1:n-1) = -1.0D+00
d(1:n) = 2.0D+00
e(1:n-1) = -1.0D+00
nrhs = 1
ldb = n
!
! Factor and solve the linear system.
!
call dgtsv ( n, nrhs, c, d, e, b, ldb, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Solution procedure failed.'
write ( *, '(a,i8)' ) ' INFO = ', info
return
end if
call r8vec_print_some ( n, b, 1, 5, ' Partial solution (Should be 1,2,3...)' )
return
end
subroutine test12 ( )
!*****************************************************************************80
!
!! TEST12 tests DPBTRF.
!
! Discussion:
!
! We want to compute the lower triangular Cholesky factor L
! of a positive definite symmetric band matrix.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 5
integer ( kind = 4 ), parameter :: nband = 1
real ( kind = 8 ) a(nband+1,n)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) j
real ( kind = 8 ) l(nband+1,n)
real ( kind = 8 ) l_row(n)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST12'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in positive definite band storage mode (PB):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DPBTRF computes'
write ( *, '(a)' ) ' the lower Cholesky factor A = L*L'' or'
write ( *, '(a)' ) ' the upper Cholesky factor A = U''*U;'
write ( *, '(a)' ) ' '
!
! Zero out the matrix.
!
a(1:nband+1,1:n) = 0.0D+00
!
! Store the diagonal of a symmetric band matrix.
!
a(1,1:n) = 2.0D+00
!
! Store the subdiagonal of a symmetric band matrix.
!
a(2,1:n-1) = -1.0D+00
!
! Get the lower triangular Cholesky factor L:
!
l(1:nband+1,1:n) = a(1:nband+1,1:n)
call dpbtrf ( 'L', n, nband, l, nband+1, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Factorization failed, INFO = ', info
return
end if
!
! Print the relevant entries of L:
!
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' The lower Cholesky factor L:'
write ( *, '(a)' ) ' '
do i = 1, n
do j = 1, n
if ( 0 <= i - j .and. i-j <= nband ) then
l_row(j) = l(i-j+1,j)
else
l_row(j) = 0.0D+00
end if
end do
write ( *, '(5f10.6)' ) l_row(1:n)
end do
return
end
subroutine test13 ( )
!*****************************************************************************80
!
!! TEST13 tests DPBTRF and DPBTRS.
!
! Discussion:
!
! 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)
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 25
integer ( kind = 4 ), parameter :: nband = 1
integer ( kind = 4 ), parameter :: lda = nband + 1
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) b(n)
integer ( kind = 4 ) info
integer ( kind = 4 ) nrhs
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST13'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in positive definite band storage mode (PB):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' For a positive definite symmetric band matrix:'
write ( *, '(a)' ) ' DPBTRF factors;'
write ( *, '(a)' ) ' DPBTRS solves linear systems.'
write ( *, '(a)' ) ' '
!
! Zero out the matrix.
!
a(1:lda,1:n) = 0.0D+00
!
! Super (and sub) diagonal.
!
a(1,2:n) = -1.0D+00
!
! Diagonal.
!
a(2,1:n) = 2.0D+00
!
! Set the right hand side.
!
b(1) = 1.0D+00
b(2:n-1) = 0.0D+00
b(n) = 1.0D+00
!
! Factor the matrix.
!
call dpbtrf ( 'u', n, nband, a, lda, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Factorization failed, INFO = ', info
return
end if
!
! Solve the linear system.
!
nrhs = 1
call dpbtrs ( 'u', n, nband, nrhs, a, lda, b, n, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Solution failed, INFO = ', info
end if
call r8vec_print_some ( n, b, 1, 5, ' Partial solution (all should be 1)' )
return
end
subroutine test14 ( )
!*****************************************************************************80
!
!! TEST14 tests DPOTRF and DPOTRI.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 5
integer ( kind = 4 ), parameter :: lda = n
real ( kind = 8 ) a(n,n)
real ( kind = 8 ) a_inv(n,n)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
real ( kind = 8 ) r(n,n)
real ( kind = 8 ) temp(n,n)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST14'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in positive definite storage mode (PO):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DPOTRF computes the Cholesky factor A = R''*R;'
write ( *, '(a)' ) ' DPOTRI computes the inverse.'
write ( *, '(a)' ) ' '
!
! Zero out the matrix.
!
a(1:n,1:n) = 0.0D+00
!
! Subdiagonal.
!
do i = 2, n
a(i,i-1) = -1.0D+00
end do
!
! Diagonal.
!
do i = 1, n
a(i,i) = 2.0D+00
end do
!
! Superdiagonal.
!
do i = 1, n - 1
a(i,i+1) = -1.0D+00
end do
call r8mat_print ( n, n, a, ' The matrix A:' )
!
! Factor the matrix.
!
r(1:n,1:n) = a(1:n,1:n)
call dpotrf ( 'u', n, r, lda, info )
if ( info /= 0 ) then
write ( *, '(a,i8)' ) ' DPOTRF returns INFO = ', info
return
end if
do i = 1, n
r(i,1:i-1) = 0.0D+00
end do
call r8mat_print ( n, n, r, ' The Cholesky factor R:' )
temp(1:n,1:n) = matmul ( transpose ( r(1:n,1:n) ) , r(1:n,1:n) )
call r8mat_print ( n, n, temp, ' The product R'' * R' )
!
! Compute the inverse matrix.
!
a_inv(1:n,1:n) = r(1:n,1:n)
call dpotri ( 'u', n, a_inv, lda, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' The inversion procedure failed, INFO = ', info
return
end if
do i = 1, n
a_inv(i,1:i-1) = a_inv(1:i-1,i)
end do
call r8mat_print ( n, n, a_inv, ' The inverse matrix B:' )
temp(1:n,1:n) = matmul ( transpose ( a_inv(1:n,1:n) ) , a(1:n,1:n) )
call r8mat_print ( n, n, temp, ' The product B * A' )
return
end
subroutine test15 ( )
!*****************************************************************************80
!
!! TEST15 tests DSBGVX.
!
! Discussion:
!
! DSBGVX deals with the generalized eigenvalue problem:
!
! A * x = lambda * B * x
!
! where A and B are symmetric and banded (and stored in LAPACK
! symmetric
! band storage mode). B is additionally assumed to be positive
! definite.
!
! This is an "expert" interface, and the user is requesting
! only some of the eigenvalues and eigenvectors. In this example,
! only the largest and smallest (in magnitude) eigenvalues will
! be requested.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 17 March 2003
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! real ( kind = 8 ) AB(LDAB,N), contains, on input, the upper or
! lower
! triangle of the symmetric band matrix A, stored in the first KA+1
! rows
! of the array AB.
! If UPLO = 'U', then
! AB(KA+1+I-J,J) = A(I,J) for max(1,J-KA) <= I <= J;
! If UPLO = 'L', then
! AB(1+I-J,J) = A(I,J) for J <= I <= min(N,J+KA).
!
! real ( kind = 8 ) ABSTOL, the absolute error tolerance for the
! eigenvalues.
! If the input value of ABSTOL is not positive, then an appropriate
! value will be determined internally and used instead.
!
! real ( kind = 8 ) BB(LDBB,N), contains, on input, the upper or
! lower
! triangle of the positive definite symmetric band matrix B, stored
! in
! the first KB+1 rows of the array BB.
! If UPLO = 'U', then
! BB(KB+1+I-J,J) = B(I,J) for max(1,J-KB) <= I <= J;
! If UPLO = 'L', then
! BB(1+I-J,J) = B(I,J) for J <= I <= min(N,J+KB).
!
! integer ( kind = 4 ) IFAIL(N), if JOBZ = 'V', then if INFO is 0,
! the first
! M elements of IFAIL have been set to zero by DSBGVX, but if INFO
! is nonzero, IFAIL contains the indices of the eigenvalues
! for which the eigenvectors failed to converge. If JOBZ = 'N',
! then IFAIL is not referenced.
!
! integer ( kind = 4 ) IL, IU, the indices of the first (smallest)
! and last
! (largest) eigenvalues to be returned. These values are only
! used if RANGE = 'I'. It must be the case that 1 <= IL <= IU <= N.
!
! Integer INFO, is 0 for a successful computation,
! negative if an input argument was illegal (the index of that
! argument is the value of -INFO), or positive, in which case,
! if 0 < INFO <= N, then INFO off diagonal elements of an
! intermediate tridiagonal form did not converge to zero, or
! if N < INFO, B is not positive definite and the computation
! could not be completed.
!
! integer ( kind = 4 ) IWORK(5*N), workspace.
!
! character JOBZ, is 'N' if only eigenvalues are desired, or 'V'
! if eigenvectors will also be required.
!
! Integer KA, the number of superdiagonals (if UPLO = 'U') or
! subdiagonals (if UPLO = 'L') of A that are nonzero.
!
! integer ( kind = 4 ) KB, the number of superdiagonals (if UPLO =
! 'U') or
! subdiagonals (if UPLO = 'L') of B that are nonzero.
!
! integer ( kind = 4 ) LDAB, the leading dimension of the array AB,
! which
! must be at least KA+1.
!
! integer ( kind = 4 ) LDBB, the leading dimension of the array BB,
! which
! must be at least KB+1.
!
! integer ( kind = 4 ) LDQ, the leading dimension of the array Q.
! If JOBZ = 'N', then Q is not used, and LDQ should be any
! positive value such as 1. If JOBZ = 'V', then LDQ must be
! at least N.
!
! integer ( kind = 4 ) LDZ, the leading dimension of the array Z.
! If JOBZ = 'N', then Z is not used, and LDZ should be any
! positive value such as 1. If JOBZ = 'V', then LDZ must be
! at least N.
!
! integer ( kind = 4 ) M, the number of eigenvalues found by DSBGVX.
!
! integer ( kind = 4 ) N, the order of the matrices A and B.
!
! real ( kind = 8 ) Q(LDQ,N), if JOBZ = 'V', the N by N matrix used
! to
! reduce the problem to standard form: "C * x = lambda * x"
! and then to reduce the matrix C to tridiagonal form. But
! if JOBZ is not 'V', Q is not referenced.
!
! character RANGE, specifies which eigenvalues are desired.
! 'A' means all, 'V' means a real interval will be specified in which
! eigenvalues are to be sought, 'I' means a range of indices will
! be specified.
!
! character UPLO, is 'U' if the upper triangles of A and B are
! stored,
! 'L' if the lower triangles are stored.
!
! real ( kind = 8 ) VL, VU, the lower and upper bounds of an interval
! to be
! searched for eigenvalues. In this case, VL must be less than VU.
! These values are used only if RANGE = 'V'.
!
! real ( kind = 8 ) W(N), the requested eigenvalues, in ascending
! order.
!
! real ( kind = 8 ) WORK(7*N), workspace.
!
! real ( kind = 8 ) Z(LDZ,N), if JOBZ = 'V', the I-th column of Z
! contains
! the eigenvector associated with the I-th eigenvalue W(I).
!
implicit none
integer ( kind = 4 ), parameter :: n = 11
integer ( kind = 4 ), parameter :: ka = 2
integer ( kind = 4 ), parameter :: kb = 1
integer ( kind = 4 ), parameter :: ldab = ka+1
integer ( kind = 4 ), parameter :: ldbb = kb+1
integer ( kind = 4 ), parameter :: ldq = 1
integer ( kind = 4 ), parameter :: ldz = 1
real ( kind = 8 ) ab(ldab,n)
real ( kind = 8 ), parameter :: abstol = 0.0D+00
real ( kind = 8 ) bb(ldbb,n)
integer ( kind = 4 ) i
integer ( kind = 4 ) ifail(n)
integer ( kind = 4 ) il
integer ( kind = 4 ) ilo
integer ( kind = 4 ) info
integer ( kind = 4 ) iu
integer ( kind = 4 ) iwork(5*n)
integer ( kind = 4 ) j
character :: jobz = 'N'
integer ( kind = 4 ) m
real ( kind = 8 ) q(ldq,n)
character :: range = 'I'
integer ( kind = 4 ) test
character :: uplo = 'U'
real ( kind = 8 ) value
real ( kind = 8 ), parameter :: vl = 0.0D+00
real ( kind = 8 ), parameter :: vu = 1.0D+00
real ( kind = 8 ) w(n)
real ( kind = 8 ) work(7*n)
real ( kind = 8 ) z(ldz,n)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST15'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in symmetric band storage mode (SB):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' For a symmetric banded NxN matrix A, and a symmetric'
write ( *, '(a)' ) ' banded positive definite NxN matrix B,'
write ( *, '(a)' ) ' DSBGVX solves the generalized eigenvalue problem'
write ( *, '(a)' ) ' A * X = LAMBDA * B * X'
write ( *, '(a)' ) ' '
do test = 1, 2
!
! Set A.
!
do j = 1, n
ilo = max ( j - ka, 1 )
do i = ilo, j
if ( j == i-2 ) then
value = -1.0D+00
else if ( j == i-1 ) then
value = -1.0D+00
else if ( j == i ) then
value = +4.0D+00
else if ( j == i+1 ) then
value = -1.0D+00
else if ( j == i+2 ) then
value = -1.0D+00
else
value = 0.0D+00
end if
ab(ka+1+i-j,j) = value
end do
end do
!
! Set B.
!
do j = 1, n
ilo = max ( j - kb, 1 )
do i = ilo, j
if ( j == i-1 ) then
value = -1.0D+00
else if ( j == i ) then
value = +2.0D+00
else if ( j == i+1 ) then
value = -1.0D+00
else
value = 0.0D+00
end if
bb(kb+1+i-j,j) = value
end do
end do
!
! Request the value of the SMALLEST or LARGEST eigenvalue:
!
if ( test == 1 ) then
il = 1
iu = 1
else if ( test == 2 ) then
il = n
iu = n
end if
call dsbgvx ( jobz, range, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, &
ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info )
if ( info < 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Illegal value for input argument ', -info
return
else if ( 0 < info ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' The eigenvalue or eigenvector iterations'
write ( *, '(a)' ) ' did not converge.'
cycle
end if
call r8vec_print ( m, w, ' Computed eigenvalues' )
end do
return
end
subroutine test16 ( )
!*****************************************************************************80
!
!! TEST16 tests DSYEV.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2006
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: n = 7
integer ( kind = 4 ), parameter :: lwork = 3 * n - 1
real ( kind = 8 ) a(n,n)
integer ( kind = 4 ) info
character jobz
real ( kind = 8 ) lambda(n)
character uplo
real ( kind = 8 ) work(lwork)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST16'
write ( *, '(a)' ) ' For a double precision real matrix (D)'
write ( *, '(a)' ) ' in symmetric storage mode (SY):'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' For a symmetric matrix in general storage,'
write ( *, '(a)' ) ' DSYEV computes eigenvalues and eigenvectors;'
write ( *, '(a)' ) ' '
!
! Set A.
!
call clement2 ( n, a )
call r8mat_print ( n, n, a, ' The matrix A:' )
!
! Compute the eigenvalues and eigenvectors.
!
jobz = 'V'
uplo = 'U'
call dsyev ( jobz, uplo, n, a, n, lambda, work, lwork, info )
if ( info /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' DSYEV returned nonzero INFO = ', info
else
call r8vec_print ( n, lambda, ' The eigenvalues:' )
if ( jobz == 'V' ) then
call r8mat_print ( n, n, a, ' The eigenvector matrix:' )
end if
end if
return
end
subroutine test17 ( )
!*********************************************************************72
!
!! TEST17 tests DGEQRF.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 04 March 2010
!
! Author:
!
! John Burkardt
!
implicit none
integer ( kind = 4 ), parameter :: m = 5
integer ( kind = 4 ), parameter :: n = 5
integer ( kind = 4 ), parameter :: nrhs = 1
integer ( kind = 4 ), parameter :: lda = m
integer ( kind = 4 ), parameter :: ldb = max ( m, n )
integer ( kind = 4 ), parameter :: lwork = n
integer ( kind = 4 ), parameter :: tau_size = n
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) b(ldb,nrhs)
integer ( kind = 4 ) i
integer ( kind = 4 ) info
integer ( kind = 4 ) j
real ( kind = 8 ) tau(tau_size)
real ( kind = 8 ) work(lwork)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST17'
write ( *, '(a)' ) ' Call the standard LAPACK routine'
write ( *, '(a)' ) ' DGEQRF to get the QR factorization of'
write ( *, '(a)' ) ' a matrix stored in GE format.'
do j = 1, n
do i = 1, m
if ( j == i - 1 ) then
a(i,j) = - 1.0D+00
else if ( j == i ) then
a(i,j) = 2.0D+00
else if ( j == i + 1 ) then
a(i,j) = -1.0D+00
else
a(i,j) = 0.0D+00
end if
end do
end do
call r8mat_print ( m, n, a, ' Input matrix:' )
call dgeqrf ( m, n, a, lda, tau, work, lwork, info )
if ( info == 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRF called successfully.'
else
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRF returned error flag.'
write ( *, '(a,i8)' ) ' INFO = ', info
return
end if
call r8mat_print ( m, n, a, ' Factored matrix:' )
call r8vec_print ( tau_size, tau, ' Tau:' );
!
! Set up and solve a linear system using DQEQRS.
!
do i = 1, n
b(i,1) = 0.0D+00
end do
b(n,1) = dble ( n + 1 )
call dgeqrs ( m, n, nrhs, a, lda, tau, b, ldb, work, lwork, info )
if ( info == 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRS called successfully.'
else
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' DGEQRS returned error flag.'
write ( *, '(a,i8)' ) ' INFO = ', info
return
end if
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Solution from DGEQRS:'
write ( *, '(a)' ) ' '
do i = 1, m
write ( *, '(2x,g14.6)' ) b(i,1)
end do
return
end
subroutine clement2 ( n, a )
!*****************************************************************************80
!
!! CLEMENT2 returns the Clement2 matrix.
!
! Formula:
!
! if ( J = I+1 )
! A(I,J) = sqrt(I*(N-I))
! else if ( I = J+1 )
! A(I,J) = sqrt(J*(N-J))
! else
! A(I,J) = 0
!
! Example:
!
! N = 5
!
! . sqrt(4) . . .
! sqrt(4) . sqrt(6) . .
! . sqrt(6) . sqrt(6) .
! . . sqrt(6) . sqrt(4)
! . . . sqrt(4) .
!
! Properties:
!
! A is tridiagonal.
!
! Because A is tridiagonal, it has property A (bipartite).
!
! A is symmetric: A' = A.
!
! Because A is symmetric, it is normal.
!
! Because A is normal, it is diagonalizable.
!
! A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
!
! The diagonal of A is zero.
!
! A is singular if N is odd.
!
! About 64 percent of the entries of the inverse of A are zero.
!
! The eigenvalues are plus and minus the numbers
!
! N-1, N-3, N-5, ..., (1 or 0).
!
! If N is even,
!
! det ( A ) = (-1)**(N/2) * (N-1) * (N+1)**(N/2)
!
! and if N is odd,
!
! det ( A ) = 0
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 15 April 1999
!
! Author:
!
! John Burkardt
!
! Reference:
!
! P A Clement,
! A class of triple-diagonal matrices for test purposes,
! SIAM Review,
! Volume 1, 1959, pages 50-52.
!
! Parameters:
!
! Input, integer N, the order of A.
!
! Output, real ( kind = 8 ) A(N,N), the matrix.
!
implicit none
integer ( kind = 4 ) n
real ( kind = 8 ) a(n,n)
integer ( kind = 4 ) i
integer ( kind = 4 ) j
do i = 1, n
do j = 1, n
if ( j == i + 1 ) then
a(i,j) = sqrt ( real ( i * ( n - i ), kind = 8 ) )
else if ( i == j + 1 ) then
a(i,j) = sqrt ( real ( j * ( n - j ), kind = 8 ) )
else
a(i,j) = 0.0D+00
end if
end do
end do
return
end
subroutine dgeqrs ( m, n, nrhs, a, lda, tau, b, ldb, work, lwork, info )
!*****************************************************************************80
!
!! DGEQRS solves a linear system factored by DGEQRF.
!
! Discussion:
!
! This routine solves the least squares problem
! min || A*X - B ||
! using the QR factorization of the M x N matrix A:
! A = Q*R
! computed by DGEQRF.
!
! Note that this routine requires that N <= M.
!
! Parameters:
!
! Input, integer ( kind = 4 ) M,
! The number of rows of the matrix A. M >= 0.
!
! Input, integer ( kind = 4 ) N,
! The number of columns of the matrix A. M >= N >= 0.
!
! Input, integer ( kind = 4 ) NRHS,
! The number of columns of B. NRHS >= 0.
!
! Input, real ( kind = 8 ) A(LDA,N),
! Details of the QR factorization of the original matrix A as
! returned by DGEQRF.
!
! Input, integer ( kind = 4 ) LDA,
! The leading dimension of the array A. LDA >= M.
!
! Input, real ( kind = 8 ) TAU(N),
! Details of the orthogonal matrix Q.
!
! Input/output, real ( kind = 8 ) B(LDA,NHRS),
! On entry, the m-by-nrhs right hand side matrix B.
! On exit, the n-by-nrhs solution matrix X.
!
! Input, integer ( kind = 4 ) LDB,
! The leading dimension of the array B. LDB >= M.
!
! Workspace, real ( kind = 8 ) WORK(LWORK).
!
! Input, integer ( kind = 4 ) LWORK,
! The length of the array WORK. LWORK must be at least NRHS,
! and should be at least NRHS*NB, where NB is the block size
! for this environment.
!
! Output, integer ( kind = 4 ) INFO,
! = 0: successful exit
! < 0: if INFO = -i, the i-th argument had an illegal value
!
implicit none
integer ( kind = 4 ) lda
integer ( kind = 4 ) ldb
integer ( kind = 4 ) lwork
integer ( kind = 4 ) n
integer ( kind = 4 ) nrhs
real ( kind = 8 ) a(lda,n)
real ( kind = 8 ) b(ldb,nrhs)
integer ( kind = 4 ) info
integer ( kind = 4 ) m
real ( kind = 8 ), parameter :: one = 1.0D+00
real ( kind = 8 ) tau(n)
real ( kind = 8 ) work(lwork)
info = 0
if ( m < 0 ) then
info = - 1
else if ( n < 0 .or. m < n ) then
info = - 2
else if ( nrhs < 0 ) then
info = - 3
else if ( lda < max( 1, m ) ) then
info = - 5
else if ( ldb < max( 1, m ) ) then
info = - 8
else if ( lwork < 1 .or. lwork < nrhs .and. 0 < m .and. 0 < n ) then
info = - 10
end if
if ( info /= 0 ) then
call xerbla ( 'DGEQRS', - info )
return
end if
!
! Quick return if possible.
!
if ( n == 0 .or. nrhs == 0 .or. m == 0 ) then
return
end if
!
! Compute B := Q' * B
!
call dormqr( 'Left', 'Transpose', m, nrhs, n, a, lda, tau, b, ldb, &
work, lwork, info )
!
! Solve R * X = B(1:n,:)
!
call dtrsm ( 'Left', 'Upper', 'No transpose', 'Non-unit', n, nrhs, &
one, a, lda, b, ldb )
return
end
function r8mat_norm_li ( m, n, a )
!*****************************************************************************80
!
!! R8MAT_NORM_LI returns the matrix L-infinity norm of an R8MAT.
!
! Discussion:
!
! An R8MAT is a two dimensional matrix of double precision real
! values.
!
! The matrix L-infinity norm is defined as:
!
! R8MAT_NORM_LI = max ( 1 <= I <= M ) sum ( 1 <= J <= N ) abs (
! A(I,J) ).
!
! The matrix L-infinity norm is derived from the vector L-infinity
! norm,
! and satisifies:
!
! r8vec_norm_li ( A * x ) <= r8mat_norm_li ( A ) * r8vec_norm_li (
! x ).
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 11 December 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer M, the number of rows in A.
!
! Input, integer N, the number of columns in A.
!
! Input, real ( kind = 8 ) A(M,N), the matrix whose L-infinity
! norm is desired.
!
! Output, real ( kind = 8 ) R8MAT_NORM_LI, the L-infinity norm of A.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) n
real ( kind = 8 ) a(m,n)
integer ( kind = 4 ) i
real ( kind = 8 ) r8mat_norm_li
r8mat_norm_li = 0.0D+00
do i = 1, m
r8mat_norm_li = max ( r8mat_norm_li, sum ( abs ( a(i,1:n) ) ) )
end do
return
end
subroutine r8mat_print ( m, n, a, title )
!*****************************************************************************80
!
!! R8MAT_PRINT prints an R8MAT.
!
! Discussion:
!
! An R8MAT is a two dimensional matrix of double precision real
! values.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 September 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer M, the number of rows in A.
!
! Input, integer N, the number of columns in A.
!
! Input, real ( kind = 8 ) A(M,N), the matrix.
!
! Input, character ( len = * ) TITLE, a title to be printed.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) n
real ( kind = 8 ) 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 an R8MAT.
!
! Discussion:
!
! An R8MAT is a two dimensional matrix of double precision real
! values.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 26 March 2005
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer M, N, the number of rows and columns.
!
! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed.
!
! Input, integer ILO, JLO, the first row and column to print.
!
! Input, integer IHI, JHI, the last row and column to print.
!
! Input, character ( len = * ) TITLE, an optional title.
!
implicit none
integer ( kind = 4 ), parameter :: incx = 5
integer ( kind = 4 ) m
integer ( kind = 4 ) n
real ( kind = 8 ) a(m,n)
character ( len = 14 ) ctemp(incx)
integer ( kind = 4 ) i
integer ( kind = 4 ) i2hi
integer ( kind = 4 ) i2lo
integer ( kind = 4 ) ihi
integer ( kind = 4 ) ilo
integer ( kind = 4 ) inc
integer ( kind = 4 ) j
integer ( kind = 4 ) j2
integer ( kind = 4 ) j2hi
integer ( kind = 4 ) j2lo
integer ( kind = 4 ) jhi
integer ( kind = 4 ) jlo
character ( len = * ) title
if ( 0 < len_trim ( title ) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) trim ( title )
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
if ( a(i,j) == real ( int ( a(i,j) ), kind = 8 ) ) then
write ( ctemp(j2), '(f8.0,6x)' ) a(i,j)
else
write ( ctemp(j2), '(g14.6)' ) a(i,j)
end if
end do
write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j), j = 1, inc )
end do
end do
write ( *, '(a)' ) ' '
return
end
subroutine r8mat_uniform_01 ( m, n, seed, r )
!*****************************************************************************80
!
!! R8MAT_UNIFORM_01 fills an R8MAT with unit pseudorandom numbers.
!
! Discussion:
!
! An R8MAT is a two dimensional matrix of double precision real
! values.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 11 August 2004
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Paul Bratley, Bennett Fox, Linus Schrage,
! A Guide to Simulation,
! Springer Verlag, pages 201-202, 1983.
!
! Bennett Fox,
! Algorithm 647:
! Implementation and Relative Efficiency of Quasirandom
! Sequence Generators,
! ACM Transactions on Mathematical Software,
! Volume 12, Number 4, pages 362-376, 1986.
!
! Peter Lewis, Allen Goodman, James Miller,
! A Pseudo-Random Number Generator for the System/360,
! IBM Systems Journal,
! Volume 8, pages 136-143, 1969.
!
! Parameters:
!
! Input, integer M, N, the number of rows and columns in the array.
!
! Input/output, integer SEED, the "seed" value, which should NOT be
! 0.
! On output, SEED has been updated.
!
! Output, real ( kind = 8 ) R(M,N), the array of pseudorandom values.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) n
integer ( kind = 4 ) i
integer ( kind = 4 ) j
integer ( kind = 4 ) k
integer ( kind = 4 ) seed
real ( kind = 8 ) r(m,n)
do j = 1, n
do i = 1, m
k = seed / 127773
seed = 16807 * ( seed - k * 127773 ) - k * 2836
if ( seed < 0 ) then
seed = seed + huge ( seed )
end if
r(i,j) = real ( seed, kind = 8 ) * 4.656612875D-10
end do
end do
return
end
subroutine r8vec_print ( n, a, title )
!*****************************************************************************80
!
!! R8VEC_PRINT prints an R8VEC.
!
! Discussion:
!
! An R8VEC is an array of double precision real values.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 22 August 2000
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of components of the vector.
!
! Input, real ( kind = 8 ) A(N), the vector to be printed.
!
! Input, character ( len = * ) TITLE, an optional title.
!
implicit none
integer ( kind = 4 ) n
real ( kind = 8 ) a(n)
integer ( kind = 4 ) i
character ( len = * ) title
if ( 0 < len_trim ( title ) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) trim ( title )
end if
write ( *, '(a)' ) ' '
do i = 1, n
write ( *, '(2x,i8,2x,g16.8)' ) i, a(i)
end do
return
end
subroutine r8vec_print_some ( n, a, i_lo, i_hi, title )
!*****************************************************************************80
!
!! R8VEC_PRINT_SOME prints "some" of an R8VEC.
!
! Discussion:
!
! An R8VEC is a vector of R8 values.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 16 October 2006
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of entries of the vector.
!
! Input, real ( kind = 8 ) A(N), the vector to be printed.
!
! Input, integer I_LO, I_HI, the first and last indices to print.
! The routine expects 1 <= I_LO <= I_HI <= N.
!
! Input, character ( len = * ) TITLE, an optional title.
!
implicit none
integer ( kind = 4 ) n
real ( kind = 8 ) a(n)
integer ( kind = 4 ) i
integer ( kind = 4 ) i_hi
integer ( kind = 4 ) i_lo
character ( len = * ) title
if ( 0 < len_trim ( title ) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) trim ( title )
end if
write ( *, '(a)' ) ' '
do i = max ( i_lo, 1 ), min ( i_hi, n )
write ( *, '(2x,i8,2x,g14.8)' ) i, a(i)
end do
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 GNU LGPL license.
!
! Modified:
!
! 18 May 2013
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! None
!
implicit none
character ( len = 8 ) ampm
integer ( kind = 4 ) d
integer ( kind = 4 ) h
integer ( kind = 4 ) m
integer ( kind = 4 ) mm
character ( len = 9 ), parameter, dimension(12) :: month = (/ &
'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', &
'September', 'October ', 'November ', 'December ' /)
integer ( kind = 4 ) n
integer ( kind = 4 ) s
integer ( kind = 4 ) values(8)
integer ( kind = 4 ) 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,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