program main

!*****************************************************************************80
!
!! FFTPACK5_PRB calls the FFTPACK5 test routines.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )

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

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

  call test10
  call test11
  call test12
  call test13
  call test14
  call test15

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

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

  stop
end
subroutine test01

!*****************************************************************************80
!
!! TEST01 tests CFFT1B, CFFT1F and CFFT1I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

  integer ( kind = 4 ), parameter :: lenwrk = 2 * n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call c4vec_uniform_01 ( n, seed, c )

  call c4vec_print_some ( n, c, 1, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n * int ( log ( real ( n, kind = 4 ) ) ) + 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 c4vec_print_some ( n, c, 1, 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 c4vec_print_some ( n, c, 1, 10, '  The retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test02

!*****************************************************************************80
!
!! TEST02 tests CFFT2B, CFFT2F and CFFT2I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: l = 32
  integer ( kind = 4 ), parameter :: m = 64

  integer ( kind = 4 ), parameter :: lenwrk = 2 * l * m

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  L = ', l
  write ( *, '(a,i8)' ) '  M = ', m
!
!  Set the data values.
!
  seed = 1973

  call c4mat_uniform_01 ( l, m, seed, c )

  call c4mat_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, kind = 4 ) ) ) &
    + 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 c4mat_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 c4mat_print_some ( l, m, c, 1, 1, 5, 5, '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test03

!*****************************************************************************80
!
!! TEST03 tests CFFTMB, CFFTMF and CFFTMI.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 32
  integer ( kind = 4 ), parameter :: lot = 6
  integer ( kind = 4 ), parameter :: lenc = n * lot

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

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of sequences is LOT =   ', lot
  write ( *, '(a,i8)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call c4mat_uniform_01 ( n, lot, seed, c )

  call c4mat_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, kind = 4 ) ) ) + 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 c4mat_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 c4mat_print_some ( n, lot, c, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test04

!*****************************************************************************80
!
!! TEST04 tests COSQ1B, COSQ1F and COSQ1I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

  integer ( kind = 4 ), parameter :: lenwrk = n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4vec_uniform_01 ( n, seed, r )

  call r4vec_print_some ( n, r, 1, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n, kind = 4 ) ) ) + 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 r4vec_print_some ( n, r, 1, 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 r4vec_print_some ( n, r, 1, 10, '  The retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test05

!*****************************************************************************80
!
!! TEST05 tests COSQMB, COSQMF and COSQMI.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 32
  integer ( kind = 4 ), parameter :: lot = 6
  integer ( kind = 4 ), parameter :: lenr = n * lot

  integer ( kind = 4 ), parameter :: lenwrk = lot * n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of sequences is LOT =   ', lot
  write ( *, '(a,i8)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4mat_uniform_01 ( n, lot, seed, r )

  call r4mat_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, kind = 4 ) ) ) + 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 r4mat_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 r4mat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test06

!*****************************************************************************80
!
!! TEST06 tests COST1B, COST1F and COST1I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

  integer ( kind = 4 ), parameter :: lenwrk = n - 1

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4vec_uniform_01 ( n, seed, r )

  call r4vec_print_some ( n, r, 1, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n, kind = 4 ) ) ) + 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 r4vec_print_some ( n, r, 1, 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 r4vec_print_some ( n, r, 1, 10, '  The retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test07

!*****************************************************************************80
!
!! TEST07 tests COSTMB, COSTMF and COSTMI.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 32
  integer ( kind = 4 ), parameter :: lot = 6
  integer ( kind = 4 ), parameter :: lenr = n * lot

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

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of sequences is LOT =   ', lot
  write ( *, '(a,i8)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4mat_uniform_01 ( n, lot, seed, r )

  call r4mat_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, kind = 4 ) ) ) + 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 r4mat_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 r4mat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test075

!*****************************************************************************80
!
!! TEST075 tests DCOSQ1B, DCOSQ1F and DCOSQ1I.
!
!  Modified:
!
!    17 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

  integer ( kind = 4 ), parameter :: lenwrk = n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST075'
  write ( *, '(a)' ) '  For double precision real fast cosine transforms, 1D,'
  write ( *, '(a)' ) '  DCOSQ1I initializes the transforms,'
  write ( *, '(a)' ) '  DCOSQ1F does a forward transforms;'
  write ( *, '(a)' ) '  DCOSQ1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r8vec_uniform_01 ( n, seed, r )

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

  allocate ( wsave(1:lensav) )

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

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

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

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

  deallocate ( wsave )

  return
end
subroutine test08

!*****************************************************************************80
!
!! TEST08 tests DFFT1B, DFFT1F and DFFT1I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

  integer ( kind = 4 ), parameter :: lenwrk = 2 * n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  For double precision real fast Fourier transforms, 1D,'
  write ( *, '(a)' ) '  DFFT1I initializes the transforms,'
  write ( *, '(a)' ) '  DFFT1F does a forward transforms;'
  write ( *, '(a)' ) '  DFFT1B does a backward transforms.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r8vec_uniform_01 ( n, seed, r )

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

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

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

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

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

  deallocate ( wsave )

  return
end
subroutine test09

!*****************************************************************************80
!
!! TEST09 tests RFFT1B, RFFT1F and RFFT1I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

  integer ( kind = 4 ), parameter :: lenwrk = 2 * n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4vec_uniform_01 ( n, seed, r )

  call r4vec_print_some ( n, r, 1, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = n * int ( log ( real ( n, kind = 4 ) ) ) + 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 r4vec_print_some ( n, r, 1, 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 r4vec_print_some ( n, r, 1, 10, '  The retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test10

!*****************************************************************************80
!
!! TEST10 tests RFFT2B, RFFT2F and RFFT2I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

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

  integer ( kind = 4 ), parameter :: lenwrk = 2 * ldim * m

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  L =    ', l
  write ( *, '(a,i8)' ) '  LDIM = ', ldim
  write ( *, '(a,i8)' ) '  M =    ', m
!
!  Set the data values.
!
  seed = 1973

  call r4mat_uniform_01 ( ldim, m, seed, r )

  call r4mat_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, kind = 4 ) ) ) &
    + 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 r4mat_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 r4mat_print_some ( ldim, m, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test11

!*****************************************************************************80
!
!! TEST11 tests RFFTMB, RFFTMF and RFFTMI.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 32
  integer ( kind = 4 ), parameter :: lot = 6
  integer ( kind = 4 ), parameter :: lenr = n * lot

  integer ( kind = 4 ), parameter :: lenwrk = lot * n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of sequences is LOT =   ', lot
  write ( *, '(a,i8)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4mat_uniform_01 ( n, lot, seed, r )

  call r4mat_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, kind = 4 ) ) ) + 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 r4mat_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 r4mat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test12

!*****************************************************************************80
!
!! TEST12 tests SINQ1B, SINQ1F and SINQ1I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

  integer ( kind = 4 ), parameter :: lenwrk = n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4vec_uniform_01 ( n, seed, r )

  call r4vec_print_some ( n, r, 1, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = 2 * n + int ( log ( real ( n, kind = 4 ) ) ) + 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 r4vec_print_some ( n, r, 1, 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 r4vec_print_some ( n, r, 1, 10, '  The retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test13

!*****************************************************************************80
!
!! TEST13 tests SINQMB, SINQMF and SINQMI.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 32
  integer ( kind = 4 ), parameter :: lot = 6
  integer ( kind = 4 ), parameter :: lenr = n * lot

  integer ( kind = 4 ), parameter :: lenwrk = lot * n

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of sequences is LOT =   ', lot
  write ( *, '(a,i8)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4mat_uniform_01 ( n, lot, seed, r )

  call r4mat_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, kind = 4 ) ) ) + 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 r4mat_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 r4mat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test14

!*****************************************************************************80
!
!! TEST14 tests SINT1B, SINT1F and SINT1I.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4096

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

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST14'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of data items is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4vec_uniform_01 ( n, seed, r )

  call r4vec_print_some ( n, r, 1, 10, '  The original data:' )
!
!  Allocate and initialize the WSAVE array.
!
  lensav = n/2 + n + int ( log ( real ( n, kind = 4 ) ) ) + 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 r4vec_print_some ( n, r, 1, 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 r4vec_print_some ( n, r, 1, 10, '  The retrieved data:' )

  deallocate ( wsave )

  return
end
subroutine test15

!*****************************************************************************80
!
!! TEST15 tests SINTMB, SINTMF and SINTMI.
!
!  Modified:
!
!    16 November 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 32
  integer ( kind = 4 ), parameter :: lot = 6
  integer ( kind = 4 ), parameter :: lenr = n * lot

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

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

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST15'
  write ( *, '(a)' ) '  For single precision 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,i8)' ) '  The number of sequences is LOT =   ', lot
  write ( *, '(a,i8)' ) '  The length of each sequence is N = ', n
!
!  Set the data values.
!
  seed = 1973

  call r4mat_uniform_01 ( n, lot, seed, r )

  call r4mat_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, kind = 4 ) ) ) + 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 r4mat_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 r4mat_print_some ( n, lot, r, 1, 1, 5, 5, &
    '  Part of the retrieved data:' )

  deallocate ( wsave )

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

!*****************************************************************************80
!
!! C4MAT_PRINT_SOME prints some of a C4MAT.
!
!  Discussion:
!
!    A C4MAT is a matrix of complex ( kind = 4 ) values.
!
!  Modified:
!
!    23 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns 
!    in the matrix.
!
!    Input, complex ( kind = 4 ) A(M,N), the matrix.
!
!    Input, integer ( kind = 4 ) 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 ( kind = 4 ), parameter :: incx = 4
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  complex ( kind = 4 ) a(m,n)
  character ( len = 20 ) ctemp(incx)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i2hi
  integer ( kind = 4 ) i2lo
  integer ( kind = 4 ) ihi
  integer ( kind = 4 ) ilo
  integer ( kind = 4 ) inc
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j2
  integer ( kind = 4 ) j2hi
  integer ( kind = 4 ) j2lo
  integer ( kind = 4 ) jhi
  integer ( kind = 4 ) jlo
  character ( len = * ) title
  complex ( kind = 4 ) zero

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

  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 ( imag ( a(i,j) ) == 0.0E+00 ) then
          write ( ctemp(j2), '(g10.3,10x)' ) real ( a(i,j), kind = 4 )
        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 c4mat_uniform_01 ( m, n, seed, c )

!*****************************************************************************80
!
!! C4MAT_UNIFORM_01 returns a unit pseudorandom C4MAT.
!
!  Discussion:
!
!    A C4MAT is a matrix of complex ( kind = 4 ) values.
!
!    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:
!
!    31 May 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, Number 2, 1969, pages 136-143.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns 
!    in the matrix.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which 
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, complex ( kind = 4 ) C(M,N), the pseudorandom complex matrix.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

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

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'C4MAT_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

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

      k = seed / 127773

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

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

      r = sqrt ( real ( seed, kind = 4 ) * 4.656612875E-10 )

      k = seed / 127773

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

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

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

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

    end do

  end do

  return
end
subroutine c4vec_print_some ( n, x, i_lo, i_hi, title )

!*****************************************************************************80
!
!! C4VEC_PRINT_SOME prints some of a C4VEC.
!
!  Discussion:
!
!    A C4VEC is a vector of complex ( kind = 4 ) values.
!
!  Modified:
!
!    18 October 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries of the vector.
!
!    Input, complex ( kind = 4 ) X(N), the vector to be printed.
!
!    Input, integer ( kind = 4 ) I_LO, I_HI, the first and last entries
!    to print.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) i_hi
  integer ( kind = 4 ) i_lo
  character ( len = * ) title
  complex ( kind = 4 ) x(n)

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

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

  do i = max ( 1, i_lo ), min ( n, i_hi )
    write ( *, '(2x,i8,2x,2g14.6)' ) i, x(i)
  end do

  return
end
subroutine c4vec_uniform_01 ( n, seed, c )

!*****************************************************************************80
!
!! C4VEC_UNIFORM_01 returns a unit pseudorandom C4VEC.
!
!  Discussion:
!
!    A C4VEC is a vector of complex ( kind = 4 ) values.
!
!    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:
!
!    31 May 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, Number 2, 1969, pages 136-143.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of values to compute.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which 
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, complex ( kind = 4 ) C(N), the pseudorandom complex vector.
!
  implicit none

  integer ( kind = 4 ) n

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

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'C4VEC_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

  do i = 1, n

    k = seed / 127773

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

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

    r = sqrt (  real ( seed, kind = 4 ) * 4.656612875E-10 )

    k = seed / 127773

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

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

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

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

  end do

  return
end
function i4_huge ( )

!*****************************************************************************80
!
!! I4_HUGE returns a "huge" I4.
!
!  Discussion:
!
!    On an IEEE 32 bit machine, I4_HUGE should be 2**31 - 1, and its
!    bit pattern should be
!
!     01111111111111111111111111111111
!
!    In this case, its numerical value is 2147483647.
!
!    Using the Dec/Compaq/HP Alpha FORTRAN compiler FORT, I could
!    use I4_HUGE() and HUGE interchangeably.
!
!    However, when using the G95, the values returned by HUGE were
!    not equal to 2147483647, apparently, and were causing severe
!    and obscure errors in my random number generator, which needs to
!    add I4_HUGE to the seed whenever the seed is negative.  So I
!    am backing away from invoking HUGE, whereas I4_HUGE is under
!    my control.
!
!    Explanation: because under G95 the default integer type is 64 bits!
!    So HUGE ( 1 ) = a very very huge integer indeed, whereas
!    I4_HUGE ( ) = the same old 32 bit big value.
!
!    An I4 is an integer ( kind = 4 ) value.
!
!  Modified:
!
!    26 January 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ( kind = 4 ) I4_HUGE, a "huge" I4.
!
  implicit none

  integer ( kind = 4 ) i4
  integer ( kind = 4 ) i4_huge

  i4_huge = 2147483647

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

!*****************************************************************************80
!
!! R4MAT_PRINT_SOME prints some of an R4MAT.
!
!  Discussion:
!
!    An R4MAT is an array of R4 values.
!
!  Modified:
!
!    26 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns.
!
!    Input, real ( kind = 4 ) A(M,N), an M by N matrix to be printed.
!
!    Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print.
!
!    Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer ( kind = 4 ), parameter :: incx = 5
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real    ( kind = 4 ) a(m,n)
  character ( len = 14 ) ctemp(incx)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i2hi
  integer ( kind = 4 ) i2lo
  integer ( kind = 4 ) ihi
  integer ( kind = 4 ) ilo
  integer ( kind = 4 ) inc
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j2
  integer ( kind = 4 ) j2hi
  integer ( kind = 4 ) j2lo
  integer ( kind = 4 ) jhi
  integer ( kind = 4 ) 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), '(i8,6x)' ) 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) ), kind = 4 ) ) 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 r4mat_uniform_01 ( m, n, seed, r )

!*****************************************************************************80
!
!! R4MAT_UNIFORM_01 returns a unit pseudorandom R4MAT.
!
!  Discussion:
!
!    An R4MAT is an array of real ( kind = 4 ) values.
!
!  Modified:
!
!    31 May 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, Number 2, 1969, pages 136-143.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns 
!    in the array.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, 
!    which should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 4 ) R(M,N), the array of pseudorandom values.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) i4_huge
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed
  real    ( kind = 4 ) r(m,n)

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R4MAT_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

  do j = 1, n

    do i = 1, m

      k = seed / 127773

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

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

      r(i,j) = real ( seed, kind = 4 ) * 4.656612875E-10

    end do
  end do

  return
end
subroutine r4vec_print_some ( n, a, i_lo, i_hi, title )

!*****************************************************************************80
!
!! R4VEC_PRINT_SOME prints "some" of an R4VEC.
!
!  Discussion:
!
!    An R4VEC is a vector of R4 values.
!
!  Modified:
!
!    16 October 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries of the vector.
!
!    Input, real ( kind = 4 ) A(N), the vector to be printed.
!
!    Input, integer ( kind = 4 ) I_LO, I_HI, the first and last indices
!    to print.  The routine expects 1 <= I_LO <= I_HI <= N.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer ( kind = 4 ) n

  real    ( kind = 4 ) a(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i_hi
  integer ( kind = 4 ) i_lo
  character ( len = * ) title

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

  write ( *, '(a)' ) ' '
  do i = max ( i_lo, 1 ), min ( i_hi, n )
    write ( *, '(2x,i8,2x,g14.8)' ) i, a(i)
  end do

  return
end
subroutine r4vec_uniform_01 ( n, seed, r )

!*****************************************************************************80
!
!! R4VEC_UNIFORM_01 returns a unit pseudorandom R4VEC.
!
!  Discussion:
!
!    An R4VEC is an array of real ( kind = 4 ) values.
!
!  Modified:
!
!    31 May 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, Number 2, 1969, pages 136-143.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the vector.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, 
!    which should NOT be 0.
!    On output, SEED has been updated.
!
!    Output, real ( kind = 4 ) R(N), the vector of pseudorandom values.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) i4_huge
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed
  real    ( kind = 4 ) r(n)

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R4VEC_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

  do i = 1, n

    k = seed / 127773

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

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

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

  end do

  return
end
subroutine r8vec_print_some ( n, a, i_lo, i_hi, title )

!*****************************************************************************80
!
!! R8VEC_PRINT_SOME prints "some" of an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8 values.
!
!  Modified:
!
!    16 October 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries of the vector.
!
!    Input, real ( kind = 8 ) A(N), the vector to be printed.
!
!    Input, integer ( kind = 4 ) I_LO, I_HI, the first and last indices
!    to print.  The routine expects 1 <= I_LO <= I_HI <= N.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer ( kind = 4 ) n

  real    ( kind = 8 ) a(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i_hi
  integer ( kind = 4 ) i_lo
  character ( len = * ) title

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

  write ( *, '(a)' ) ' '
  do i = max ( i_lo, 1 ), min ( i_hi, n )
    write ( *, '(2x,i8,2x,g14.8)' ) i, a(i)
  end do

  return
end
subroutine r8vec_uniform_01 ( n, seed, r )

!*****************************************************************************80
!
!! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of real ( kind = 8 ) values.
!
!    For now, the input quantity SEED is an integer variable.
!
!  Modified:
!
!    31 May 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, Number 2, 1969, pages 136-143.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the vector.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which 
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 8 ) R(N), the vector of pseudorandom values.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) i4_huge
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed
  real    ( kind = 8 ) r(n)

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

  do i = 1, n

    k = seed / 127773

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

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

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

  end do

  return
end