program main !*****************************************************************************80 ! !! fftpack5_test() tests fftpack5(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 April 2020 ! ! Author: ! ! John Burkardt ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FFTPACK5_TEST' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Test FFTPACK5.' call cfft1_test ( ) call cfft2_test ( ) call cfftm_test ( ) call cosq1_test ( ) call cosqm_test ( ) call cost1_test ( ) call costm_test ( ) call dcosq1_test ( ) call dcost1_test ( ) call dfft1_test ( ) call dsint1_test ( ) call rfft1_test ( ) call rfft2_test ( ) call rfftm_test ( ) call sinq1_test ( ) call sinqm_test ( ) call sint1_test ( ) call sintm_test ( ) call zfft1_test ( ) call zfft2_test ( ) call zfftm_test ( ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FFTPACK5_TEST' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine cfft1_test ( ) !*****************************************************************************80 ! !! cfft1_test tests CFFT1B, CFFT1F and CFFT1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 4096 complex c(n) integer ier integer inc integer lenc integer lensav integer lenwrk integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'cfft1_test' write ( *, '(a)' ) ' For complex single precision fast Fourier transforms, 1D,' write ( *, '(a)' ) ' CFFT1I initializes the transform,' write ( *, '(a)' ) ' CFFT1F does a forward transform;' write ( *, '(a)' ) ' CFFT1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Allocate the work arrays. ! lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 lenwrk = 2 * n write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call cfft1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call c4vec_uniform_01 ( n, seed, c ) call c4vec_print_part ( n, c, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenc = n call cfft1f ( n, inc, c, lenc, wsave, lensav, work, lenwrk, ier ) call c4vec_print_part ( 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 c4vec_print_part ( n, c, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine cfft2_test ( ) !*****************************************************************************80 ! !! cfft2_test tests CFFT2B, CFFT2F and CFFT2I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: l = 32 integer, parameter :: m = 64 complex c(l,m) integer ier integer ldim integer lensav integer lenwrk integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'cfft2_test' write ( *, '(a)' ) ' For complex single precision fast Fourier transforms, 2D,' write ( *, '(a)' ) ' CFFT2I initializes the transform,' write ( *, '(a)' ) ' CFFT2F does a forward transform;' write ( *, '(a)' ) ' CFFT2B does a backward transform.' 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 ! ! Allocate work arrays. ! lenwrk = 2 * l * m lensav = 2 * l + int ( log ( real ( l ) ) / log ( 2.0E+00 ) ) & + 2 * m + int ( log ( real ( m ) ) / log ( 2.0E+00 ) ) & + 8 write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call cfft2i ( l, m, wsave, lensav, ier ) ! ! 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:' ) ! ! 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 ( work ) deallocate ( wsave ) return end subroutine cfftm_test ( ) !*****************************************************************************80 ! !! cfftm_test tests CFFTMB, CFFTMF and CFFTMI. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 32 integer, parameter :: lot = 6 complex, allocatable, dimension ( : ) :: c integer ier integer inc integer jump integer lenc integer lensav integer lenwrk integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'cfftm_test' write ( *, '(a)' ) ' For complex single precision fast Fourier transforms, 1D, multiple' write ( *, '(a)' ) ' CFFTMI initializes the transform,' write ( *, '(a)' ) ' CFFTMF does a forward transform;' write ( *, '(a)' ) ' CFFTMB does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of sequences is LOT = ', lot write ( *, '(a,i8)' ) ' The length of each sequence is N = ', n ! ! Set work vectors. ! lenc = n * lot lenwrk = 2 * lot * n lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENC = ', lenc write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( c(1:lenc) ) allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call cfftmi ( n, wsave, lensav, ier ) ! ! 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:' ) ! ! 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 ( c ) deallocate ( wsave ) deallocate ( work ) return end subroutine cosq1_test ( ) !*****************************************************************************80 ! !! cosq1_test tests COSQ1B, COSQ1F and COSQ1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 November 2007 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 4096 integer ier integer inc integer lenr integer lensav integer lenwrk real r(n) integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'cosq1_test' write ( *, '(a)' ) ' For real single precision fast cosine transforms, 1D,' write ( *, '(a)' ) ' COSQ1I initializes the transform,' write ( *, '(a)' ) ' COSQ1F does a forward transform;' write ( *, '(a)' ) ' COSQ1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work vectors. ! lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 lenwrk = n write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call cosq1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r4vec_print_part ( n, r, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call cosq1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r4vec_print_part ( 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 r4vec_print_part ( n, r, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine cosqm_test ( ) !*****************************************************************************80 ! !! cosqm_test tests COSQMB, COSQMF and COSQMI. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 32 integer, parameter :: lot = 6 integer ier integer inc integer jump integer lenr integer lensav integer lenwrk real, allocatable, dimension ( :, : ) :: r integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'cosqm_test():' write ( *, '(a)' ) ' For real single precision fast cosine transform, 1D, multiple' write ( *, '(a)' ) ' COSQMI initializes the transform,' write ( *, '(a)' ) ' COSQMF does a forward transform;' write ( *, '(a)' ) ' COSQMB does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of sequences is LOT = ', lot write ( *, '(a,i8)' ) ' The length of each sequence is N = ', n ! ! Set work arrays. ! lenr = n * lot lenwrk = lot * n lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENR = ', lenr write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( r(1:n,1:lot) ) allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call cosqmi ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n,1:lot) ) call r4mat_print_some ( n, lot, r, 1, 1, 5, 5, ' Part of the original data:' ) ! ! 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 ( r ) deallocate ( work ) deallocate ( wsave ) return end subroutine cost1_test ( ) !*****************************************************************************80 ! !! cost1_test tests COST1B, COST1F and COST1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 4096 integer ier integer inc integer lenr integer lensav integer lenwrk real r(n) integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'cost1_test' write ( *, '(a)' ) ' For real single precision fast cosine transforms, 1D,' write ( *, '(a)' ) ' COST1I initializes the transform,' write ( *, '(a)' ) ' COST1F does a forward transform;' write ( *, '(a)' ) ' COST1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work arrays. ! lenwrk = n - 1 lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call cost1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r4vec_print_part ( n, r, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call cost1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r4vec_print_part ( n, r, 10, ' The FFT coefficients:' ) ! ! Compute inverse FFT of coefficients. ! call cost1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r4vec_print_part ( n, r, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine costm_test ( ) !*****************************************************************************80 ! !! costm_test tests COSTMB, COSTMF and COSTMI. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 32 integer, parameter :: lot = 6 integer ier integer inc integer jump integer lenr integer lensav integer lenwrk real, allocatable, dimension ( :, : ) :: r integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'costm_test():' write ( *, '(a)' ) ' For real single precision fast cosine transforms, 1D, multiple' write ( *, '(a)' ) ' COSTMI initializes the transform,' write ( *, '(a)' ) ' COSTMF does a forward transform;' write ( *, '(a)' ) ' COSTMB does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of sequences is LOT = ', lot write ( *, '(a,i8)' ) ' The length of each sequence is N = ', n ! ! Set work arrays. ! lenr = n * lot lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 lenwrk = lot * ( n + 1 ) write ( *, '(a,i8)' ) ' LENR = ', lenr write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( r(1:n,1:lot) ) allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call costmi ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n,1:lot) ) call r4mat_print_some ( n, lot, r, 1, 1, 5, 5, ' Part of the original data:' ) ! ! 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 ( r ) deallocate ( work ) deallocate ( wsave ) return end subroutine dcosq1_test ( ) !*****************************************************************************80 ! !! dcosq1_test tests DCOSQ1B, DCOSQ1F and DCOSQ1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 November 2007 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 4096 integer, parameter :: lenwrk = n integer ier integer inc integer lenr integer lensav real ( kind = rk ) r(n) integer seed real ( kind = rk ) work(lenwrk) real ( kind = rk ), allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'dcosq1_test' write ( *, '(a)' ) ' For real double precision fast cosine transform, 1D,' write ( *, '(a)' ) ' DCOSQ1I initializes the transform,' write ( *, '(a)' ) ' DCOSQ1F does a forward transform;' write ( *, '(a)' ) ' DCOSQ1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work arrays. ! lensav = 2 * n + int ( log ( real ( n ) ) ) + 4 allocate ( wsave(1:lensav) ) call dcosq1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r8vec_print_part ( n, r, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call dcosq1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r8vec_print_part ( n, r, 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_part ( n, r, 10, ' The retrieved data:' ) deallocate ( wsave ) return end subroutine dcost1_test ( ) !*****************************************************************************80 ! !! dcost1_test tests DCOST1B, DCOST1F, DCOST1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 November 2007 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 4096 integer, parameter :: lenwrk = n - 1 integer, parameter :: lensav = 8204 integer ier integer inc integer lenr real ( kind = rk ) r(n) integer seed real ( kind = rk ) work(lenwrk) real ( kind = rk ) wsave(lensav) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'dcost1_test' write ( *, '(a)' ) ' For double precision fast cosine transforms, 1D,' write ( *, '(a)' ) ' DCOST1I initializes the transforms,' write ( *, '(a)' ) ' DCOST1F does a forward transforms;' write ( *, '(a)' ) ' DCOST1B 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_number ( harvest = r(1:n) ) call r8vec_print_part ( n, r, 10, ' The original data:' ) ! ! Allocate and initialize the WSAVE array. ! call dcost1i ( n, wsave, lensav, ier ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call dcost1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r8vec_print_part ( n, r, 10, ' The FFT coefficients:' ) ! ! Compute inverse FFT of coefficients. Should get back the ! original data. ! call dcost1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r8vec_print_part ( n, r, 10, ' The retrieved data:' ) return end subroutine dfft1_test ( ) !*****************************************************************************80 ! !! dfft1_test tests DFFT1B, DFFT1F and DFFT1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 November 2007 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 4096 integer ier integer inc integer lenr integer lensav integer lenwrk real ( kind = rk ) r(n) integer seed real ( kind = rk ), allocatable, dimension ( : ) :: work real ( kind = rk ), allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'dfft1_test' write ( *, '(a)' ) ' For real double precision fast Fourier transform, 1D,' write ( *, '(a)' ) ' DFFT1I initializes the transform,' write ( *, '(a)' ) ' DFFT1F does a forward transform;' write ( *, '(a)' ) ' DFFT1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work vectors. ! lensav = n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 lenwrk = n write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call dfft1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r8vec_print_part ( n, r, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call dfft1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r8vec_print_part ( n, r, 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_part ( n, r, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine dsint1_test ( ) !*****************************************************************************80 ! !! dsint1_test tests DSINT1B, DSINT1F, DSINT1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 November 2007 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 4096 integer, parameter :: lenwrk = 8194 integer, parameter :: lensav = 6156 integer ier integer inc integer lenr real ( kind = rk ) r(n) integer seed real ( kind = rk ) work(lenwrk) real ( kind = rk ) wsave(lensav) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DSINT1_TEST' write ( *, '(a)' ) ' For double precision fast sine transforms, 1D,' write ( *, '(a)' ) ' DSINT1I initializes the transforms,' write ( *, '(a)' ) ' DSINT1F does a forward transforms;' write ( *, '(a)' ) ' DSINT1B does a backward transforms.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r8vec_print_part ( n, r, 10, ' The original data:' ) ! ! Allocate and initialize the WSAVE array. ! call dsint1i ( n, wsave, lensav, ier ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call dsint1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r8vec_print_part ( n, r, 10, ' The FFT coefficients:' ) ! ! Compute inverse FFT of coefficients. Should get back the ! original data. ! call dsint1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r8vec_print_part ( n, r, 10, ' The retrieved data:' ) return end subroutine rfft1_test ( ) !*****************************************************************************80 ! !! rfft1_test tests RFFT1B, RFFT1F and RFFT1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 4096 integer ier integer inc integer lenr integer lensav integer lenwrk real r(n) integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'rfft1_test' write ( *, '(a)' ) ' For real single precision fast Fourier transforms, 1D,' write ( *, '(a)' ) ' RFFT1I initializes the transform,' write ( *, '(a)' ) ' RFFT1F does a forward transform;' write ( *, '(a)' ) ' RFFT1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work vectors. ! lensav = n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 lenwrk = n write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call rfft1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r4vec_print_part ( n, r, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call rfft1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r4vec_print_part ( 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 r4vec_print_part ( n, r, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine rfft2_test ( ) !*****************************************************************************80 ! !! rfft2_test tests RFFT2B, RFFT2F and RFFT2I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: l = 32 integer, parameter :: m = 64 integer, parameter :: ldim = 2 * ( l / 2 + 1 ) integer ier integer lensav integer lenwrk real r(ldim,m) integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'rfft2_test' write ( *, '(a)' ) ' For real single precision fast Fourier transform, 2D,' write ( *, '(a)' ) ' RFFT2I initializes the transform,' write ( *, '(a)' ) ' RFFT2F does a forward transform;' write ( *, '(a)' ) ' RFFT2B does a backward transform.' 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 work arrays. ! lenwrk = 2 * ldim * m lensav = 2 * ( l + m ) + int ( log ( real ( l ) ) ) & + int ( log ( real ( m ) ) ) + 8 write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call rfft2i ( l, m, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:ldim,1:m) ) call r4mat_print_some ( ldim, m, r, 1, 1, 5, 5, & ' Part of the original data:' ) ! ! 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 ( work ) deallocate ( wsave ) return end subroutine rfftm_test ( ) !*****************************************************************************80 ! !! rfftm_test tests RFFTMB, RFFTMF and RFFTMI. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 32 integer, parameter :: lot = 6 integer ier integer inc integer jump integer lenr integer lensav integer lenwrk real, allocatable, dimension ( :, : ) :: r integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'rfftm_test()' write ( *, '(a)' ) ' For real single precision fast Fourier transform, 1D, multiple' write ( *, '(a)' ) ' RFFTMI initializes the transform,' write ( *, '(a)' ) ' RFFTMF does a forward transform;' write ( *, '(a)' ) ' RFFTMB does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of sequences is LOT = ', lot write ( *, '(a,i8)' ) ' The length of each sequence is N = ', n ! ! Set work arrays. ! lenr = n * lot lenwrk = lot * n lensav = n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENR = ', lenr write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( r(1:n,1:lot) ) allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call rfftmi ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n,1:lot) ) call r4mat_print_some ( n, lot, r, 1, 1, 5, 5, ' Part of the original data:' ) ! ! 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 ( r ) deallocate ( work ) deallocate ( wsave ) return end subroutine sinq1_test ( ) !*****************************************************************************80 ! !! sinq1_test tests SINQ1B, SINQ1F and SINQ1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 4096 integer ier integer inc integer lenr integer lensav integer lenwrk real r(n) integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'sinq1_test' write ( *, '(a)' ) ' For real single precision fast sine transforms, 1D,' write ( *, '(a)' ) ' SINQ1I initializes the transform,' write ( *, '(a)' ) ' SINQ1F does a forward transform;' write ( *, '(a)' ) ' SINQ1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work arrays. ! lenwrk = n lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call sinq1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r4vec_print_part ( n, r, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call sinq1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r4vec_print_part ( 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 r4vec_print_part ( n, r, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine sinqm_test ( ) !*****************************************************************************80 ! !! sinqm_test tests SINQMB, SINQMF and SINQMI. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 August 2022 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 32 integer, parameter :: lot = 6 integer ier integer inc integer jump integer lenr integer lensav integer lenwrk real, allocatable, dimension ( :, : ) :: r integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'sinqm_test()' write ( *, '(a)' ) ' For real single precision fast sine transforms, 1D, multiple' write ( *, '(a)' ) ' SINQMI initializes the transform,' write ( *, '(a)' ) ' SINQMF does a forward transform;' write ( *, '(a)' ) ' SINQMB does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of sequences is LOT = ', lot write ( *, '(a,i8)' ) ' The length of each sequence is N = ', n ! ! Set work arrays. ! lenr = n * lot lenwrk = lot * n lensav = 2 * n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENR = ', lenr write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( r(1:n,1:lot) ) allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call sinqmi ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n,1:lot) ) call r4mat_print_some ( n, lot, r, 1, 1, 5, 5, ' Part of the original data:' ) ! ! 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 ( r ) deallocate ( work ) deallocate ( wsave ) return end subroutine sint1_test ( ) !*****************************************************************************80 ! !! sint1_test tests SINT1B, SINT1F and SINT1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 4096 integer ier integer inc integer lenr integer lensav integer lenwrk real r(n) integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'sint1_test' write ( *, '(a)' ) ' For real single precision fast sine transforms, 1D,' write ( *, '(a)' ) ' SINT1I initializes the transform,' write ( *, '(a)' ) ' SINT1F does a forward transform;' write ( *, '(a)' ) ' SINT1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work arrays. ! lenwrk = 2 * ( n + 1 ) lensav = n / 2 + n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call sint1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n) ) call r4vec_print_part ( n, r, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenr = n call sint1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) call r4vec_print_part ( 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 r4vec_print_part ( n, r, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine sintm_test ( ) !*****************************************************************************80 ! !! sintm_test tests SINTMB, SINTMF and SINTMI. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: n = 32 integer, parameter :: lot = 6 integer ier integer inc integer jump integer lenr integer lensav integer lenwrk real, allocatable, dimension ( :, : ) :: r integer seed real, allocatable, dimension ( : ) :: work real, allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'sintm_test' write ( *, '(a)' ) ' For real single precision fast sine transforms, 1D, multiple' write ( *, '(a)' ) ' SINTMI initializes the transform,' write ( *, '(a)' ) ' SINTMF does a forward transform;' write ( *, '(a)' ) ' SINTMB does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of sequences is LOT = ', lot write ( *, '(a,i8)' ) ' The length of each sequence is N = ', n ! ! Set work arrays. ! lenr = n * lot lenwrk = lot * 2 * ( n + 2 ) lensav = n / 2 + n + int ( log ( real ( n ) ) / log ( 2.0E+00 ) ) + 4 write ( *, '(a,i8)' ) ' LENR = ', lenr write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( r(1:n,1:lot) ) allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call sintmi ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call random_number ( harvest = r(1:n,1:lot) ) call r4mat_print_some ( n, lot, r, 1, 1, 5, 5, ' Part of the original data:' ) ! ! 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 ( r ) deallocate ( work ) deallocate ( wsave ) return end subroutine zfft1_test ( ) !*****************************************************************************80 ! !! zfft1_test tests ZFFT1B, ZFFT1F and ZFFT1I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 4096 complex ( kind = ck ) c(n) integer ier integer inc integer lenc integer lensav integer lenwrk integer seed real ( kind = rk ), allocatable, dimension ( : ) :: work real ( kind = rk ), allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'zfft1_test' write ( *, '(a)' ) ' For complex double precision fast Fourier transforms, 1D,' write ( *, '(a)' ) ' ZFFT1I initializes the transform,' write ( *, '(a)' ) ' ZFFT1F does a forward transform;' write ( *, '(a)' ) ' ZFFT1B does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of data items is N = ', n ! ! Set work arrays. ! lensav = 2 * n + int ( log ( real ( n, kind = rk ) ) ) + 4 lenwrk = 2 * n write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call zfft1i ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call c8vec_uniform_01 ( n, seed, c ) call c8vec_print_part ( n, c, 10, ' The original data:' ) ! ! Compute the FFT coefficients. ! inc = 1 lenc = n call zfft1f ( n, inc, c, lenc, wsave, lensav, work, lenwrk, ier ) call c8vec_print_part ( n, c, 10, ' The FFT coefficients:' ) ! ! Compute inverse FFT of coefficients. Should get back the ! original data. ! call zfft1b ( n, inc, c, lenc, wsave, lensav, work, lenwrk, ier ) call c8vec_print_part ( n, c, 10, ' The retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine zfft2_test ( ) !*****************************************************************************80 ! !! zfft2_test tests ZFFT2B, ZFFT2F and ZFFT2I. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: l = 32 integer, parameter :: m = 64 complex ( kind = ck ) c(l,m) integer ier integer ldim integer lensav integer lenwrk integer seed real ( kind = rk ), allocatable, dimension ( : ) :: work real ( kind = rk ), allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'zfft2_test' write ( *, '(a)' ) ' For complex double precision fast Fourier transforms, 2D,' write ( *, '(a)' ) ' ZFFT2I initializes the transform,' write ( *, '(a)' ) ' ZCFFT2F does a forward transform;' write ( *, '(a)' ) ' ZFFT2B does a backward transform.' 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 work arrays. ! lensav = 2 * ( l + m ) + int ( log ( real ( l, kind = rk ) ) ) & + int ( log ( real ( m, kind = rk ) ) ) + 8 lenwrk = 2 * l * m write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( wsave(1:lensav) ) allocate ( work(1:lenwrk) ) call zfft2i ( l, m, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call c8mat_uniform_01 ( l, m, seed, c ) call c8mat_print_some ( l, m, c, 1, 1, 5, 5, ' Part of the original data:' ) ! ! Compute the FFT coefficients. ! ldim = l call zfft2f ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier ) call c8mat_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 zfft2b ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier ) call c8mat_print_some ( l, m, c, 1, 1, 5, 5, ' Part of the retrieved data:' ) deallocate ( work ) deallocate ( wsave ) return end subroutine zfftm_test ( ) !*****************************************************************************80 ! !! zfftm_test tests ZFFTMB, ZFFTMF and ZFFTMI. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 August 2011 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 32 integer, parameter :: lot = 6 complex ( kind = ck ), allocatable, dimension ( : ) :: c integer ier integer inc integer jump integer lenc integer lensav integer lenwrk integer seed real ( kind = rk ), allocatable, dimension ( : ) :: work real ( kind = rk ), allocatable, dimension ( : ) :: wsave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'zfftm_test' write ( *, '(a)' ) ' For complex double precision fast Fourier transforms, 1D, multiple' write ( *, '(a)' ) ' ZFFTMI initializes the transform,' write ( *, '(a)' ) ' ZFFTMF does a forward transform;' write ( *, '(a)' ) ' ZFFTMB does a backward transform.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of sequences is LOT = ', lot write ( *, '(a,i8)' ) ' The length of each sequence is N = ', n ! ! Set work arrays. ! lenc = n * lot lensav = 2 * n + int ( log ( real ( n, kind = rk ) ) ) + 4 lenwrk = 2 * lot * n write ( *, '(a,i8)' ) ' LENC = ', lenc write ( *, '(a,i8)' ) ' LENSAV = ', lensav write ( *, '(a,i8)' ) ' LENWRK = ', lenwrk allocate ( c(1:lenc) ) allocate ( work(1:lenwrk) ) allocate ( wsave(1:lensav) ) call zfftmi ( n, wsave, lensav, ier ) ! ! Set the data values. ! seed = 1973 call c8mat_uniform_01 ( n, lot, seed, c ) call c8mat_print_some ( n, lot, c, 1, 1, 5, 5, ' Part of the original data:' ) ! ! Compute the FFT coefficients. ! jump = n inc = 1 call zfftmf ( lot, jump, n, inc, c, lenc, wsave, lensav, work, lenwrk, ier ) call c8mat_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 zfftmb ( lot, jump, n, inc, c, lenc, wsave, lensav, work, lenwrk, ier ) call c8mat_print_some ( n, lot, c, 1, 1, 5, 5, & ' Part of the retrieved data:' ) deallocate ( c ) deallocate ( work ) 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 C4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! 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. ! 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 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) ! ! 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) ) 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 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 C4's. ! ! 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. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! 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 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 i4_huge integer j real r integer k real, parameter :: pi = 3.1415926E+00 integer seed real 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 ) * 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 ) * 4.656612875E-10 ) c(i,j) = r * cmplx ( cos ( theta ), sin ( theta ) ) end do end do return end subroutine c4vec_print_part ( n, a, max_print, title ) !*****************************************************************************80 ! !! C4VEC_PRINT_PART prints "part" of a C4VEC. ! ! Discussion: ! ! The user specifies MAX_PRINT, the maximum number of lines to print. ! ! If N, the size of the vector, is no more than MAX_PRINT, then ! the entire vector is printed, one entry per line. ! ! Otherwise, if possible, the first MAX_PRINT-2 entries are printed, ! followed by a line of periods suggesting an omission, ! and the last entry. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 June 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries of the vector. ! ! Input, complex A(N), the vector to be printed. ! ! Input, integer MAX_PRINT, the maximum number of lines ! to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer n complex 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' if ( n <= max_print ) then do i = 1, n write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) end do else if ( 3 <= max_print ) then do i = 1, max_print - 2 write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) end do write ( *, '(a)' ) ' ........ .............. ..............' i = n write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) else do i = 1, max_print - 1 write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) end do i = max_print write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6,2x,a)' ) i, ':', a(i), & '...more entries...' end if 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 C4's. ! ! 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. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! 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 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 integer i4_huge integer k real, parameter :: pi = 3.1415926E+00 real r integer seed real 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 ) * 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 ) * 4.656612875E-10 ) c(i) = r * cmplx ( cos ( theta ), sin ( theta ) ) end do return end subroutine c8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! C8MAT_PRINT_SOME prints some of a C8MAT. ! ! Discussion: ! ! A C8MAT is a matrix of C8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns ! in the matrix. ! ! Input, complex ( kind = ck ) 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. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 4 integer m integer n complex ( kind = ck ) 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 ( kind = ck ) zero zero = cmplx ( 0.0D+00, 0.0D+00, kind = ck ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) ! ! 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.0D+00 ) then write ( ctemp(j2), '(g10.3,10x)' ) real ( a(i,j), kind = rk ) 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 return end subroutine c8mat_uniform_01 ( m, n, seed, c ) !*****************************************************************************80 ! !! C8MAT_UNIFORM_01 returns a unit pseudorandom C8MAT. ! ! Discussion: ! ! A C8MAT is a matrix of C8's. ! ! For now, the input quantity SEED is an integer variable. ! ! 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. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! 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 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 ( kind = ck ) C(M,N), the pseudorandom complex matrix. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n complex ( kind = ck ) c(m,n) integer i integer, parameter :: i4_huge = 2147483647 integer j real ( kind = rk ) r integer k real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 integer seed real ( kind = rk ) theta if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'C8MAT_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 = rk ) * 4.656612875D-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 = rk ) * 4.656612875D-10 ) c(i,j) = r * cmplx ( cos ( theta ), sin ( theta ), kind = ck ) end do end do return end subroutine c8vec_print_part ( n, a, max_print, title ) !*****************************************************************************80 ! !! C8VEC_PRINT_PART prints "part" of a C8VEC. ! ! 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. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 June 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries of the vector. ! ! Input, complex ( kind = ck ) A(N), the vector to be printed. ! ! Input, integer MAX_PRINT, the maximum number of lines ! to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer n complex ( kind = ck ) 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' if ( n <= max_print ) then do i = 1, n write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) end do else if ( 3 <= max_print ) then do i = 1, max_print - 2 write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) end do write ( *, '(a)' ) ' ........ .............. ..............' i = n write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) else do i = 1, max_print - 1 write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(i) end do i = max_print write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6,2x,a)' ) i, ':', a(i), & '...more entries...' end if return end subroutine c8vec_uniform_01 ( n, seed, c ) !*****************************************************************************80 ! !! C8VEC_UNIFORM_01 returns a unit pseudorandom C8VEC. ! ! Discussion: ! ! A C8VEC is a vector of C8's. ! ! For now, the input quantity SEED is an integer variable. ! ! 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. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! 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 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 ( kind = ck ) C(N), the pseudorandom complex vector. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer n complex ( kind = ck ) c(n) integer i integer, parameter :: i4_huge = 2147483647 real ( kind = rk ) r integer k real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 integer seed real ( kind = rk ) theta if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'C8VEC_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 = rk ) * 4.656612875D-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 = rk ) * 4.656612875D-10 ) c(i) = r * cmplx ( cos ( theta ), sin ( theta ), kind = ck ) 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 value. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 January 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer I4_HUGE, a "huge" I4. ! implicit none integer 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's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 March 2005 ! ! 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, a 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) 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) ) ) ) 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 return end subroutine r4vec_print_part ( n, a, max_print, title ) !*****************************************************************************80 ! !! R4VEC_PRINT_PART prints "part" of an R4VEC. ! ! Discussion: ! ! The user specifies MAX_PRINT, the maximum number of lines to print. ! ! If N, the size of the vector, is no more than MAX_PRINT, then ! the entire vector is printed, one entry per line. ! ! Otherwise, if possible, the first MAX_PRINT-2 entries are printed, ! followed by a line of periods suggesting an omission, ! and the last entry. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 June 2010 ! ! 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, a 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' if ( n <= max_print ) then do i = 1, n write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) end do else if ( 3 <= max_print ) then do i = 1, max_print - 2 write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) end do write ( *, '(a)' ) ' ........ ..............' i = n write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) else do i = 1, max_print - 1 write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) end do i = max_print write ( *, '(2x,i8,a,1x,g14.6,2x,a)' ) i, ':', a(i), '...more entries...' end if return end subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_PRINT_SOME prints some of an R8MAT. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) 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, a title. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 5 integer m integer n real ( kind = rk ) 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) 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 = rk ) ) 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 return end subroutine r8vec_print_part ( n, a, max_print, title ) !*****************************************************************************80 ! !! R8VEC_PRINT_PART prints "part" of an R8VEC. ! ! 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. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries of the vector. ! ! Input, real ( kind = rk ) A(N), the vector to be printed. ! ! Input, integer MAX_PRINT, the maximum number of lines ! to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' if ( n <= max_print ) then do i = 1, n write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) end do else if ( 3 <= max_print ) then do i = 1, max_print - 2 write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) end do write ( *, '(a)' ) ' ........ ..............' i = n write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) else do i = 1, max_print - 1 write ( *, '(2x,i8,a,1x,g14.6)' ) i, ':', a(i) end do i = max_print write ( *, '(2x,i8,a,1x,g14.6,2x,a)' ) i, ':', a(i), '...more entries...' end if return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 August 2005 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end