program fftpack_prb
!
!*******************************************************************************
!
!! FFTPACK_PRB calls the FFTPACK test routines.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FFTPACK_PRB'
  write ( *, '(a)' ) '  A set of tests for FFTPACK.'

  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
  call test14
  call test15
  call test16
  call test17
  call test18
  call test19
  call test20

  call test21
  call test22
  call test23

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FFTPACK_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests CFFTB.
!! TEST01 tests CFFTF.
!! TEST01 tests CFFTI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  integer seed
  real wsave(4*n+15)
  complex x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  For complex fast Fourier transforms,'
  write ( *, '(a)' ) '  CFFTI initializes the transforms,'
  write ( *, '(a)' ) '  CFFTF does a forward transforms;'
  write ( *, '(a)' ) '  CFFTB does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call cvec_random ( alo, ahi, n, x )

  call cvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call cffti ( n, wsave )
!
!  Compute the FFT coefficients.
!
  call cfftf ( n, x, wsave )

  call cvec_print_some ( n, x, 10, '  The FFT coefficients:' )
!
!  Now compute inverse FFT of coefficients.  Should get back the
!  original data.

  call cfftb ( n, x, wsave )

  x(1:n) = x(1:n) / real ( n )

  call cvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests CFFTB_2D.
!! TEST02 tests CFFTF_2D.
!
!  Plot the image and transform of an 8 by 8 unit source
!  in a 64 by 64 array.
!
  implicit none
!
  integer, parameter :: n = 64
  integer, parameter :: lda = n
!
  real a(n,n)
  real dat
  real ermax
  real err
  integer i
  complex image(lda,n)
  complex image2(lda,n)
  integer j
  real wsave(4*n+15)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  For complex fast Fourier transforms of 2D data:'
  write ( *, '(a)' ) '  CFFTF_2D computes the forward transform;'
  write ( *, '(a)' ) '  CFFTB_2D computes the backward transform.'
  write ( *, '(a)' ) ' '
!
!  Initialize the WSAVE array.
!
  call cffti ( n, wsave )
!
!  Set up the data.
!
!    IMAGE is the original data.
!
!    IMAGE2 is IMAGE scaled by (-1)**(I+J), to place the Fourier coefficients
!    in the correct place for viewing.
!
  do i = 1, n
    do j = 1, n

      if ( i >=(n/2-4) .and.  i<=(n/2+4) .and. &
           j>=(n/2-4) .and. j<=(n/2+4) ) then
        a(i,j) = 1.0E+00
      else
        a(i,j) = 0.0E+00
      end if

      image(i,j) = cmplx ( a(i,j), 0.0E+00 )
      image2(i,j) = image(i,j) * (-1)**(i-1+j-1)

    end do
  end do
!
!  Compute the forward Fourier transform of IMAGE and IMAGE2.
!
  call cfftf_2d ( lda, n, image, wsave )

  call cfftf_2d ( lda, n, image2, wsave )
!
!  Compute the magnitude of the components of transforms.
!  The actual transforms are unscaled and need to be divided by N*N
!  to be correct.
!
  a(1:n,1:n) = abs ( image(1:n,1:n) )
!
!  Compute the inverse Fourier transform of IMAGE.
!
  call cfftb_2d ( lda, n, image, wsave )
!
!  The transforms need to be divided by N*N to be correct.
!
  image(1:n,1:n) = image(1:n,1:n) / real ( n**2 )
!
!  See if the result agrees with the original data.
!
  ermax = 0.0E+00

  do i = 1, n
    do j = 1, n

      if ( i >=(n/2-4) .and. i <=(n/2+4) .and. &
           j >=(n/2-4) .and. j <=(n/2+4)) then
        dat = 1.0E+00
      else
        dat = 0.0E+00
      end if

      err = abs ( dat - abs ( image(i,j) ) )

      ermax = max ( ermax, err )

    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Maximum error in CFFT2D calculation:'
  write ( *, '(a)' ) ' '
  write ( *, '(g14.6)' ) ermax

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests COSQB.
!! TEST03 tests COSQF.
!! TEST03 tests COSQI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  integer seed
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  For real fast cosine quarter wave transforms,'
  write ( *, '(a)' ) '  COSQI initializes the transforms.'
  write ( *, '(a)' ) '  COSQF does a forward transform;'
  write ( *, '(a)' ) '  COSQB does a backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call cosqi ( n, wsave )
!
!  Compute the coefficients.
!
  call cosqf ( n, x, wsave )

  call rvec_print_some ( n, x, 10, '  The cosine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.
!
  call cosqb ( n, x, wsave )

  x(1:n) = x(1:n) / real ( 4 * n )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests RSQSTB.
!! TEST04 tests RSQSTF.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  integer seed
  real x(n)
  real y(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  For real slow quarter wave cosine transforms,'
  write ( *, '(a)' ) '  RSQCTF does a forward transform;'
  write ( *, '(a)' ) '  RSQCTB does a backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Compute the coefficients.
!
  call rsqctf ( n, x, y )

  call rvec_print_some ( n, y, 10, '  The cosine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call rsqctb ( n, y, x )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests EZFFTB.
!! TEST05 tests EZFFTF.
!! TEST05 tests EZFFTI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real a(n/2)
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real azero
  real b(n/2)
  real r_pi
  integer seed
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  For real fast Fourier transforms,'
  write ( *, '(a)' ) '  EZFFTI initializes the transforms.'
  write ( *, '(a)' ) '  EZFFTF does a forward transform;'
  write ( *, '(a)' ) '  EZFFTB does a backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call ezffti ( n, wsave )
!
!  Compute FFT
!
  call ezfftf ( n, x, azero, a, b, wsave )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The A0 coefficient:'
  write ( *, '(a)' ) ' '
  write ( *, '(g14.6)' ) azero

  call rvec_print_some ( n/2, a, 10, '  The A coefficients:' )

  call rvec_print_some ( n/2, b, 10, '  The B coefficients:' )
!
!  Now compute inverse FFT of coefficients.  Should get back the
!  original data.  First destroy original data so we're sure
!  that the routine had to recreate them!
!
  x(1:n) = 0.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Retrieve data from FFT coeficients.'

  call ezfftb ( n, x, azero, a, b, wsave )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )
  
  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests EZFFTB.
!! TEST06 tests EZFFTF.
!! TEST06 tests EZFFTI.
!
  implicit none
!
  integer, parameter :: n = 3087
!
  real a(n/2)
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real azero
  real b(n/2)
  real r_pi
  integer seed
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  For real fast Fourier transforms,'
  write ( *, '(a)' ) '  EZFFTI initializes the transforms.'
  write ( *, '(a)' ) '  EZFFTF does a forward transform;'
  write ( *, '(a)' ) '  EZFFTB does a backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
  write ( *, '(a)' ) '  which is not a multiple of two!'
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call ezffti ( n, wsave )
!
!  Compute FFT
!
  call ezfftf ( n, x, azero, a, b, wsave )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The A0 coefficient:'
  write ( *, '(a)' ) ' '
  write ( *, '(g14.6)' ) azero

  call rvec_print_some ( n/2, a, 10, '  The A coefficients:' )

  call rvec_print_some ( n/2, b, 10, '  The B coefficients:' )
!
!  Now compute inverse FFT of coefficients.  Should get back the
!  original data.  First destroy original data so we're sure
!  that the routine had to recreate them!
!
  x(1:n) = 0.0E+00

  call ezfftb ( n, x, azero, a, b, wsave )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test07
!
!*******************************************************************************
!
!! TEST07 tests EZFFTF.
!
  implicit none
!
  integer, parameter :: n = 36
!
  real a(n/2)
  real azero
  real b(n/2)
  real error
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  EZFFTF can take the Fourier transform of a real vector'
  write ( *, '(a)' ) '  of data.  In this case, the vector is'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (1,1,1,...,1)'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  and the transform should be'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (N,0,0,...,0),'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  where N is the number of entries, ', n

  x(1:n) = 1.0E+00
!
!  Initialize the WSAVE array.
!
  call ezffti ( n, wsave )
!
!  Compute the Fourier transform of the data.
!
  call ezfftf ( n, x, azero, a, b, wsave )
!
!  Compute the maximum error.
!
  error = max ( &
             abs ( 1.0E+00 - azero ), &
    maxval ( abs ( a(1:n/2) ) ), &
    maxval ( abs ( b(1:n/2) ) ) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The maximum error in computation is ', error

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08 tests EZFFTF.
!! TEST08 tests RSFTF.
!
  implicit none
!
  integer, parameter :: n = 36
!
  real a_f(n/2)
  real a_s(n/2)
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real azero_f
  real azero_s
  real b_f(n/2)
  real b_s(n/2)
  integer i
  integer seed
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  For real data,'
  write ( *, '(a)' ) '  EZFFTF takes the fast Fourier transform forward.'
  write ( *, '(a)' ) '  RSFTF computes the "slow" Fourier transform forward.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data values, N = ', n

  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )
!
!  Initialize the WSAVE array.
!
  call ezffti ( n, wsave )
!
!  Compute the fast Fourier transform of the data.
!
  call ezfftf ( n, x, azero_f, a_f, b_f, wsave )
!
!  Compute the slow Fourier transform of the data.
!
  call rsftf ( n, x, azero_s, a_s, b_s )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Fast    Slow'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A coefficients:'
  write ( *, '(a)' ) ' '

  write ( *, '(i3,2g14.6)' ) 0, azero_f, azero_s

  do i = 1, n/2
    write ( *, '(i3,2g14.6)' ) i, a_f(i), a_s(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  B coefficients:'
  write ( *, '(a)' ) ' '

  do i = 1, n/2
    write ( *, '(i3,2g14.6)' ) i, b_f(i), b_s(i)
  end do

  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09 tests EZFFTB.
!
!  The input (1,0,0,...,0) should produce the output (1,0,0,...,0).
!
  implicit none
!
  integer, parameter :: n = 36
!
  real a(n/2)
  real azero
  real b(n/2)
  real error
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  EZFFTB can be used to recover a real data vector'
  write ( *, '(a)' ) '  from a Fourier coefficient vector.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, the Fourier coefficient vector is:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (1,0,0,...,0)'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  and the recovered data vector should be'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (1,1,1,...,1).'

  azero = 1.0E+00
  a(1:n/2) = 0.0E+00
  b(1:n/2) = 0.0E+00
!
!  Initialize the WSAVE array.
!
  call ezffti ( n, wsave )
!
!  Compute the inverse Fourier transform.
!
  call ezfftb ( n, x, azero, a, b, wsave )
!
!  Compute the maximum error.
!
  error = maxval ( abs ( x(1:n) - 1.0E+00 ) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The maximum error in the computation was ', error

  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 tests RFFTF.
!
!  The input vector is (1,1,1,...,1) which should produce
!  the output (N,0,0,...,0).
!
  implicit none
!
  integer, parameter :: n = 36
!
  real error
  real wsave(2*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  RFFTF can compute the Fourier transform of a real'
  write ( *, '(a)' ) '  vector of data.  In this case, the vector is'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (1,1,1,...,1)'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  and the transform should be'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (N,0,0,...,0), '
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  where N is the number of entries, N = ', n

  x(1:n) = 1.0E+00
!
!  Initialize the WSAVE array.
!
  call rffti ( n, wsave )
!
!  Compute the Fourier transform.
!
  call rfftf ( n, x, wsave )
!
!  Test results.
!
  error = max ( &
    abs ( real ( n ) - x(1) ), & 
    maxval ( abs ( x(2:n) ) ) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The maximum error in computation is ', error

  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11 tests RFFTB.
!
!  The input vector is (1,0,0,...,0) which should produce
!  the output (1,0,0,...,0).
!
  implicit none
!
  integer, parameter :: n = 36
!
  real error
  real wsave(2*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  RFFTB can recover a real vector of data from Fourier'
  write ( *, '(a)' ) '  coefficients.  In this case, the coefficients are:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (1,0,0,...,0)'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  and the data should be:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    (1,1,1,...,1).'

  x(1) = 1.0E+00
  x(2:n) = 0.0E+00
!
!  Initialize the WSAVE array.
!
  call rffti ( n, wsave )
!
!  Compute the inverse Fourier transform.
!
  call rfftb ( n, x, wsave )
!
!  Compute the maximum error.
!
  error = maxval ( abs ( x(1:n) - 1.0E+00 ) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The maximum error in computation is ', error

  return
end
subroutine test12
!
!*******************************************************************************
!
!! TEST12 tests RSFTB.
!! TEST12 tests RSFTF.
!
  implicit none
!
  integer, parameter :: n = 36
!
  real a(n/2)
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real azero
  real b(n/2)
  integer i
  integer seed
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  For real slow Fourier transforms,'
  write ( *, '(a)' ) '  RSFTF computes the forward transform.'
  write ( *, '(a)' ) '  RSFTB computes the backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data values, N = ', n

  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Compute the slow Fourier transform of the data.
!
  call rsftf ( n, x, azero, a, b )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A coefficients:'
  write ( *, '(a)' ) ' '

  write ( *, '(i3,g14.6)' ) 0, azero

  do i = 1, n/2
    write ( *, '(i3,g14.6)' ) i, a(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  B coefficients:'
  write ( *, '(a)' ) ' '

  do i = 1, n/2
    write ( *, '(i3,g14.6)' ) i, b(i)
  end do

  call rsftb ( n, x, azero, a, b )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test13
!
!*******************************************************************************
!
!! TEST13 tests SINQB.
!! TEST13 tests SINQF.
!! TEST13 tests SINQI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  integer seed
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  For real fast sine quarter wave transforms,'
  write ( *, '(a)' ) '  SINQI initializes the transforms;'
  write ( *, '(a)' ) '  SINQF does a forward transform;'
  write ( *, '(a)' ) '  SINQB does a backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call sinqi ( n, wsave )
!
!  Compute the coefficients.
!
  call sinqf ( n, x, wsave )

  call rvec_print_some ( n, x, 10, '  The sine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call sinqb ( n, x, wsave )

  x(1:n) = x(1:n) / real ( 4 * n )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test14
!
!*******************************************************************************
!
!! TEST14 tests RSQSTB.
!! TEST14 tests RSQSTF.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  integer seed
  real x(n)
  real y(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST14'
  write ( *, '(a)' ) '  For real slow quarter wave sine transforms,'
  write ( *, '(a)' ) '  RSQSTF does a forward transform;'
  write ( *, '(a)' ) '  RSQSTB does a backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Compute the coefficients.
!
  call rsqstf ( n, x, y )

  call rvec_print_some ( n, y, 10, '  The sine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call rsqstb ( n, y, x )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test15
!
!*******************************************************************************
!
!! TEST15 tests RSINT.
!! TEST15 tests RSINTI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  integer seed
  real wsave((5*n+30)/2)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST15'
  write ( *, '(a)' ) '  For real fast sine transforms,'
  write ( *, '(a)' ) '  RSINTI initializes the transforms.'
  write ( *, '(a)' ) '  RSINT does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call rsinti ( n, wsave )
!
!  Compute the coefficients.
!
  call rsint ( n, x, wsave )

  call rvec_print_some ( n, x, 10, '  The sine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call rsint ( n, x, wsave )

  x(1:n) = x(1:n) / real ( 2 * ( n + 1 ) )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test16
!
!*******************************************************************************
!
!! TEST16 tests DSINT.
!! TEST16 tests DSINTI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  double precision, parameter :: ahi = 5.0D+00
  double precision, parameter :: alo = 0.0D+00
  integer seed
  double precision wsave((5*n+30)/2)
  double precision x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST16'
  write ( *, '(a)' ) '  For double precision fast sine transforms,'
  write ( *, '(a)' ) '  DSINTI initializes the transforms.'
  write ( *, '(a)' ) '  DSINT does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call dvec_random ( alo, ahi, n, x )

  call dvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call dsinti ( n, wsave )
!
!  Compute the coefficients.
!
  call dsint ( n, x, wsave )

  call dvec_print_some ( n, x, 10, '  The sine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call dsint ( n, x, wsave )

  x(1:n) = x(1:n) / dble ( 2 * ( n + 1 ) )

  call dvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test17
!
!*******************************************************************************
!
!! TEST17 tests RSST.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real c(n)
  real d(n)
  real e(n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST17'
  write ( *, '(a)' ) '  For real slow sine transforms,'
  write ( *, '(a)' ) '  RSST does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, c )

  call rvec_print_some ( n, c, 10, '  The original data:' )
!
!  Compute the coefficients.
!
  call rsst ( n, c, d )

  call rvec_print_some ( n, d, 10, '  The sine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call rsst ( n, d, e )

  e(1:n) = e(1:n) / real ( 2 * ( n + 1 ) )

  call rvec_print_some ( n, e, 10, '  The retrieved data:' )

  return
end
subroutine test18
!
!*******************************************************************************
!
!! TEST18 tests DSST.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  double precision, parameter :: ahi = 5.0D+00
  double precision, parameter :: alo = 0.0D+00
  double precision c(n)
  double precision d(n)
  double precision e(n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST18'
  write ( *, '(a)' ) '  For double precision slow sine transforms,'
  write ( *, '(a)' ) '  DSST does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call dvec_random ( alo, ahi, n, c )

  call dvec_print_some ( n, c, 10, '  The original data:' )
!
!  Compute the coefficients.
!
  call dsst ( n, c, d )

  call dvec_print_some ( n, d, 10, '  The sine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call dsst ( n, d, e )

  e(1:n) = e(1:n) / dble ( 2 * ( n + 1 ) )

  call dvec_print_some ( n, e, 10, '  The retrieved data:' )

  return
end
subroutine test19
!
!*******************************************************************************
!
!! TEST19 tests RCOST.
!! TEST19 tests RCOSTI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  integer seed
  real wsave(3*n+15)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST19'
  write ( *, '(a)' ) '  For real fast cosine transforms,'
  write ( *, '(a)' ) '  RCOSTI initializes the transforms.'
  write ( *, '(a)' ) '  RCOST does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, x )

  call rvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call rcosti ( n, wsave )
!
!  Compute the coefficients.
!
  call rcost ( n, x, wsave )

  call rvec_print_some ( n, x, 10, '  The cosine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.
!
  call rcost ( n, x, wsave )

  x(1:n) = x(1:n) / real ( 2 * ( n - 1 ) )

  call rvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test20
!
!*******************************************************************************
!
!! TEST20 tests DCOST.
!! TEST20 tests DCOSTI.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  double precision, parameter :: ahi = 5.0D+00
  double precision, parameter :: alo = 0.0D+00
  integer seed
  double precision wsave(3*n+15)
  double precision x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST20'
  write ( *, '(a)' ) '  For double precision fast cosine transforms,'
  write ( *, '(a)' ) '  DCOSTI initializes the transforms.'
  write ( *, '(a)' ) '  DCOST does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call dvec_random ( alo, ahi, n, x )

  call dvec_print_some ( n, x, 10, '  The original data:' )
!
!  Initialize the WSAVE array.
!
  call dcosti ( n, wsave )
!
!  Compute the coefficients.
!
  call dcost ( n, x, wsave )

  call dvec_print_some ( n, x, 10, '  The cosine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.
!
  call dcost ( n, x, wsave )

  x(1:n) = x(1:n) / dble ( 2 * ( n - 1 ) )

  call dvec_print_some ( n, x, 10, '  The retrieved data:' )

  return
end
subroutine test21
!
!*******************************************************************************
!
!! TEST21 tests RSCT.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real c(n)
  real d(n)
  real e(n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST21'
  write ( *, '(a)' ) '  For real slow cosine transforms,'
  write ( *, '(a)' ) '  RSCT does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, c )

  call rvec_print_some ( n, c, 10, '  The original data:' )
!
!  Compute the coefficients.
!
  call rsct ( n, c, d )

  call rvec_print_some ( n, d, 10, '  The cosine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call rsct ( n, d, e )

  e(1:n) = e(1:n) / real ( 2 * ( n + 1 ) )

  call rvec_print_some ( n, e, 10, '  The retrieved data:' )

  return
end
subroutine test22
!
!*******************************************************************************
!
!! TEST22 tests DSCT.
!
  implicit none
!
  integer, parameter :: n = 4096
!
  double precision, parameter :: ahi = 5.0D+00
  double precision, parameter :: alo = 0.0D+00
  double precision c(n)
  double precision d(n)
  double precision e(n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST22'
  write ( *, '(a)' ) '  For double precision slow cosine transforms,'
  write ( *, '(a)' ) '  DSCT does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call dvec_random ( alo, ahi, n, c )

  call dvec_print_some ( n, c, 10, '  The original data:' )
!
!  Compute the coefficients.
!
  call dsct ( n, c, d )

  call dvec_print_some ( n, d, 10, '  The cosine coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call dsct ( n, d, e )

  e(1:n) = e(1:n) / real ( 2 * ( n + 1 ) )

  call dvec_print_some ( n, e, 10, '  The retrieved data:' )

  return
end
subroutine test23
!
!*******************************************************************************
!
!! TEST23 tests RSHT.
!
  implicit none
!
  integer, parameter :: n = 17
!
  real, parameter :: ahi = 5.0E+00
  real, parameter :: alo = 0.0E+00
  real c(n)
  real d(n)
  real e(n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST23'
  write ( *, '(a)' ) '  For real slow Hartley transforms,'
  write ( *, '(a)' ) '  RSHT does a forward or backward transform.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973
  call random_initialize ( seed )

  call rvec_random ( alo, ahi, n, c )

  call rvec_print_some ( n, c, 10, '  The original data:' )
!
!  Compute the coefficients.
!
  call rsht ( n, c, d )

  call rvec_print_some ( n, d, 10, '  The Hartley coefficients:' )
!
!  Now compute inverse transform of coefficients.  Should get back the
!  original data.

  call rsht ( n, d, e )

  call rvec_print_some ( n, e, 10, '  The retrieved data:' )

  return
end