program main

!*****************************************************************************80
!
!! MAIN is the main program for ESSL_PRB.
!
!  Discussion:
!
!    ESSL_PRB calls some of the ESSL routines.
!
!  Modified:
!
!    04 May 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ESSL_PRB'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Tests for the ESSL library.'

  call test005
  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09

  call test10
  call test11
  call test12
  call test13
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ESSL_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine test005

!*****************************************************************************80
!
!! TEST005 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 ).
!
  implicit none

  integer, parameter :: m = 8
  integer, parameter :: n = 6

  integer, parameter :: k = min ( m, n )
  integer, parameter :: lwork = n

  real ( kind = 8 ) a(m,n)
  integer i
  integer info
  integer lda
  real ( kind = 8 ) q(m,k)
  real ( kind = 8 ) r(k,n)
  integer seed
  real ( kind = 8 ) tau(k)
  real ( kind = 8 ) work(lwork)

  seed = 123456789

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST005'
  write ( *, '(a)' ) '  For a double precision 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,i6)' ) '  M = ', m
  write ( *, '(a,i6)' ) '  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,i6)' ) '  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,i6)' ) '  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 test01

!*****************************************************************************80
!
!! TEST01 tests DSRIS.
!
!  Modified:
!
!    21 September 2005
!
  implicit none

  integer, parameter :: n = 9
  integer, parameter :: naux1 = 122
  integer, parameter :: naux2 = 63

  double precision, dimension ( 22 ) :: ar = (/ &
     2.0D+00,  2.0D+00, -1.0D+00,  1.0D+00,  2.0D+00, &
     1.0D+00,  2.0D+00, -1.0D+00,  1.0D+00,  2.0D+00, &
    -1.0D+00,  1.0D+00,  2.0D+00, -1.0D+00,  1.0D+00, &
     2.0D+00, -1.0D+00,  1.0D+00,  2.0D+00, -1.0D+00, &
     1.0D+00,  2.0D+00 /)
  double precision, dimension ( naux1 ) :: aux1
  double precision, dimension ( naux2 ) :: aux2
  double precision, dimension ( n ) :: b = (/ &
     2.0D+00, 1.0D+00, 3.0D+00, 2.0D+00, 2.0D+00, &
     2.0D+00, 2.0D+00, 2.0D+00, 3.0D+00 /)
  integer i
  integer, dimension ( n+1 ) :: ia = (/ &
       1, 2, 4, 6, 9, 12, 15, 18, 21, 23 /)
  character :: init = 'I'
  integer, dimension ( 6 ) :: iparm = (/ 20, 4, 0, 1, 10, 0 /)
  integer, dimension ( 22 ) :: ja = (/ &
       1,   2,    3,    2,   3,   1,    4,    5,   4,    5,  &
       6,   5,    6,    7,   6,   7,    8,    7,   8,    9,  &
       8,   9 /)
  character :: stor = 'G'
  double precision, dimension ( 3 ) :: rparm = (/ 1.0D-7, 0.0D+00, 1.0D+00 /)
  double precision, dimension ( n ) :: x = (/ &
     0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
     0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /)
  double precision, dimension ( n ) :: x_exact = (/ &
     1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
     1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  DSRIS is an iterative solver for a linear system.'

  call dsris ( stor, init, n, ar, ja, ia, b, x, iparm, rparm, aux1, naux1, &
    aux2, naux2 )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of iterations taken = ', iparm(6)
  write ( *, '(a,g14.6)' ) '  Relative accuracy = ', rparm(2)
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Computed  Exact      Error'
  write ( *, '(a)' ) '  Solution  Solution'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(3g14.6)' ) x(i), x_exact(i), abs ( x(i) - x_exact(i) )
  end do

  return
end
subroutine test02

!*****************************************************************************80
!
!! TEST02 tests SCSINT.
!
  implicit none

  integer, parameter :: m = 20
  integer, parameter :: n = 7

  real c(n,4)
  integer i
  integer incx
  integer init
  real runge
  real s(m)
  real s2(m)
  real t(m)
  real x(n)
  real y(n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  SCSINT sets up cubic spline interpolation.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  call random_number ( harvest = x(1:n) )
  x(1:n) = -5.0E+00 * ( 1.0E+00 - x(1:n) ) + 5.0E+00 * x(1:n)
  incx = 1
  call ssort ( x, incx, n )

  do i = 1, n
    y(i) = runge ( x(i) )
  end do

  call r4vec2_print ( n, x, y, '  The (X,Y) data:' )
!
!  Compute the cubic spline interpolant.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SCSINT.'

  call random_number ( harvest = t(1:m) )
  t(1:m) = -5.0E+00 * ( 1.0E+00 - t(1:m) ) + 5.0E+00 * t(1:m)
  incx = 1
  call ssort ( t, incx, m )

  init = 0
  call scsint ( x, y, c, n, init, t, s, m )

  do i = 1, m
    s2(i) = runge ( t(i) )
  end do

  call r4vec3_print ( m, t, s, s2, '  X, interpolated Y, Runge(X):' )
  
  return
end
subroutine test03

!*****************************************************************************80
!
!! TEST03 tests SGEEV.
!
!  Discussion:
!
!    The problem is just an enlarged version of the
!    problem for n = 5, which is:
!
!    Matrix A is ( 1  C  C  C  C)
!                ( C  1  C  C  C)
!                ( C  C  1  C  C)
!                ( C  C  C  1  C)
!                ( C  C  C  C  1)
!
!    where C is a constant.
!
!    There are N-1 eigenvalues equal to 1-C, and 1 eigenvalue
!    equal to 1+C*(N-1).
!
  implicit none

  integer, parameter :: n = 25

  real a(n,n)
  real aux(1)
  integer i
  integer iopt
  logical select(n)
  complex w(n)
  complex z(n,n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  SGEEV computes the eigenvalues and eigenvectors'
  write ( *, '(a)' ) '  of a real, general matrix.'
  write ( *, '(a)' ) ' '
!
!  Assign values to the matrix.
!
  a(1:n,1:n) = 5.0E+00

  do i = 1, n
    a(i,i) = 1.0E+00
  end do
!
!  IOPT = 
!    0 for eigenvalues only, 
!    1 for all eigenvalues and eigenvectors,
!    2 for selected eigenvalues and eigenvectors.
!
  iopt = 1

  call sgeev ( iopt, a, n, w, z, n, select, n, aux, 0 )
 
  call c4vec_print_some ( n, w, 10, '  Computed eigenvalues:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Eigenvectors #1 and #N:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,i3,2g14.6,5x,2g14.6)' ) i, z(i,1), z(i,n)
  end do

  return
end
subroutine test04

!*****************************************************************************80
!
!! TEST04 tests SGEF and SGES.
!
  implicit none

  integer, parameter :: n = 3

  real a(n,n)
  real b(n)
  integer i
  integer iopt
  integer ipvt(n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  SGEF factors a general matrix;'
  write ( *, '(a)' ) '  SGES solves a linear system;'
!
!  Set the matrix.
!
  a(1,1) = 1.0E+00
  a(1,2) = 2.0E+00
  a(1,3) = 3.0E+00

  a(2,1) = 4.0E+00
  a(2,2) = 5.0E+00
  a(2,3) = 6.0E+00

  a(3,1) = 7.0E+00
  a(3,2) = 8.0E+00
  a(3,3) = 0.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  System matrix A:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,5g14.6)' ) a(i,1:n)
  end do
!
!  Set the right hand side.
!
  b(1:3) = (/ 14.0E+00, 32.0E+00, 23.0E+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Right hand side b:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do
!
!  Factor the matrix.
!
  call sgef ( a, n, n, ipvt )
!
!  Solve the linear system.
!
  iopt = 0
  call sges ( a, n, n, ipvt, b, iopt )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Solution X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(g14.6)' ) b(i)
  end do

  return
end
subroutine test05

!*****************************************************************************80
!
!! TEST05 demonstrates SGEMM.
!
  implicit none

  integer, parameter :: nn = 5

  real a(nn,nn)
  real alpha
  real b(nn,nn)
  real beta
  real c(nn,nn)
  integer i
  integer k
  integer lda
  integer ldb
  integer ldc
  integer m
  integer n

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  SGEMM multiplies two real general matrices.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Matrix order N = ', nn

  call helmert ( nn, a )

  b(1:nn,1:nn) = a(1:nn,1:nn)

  m = nn
  n = nn
  k = nn
  alpha = 1.0E+00
  beta = 0.0E+00
  lda = nn
  ldb = nn
  ldc = nn

  call sgemm ( 'Transpose', 'No transpose', m, n, k, alpha, a, lda, &
    b, ldb, beta, c, ldc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Product matrix C'
  write ( *, '(a)' ) '  (Should be the identity matrix)'
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(5f12.4)' ) c(i,1:n)
  end do

  return
end
subroutine test06

!*****************************************************************************80
!
!! TEST06 tests SGETRF and SGETRS.
!
!  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.
!
  implicit none

  integer, parameter :: n = 25
  integer, parameter :: lda = n

  real a(lda,n)
  real b(n)
  integer i
  integer info
  integer ipvt(n)
  integer j

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  SGETRF factors a general matrix;'
  write ( *, '(a)' ) '  SGETRS 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) = n
      else
        a(i,j) = -1.0E+00
      end if
    end do
  end do

  b(1:n) = 1.0E+00
!
!  Factor the matrix.
!
  call sgetrf ( n, n, a, lda, ipvt, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Matrix is singular, INFO = ', info
    return
  end if
!
!  Solve the linear system.
!
  call sgetrs ( 'n', n, 1, a, lda, ipvt, b, n, info )
 
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Solution procedure failed, INFO = ', info
    return
  end if

  call r4vec_print_some ( n, b, 10, '  Computed solution:' )
 
  return
end
subroutine test07

!*****************************************************************************80
!
!! TEST07 tests SGLNQ.
!
  implicit none

  integer, parameter :: n_size = 20

  real, parameter :: a = 0.0E+00
  real, parameter :: b = 1.0E+00
  real error
  real, parameter :: exact = -4.0E+00 
  integer i
  integer n
  integer, dimension ( n_size ) :: n_array = (/ &
    1, 2, 3, 4, 5, 6, 8, 10, 12, 14, 16, 20, 24, 32, 40, 48, 64, 96, 128, 256 /)
  real result
  real sglnq
  external test07_f

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  SGLNQ integrates a function F(X) from A to B.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We will integrate log(x)/sqrt(x) from 0 to 1.'
  write ( *, '(a)' ) '  The exact answer is -4.'
  write ( *, '(a)' ) '  We will increase N, the quadrature order.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     N  Result   Abs(Result+4)'
  write ( *, '(a)' ) ' '

  do i = 1, n_size
    n = n_array(i)
    result = sglnq ( test07_f, a, b, n )
    error = abs ( result - exact )
    write ( *, '(i6,f10.6,g14.6)' ) n, result, error
  end do

  return
end
subroutine test07_f ( x, y, n )

!*****************************************************************************80
!
!! TEST07_F evaluates the function to be integrated.
!
  implicit none

  integer n

  real x(n)
  real y(n)

  y(1:n) = log ( x(1:n) ) / sqrt ( x(1:n) )

  return
end
subroutine test08

!*****************************************************************************80
!
!! TEST08 tests SNRAND.
!
  implicit none

  integer, parameter :: n_max = 1000

  real aux(1)
  integer n
  integer, parameter :: naux = 0
  double precision seed
  real x(n_max)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  SNRAND generates normal random numbers.'

  seed = 123456789.0D+00
  n = n_max / 2

  write ( *, '(a)' ) ' '
  write ( *, '(a,g17.9)' ) '  Initial seed = ', seed
  write ( *, '(a,i6)' ) '  Number of values to generate = ', n

  call snrand ( seed, n, x, aux, naux )
  call r4vec_print_some ( n, x, 10, '  Random vector:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g17.9)' ) '  Current seed = ', seed
  write ( *, '(a,i6)' ) '  Number of values to generate = ', n

  call snrand ( seed, n, x, aux, naux )
  call r4vec_print_some ( n, x, 10, '  Random vector:' )

  write ( *, * ) ' '
  write ( *, * ) '  Now restore seed to initial value,'
  write ( *, * ) '  and recompute the two sets of data at one call.'

  seed = 123456789.0D+00
  n = n_max
  write ( *, '(a)' ) ' '
  write ( *, '(a,g17.9)' ) '  Current seed = ', seed
  write ( *, '(a,i6)' ) '  Number of values to generate = ', n

  call snrand ( seed, n, x, aux, naux )
  call r4vec_print_some ( n, x, 10, '  Random vector:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Maximum value = ', maxval ( x(1:n) )
  write ( *, '(a,g14.6)' ) '  Minimum value = ', minval ( x(1:n) )
  write ( *, '(a,g14.6)' ) '  Average value = ', sum ( x(1:n) ) / real ( n )

  return
end
subroutine test09

!*****************************************************************************80
!
!! TEST09 tests SPINT.
!
  implicit none

  integer, parameter :: m = 20
  integer, parameter :: n = 7

  real c(n)
  integer incx
  integer ninit
  real s(m)
  real s2(m)
  real t(m)
  real x(n)
  real y(n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  SPINT sets up polynomial interpolation.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  call random_number ( harvest = x(1:n) )
  x(1:n) = 10.0E+00 * x(1:n)
  incx = 1
  call ssort ( x, incx, n )

  y(1:n) = sin ( x(1:n) )

  call r4vec2_print ( n, x, y, '  The (X,Y) data:' )
!
!  Compute the polynomial interpolant.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SPINT.'

  call random_number ( harvest = t(1:m) )
  t(1:m) = 10.0E+00 * t(1:m)
  incx = 1
  call ssort ( t, incx, m )

  ninit = 0
  call spint ( x, y, n, c, ninit, t, s, m )

  s2(1:m) = sin ( t(1:m) )
  call r4vec3_print ( m, t, s, s2, '  X, interpolated Y, SIN(X):' )
  
  return
end
subroutine test10

!*****************************************************************************80
!
!! TEST10 tests SRCFT.
!
!  Discussion:
!
!    I found the interface to SRCFT to be awkward and poorly documented.
!    The separate initialization for forward AND reverse calculations
!    is confusing.  A formula for the size of the first auxilliary 
!    storage is not given; the suggestion is that maybe 25,000 would
!    be suitable for problems with n <= 16384!  All in all, a clumsy, 
!    unfriendly and error-prone tool, JVB 13 December 2002.
!
  implicit none

  integer, parameter :: m = 1
  integer, parameter :: n = 10
  integer, parameter :: naux1 = 100

  real aux1(naux1)
  real aux2(1)
  real aux3(1)
  integer inc2x
  integer inc2y
  integer init
  integer isgn
  integer naux2
  integer naux3
  real scale
  real x(n,m)
  complex y(1:(n/2)+1,m)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  For real fast Fourier transforms,'
  write ( *, '(a)' ) '  SRCFT computes the complex FFT of real data.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the X values as random numbers.
!
  call random_number ( harvest = x(1:n,1:m) )

  call r4vec_print_some ( n, x, 10, '  The original data:' )
!
!  Set some values that don't change.
!
  inc2x = n
  inc2y = ( n / 2 ) + 1
  naux2 = 0
  naux3 = 0
!
!  Initialize the SRCFT routine.
!  (Apparently, you MUST specify a value for ISGN, even though
!  it is NOT needed yet!)
!  (And you have to specify SCALE here...)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Initialize the FFT routine.'

  init = 1
  isgn = +1
  scale = 1.0E+00

  call srcft ( init, x, inc2x, y, inc2y, n, m, isgn, scale, aux1, naux1, &
    aux2, naux2, aux3, naux3 )
!
!  Compute the FFT of the real data.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Compute the FFT coeficients.'

  init = 0
  isgn = +1

  call srcft ( init, x, inc2x, y, inc2y, n, m, isgn, scale, aux1, naux1, &
    aux2, naux2, aux3, naux3 )

  call c4vec_print_some ( (n/2)+1, y, 10, '  The Fourier coefficients:' )
!
!  Reinitialize the SRCFT routine for inverse operations.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Reinitialize the FFT routine for inverse operations.'

  init = 1
  isgn = -1
  scale = 1.0E+00

  call srcft ( init, x, inc2x, y, inc2y, n, m, isgn, scale, aux1, naux1, &
    aux2, naux2, aux3, naux3 )
!
!  Now compute the inverse FFT of coefficients.  Should get back the
!  original data.  
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Recover the data from the FFT coeficients.'

  init = 0
  isgn = -1

  call srcft ( init, x, inc2x, y, inc2y, n, m, isgn, scale, aux1, naux1, &
    aux2, naux2, aux3, naux3 )

  call r4vec_print_some ( n, x, 10, '  The retrieved data:' )
  
  return
end
subroutine test11

!*****************************************************************************80
!
!! TEST11 tests SSORT.
!
  implicit none

  integer, parameter :: n = 1000

  integer incx
  real x_vec(n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  SSORT sorts a vector of real data.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the X values as random numbers.
!
  call random_number ( harvest = x_vec(1:n) )

  call r4vec_print_some ( n, x_vec, 10, '  The original data:' )
!
!  Sort the data.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call SSORT.'

  incx = 1
  call ssort ( x_vec, incx, n )

  call r4vec_print_some ( n, x_vec, 10, '  The sorted data:' )
  
  return
end
subroutine test12

!*****************************************************************************80
!
!! TEST12 tests SURAND.
!
  implicit none

  integer, parameter :: n_max = 1000

  integer n
  double precision seed
  real x(n_max)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  SURAND generates real random numbers.'

  seed = 123456789.0D+00
  n = n_max / 2

  write ( *, '(a)' ) ' '
  write ( *, '(a,g17.9)' ) '  Initial seed = ', seed
  write ( *, '(a,i6)' ) '  Number of values to generate = ', n

  call surand ( seed, n, x )
  call r4vec_print_some ( n, x, 10, '  Random vector:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g17.9)' ) '  Current seed = ', seed
  write ( *, '(a,i6)' ) '  Number of values to generate = ', n

  call surand ( seed, n, x )
  call r4vec_print_some ( n, x, 10, '  Random vector:' )

  write ( *, * ) ' '
  write ( *, * ) '  Now restore seed to initial value,'
  write ( *, * ) '  and recompute the two sets of data at one call.'

  seed = 123456789.0D+00
  n = n_max
  write ( *, '(a)' ) ' '
  write ( *, '(a,g17.9)' ) '  Current seed = ', seed
  write ( *, '(a,i6)' ) '  Number of values to generate = ', n

  call surand ( seed, n, x )
  call r4vec_print_some ( n, x, 10, '  Random vector:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Maximum value = ', maxval ( x(1:n) )
  write ( *, '(a,g14.6)' ) '  Minimum value = ', minval ( x(1:n) )
  write ( *, '(a,g14.6)' ) '  Average value = ', sum ( x(1:n) ) / real ( n )

  return
end
subroutine test13

!*****************************************************************************80
!
!! TEST13 tests DGBF and DGBS.
!
!  Discussion:
!
!    This routine demonstrates how to set up and solve a banded
!    linear system with the IBM ESSL library.
!
!    Note that the matrix format is similar, but the names and
!    argument lists are DIFFERENT from the ones used in LAPACK.
!    The LAPACK routines SGBTRF/DGBTFR and SGBTRS/DGBTRS are
!    not available in ESSL!
!
  implicit none

  integer, parameter :: n = 25
  integer, parameter :: ml = 3
  integer, parameter :: mu = 3
!
!  Note the value of LDA here, which is NOT the LINPACK/LAPACK standard!
!
  integer, parameter :: lda = 2 * ml + mu + 16

  real ( kind = 8 ) a(lda,n)
  real ( kind = 8 ) b(n)
  integer i
  integer ihi
  integer ilo
  integer info
  integer ipiv(n)
  integer j
  integer m
  real ( kind = 8 ) temp

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  For a double precision matrix (D)'
  write ( *, '(a)' ) '  in general band storage mode (GB):'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  DGBF factors a general band matrix.'
  write ( *, '(a)' ) '  DGBS 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,i6)' ) '  Bandwidth is ', m

  call dgbf ( a, lda, n, ml, mu, ipiv )
!
!   Solve the linear system.
!
  call dgbs ( a, lda, n, ml, mu, ipiv, b )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First five entries of solution.'
  write ( *, '(a)' ) '  (Should all be 1.)'
  write ( *, '(a)' ) ' '

  do i = 1, 5
    write ( *, '(2x,i6,2x,g14.6)' ) i, b(i)
  end do

  return
end
subroutine c4vec_print_some ( n, x, max_print, title )

!*****************************************************************************80
!
!! C4VEC_PRINT_SOME prints some of a C4VEC.
!
!  Discussion:
!
!    The user specifies MAX_PRINT, the maximum number of lines to print.
!
!    If N, the size of the vector, is no more than MAX_PRINT, then
!    the entire vector is printed, one entry per line.
!
!    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
!    followed by a line of periods suggesting an omission,
!    and the last entry.
!
!  Modified:
!
!    17 December 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries of the vector.
!
!    Input, complex X(N), the vector to be printed.
!
!    Input, integer MAX_PRINT, the maximum number of lines to print.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer n

  integer i
  integer max_print
  character ( len = * ) title
  complex x(n)

  if ( max_print <= 0 ) then
    return
  end if

  if ( n <= 0 ) then
    return
  end if

  if ( len_trim ( title ) > 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
    write ( *, '(a)' ) ' '
  end if

  if ( n <= max_print ) then

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

  else if ( max_print >= 3 ) then

    do i = 1, max_print-2
      write ( *, '(i6,2x,2g14.6)' ) i, x(i)
    end do
    write ( *, '(a)' ) '......  ..............'
    i = n
    write ( *, '(i6,2x,2g14.6)' ) i, x(i)

  else

    do i = 1, max_print - 1
      write ( *, '(i6,2x,2g14.6)' ) i, x(i)
    end do
    i = max_print
    write ( *, '(i6,2x,2g14.6,2x,a)' ) i, x(i), '...more entries...'

  end if

  return
end
subroutine helmert ( n, a )

!*****************************************************************************80
!
!! HELMERT returns the classic Helmert matrix.
!
!  Formula:
!
!    If I = 1 then
!      A(I,J) = 1 / sqrt ( N )
!    else if J < I then
!      A(I,J) = 1 / sqrt ( I * ( I - 1 ) )
!    else if J = I then
!      A(I,J) = (1-I) / sqrt ( I * ( I - 1 ) )
!    else
!      A(I,J) = 0
!
!  Discussion:
!
!    The matrix given above by Helmert is the classic example of
!    a family of matrices which are now called Helmertian or
!    Helmert matrices.
!
!    A matrix is a (standard) Helmert matrix if it is orthogonal,
!    and the elements which are above the diagonal and below the
!    first row are zero.
!
!    If the elements of the first row are all strictly positive,
!    then the matrix is a strictly Helmertian matrix.
!
!    It is possible to require in addition that all elements below
!    the diagonal be strictly positive, but in the reference, this
!    condition is discarded as being cumbersome and not useful.
!
!    A Helmert matrix can be regarded as a change of basis matrix
!    between a pair of orthonormal coordinate bases.  The first row
!    gives the coordinates of the first new basis vector in the old
!    basis.  Each later row describes combinations of (an increasingly
!    extensive set of) old basis vectors that constitute the new
!    basis vectors.
!
!    Helmert matrices have important applications in statistics.
!
!  Example:
!
!    N = 5
!
!    0.4472    0.4472    0.4472    0.4472    0.4472
!    0.7071   -0.7071         0         0         0
!    0.4082    0.4082   -0.8165         0         0
!    0.2887    0.2887    0.2887   -0.8660         0
!    0.2236    0.2236    0.2236    0.2236   -0.8944
!
!  Properties:
!
!    A is generally not symmetric: A' /= A.
!
!    A is orthogonal: A' * A = A * A' = I.
!
!    Because A is orthogonal, it is normal: A' * A = A * A'.
!
!    A is not symmetric: A' /= A.
!
!  Reference:
!
!    H O Lancaster,
!    The Helmert Matrices,
!    American Mathematical Monthly,
!    Volume 72, 1965, pages 4-12.
!
!  Modified:
!
!    15 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the order of A.
!
!    Output, real A(N,N), the matrix.
!
  implicit none

  integer n

  real a(n,n)
  integer i
  integer j
!
!  A begins with the first row, diagonal, and lower triangle set to 1.
!
  do i = 1, n
    do j = 1, n

      if ( i == 1 ) then
        a(i,j) = 1.0E+00 / sqrt ( real ( n ) )
      else if ( j < i ) then
        a(i,j) = 1.0E+00 / sqrt ( real ( i * ( i - 1 ) ) )
      else if ( i == j ) then
        a(i,j) = real ( 1 - i ) / sqrt ( real ( i * ( i - 1 ) ) )
      else
        a(i,j) = 0.0E+00
      end if

    end do
  end do

  return
end
subroutine r8mat_print ( m, n, a, title )

!*****************************************************************************80
!
!! R8MAT_PRINT prints an R8MAT.
!
!  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 m
  integer 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.
!
!  Modified:
!
!    12 September 2004
!
!  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, parameter :: incx = 5
  integer m
  integer n

  real ( kind = 8 ) a(m,n)
  character ( len = 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 ( 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), '(i7,7x)') 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

        write ( ctemp(j2), '(g14.6)' ) a(i,j)

      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.
!
!  Modified:
!
!    11 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, L E 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.
!
!    P A Lewis, A S Goodman, J M 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 m
  integer n

  integer i
  integer j
  integer k
  integer 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 + 2147483647
      end if

      r(i,j) = real ( seed, kind = 8 ) * 4.656612875D-10

    end do
  end do

  return
end
function runge ( x )

!*****************************************************************************80
!
!! RUNGE evaluates Runge's function.
!
!  Modified:
!
!    11 November 1999
!
!  Formula:
!
!    RUNGE(X) = 1 / ( 1 + 25 * X**2 )
!
!  Discussion:
!
!    Runge's function is a common test case for interpolation
!    and approximation, over the interval [-1,1].
!
!  Parameters:
!
!    Input, real X, the argument of Runge's function.
!
!    Output, real RUNGE, the value of Runge's function.
!
  implicit none

  real runge
  real x

  runge = 1.0E+00 / ( 1.0E+00 + 25.0E+00 * x**2 )

  return
end
subroutine r4vec_print ( n, a, title )

!*****************************************************************************80
!
!! R4VEC_PRINT prints an R4VEC.
!
!  Discussion:
!
!    If all the entries are integers, the data if printed
!    in integer format.
!
!  Modified:
!
!    19 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none

  integer n

  real a(n)
  integer i
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '

  if ( all ( a(1:n) == aint ( a(1:n) ) ) ) then
    do i = 1, n
      write ( *, '(i6,i6)' ) i, int ( a(i) )
    end do
  else if ( all ( abs ( a(1:n) ) < 1000000.0E+00 ) ) then
    do i = 1, n
      write ( *, '(i6,f14.6)' ) i, a(i)
    end do
  else
    do i = 1, n
      write ( *, '(i6,g14.6)' ) i, a(i)
    end do
  end if

  return
end
subroutine r4vec_print_some ( n, a, max_print, title )

!*****************************************************************************80
!
!! R4VEC_PRINT_SOME prints "some" of an R4VEC.
!
!  Discussion:
!
!    The user specifies MAX_PRINT, the maximum number of lines to print.
!
!    If N, the size of the vector, is no more than MAX_PRINT, then
!    the entire vector is printed, one entry per line.
!
!    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
!    followed by a line of periods suggesting an omission,
!    and the last entry.
!
!  Modified:
!
!    13 December 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries of the vector.
!
!    Input, real A(N), the vector to be printed.
!
!    Input, integer MAX_PRINT, the maximum number of lines to print.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer n

  real a(n)
  integer i
  integer max_print
  character ( len = * ) title

  if ( max_print <= 0 ) then
    return
  end if

  if ( n <= 0 ) then
    return
  end if

  if ( len_trim ( title ) > 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
    write ( *, '(a)' ) ' '
  end if

  if ( n <= max_print ) then

    if ( all ( a(1:n) == aint ( a(1:n) ) ) ) then
      do i = 1, n
        write ( *, '(i6,2x,i6)' ) i, int ( a(i) )
      end do
    else if ( all ( abs ( a(1:n) ) < 1000000.0E+00 ) ) then
      do i = 1, n
        write ( *, '(i6,2x,f14.6)' ) i, a(i)
      end do
    else
      do i = 1, n
        write ( *, '(i6,2x,g14.6)' ) i, a(i)
      end do
    end if

  else if ( max_print >= 3 ) then

    if ( all ( a(1:max_print-2) == aint ( a(1:max_print-2) ) ) ) then
      do i = 1, max_print-2
        write ( *, '(i6,2x,i6)' ) i, int ( a(i) )
      end do
    else if ( all ( abs ( a(1:max_print-2) ) < 1000000.0E+00 ) ) then
      do i = 1, max_print-2
        write ( *, '(i6,2x,f14.6)' ) i, a(i)
      end do
    else
      do i = 1, max_print-2
        write ( *, '(i6,2x,g14.6)' ) i, a(i)
      end do
    end if

    write ( *, '(a)' ) '......  ..............'
    i = n

    if ( a(n) == aint ( a(n) ) ) then
      write ( *, '(i6,2x,i6)' ) i, int ( a(i) )
    else if (  abs ( a(n) ) < 1000000.0E+00 ) then
      write ( *, '(i6,2x,f14.6)' ) i, a(i)
    else
      write ( *, '(i6,2x,g14.6)' ) i, a(i)
    end if

  else

    if ( all ( a(1:max_print-1) == aint ( a(1:max_print-1) ) ) ) then
      do i = 1, max_print-1
        write ( *, '(i6,2x,i6)' ) i, int ( a(i) )
      end do
    else if ( all ( abs ( a(1:max_print-1) ) < 1000000.0E+00 ) ) then
      do i = 1, max_print-1
        write ( *, '(i6,2x,f14.6)' ) i, a(i)
      end do
    else
      do i = 1, max_print-1
        write ( *, '(i6,2x,g14.6)' ) i, a(i)
      end do
    end if

    i = max_print

    if ( a(n) == aint ( a(n) ) ) then
      write ( *, '(i6,2x,i6,a)' ) i, int ( a(i) ), '...more entries...'
    else if (  abs ( a(n) ) < 1000000.0E+00 ) then
      write ( *, '(i6,2x,f14.6,a)' ) i, a(i), '...more entries...'
    else
      write ( *, '(i6,2x,g14.6,a)' ) i, a(i), '...more entries...'
    end if

  end if

  return
end
subroutine r4vec2_print ( n, a1, a2, title )

!*****************************************************************************80
!
!! R4VEC2_PRINT prints a pair of R4VEC's.
!
!  Modified:
!
!    13 December 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real A1(N), A2(N), the vectors to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none

  integer n

  real a1(n)
  real a2(n)
  integer i
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '

  if ( all ( a1(1:n) == aint ( a1(1:n) ) ) .and. &
       all ( a2(1:n) == aint ( a2(1:n) ) ) ) then
    do i = 1, n
      write ( *, '(i6,2i6)' ) i, int ( a1(i) ), int ( a2(i) )
    end do
  else if ( all ( abs ( a1(1:n) ) < 1000000.0E+00 ) .and. &
            all ( abs ( a2(1:n) ) < 1000000.0E+00 ) ) then
    do i = 1, n
      write ( *, '(i6,2f14.6)' ) i, a1(i), a2(i)
    end do
  else
    do i = 1, n
      write ( *, '(i6,2g14.6)' ) i, a1(i), a2(i)
    end do
  end if

  return
end
subroutine r4vec3_print ( n, a1, a2, a3, title )

!*****************************************************************************80
!
!! R4VEC3_PRINT prints a trio of R4VEC's.
!
!  Modified:
!
!    13 December 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real A1(N), A2(N), A3(N), the vectors to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none

  integer n

  real a1(n)
  real a2(n)
  real a3(n)
  integer i
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '

  if ( all ( a1(1:n) == aint ( a1(1:n) ) ) .and. &
       all ( a2(1:n) == aint ( a2(1:n) ) ) .and. &
       all ( a3(1:n) == aint ( a3(1:n) ) )) then
    do i = 1, n
      write ( *, '(i6,3i6)' ) i, int ( a1(i) ), int ( a2(i) ), int ( a3(i) )
    end do
  else if ( all ( abs ( a1(1:n) ) < 1000000.0E+00 )  .and. &
            all ( abs ( a2(1:n) ) < 1000000.0E+00 )  .and. &
            all ( abs ( a3(1:n) ) < 1000000.0E+00 ) ) then
    do i = 1, n
      write ( *, '(i6,3f14.6)' ) i, a1(i), a2(i), a3(i)
    end do
  else
    do i = 1, n
      write ( *, '(i6,3g14.6)' ) i, a1(i), a2(i), a3(i)
    end do
  end if

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Modified:
!
!    06 August 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer 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