program main

c*********************************************************************72
c
cc blas_test() tests blas().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    18 September 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'blas_test():'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test blas().'

      call test01 ( )
      call test02 ( )
      call test03 ( )
      call test04 ( )

      call test05 ( )
      call test06 ( )
      call test07 ( )
      call test08 ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'blas_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      subroutine test01 ( )

c*********************************************************************72
c
cc TEST01 tests CAXPY.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 February 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 5 )

      integer i
      complex s
      complex x(n)
      complex y(n)

      x(1) = (  2.0E+00, -1.0E+00 )
      x(2) = ( -4.0E+00, -2.0E+00 )
      x(3) = (  3.0E+00,  1.0E+00 )
      x(4) = (  2.0E+00,  2.0E+00 )
      x(5) = ( -1.0E+00, -1.0E+00 )

      y(1) = ( -1.0E+00,  0.0E+00 )
      y(2) = (  0.0E+00, -3.0E+00 )
      y(3) = (  4.0E+00,  0.0E+00 )
      y(4) = ( -3.0E+00,  4.0E+00 )
      y(5) = ( -2.0E+00,  0.0E+00 )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01'
      write ( *, '(a)' ) '  CAXPY adds a multiple of one'
      write ( *, '(a)' ) '  single precision complex vector to another.'

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  X = '
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i6,2x,f6.1,2x,f6.1)' ) i, x(i)
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Y = '
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i6,2x,f6.1,2x,f6.1)' ) i, y(i)
      end do

      s = ( 0.50E+00, -1.00E+00 )

      write ( *, '(a)' ) ' '
      write ( *, '(a,2g14.6)' ) '  The scalar multiplier is: ', s

      call caxpy ( n, s, x, 1, y, 1 )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  A * X + Y = '
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i6,2f10.6)' ) i, y(i)
      end do

      return
      end
      subroutine test02 ( )

c*********************************************************************72
c
cc test02 tests DASUM.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 March 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer lda
      integer ma
      integer na
      integer nx

      parameter ( lda = 6 )
      parameter ( ma = 5 )
      parameter ( na = 4 )
      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 * dble ( 2 * i )
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test02'
      write ( *, '(a)' ) '  DASUM adds the absolute values of '
      write ( *, '(a)' ) '  elements of a double precision vector.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  X = '
      write ( *, '(a)' ) ' '
      do i = 1, nx
        write ( *, '(2x,i6,g14.6)' ) i, x(i)
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  DASUM ( NX, X, 1 ) =   ',
     &  dasum ( nx, x, 1 )
      write ( *, '(a,g14.6)' ) '  DASUM ( NX/2, X, 2 ) = ',
     &  dasum ( nx/2, x, 2 )
      write ( *, '(a,g14.6)' ) '  DASUM ( 2, X, NX/2 ) = ',
     &  dasum ( 2, x, nx/2 )

      do i = 1, lda
        do j = 1, na
          a(i,j) = 0.0D+00
        end do
      end do

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

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Demonstrate with a matrix A:'
      write ( *, '(a)' ) ' '
      do i = 1, ma
        write ( *, '(2x,5g14.6)' ) ( a(i,j), j = 1, na )
      end do

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

      return
      end
      subroutine test03 ( )

c*********************************************************************72
c
cc TEST03 tests SDOT.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    22 February 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      integer lda
      integer ldb
      integer ldc

      parameter ( n = 5 )
      parameter ( lda = 10 )
      parameter ( ldb = 7 )
      parameter ( ldc = 6 )

      real a(lda,lda)
      real b(ldb,ldb)
      real c(ldc,ldc)
      integer i
      integer j
      real sdot
      real sum1
      real x(n)
      real y(n)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST03'
      write ( *, '(a)' ) '  SDOT computes the dot product of vectors.'

      do i = 1, n
        x(i) = real ( i )
      end do

      do i = 1, n
        y(i) = - real ( i )
      end do

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

      do i = 1, n
        do j = 1, n
          b(i,j) = real ( i - j )
        end do
      end do
c
c  To compute a simple dot product of two vectors, use a
c  call like this:
c
      sum1 = sdot ( n, x, 1, y, 1 )

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  Dot product of X and Y is ', sum1
c
c  To multiply a ROW of a matrix A times a vector X, we need to
c  specify the increment between successive entries of the row of A:
c
      sum1 = sdot ( n, a(2,1), lda, x, 1 )

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  Product of row 2 of A and X is ', sum1
c
c  Product of a column of A and a vector is simpler:
c
      sum1 = sdot ( n, a(1,2), 1, x, 1 )

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' )
     &  '  Product of column 2 of A and X is ', sum1
c
c  Here's how matrix multiplication, c = a*b, could be done
c  with SDOT:
c
      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 ( *, '(2x,5g14.6)' ) ( c(i,j), j = 1, n )
      end do

      return
      end
      subroutine test04 ( )

c*********************************************************************72
c
cc TEST04 tests ZSCAL.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    02 March 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 6 )

      double complex da
      integer i
      double complex x(n)

      do i = 1, n
        x(i) = cmplx ( 10 * i, i )
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST04'
      write ( *, '(a)' ) '  ZSCAL multiplies a double complex scalar '
      write ( *, '(a)' ) '  times a vector.'

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  X = '
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i6,2x,f6.1,2x,f6.1)' ) i, x(i)
      end do

      da = ( 5.0D+00, 0.0D+00 )
      call zscal ( n, da, x, 1 )
      write ( *, '(a)' ) ' '
      write ( *, '(a,2f8.4,a)' ) '  ZSCAL ( N, (', da, '), X, 1 )'
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i6,2x,f6.1,2x,f6.1)' ) i, x(i)
      end do

      do i = 1, n
        x(i) = cmplx ( 10 * i, i )
      end do

      da = ( -2.0D+00, 1.0D+00 )
      call zscal ( 3, da, x, 2 )
      write ( *, '(a)' ) ' '
      write ( *, '(a,2f8.4,a)' ) '  ZSCAL ( 3, (', da, '), X, 2 )'
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i6,2x,f6.1,2x,f6.1)' ) i, x(i)
      end do

      return
      end
      subroutine test05 ( )

c*********************************************************************72
c
cc TEST05 tests CGEMV.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    17 January 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer m
      parameter ( m = 5 )
      integer n
      parameter ( n = 5 )

      complex a(m,n)
      complex alpha
      complex beta
      integer i
      integer incx
      integer incy
      integer j
      integer lda
      character trans
      complex x(n)
      real x1
      real x2
      complex y(m)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST05'
      write ( *, '(a)' ) '  For a general matrix A,'
      write ( *, '(a)' ) 
     &  '  CGEMV computes y := alpha * A * x + beta * y'
    
      trans = 'N'
      alpha = ( 10.0E+00, 1.0E+00 ) 
      lda = m
      incx = 1
      beta = 3.0E+00
      incy = 1

      do i = 1, m
        do j = 1, n
          if ( i .eq. j ) then
            a(i,j) = ( 2.0E+00, 0.0E+00 )
          else if ( i .eq. j - 1 .or. i .eq. j + 1 ) then
            a(i,j) = ( -1.0E+00, 0.0E+00 )
          else
            a(i,j) = ( 0.0E+00, 0.0E+00 )
          end if
        end do
      end do

      x1 = 0.0E+00
      x2 = real ( n )
      do i = 1, n
        x(i) = cmplx ( x1, x2 )
        x1 = x1 + 1.0E+00
        x2 = x2 - 2.0E+00
      end do

      do i = 1, m
        y(i) = ( 100.0E+00, 1.0E+00 )
      end do

      call cgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Result vector Y = '
      write ( *, '(a)' ) ' '

      do i = 1, m
        write ( *, '(2x,2g14.6)' ) y(i)
      end do

      return
      end
      subroutine test06 ( )

c*********************************************************************72
c
cc TEST06 tests DSYMV.
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 )

      double precision a(lda,n)
      double precision alpha
      double precision beta
      integer i
      integer incx
      integer incy
      integer j
      character uplo
      double precision x(n)
      double precision y(n)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST06'
      write ( *, '(a)' ) '  For a general symmetric matrix A,'
      write ( *, '(a)' ) '  DSYMV computes '
      write ( *, '(a)' ) '  y := alpha * A * x + beta * y'

      uplo = 'U'
      alpha = 2.0D+00
      incx = 1
      beta = 3.0D+00
      incy = 1

      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
            a(i,j) = 0.0D+00
          end if
        end do
      end do

      do i = 1, n
        x(i) = dble ( i )
      end do

      do i = 1, n
        y(i) = dble ( 10 * i )
      end do

      call dsymv ( uplo, n, alpha, a, lda, x, incx, beta, y, incy )

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Result vector y:'
      write ( *, '(a)' ) ''
      do i = 1, n
        write ( *, '(g14.6)' ) y(i)
      end do

      return
      end
      subroutine test07 ( )

c*********************************************************************72
c
cc TEST07 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)' ) 'TEST07'
      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 )

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Result vector y:'
      write ( *, '(a)' ) ''
      do i = 1, m
        write ( *, '(g14.6)' ) y(i)
      end do

      return
      end
      subroutine test08 ( )

c*********************************************************************72
c
cc TEST08 tests ZGEMV.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    18 January 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer m
      parameter ( m = 5 )
      integer n
      parameter ( n = 5 )

      double complex a(m,n)
      double complex alpha
      double complex beta
      integer i
      integer incx
      integer incy
      integer j
      integer lda
      character trans
      double complex x(n)
      double precision x1
      double precision x2
      double complex y(m)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST08'
      write ( *, '(a)' ) '  For a general matrix A,'
      write ( *, '(a)' ) 
     &  '  ZGEMV computes y := alpha * A * x + beta * y'
    
      trans = 'N'
      alpha = ( 10.0D+00, 1.0D+00 ) 
      lda = m
      incx = 1
      beta = 3.0D+00
      incy = 1

      do i = 1, m
        do j = 1, n
          if ( i .eq. j ) then
            a(i,j) = ( 2.0D+00, 0.0D+00 )
          else if ( i .eq. j - 1 .or. i .eq. j + 1 ) then
            a(i,j) = ( -1.0D+00, 0.0D+00 )
          else
            a(i,j) = ( 0.0D+00, 0.0D+00 )
          end if
        end do
      end do

      x1 = 0.0D+00
      x2 = dble ( n )
      do i = 1, n
        x(i) = dcmplx ( x1, x2 )
        x1 = x1 + 1.0D+00
        x2 = x2 - 2.0D+00
      end do

      do i = 1, m
        y(i) = ( 100.0D+00, 1.0d+00 )
      end do

      call zgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Result vector Y = '
      write ( *, '(a)' ) ' '

      do i = 1, m
        write ( *, '(2x,2g14.6)' ) y(i)
      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 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