program main

!*******************************************************************************
!
!! FFTPACK5_PRB calls the FFTPACK5 test routines.
!
  implicit none

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FFTPACK5_PRB'
  write ( *, '(a)' ) '  A set of tests for FFTPACK5.'

  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

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

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

  stop
end
subroutine test01

!*******************************************************************************
!
!! TEST01 tests CFFT1B, CFFT1F and CFFT1I.
!
  implicit none

  integer, parameter :: n = 4096

  integer, parameter :: lenwrk = 2 * n

  complex c(n)
  integer ier
  integer inc
  integer lenc
  integer lensav
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  For complex fast Fourier transforms, 1D,'
  write ( *, '(a)' ) '  CFFT1I initializes the transforms,'
  write ( *, '(a)' ) '  CFFT1F does a forward transforms;'
  write ( *, '(a)' ) '  CFFT1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call cvec_uniform_01 ( n, seed, c )

  call cvec_print_some ( n, c, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n * int ( log ( real ( n ) ) ) + 4
  allocate ( wsave(1:lensav) )

  call cfft1i ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  inc = 1
  lenc = n

  call cfft1f ( n, inc, c, lenc, wsave, lensav, work, lenwrk, ier )

  call cvec_print_some ( n, c, 10, '  The FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call cfft1b ( n, inc, c, lenc, wsave, lensav, work, lenwrk, ier )

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

  deallocate ( wsave )

  return
end
subroutine test02

!*******************************************************************************
!
!! TEST02 tests CFFT2B, CFFT2F and CFFT2I.
!
  implicit none

  integer, parameter :: l = 32
  integer, parameter :: m = 64

  integer, parameter :: lenwrk = 2 * l * m

  complex c(l,m)
  integer ier
  integer ldim
  integer lensav
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  For complex fast Fourier transforms, 2D,'
  write ( *, '(a)' ) '  CFFT2I initializes the transforms,'
  write ( *, '(a)' ) '  CFFT2F does a forward transforms;'
  write ( *, '(a)' ) '  CFFT2B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The data is stored in an L by M array, with'
  write ( *, '(a,i4)' ) '  L = ', l
  write ( *, '(a,i4)' ) '  M = ', m
!
!  Set the data values.
!
  seed = 1973

  call cmat_uniform_01 ( l, m, seed, c )

  call cmat_print_some ( l, m, c, 1, 1, 5, 5, '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * ( l + m ) + int ( log ( real ( l ) ) ) &
    + int ( log ( real ( m ) ) ) + 8

  allocate ( wsave(1:lensav) )

  call cfft2i ( l, m, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  ldim = l

  call cfft2f ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier )

  call cmat_print_some ( l, m, c, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call cfft2b ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier )

  call cmat_print_some ( l, m, c, 1, 1, 5, 5, '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test03

!*******************************************************************************
!
!! TEST03 tests CFFTMB, CFFTMF and CFFTMI.
!
  implicit none

  integer, parameter :: n = 32
  integer, parameter :: lot = 6
  integer, parameter :: lenc = n * lot

  integer, parameter :: lenwrk = 2 * lot * n

  complex c(lenc)
  integer ier
  integer inc
  integer jump
  integer lensav
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  For complex fast Fourier transforms, 1D, multiple'
  write ( *, '(a)' ) '  CFFTMI initializes the transforms,'
  write ( *, '(a)' ) '  CFFTMF does a forward transforms;'
  write ( *, '(a)' ) '  CFFTMB does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of sequences is LOT = ', lot
  write ( *, '(a,i4)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call cmat_uniform_01 ( n, lot, seed, c )

  call cmat_print_some ( n, lot, c, 1, 1, 5, 5, '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call cfftmi ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  jump = n
  inc = 1

  call cfftmf ( lot, jump, n, inc, c, lenc, wsave, lensav, work, lenwrk, ier )

  call cmat_print_some ( n, lot, c, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call cfftmb ( lot, jump, n, inc, c, lenc, wsave, lensav, work, lenwrk, ier )

  call cmat_print_some ( n, lot, c, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test04

!*******************************************************************************
!
!! TEST04 tests RFFT1B, RFFT1F and RFFT1I.
!
  implicit none

  integer, parameter :: n = 4096

  integer, parameter :: lenwrk = 2 * n

  integer ier
  integer inc
  integer lenr
  integer lensav
  real r(n)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  For real fast Fourier transforms, 1D,'
  write ( *, '(a)' ) '  RFFT1I initializes the transforms,'
  write ( *, '(a)' ) '  RFFT1F does a forward transforms;'
  write ( *, '(a)' ) '  RFFT1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rvec_uniform_01 ( n, seed, r )

  call rvec_print_some ( n, r, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = n * int ( log ( real ( n ) ) ) + 4
  allocate ( wsave(1:lensav) )

  call rfft1i ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  inc = 1
  lenr = n

  call rfft1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rvec_print_some ( n, r, 10, '  The FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call rfft1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

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

  deallocate ( wsave )

  return
end
subroutine test05

!*******************************************************************************
!
!! TEST05 tests RFFT2B, RFFT2F and RFFT2I.
!
  implicit none

  integer, parameter :: l = 32
  integer, parameter :: ldim = 2 * ( l / 2 + 1 )
  integer, parameter :: m = 64

  integer, parameter :: lenwrk = 2 * ldim * m

  integer ier
  integer lensav
  real r(ldim,m)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  For real fast Fourier transforms, 2D,'
  write ( *, '(a)' ) '  RFFT2I initializes the transforms,'
  write ( *, '(a)' ) '  RFFT2F does a forward transforms;'
  write ( *, '(a)' ) '  RFFT2B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The L by M data is stored in an LDIM by M array, with'
  write ( *, '(a,i4)' ) '  L = ', l
  write ( *, '(a,i4)' ) '  LDIM = ', ldim
  write ( *, '(a,i4)' ) '  M = ', m
!
!  Set the data values.
!
  seed = 1973

  call rmat_uniform_01 ( ldim, m, seed, r )

  call rmat_print_some ( ldim, m, r, 1, 1, 5, 5, &
    '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * ( l + m ) + int ( log ( real ( l ) ) ) &
    + int ( log ( real ( m ) ) ) + 8

  allocate ( wsave(1:lensav) )

  call rfft2i ( l, m, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  call rfft2f ( ldim, l, m, r, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( ldim, m, r, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call rfft2b ( ldim, l, m, r, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( ldim, m, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test06

!*******************************************************************************
!
!! TEST06 tests RFFTMB, RFFTMF and RFFTMI.
!
  implicit none

  integer, parameter :: n = 32
  integer, parameter :: lot = 6
  integer, parameter :: lenr = n * lot

  integer, parameter :: lenwrk = lot * n

  integer ier
  integer inc
  integer jump
  integer lensav
  real r(lenr)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  For real fast Fourier transforms, 1D, multiple'
  write ( *, '(a)' ) '  RFFTMI initializes the transforms,'
  write ( *, '(a)' ) '  RFFTMF does a forward transforms;'
  write ( *, '(a)' ) '  RFFTMB does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of sequences is LOT = ', lot
  write ( *, '(a,i4)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rmat_uniform_01 ( n, lot, seed, r )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call rfftmi ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  jump = n
  inc = 1

  call rfftmf ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call rfftmb ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test07

!*******************************************************************************
!
!! TEST07 tests COST1B, COST1F and COST1I.
!
  implicit none

  integer, parameter :: n = 4096

  integer, parameter :: lenwrk = n - 1

  integer ier
  integer inc
  integer lenr
  integer lensav
  real r(n)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  For real fast cosine transforms, 1D,'
  write ( *, '(a)' ) '  COST1I initializes the transforms,'
  write ( *, '(a)' ) '  COST1F does a forward transforms;'
  write ( *, '(a)' ) '  COST1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rvec_uniform_01 ( n, seed, r )

  call rvec_print_some ( n, r, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call cost1i ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  inc = 1
  lenr = n

  call cost1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rvec_print_some ( n, r, 10, '  The FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call cost1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

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

  deallocate ( wsave )

  return
end
subroutine test08

!*******************************************************************************
!
!! TEST08 tests COSTMB, COSTMF and COSTMI.
!
  implicit none

  integer, parameter :: n = 32
  integer, parameter :: lot = 6
  integer, parameter :: lenr = n * lot

  integer, parameter :: lenwrk = lot * ( n + 1 )

  integer ier
  integer inc
  integer jump
  integer lensav
  real r(lenr)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  For real fast cosine transforms, 1D, multiple'
  write ( *, '(a)' ) '  COSTMI initializes the transforms,'
  write ( *, '(a)' ) '  COSTMF does a forward transforms;'
  write ( *, '(a)' ) '  COSTMB does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of sequences is LOT = ', lot
  write ( *, '(a,i4)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rmat_uniform_01 ( n, lot, seed, r )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call costmi ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  jump = n
  inc = 1

  call costmf ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call costmb ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test09

!*******************************************************************************
!
!! TEST09 tests SINT1B, SINT1F and SINT1I.
!
  implicit none

  integer, parameter :: n = 4096

  integer, parameter :: lenwrk = 2 * ( n + 1 )

  integer ier
  integer inc
  integer lenr
  integer lensav
  real r(n)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  For real fast sine transforms, 1D,'
  write ( *, '(a)' ) '  SINT1I initializes the transforms,'
  write ( *, '(a)' ) '  SINT1F does a forward transforms;'
  write ( *, '(a)' ) '  SINT1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rvec_uniform_01 ( n, seed, r )

  call rvec_print_some ( n, r, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = n/2 + n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call sint1i ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  inc = 1
  lenr = n

  call sint1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rvec_print_some ( n, r, 10, '  The FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call sint1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

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

  deallocate ( wsave )

  return
end
subroutine test10

!*******************************************************************************
!
!! TEST10 tests SINTMB, SINTMF and SINTMI.
!
  implicit none

  integer, parameter :: n = 32
  integer, parameter :: lot = 6
  integer, parameter :: lenr = n * lot

  integer, parameter :: lenwrk = lot * 2 * ( n + 2 )

  integer ier
  integer inc
  integer jump
  integer lensav
  real r(lenr)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  For real fast sine transforms, 1D, multiple'
  write ( *, '(a)' ) '  SINTMI initializes the transforms,'
  write ( *, '(a)' ) '  SINTMF does a forward transforms;'
  write ( *, '(a)' ) '  SINTMB does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of sequences is LOT = ', lot
  write ( *, '(a,i4)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rmat_uniform_01 ( n, lot, seed, r )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = n / 2 + n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call sintmi ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  jump = n
  inc = 1

  call sintmf ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call sintmb ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test11

!*******************************************************************************
!
!! TEST11 tests COSQ1B, COSQ1F and COSQ1I.
!
  implicit none

  integer, parameter :: n = 4096

  integer, parameter :: lenwrk = n

  integer ier
  integer inc
  integer lenr
  integer lensav
  real r(n)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  For real fast cosine transforms, 1D,'
  write ( *, '(a)' ) '  COSQ1I initializes the transforms,'
  write ( *, '(a)' ) '  COSQ1F does a forward transforms;'
  write ( *, '(a)' ) '  COSQ1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rvec_uniform_01 ( n, seed, r )

  call rvec_print_some ( n, r, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call cosq1i ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  inc = 1
  lenr = n

  call cosq1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rvec_print_some ( n, r, 10, '  The FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call cosq1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

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

  deallocate ( wsave )

  return
end
subroutine test12

!*******************************************************************************
!
!! TEST12 tests COSQMB, COSQMF and COSQMI.
!
  implicit none

  integer, parameter :: n = 32
  integer, parameter :: lot = 6
  integer, parameter :: lenr = n * lot

  integer, parameter :: lenwrk = lot * n

  integer ier
  integer inc
  integer jump
  integer lensav
  real r(lenr)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  For real fast cosine transforms, 1D, multiple'
  write ( *, '(a)' ) '  COSQMI initializes the transforms,'
  write ( *, '(a)' ) '  COSQMF does a forward transforms;'
  write ( *, '(a)' ) '  COSQMB does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of sequences is LOT = ', lot
  write ( *, '(a,i4)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rmat_uniform_01 ( n, lot, seed, r )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call cosqmi ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  jump = n
  inc = 1

  call cosqmf ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call cosqmb ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test13

!*******************************************************************************
!
!! TEST13 tests SINQ1B, SINQ1F and SINQ1I.
!
  implicit none

  integer, parameter :: n = 4096

  integer, parameter :: lenwrk = n

  integer ier
  integer inc
  integer lenr
  integer lensav
  real r(n)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  For real fast sine transforms, 1D,'
  write ( *, '(a)' ) '  SINQ1I initializes the transforms,'
  write ( *, '(a)' ) '  SINQ1F does a forward transforms;'
  write ( *, '(a)' ) '  SINQ1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rvec_uniform_01 ( n, seed, r )

  call rvec_print_some ( n, r, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call sinq1i ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  inc = 1
  lenr = n

  call sinq1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rvec_print_some ( n, r, 10, '  The FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call sinq1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

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

  deallocate ( wsave )

  return
end
subroutine test14

!*******************************************************************************
!
!! TEST14 tests SINQMB, SINQMF and SINQMI.
!
  implicit none

  integer, parameter :: n = 32
  integer, parameter :: lot = 6
  integer, parameter :: lenr = n * lot

  integer, parameter :: lenwrk = lot * n

  integer ier
  integer inc
  integer jump
  integer lensav
  real r(lenr)
  integer seed
  real work(lenwrk)
  real, allocatable, dimension ( : ) :: wsave

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST14'
  write ( *, '(a)' ) '  For real fast sine transforms, 1D, multiple'
  write ( *, '(a)' ) '  SINQMI initializes the transforms,'
  write ( *, '(a)' ) '  SINQMF does a forward transforms;'
  write ( *, '(a)' ) '  SINQMB does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The number of sequences is LOT = ', lot
  write ( *, '(a,i4)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call rmat_uniform_01 ( n, lot, seed, r )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, '  Part of the original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n ) ) ) + 4

  allocate ( wsave(1:lensav) )

  call sinqmi ( n, wsave, lensav, ier )
!
!  Compute the FFT coefficients.
!
  jump = n
  inc = 1

  call sinqmf ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the FFT coefficients:' )
!
!  Compute inverse FFT of coefficients.  Should get back the
!  original data.
!
  call sinqmb ( lot, jump, n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

  call rmat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

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

!*******************************************************************************
!
!! CVEC_PRINT_SOME prints some of a complex vector.
!
!  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 ( 0 < len_trim ( title ) ) 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 ( 3 <= max_print ) 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 cvec_uniform_01 ( n, seed, c )

!*******************************************************************************
!
!! CVEC_UNIFORM_01 returns a unit double precision complex pseudorandom vector.
!
!  Discussion:
!
!    The angles should be uniformly distributed between 0 and 2 * PI,
!    the square roots of the radius uniformly distributed between 0 and 1.
!
!    This results in a uniform distribution of values in the unit circle.
!
!  Modified:
!
!    15 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of values to compute.
!
!    Input/output, integer SEED, the "seed" value, which should NOT be 0.
!    On output, SEED has been updated.
!
!    Output, complex C(N), the pseudorandom complex vector.
!
  implicit none

  integer n

  complex c(n)
  integer i
  real r
  integer k
  real, parameter :: pi = 3.1415926E+00
  integer seed
  real theta

  do i = 1, n

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
      seed = seed + 2147483647
    end if

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

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
      seed = seed + 2147483647
    end if

    theta = 2.0E+00 * pi * real ( real ( seed, kind = 8 ) * 4.656612875D-10 )

    c(i) = r * cmplx ( cos ( theta ), sin ( theta ) )

  end do

  return
end
subroutine timestamp ( )

!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    26 February 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

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

!*******************************************************************************
!
!! CMAT_PRINT prints a CMAT.
!
!  Modified:
!
!    23 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, N, the number of rows and columns in the matrix.
!
!    Input, complex A(M,N), the matrix.
!
!    Input, character ( len = * ) TITLE, a title to be printed.
!
  implicit none

  integer m
  integer n

  complex a(m,n)
  character ( len = * ) title

  call cmat_print_some ( m, n, a, 1, 1, m, n, title )

  return
end
subroutine cmat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

!*******************************************************************************
!
!! CMAT_PRINT_SOME prints some of a CMAT.
!
!  Modified:
!
!    23 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, N, the number of rows and columns in the matrix.
!
!    Input, complex A(M,N), the matrix.
!
!    Input, integer ILO, JLO, IHI, JHI, the first row and
!    column, and the last row and column to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed.
!
  implicit none

  integer, parameter :: incx = 4
  integer m
  integer n

  complex a(m,n)
  character ( len = 20 ) 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
  complex zero

  zero = cmplx ( 0.0E+00, 0.0E+00 )

  if ( 0 < len_trim ( title ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if
!
!  Print the columns of the matrix, in strips of INCX.
!
  do j2lo = jlo, jhi, 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), '(i10,10x)' ) j
    end do

    write ( *, '(a,4a20)' ) '  Col: ', ( ctemp(j2), j2 = 1, inc )
    write ( *, '(a)' ) '  Row'
    write ( *, '(a)' ) '  ---'
!
!  Determine the range of the rows in this strip.
!
    i2lo = max ( ilo, 1 )
    i2hi = min ( ihi, m )

    do i = i2lo, i2hi
!
!  Print out (up to) INCX entries in row I, that lie in the current strip.
!
      do j2 = 1, inc

        j = j2lo - 1 + j2

        if ( a(i,j) == zero ) then
          ctemp(j2) = '    0.0'
        else if ( aimag ( a(i,j) ) == 0.0E+00 ) then
          write ( ctemp(j2), '(g10.3,10x)' ) real ( a(i,j) )
        else
          write ( ctemp(j2), '(2g10.3)' ) a(i,j)
        end if

      end do

      write ( *, '(i5,1x,4a20)' ) i, ( ctemp(j2), j2 = 1, inc )

    end do

  end do

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

  return
end
subroutine cmat_uniform_01 ( m, n, seed, c )

!*******************************************************************************
!
!! CMAT_UNIFORM_01 returns a unit complex pseudorandom matrix.
!
!  Discussion:
!
!    The angles should be uniformly distributed between 0 and 2 * PI,
!    the square roots of the radius uniformly distributed between 0 and 1.
!
!    This results in a uniform distribution of values in the unit circle.
!
!  Modified:
!
!    15 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, N, the number of rows and columns in the matrix.
!
!    Input/output, integer SEED, the "seed" value, which should NOT be 0.
!    On output, SEED has been updated.
!
!    Output, complex C(M,N), the pseudorandom complex matrix.
!
  implicit none

  integer m
  integer n

  complex c(m,n)
  integer i
  integer j
  real r
  integer k
  real, parameter :: pi = 3.1415926E+00
  integer seed
  real theta

  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 = sqrt ( real ( real ( seed, kind = 8 ) * 4.656612875D-10 ) )

      k = seed / 127773

      seed = 16807 * ( seed - k * 127773 ) - k * 2836

      if ( seed < 0 ) then
        seed = seed + 2147483647
      end if

      theta = 2.0D+00 * pi * real ( real ( seed, kind = 8 ) * 4.656612875D-10 )

      c(i,j) = r * cmplx ( cos ( theta ), sin ( theta ) )

    end do

  end do

  return
end
subroutine rvec_uniform_01 ( n, seed, r )

!*******************************************************************************
!
!! RVEC_UNIFORM_01 sets a real vector to unit pseudorandom numbers.
!
!  Modified:
!
!    22 March 2005
!
!  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 N, the number of entries in the vector.
!
!    Input/output, integer SEED, the "seed" value, which should NOT be 0.
!    On output, SEED has been updated.
!
!    Output, real R(N), the vector of pseudorandom values.
!
  implicit none

  integer n

  integer i
  integer k
  integer seed
  real r(n)

  do i = 1, n

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
      seed = seed + 2147483647
    end if

    r(i) = real ( seed ) * 4.656612875E-10

  end do

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

!*******************************************************************************
!
!! RVEC_PRINT_SOME prints "some" of a DVEC.
!
!  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:
!
!    19 December 2001
!
!  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 ( 0 < len_trim ( title ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
    write ( *, '(a)' ) ' '
  end if

  if ( n <= max_print ) then

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

  else if ( 3 <= max_print ) then

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

  else

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

  end if

  return
end
subroutine rmat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

!*******************************************************************************
!
!! RMAT_PRINT_SOME prints some of a real matrix.
!
!  Modified:
!
!    12 September 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, N, the number of rows and columns.
!
!    Input, real 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 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

        if ( a(i,j) == real ( int ( a(i,j) ) ) ) then
          write ( ctemp(j2), '(f8.0,6x)' ) a(i,j)
        else
          write ( ctemp(j2), '(g14.6)' ) a(i,j)
        end if

      end do

      write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j), j = 1, inc )

    end do

  end do

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

  return
end
subroutine rmat_uniform_01 ( m, n, seed, r )

!*******************************************************************************
!
!! RMAT_UNIFORM_01 fills an array with unit pseudorandom numbers.
!
!  Modified:
!
!    26 March 2005
!
!  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 R(M,N), the array of pseudorandom values.
!
  implicit none

  integer m
  integer n

  integer i
  integer j
  integer k
  integer seed
  real 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 ) * 4.656612875E-10

    end do
  end do

  return
end