program main c*********************************************************************72 c cc blas2_s_test() tests blas(). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 04 April 2014 c c Author: c c John Burkardt c implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'blas2_s_test():' write ( *, '(a)' ) ' FORTRAN77 version' write ( *, '(a)' ) ' Test blas2_s().' call test01 ( ) call test02 ( ) call test03 ( ) call test04 ( ) call test05 ( ) call test06 ( ) c c Terminate c write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'blas2_s_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine test01 ( ) c*********************************************************************72 c cc TEST01 tests SGEMV. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 February 2014 c c Author: c c John Burkardt c implicit none integer mn_max parameter ( mn_max = 5 ) real a(mn_max,mn_max) real alpha real beta integer i integer incx integer incy integer lda integer m integer n character trans real x(mn_max) real y(mn_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' For a general matrix A,' write ( *, '(a)' ) & ' SGEMV computes y := alpha * A * x + beta * y' write ( *, '(a)' ) & ' or y := alpha * A'' * x + beta * y.' c c y = alpha * A * x + beta * y c trans = 'N' m = 5 n = 4 alpha = 2.0E+00 lda = m call r4mat_test ( trans, lda, m, n, a ) do i = 1, n x(i) = real ( i ) end do incx = 1 beta = 3.0E+00 do i = 1, m y(i) = real ( 10 * i ) end do incy = 1 call r4mat_print ( m, n, a, ' Matrix A:' ) call r4vec_print ( n, x, ' Vector X:' ) call r4vec_print ( m, y, ' Vector Y:' ) call sgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy ) call r4vec_print ( m, y, & ' Result Y = alpha * A * x + beta * y' ) c c y = alpha * A' * x + beta * y c trans = 'T' m = 5 n = 4 alpha = 2.0E+00 lda = m call r4mat_test ( trans, lda, n, m, a ) do i = 1, m x(i) = real ( i ) end do incx = 1 beta = 3.0E+00 do i = 1, n y(i) = real ( 10 * i ) end do incy = 1 call r4mat_print ( m, n, a, ' Matrix A:' ) call r4vec_print ( m, x, ' Vector X:' ) call r4vec_print ( n, y, ' Vector Y:' ) call sgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy ) call r4vec_print ( n, y, & ' Result Y = alpha * A'' * x + beta * y' ) return end subroutine test02 ( ) c*********************************************************************72 c cc TEST02 tests SGBMV. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 16 September 2005 c c Author: c c John Burkardt c implicit none integer m integer n integer kl integer ku integer lda parameter ( m = 5 ) parameter ( n = 5 ) parameter ( kl = 1 ) parameter ( ku = 1 ) parameter ( lda = kl + 1 + ku ) real a(lda,n) real alpha real beta integer i integer incx integer incy integer j integer jhi integer jlo character trans real x(n) real y(m) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) ' For a general band matrix A,' write ( *, '(a)' ) ' SGBMV computes ' write ( *, '(a)' ) ' y := alpha * A * x + beta * y' trans = 'N' alpha = 2.0E+00 incx = 1 beta = 3.0E+00 incy = 1 do i = 1, m jlo = max ( 1, i - kl ) jhi = min ( n, i + ku ) do j = jlo, jhi if ( i == j ) then a(ku+1+i-j,j) = 2.0E+00 else if ( i == j - 1 .or. i == j + 1 ) then a(ku+1+i-j,j) = -1.0E+00 else a(ku+1+i-j,j) = 0.0E+00 end if end do end do do i = 1, n x(i) = real ( i ) end do do i = 1, m y(i) = real ( 10 * i ) end do call sgbmv ( trans, m, n, kl, ku, alpha, a, lda, x, & incx, beta, y, incy ) call r4vec_print ( m, y, ' Result vector Y' ) return end subroutine test03 ( ) c*********************************************************************72 c cc TEST03 tests SSYMV. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 16 September 2005 c c Author: c c John Burkardt c implicit none integer n integer lda parameter ( n = 5 ) parameter ( lda = n ) real a(lda,n) real alpha real beta integer i integer incx integer incy integer j character uplo real x(n) real y(n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03' write ( *, '(a)' ) ' For a general symmetric matrix A,' write ( *, '(a)' ) ' SSYMV computes ' write ( *, '(a)' ) ' y := alpha * A * x + beta * y' uplo = 'U' alpha = 2.0E+00 incx = 1 beta = 3.0E+00 incy = 1 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 a(i,j) = 0.0E+00 end if end do end do do i = 1, n x(i) = real ( i ) end do do i = 1, n y(i) = real ( 10 * i ) end do call ssymv ( uplo, n, alpha, a, lda, x, incx, beta, y, incy ) call r4vec_print ( n, y, ' Result vector Y' ) return end subroutine test04 ( ) c*********************************************************************72 c cc TEST04 tests SSBMV. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 16 September 2005 c c Author: c c John Burkardt c implicit none integer m integer n integer k integer lda parameter ( m = 5 ) parameter ( n = 5 ) parameter ( k = 1 ) parameter ( lda = k + 1 ) real a(lda,n) real alpha real beta integer i integer incx integer incy integer j integer jhi character uplo real x(n) real y(m) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST04' write ( *, '(a)' ) ' For a symmetric band matrix A,' write ( *, '(a)' ) ' SSBMV computes ' write ( *, '(a)' ) ' y := alpha * A * x + beta * y' uplo = 'U' alpha = 2.0E+00 incx = 1 beta = 3.0E+00 incy = 1 do i = 1, m jhi = min ( n, i + k ) do j = i, jhi if ( i == j ) then a(k+1+i-j,j) = 2.0E+00 else if ( i == j - 1 ) then a(k+1+i-j,j) = -1.0E+00 else a(k+1+i-j,j) = 0.0E+00 end if end do end do do i = 1, n x(i) = real ( i ) end do do i = 1, m y(i) = real ( 10 * i ) end do call ssbmv ( uplo, n, k, alpha, a, lda, x, incx, beta, y, incy ) call r4vec_print ( m, y, ' Result vector Y' ) return end subroutine test05 ( ) c*********************************************************************72 c cc TEST05 tests SGER. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 04 April 2014 c c Author: c c John Burkardt c implicit none integer m parameter ( m = 5 ) integer n parameter ( n = 4 ) integer lda parameter ( lda = m ) real a(lda,n) real alpha integer i integer incx integer incy character trans real x(m) real y(n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05' write ( *, '(a)' ) ' For a general matrix A,' write ( *, '(a)' ) ' SGER computes A := A + alpha * x * y''' alpha = 2.0E+00 trans = 'N' call r4mat_test ( trans, lda, m, n, a ) do i = 1, m x(i) = real ( i ) end do incx = 1 do i = 1, n y(i) = real ( 10 * i ) end do incy = 1 call r4mat_print ( m, n, a, ' Matrix A:' ) call r4vec_print ( m, x, ' Vector X:' ) call r4vec_print ( n, y, ' Vector Y:' ) call sger ( m, n, alpha, x, incx, y, incy, a, lda ) call r4mat_print ( m, n, a, ' Result A = A + alpha * x * y' ) return end subroutine test06 ( ) c*********************************************************************72 c cc TEST06 tests STRMV. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 04 April 2014 c c Author: c c John Burkardt c implicit none integer m parameter ( m = 5 ) integer lda parameter ( lda = m ) integer n parameter ( n = m ) real a(lda,n) character diag integer i integer incx integer j integer test character trans character uplo real x(n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06' write ( *, '(a)' ) ' For a triangular matrix A,' write ( *, '(a)' ) ' STRMV computes y := A * x or y := A'' * x' do test = 1, 2 uplo = 'U' if ( test .eq. 1 ) then trans = 'N' else trans = 'T' end if diag = 'N' do j = 1, n do i = 1, j a(i,j) = real ( i + j ) end do do i = j + 1, m a(i,j) = 0.0E+00 end do end do incx = 1 do i = 1, n x(i) = real ( i ) end do call strmv ( uplo, trans, diag, n, a, lda, x, incx ) if ( trans .eq. 'N' ) then call r4vec_print ( n, x, ' Result y = A * x' ); else call r4vec_print ( n, x, ' Result y = A'' * x' ); end if end do return end function lsame ( ca, cb ) c*********************************************************************72 c cc lsame() returns TRUE if CA is the same letter as CB regardless of case. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 February 2014 c c Author: c c Original FORTRAN77 version by Jack Dongarra. c This version by John Burkardt. c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, character CA, CB, the character to compare. c c Output, logical LSAME, is TRUE if the characters are equal, c disregarding case. c implicit none character ca character cb intrinsic ichar integer inta integer intb logical lsame integer zcode c c Test if the characters are equal. c lsame = ca .eq. cb if ( lsame ) then return end if c c Now test for equivalence if both characters are alphabetic. c zcode = ichar ( 'Z' ) c c Use 'Z' rather than 'A' so that ASCII can be detected on Prime c machines, on which ICHAR returns a value with bit 8 set. c ICHAR('A') on Prime machines returns 193 which is the same as c ICHAR('A') on an EBCDIC machine. c inta = ichar ( ca ) intb = ichar ( cb ) c c ASCII is assumed. c ZCODE is the ASCII code of either lower or upper case 'Z'. c if ( zcode .eq. 90 .or. zcode .eq. 122 ) then if ( inta .ge. 97 .and. inta .le. 122 ) then inta = inta - 32 end if if ( intb .ge. 97 .and. intb .le. 122 ) then intb = intb - 32 end if c c EBCDIC is assumed. c ZCODE is the EBCDIC code of either lower or upper case 'Z'. c else if ( zcode .eq. 233 .or. zcode .eq. 169 ) then if ( inta .ge. 129 .and. inta .le. 137 .or. & inta .ge. 145 .and. inta .le. 153 .or. & inta .ge. 162 .and. inta .le. 169 ) then inta = inta + 64 end if if ( intb .ge. 129 .and. intb .le. 137 .or. & intb .ge. 145 .and. intb .le. 153 .or. & intb .ge. 162 .and. intb .le. 169 ) then intb = intb + 64 end if c c ASCII is assumed, on Prime machines. c ZCODE is the ASCII code plus 128 of either lower or upper case 'Z'. c else if ( zcode .eq. 218 .or. zcode .eq. 250 ) then if ( inta .ge. 225 .and. inta .le. 250 ) then inta = inta - 32 end if if ( intb .ge. 225 .and. intb .le. 250 ) then intb = intb - 32 end if end if lsame = inta .eq. intb return end subroutine r4mat_print ( m, n, a, title ) c*********************************************************************72 c cc r4mat_print() prints an R4MAT. c c Discussion: c c An R4MAT is an array of R4 values. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 23 August 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, the number of rows in A. c c Input, integer N, the number of columns in A. c c Input, real A(M,N), the matrix. c c Input, character*(*) TITLE, a title. c implicit none integer m integer n real a(m,n) character*(*) title call r4mat_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, & title ) c*********************************************************************72 c cc r4mat_print_some() prints some of an R4MAT. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 25 January 2007 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, N, the number of rows and columns. c c Input, real A(M,N), an M by N matrix to be printed. c c Input, integer ILO, JLO, the first row and column to print. c c Input, integer IHI, JHI, the last row and column to print. c c Input, character ( len = * ) TITLE, a title. c implicit none integer incx parameter ( incx = 5 ) integer m integer n real a(m,n) character * ( 14 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character * ( * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m .le. 0 .or. n .le. 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return 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), '(i7,7x)') j end do write ( *, '('' Col '',5a14)' ) ( ctemp(j), j = 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 write ( ctemp(j2), '(g14.6)' ) a(i,j) end do write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc ) end do end do return end subroutine r4mat_test ( trans, lda, m, n, a ) c*********************************************************************72 c cc r4mat_test() sets up a test matrix. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 10 February 2014 c c Author: c c John Burkardt. c c Parameters: c c Input, character * ( 1 ) TRANS, indicates whether matrix is to be c transposed. c 'N', no transpose. c 'T', transpose the matrix. c c Input, integer LDA, the leading dimension of the matrix. c c Input, integer M, N, the number of rows and columns of c the matrix. c c Output, real A(LDA,*), the matrix. c if TRANS is 'N', then the matrix is stored in LDA*N entries, c as an M x N matrix; c if TRANS is 'T', then the matrix is stored in LDA*M entries, c as an N x M matrix. c implicit none integer lda real a(lda,*) integer i integer j integer m integer n character * ( 1 ) trans if ( trans .eq. 'N' ) then do j = 1, n do i = 1, m a(i,j) = real ( 10 * i + j ) end do end do else do j = 1, n do i = 1, m a(j,i) = real ( 10 * i + j ) end do end do end if return end subroutine r4vec_print ( n, a, title ) c*********************************************************************72 c cc r4vec_print() prints an R4VEC. c c Discussion: c c An R4VEC is a vector of R4's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of components of the vector. c c Input, real A(N), the vector to be printed. c c Input, character * ( * ) TITLE, a title. c implicit none integer n real a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i) end do return end subroutine timestamp ( ) c*********************************************************************72 c cc timestamp() prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 16 September 2005 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character ( len = 8 ) date character ( len = 10 ) time call date_and_time ( date, time ) write ( *, '(a8,2x,a10)' ) date, time return end subroutine xerbla ( srname, info ) c*********************************************************************72 c cc xerbla() is an error handler. c c Discussion: c c XERBLA is called if an input parameter has an invalid value. c A message is printed and execution stops. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 February 2014 c c Author: c c Original FORTRAN77 version by Jack Dongarra. c This version by John Burkardt. c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, character * ( * ) SRNAME, the name of the routine c which called XERBLA. c c Input, integer INFO, the position of the invalid parameter in the c parameter list of the calling routine. c implicit none integer info character * ( * ) srname write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'XERBLA - Fatal error!' write ( *, '(a)' ) ' On entry to routine ' // trim ( srname ) write ( *, '(a)' ) ' an illegal value was detected for' write ( *, '(a,i6)' ) ' parameter number ', info stop 1 end