diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index bee0a4c..6503cd6 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -30,3 +30,4 @@ ADDTEST(fftpack_fft) ADDTEST(fftpack_rfft) ADDTEST(fftpack_dct) ADDTEST(fftpack_utils) +ADDTEST(fftpack_original) diff --git a/test/test_fftpack_original.f90 b/test/test_fftpack_original.f90 new file mode 100644 index 0000000..1cafa1e --- /dev/null +++ b/test/test_fftpack_original.f90 @@ -0,0 +1,463 @@ +module test_fftpack_original + use fftpack + use testdrive, only: new_unittest, unittest_type, error_type, check + implicit none(type, external) + private + + public :: collect_original + + real(rk), parameter :: pi = 4.0_rk*atan(1.0_rk) + real(rk), parameter :: atol = epsilon(1.0_rk) + real(rk), parameter :: rtol = sqrt(atol) + integer, parameter :: nd(1:7) = [120, 54, 49, 32, 4, 3, 2] + +contains + + subroutine collect_original(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [new_unittest("dfft", test_dfft)] + testsuite = [testsuite, new_unittest("zfft", test_zfft)] + testsuite = [testsuite, new_unittest("sint", test_sint)] + testsuite = [testsuite, new_unittest("cost", test_cost)] + testsuite = [testsuite, new_unittest("cosqt", test_cosqt)] + testsuite = [testsuite, new_unittest("dzfft", test_dzfft)] + end subroutine collect_original + + subroutine test_dfft(error) + type(error_type), allocatable, intent(out) :: error + real(rk) :: x(200), y(200), xh(200), w(2000) + integer :: i, j, k, n, np1, nm1, ns2, nz, modn + real(rk) :: dt, sum1, sum2, arg, arg1 + real(rk) :: mismatch, cf + + do nz = 1, size(nd) + !> Create multisine signal. + n = nd(nz) + modn = mod(n, 2) + np1 = n + 1; nm1 = n - 1 + do concurrent(j=1:np1) + x(j) = sin(j*sqrt(2.0_rk)) + end do + y = x; xh = x + + !> Discrete Fourier Transform. + dt = 2*pi/n + ns2 = (n + 1)/2 + if (ns2 >= 2) then + do k = 2, ns2 + sum1 = 0.0_rk; sum2 = 0.0_rk + arg = (k - 1)*dt + do i = 1, n + arg1 = (i - 1)*arg + sum1 = sum1 + x(i)*cos(arg1) + sum2 = sum2 + x(i)*sin(arg1) + end do + y(2*k - 2) = sum1 + y(2*k - 1) = -sum2 + end do + end if + sum1 = sum(x(1:nm1:2)) + sum2 = sum(x(2:nm1 + 1:2)) + if (modn == 1) sum1 = sum1 + x(n) + y(1) = sum1 + sum2 + if (modn == 0) y(n) = sum1 - sum2 + + !> Perform (real) Fourier Transform. + call dffti(n, w) + call dfftf(n, x, w) + + !> Check error. + mismatch = maxval(abs(x(:n) - y(:n)))/n + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Inverse Discrete Fourier Transform. + x(:n) = xh(:n) ! Restore signal. + do i = 1, n + sum1 = 0.5_rk*x(1) + arg = (i - 1)*dt + if (ns2 >= 2) then + do k = 2, ns2 + arg1 = (k - 1)*arg + sum1 = sum1 + x(2*k - 2)*cos(arg1) - x(2*k - 1)*sin(arg1) + end do + end if + if (modn == 0) sum1 = sum1 + 0.5_rk*(-1)**(i - 1)*x(n) + y(i) = 2*sum1 + end do + + !> Perform (real) inverse Fourier Transform. + call dfftb(n, x, w) + + !> Check error. + mismatch = maxval(abs(x(:n) - y(:n))) + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Chain direct and inverse Fourier transforms. + x(:n) = xh(:n); y(:n) = xh(:n) ! Restore signal. + call dfftb(n, y, w) + call dfftf(n, y, w) + + !> Check error. + cf = 1.0_rk/n + mismatch = maxval(abs(cf*y(:n) - x(:n))) + call check(error, mismatch < rtol) + if (allocated(error)) return + end do + end subroutine test_dfft + + subroutine test_zfft(error) + type(error_type), allocatable, intent(out) :: error + integer :: i, j, k, n, nz + complex(rk) :: cx(200), cy(200) + real(rk) :: w(2000), dt, arg1, arg2, mismatch, cf + + do nz = 1, size(nd) + !> Create signal. + n = nd(nz) + do concurrent(i=1:n) + cx(i) = cmplx(cos(sqrt(2.0_rk)*i), sin(sqrt(2.0_rk)*i**2), kind=rk) + end do + + !> Discrete Fourier Transform. + dt = 2*pi/n + do i = 1, n + arg1 = -(i - 1)*dt + cy(i) = cmplx(0.0_rk, 0.0_rk, kind=rk) + do k = 1, n + arg2 = (k - 1)*arg1 + cy(i) = cy(i) + cmplx(cos(arg2), sin(arg2), kind=rk)*cx(k) + end do + end do + + !> Fast Fourier Transform. + call zffti(n, w) + call zfftf(n, cx, w) + + !> Check error. + mismatch = maxval(abs(cx(:n) - cy(:n)))/n + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Inverse Discrete Fourier Transform. + cx(:n) = cx(:n)/n ! Scale signal. + do i = 1, n + arg1 = (i - 1)*dt + cy(i) = cmplx(0.0_rk, 0.0_rk, kind=rk) + do k = 1, n + arg2 = (k - 1)*arg1 + cy(i) = cy(i) + cmplx(cos(arg2), sin(arg2), kind=rk)*cx(k) + end do + end do + + !> Inverse Fast Fourier Transform. + call zfftb(n, cx, w) + + !> Check error. + mismatch = maxval(abs(cx(:n) - cy(:n))) + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Chain direct and inverse Fourier transforms. + cx(:n) = cy(:n) ! Restore signal. + call zfftf(n, cx, w) + call zfftb(n, cx, w) + + !> Check error. + cf = 1.0_rk/n + mismatch = maxval(abs(cf*cx(:n) - cy(:n))) + call check(error, mismatch < rtol) + if (allocated(error)) return + end do + end subroutine test_zfft + + subroutine test_sint(error) + type(error_type), allocatable, intent(out) :: error + real(rk) :: x(200), y(200), xh(200), w(2000) + integer :: i, j, k, n, np1, nm1, ns2, nz + real(rk) :: dt, sum1, sum2, arg, arg1 + real(rk) :: mismatch, cf + + do nz = 1, size(nd) + !> Create multisine signal. + n = nd(nz) + np1 = n + 1; nm1 = n - 1 + do concurrent(j=1:np1) + x(j) = sin(j*sqrt(2.0_rk)) + end do + y = x; xh = x + + !> Discrete sine transform. + dt = pi/n + + do i = 1, nm1 + y(i) = 0.0_rk + arg1 = i*dt + do k = 1, nm1 + y(i) = y(i) + x(k)*sin(k*arg1) + end do + y(i) = 2*y(i) + end do + + !> Fast Sine Transform. + call dsinti(nm1, w) + call dsint(nm1, x, w) + + !> Check error. + cf = 0.5_rk/n + mismatch = maxval(abs(x(:nm1) - y(:nm1)))*cf + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Chain direct and inverse sine transform. + x(:nm1) = xh(:nm1); y(:nm1) = x(:nm1) ! Restore signals. + call dsint(nm1, x, w) + call dsint(nm1, x, w) + + !> Check error. + mismatch = maxval(abs(cf*x(:nm1) - y(:nm1))) + call check(error, mismatch < rtol) + if (allocated(error)) return + end do + end subroutine test_sint + + subroutine test_cost(error) + type(error_type), allocatable, intent(out) :: error + real(rk) :: x(200), y(200), xh(200), w(2000) + integer :: i, j, k, n, np1, nm1, ns2, nz + real(rk) :: dt, sum1, sum2, arg, arg1 + real(rk) :: mismatch, cf + + do nz = 1, size(nd) + !> Create multisine signal. + n = nd(nz) + np1 = n + 1; nm1 = n - 1 + do concurrent(j=1:np1) + x(j) = sin(j*sqrt(2.0_rk)) + end do + y = x; xh = x + + !> Discrete sine transform. + dt = pi/n + + do i = 1, np1 + y(i) = 0.5_rk*(x(1) + (-1)**(i + 1)*x(n + 1)) + arg = (i - 1)*dt + do k = 2, n + y(i) = y(i) + cos((k - 1)*arg)*x(k) + end do + y(i) = 2*y(i) + end do + + !> Fast Sine Transform. + call dcosti(np1, w) + call dcost(np1, x, w) + + !> Check error. + cf = 0.5_rk/n + mismatch = maxval(abs(x(:np1) - y(:np1)))*cf + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Chain direct and inverse sine transform. + x(:np1) = xh(:np1); y(:np1) = x(:np1) ! Restore signals. + call dcost(np1, x, w) + call dcost(np1, x, w) + + !> Check error. + mismatch = maxval(abs(cf*x(:nm1) - y(:nm1))) + call check(error, mismatch < rtol) + if (allocated(error)) return + end do + end subroutine test_cost + + subroutine test_cosqt(error) + type(error_type), allocatable, intent(out) :: error + real(rk) :: x(200), y(200), xh(200), w(2000) + integer :: i, j, k, n, np1, nm1, ns2, nz + real(rk) :: dt, sum1, sum2, arg, arg1 + real(rk) :: mismatch, cf + + do nz = 1, size(nd) + !> Create multisine signal. + n = nd(nz) + np1 = n + 1; nm1 = n - 1 + do concurrent(j=1:np1) + x(j) = sin(j*sqrt(2.0_rk)) + end do + y = x; xh = x + + !> Discrete quater-cos transform. + dt = pi/(2*n) + + do i = 1, n + x(i) = 0.0_rk + arg = (i - 1)*dt + do k = 1, n + x(i) = x(i) + y(k)*cos((2*k - 1)*arg) + end do + x(i) = 4*x(i) + end do + + !> Fast Quarter-cos Transform. + call dcosqi(n, w) + call dcosqb(n, y, w) + + !> Check error. + cf = 0.25_rk/n + mismatch = maxval(abs(x(:n) - y(:n)))*cf + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Discrete inverse quarter-cos transform. + x(:n) = xh(:n) + do i = 1, n + y(i) = 0.5_rk*x(1) + arg = (2*i - 1)*dt + do k = 2, n + y(i) = y(i) + x(k)*cos((k - 1)*arg) + end do + y(i) = 2*y(i) + end do + + !> Fast inverse quarter-cos transform. + call dcosqf(n, x, w) + + !> Check error. + mismatch = maxval(abs(y(:n) - x(:n)))*cf + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Chain direct and inverse quarter-cos transforms. + x(:n) = xh(:n); y(:n) = xh(:n) ! Restore signals. + call dcosqb(n, x, w) + call dcosqf(n, x, w) + + !> Check error. + mismatch = maxval(abs(cf*x(:n) - y(:n))) + call check(error, mismatch < rtol) + if (allocated(error)) return + + end do + end subroutine test_cosqt + + subroutine test_dzfft(error) + type(error_type), allocatable, intent(out) :: error + real(rk) :: x(200), y(200), xh(200), w(2000) + real(rk) :: a(100), b(100), ah(100), bh(100) + real(rk) :: azero, azeroh + integer :: i, j, k, n, np1, nm1, ns2, ns2m, nz, modn + real(rk) :: dt, sum1, sum2, arg, arg1, arg2 + real(rk) :: mismatch, cf + + do nz = 1, size(nd) + !> Create multisine signal. + n = nd(nz) + modn = mod(n, 2) + np1 = n + 1; nm1 = n - 1 + do concurrent(j=1:np1) + x(j) = sin(j*sqrt(2.0_rk)) + end do + y = x; xh = x + + !> Discrete Fourier Transform. + dt = 2*pi/n + ns2 = (n + 1)/2 + ns2m = ns2 - 1 + cf = 2.0_rk/n + if (ns2m > 0) then + do k = 1, ns2m + sum1 = 0.0_rk; sum2 = 0.0_rk + arg = k*dt + do i = 1, n + arg1 = (i - 1)*arg + sum1 = sum1 + x(i)*cos(arg1) + sum2 = sum2 + x(i)*sin(arg1) + end do + a(k) = cf*sum1 + b(k) = cf*sum2 + end do + end if + nm1 = n - 1 + sum1 = sum(x(1:nm1:2)) + sum2 = sum(x(2:nm1 + 1:2)) + if (modn == 1) sum1 = sum1 + x(n) + azero = 0.5_rk*cf*(sum1 + sum2) + if (modn == 0) a(ns2) = 0.5_rk*cf*(sum1 - sum2) + + !> Fast Fourier Transform. + call dzffti(n, w) + call dzfftf(n, x, azeroh, ah, bh, w) + + !> Check error. + mismatch = abs(azeroh - azero) + if (modn == 0) mismatch = max(mismatch, abs(a(ns2) - ah(ns2))) + if (ns2m > 0) then + do i = 1, ns2m + mismatch = max(mismatch, abs(ah(i) - a(i)), abs(bh(i) - b(i))) + end do + end if + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Inverse Discrete Fourier Transform + ns2 = n/2 + if (modn == 0) b(ns2) = 0.0_rk + do i = 1, n + sum1 = azero + arg1 = (i - 1)*dt + do k = 1, ns2 + arg2 = k*arg1 + sum1 = sum1 + a(k)*cos(arg2) + b(k)*sin(arg2) + end do + x(i) = sum1 + end do + + !> Fast Inverse Fourier Transform. + call dzfftb(n, y, azero, a, b, w) + + !> Check error. + mismatch = maxval(abs(x(:n) - y(:n))) + call check(error, mismatch < rtol) + if (allocated(error)) return + + !> Chain direct and inverse Fourier transforms. + x(:n) = xh(:n) + call dzfftf(n, x, azero, a, b, w) + call dzfftb(n, x, azero, a, b, w) + + !> Check error. + mismatch = maxval(abs(x(:n) - y(:n))) + call check(error, mismatch < rtol) + if (allocated(error)) return + + end do + end subroutine test_dzfft +end module test_fftpack_original + +program test_original + use, intrinsic :: iso_fortran_env, only: error_unit + use testdrive, only: run_testsuite, new_testsuite, testsuite_type + use test_fftpack_original, only: collect_original + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("Original test suite", collect_original) & + ] + + do is = 1, size(testsuites) + write (error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write (error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program diff --git a/test/tstfft.f b/test/tstfft.f deleted file mode 100644 index 36a9a3f..0000000 --- a/test/tstfft.f +++ /dev/null @@ -1,411 +0,0 @@ - PROGRAM TSTFFT -C -C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * -C -C VERSION 4 APRIL 1985 -C -C A TEST DRIVER FOR -C A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE FAST FOURIER -C TRANSFORM OF PERIODIC AND OTHER SYMMETRIC SEQUENCES -C -C BY -C -C PAUL N SWARZTRAUBER -C -C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 -C -C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION -C -C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * -C -C -C THIS PROGRAM TESTS THE PACKAGE OF FAST FOURIER -C TRANSFORMS FOR BOTH COMPLEX AND REAL PERIODIC SEQUENCES AND -C CERTIAN OTHER SYMMETRIC SEQUENCES THAT ARE LISTED BELOW. -C -C 1. RFFTI INITIALIZE RFFTF AND RFFTB -C 2. RFFTF FORWARD TRANSFORM OF A REAL PERIODIC SEQUENCE -C 3. RFFTB BACKWARD TRANSFORM OF A REAL COEFFICIENT ARRAY -C -C 4. EZFFTI INITIALIZE EZFFTF AND EZFFTB -C 5. EZFFTF A SIMPLIFIED REAL PERIODIC FORWARD TRANSFORM -C 6. EZFFTB A SIMPLIFIED REAL PERIODIC BACKWARD TRANSFORM -C -C 7. SINTI INITIALIZE SINT -C 8. SINT SINE TRANSFORM OF A REAL ODD SEQUENCE -C -C 9. COSTI INITIALIZE COST -C 10. COST COSINE TRANSFORM OF A REAL EVEN SEQUENCE -C -C 11. SINQI INITIALIZE SINQF AND SINQB -C 12. SINQF FORWARD SINE TRANSFORM WITH ODD WAVE NUMBERS -C 13. SINQB UNNORMALIZED INVERSE OF SINQF -C -C 14. COSQI INITIALIZE COSQF AND COSQB -C 15. COSQF FORWARD COSINE TRANSFORM WITH ODD WAVE NUMBERS -C 16. COSQB UNNORMALIZED INVERSE OF COSQF -C -C 17. CFFTI INITIALIZE CFFTF AND CFFTB -C 18. CFFTF FORWARD TRANSFORM OF A COMPLEX PERIODIC SEQUENCE -C 19. CFFTB UNNORMALIZED INVERSE OF CFFTF - -C *** HACKED BY HCP FOR THE DOUBLE PREC. VERSION NOVEMEMBER 1999 - - -C - USE fftpack_kind - IMPLICIT REAL(RK) (A-H,O-Z) - DIMENSION ND(10) ,X(200) ,Y(200) ,W(2000) , - 1 A(100) ,B(100) ,AH(100) ,BH(100) , - 2 XH(200) ,CX(200) ,CY(200) - COMPLEX(RK) CX ,CY - DATA ND(1),ND(2),ND(3),ND(4),ND(5),ND(6),ND(7)/120,54,49,32,4,3,2/ - SQRT2 = SQRT(2.0D0) - NNS = 7 - DO 157 NZ=1,NNS - N = ND(NZ) - MODN = MOD(N,2) - FN = REAL(N,RK) - TFN = FN+FN - NP1 = N+1 - NM1 = N-1 - DO 101 J=1,NP1 - X(J) = SIN(REAL(J,RK)*SQRT2) - Y(J) = X(J) - XH(J) = X(J) - 101 CONTINUE -C -C TEST SUBROUTINES RFFTI,RFFTF AND RFFTB -C - CALL DFFTI (N,W) - PI = 3.14159265358979323846D0 - DT = (PI+PI)/FN - NS2 = (N+1)/2 - IF (NS2 .LT. 2) GO TO 104 - DO 103 K=2,NS2 - SUM1 = 0.0D0 - SUM2 = 0.0D0 - ARG = REAL(K-1,RK)*DT - DO 102 I=1,N - ARG1 = REAL(I-1,RK)*ARG - SUM1 = SUM1+X(I)*COS(ARG1) - SUM2 = SUM2+X(I)*SIN(ARG1) - 102 CONTINUE - Y(2*K-2) = SUM1 - Y(2*K-1) = -SUM2 - 103 CONTINUE - 104 SUM1 = 0.0D0 - SUM2 = 0.0D0 - DO 105 I=1,NM1,2 - SUM1 = SUM1+X(I) - SUM2 = SUM2+X(I+1) - 105 CONTINUE - IF (MODN .EQ. 1) SUM1 = SUM1+X(N) - Y(1) = SUM1+SUM2 - IF (MODN .EQ. 0) Y(N) = SUM1-SUM2 - CALL DFFTF (N,X,W) - RFTF = 0.0D0 - DO 106 I=1,N - RFTF = DMAX1(RFTF,ABS(X(I)-Y(I))) - X(I) = XH(I) - 106 CONTINUE - RFTF = RFTF/FN - DO 109 I=1,N - SUM = 0.5D0*X(1) - ARG = REAL(I-1,RK)*DT - IF (NS2 .LT. 2) GO TO 108 - DO 107 K=2,NS2 - ARG1 = REAL(K-1,RK)*ARG - SUM = SUM+X(2*K-2)*COS(ARG1)-X(2*K-1)*SIN(ARG1) - 107 CONTINUE - 108 IF (MODN .EQ. 0) SUM = SUM+.5*REAL((-1)**(I-1),RK)*X(N) - Y(I) = SUM+SUM - 109 CONTINUE - CALL DFFTB (N,X,W) - RFTB = 0.0D0 - DO 110 I=1,N - RFTB = DMAX1(RFTB,ABS(X(I)-Y(I))) - X(I) = XH(I) - Y(I) = XH(I) - 110 CONTINUE - CALL DFFTB (N,Y,W) - CALL DFFTF (N,Y,W) - CF = 1.0D0/FN - RFTFB = 0. - DO 111 I=1,N - RFTFB = DMAX1(RFTFB,ABS(CF*Y(I)-X(I))) - 111 CONTINUE -C -C TEST SUBROUTINES DSINTI AND DSINT -C - DT = PI/FN - DO 112 I=1,NM1 - X(I) = XH(I) - 112 CONTINUE - DO 114 I=1,NM1 - Y(I) = 0.0D0 - ARG1 = REAL(I,RK)*DT - DO 113 K=1,NM1 - Y(I) = Y(I)+X(K)*SIN(REAL(K,RK)*ARG1) - 113 CONTINUE - Y(I) = Y(I)+Y(I) - 114 CONTINUE - CALL DSINTI (NM1,W) - CALL DSINT (NM1,X,W) - CF = 0.5D0/FN - SINTT = 0.0D0 - DO 115 I=1,NM1 - SINTT = DMAX1(SINTT,ABS(X(I)-Y(I))) - X(I) = XH(I) - Y(I) = X(I) - 115 CONTINUE - SINTT = CF*SINTT - CALL DSINT (NM1,X,W) - CALL DSINT (NM1,X,W) - SINTFB = 0.0D0 - DO 116 I=1,NM1 - SINTFB = DMAX1(SINTFB,ABS(CF*X(I)-Y(I))) - 116 CONTINUE -C -C TEST SUBROUTINES COSTI AND COST -C - DO 117 I=1,NP1 - X(I) = XH(I) - 117 CONTINUE - DO 119 I=1,NP1 - Y(I) = 0.5D0*(X(1)+REAL((-1)**(I+1),RK)*X(N+1)) - ARG = REAL(I-1,RK)*DT - DO 118 K=2,N - Y(I) = Y(I)+X(K)*COS(REAL(K-1,RK)*ARG) - 118 CONTINUE - Y(I) = Y(I)+Y(I) - 119 CONTINUE - CALL DCOSTI (NP1,W) - CALL DCOST (NP1,X,W) - COSTT = 0.0D0 - DO 120 I=1,NP1 - COSTT = DMAX1(COSTT,ABS(X(I)-Y(I))) - X(I) = XH(I) - Y(I) = XH(I) - 120 CONTINUE - COSTT = CF*COSTT - CALL DCOST (NP1,X,W) - CALL DCOST (NP1,X,W) - COSTFB = 0.0D0 - DO 121 I=1,NP1 - COSTFB = DMAX1(COSTFB,ABS(CF*X(I)-Y(I))) - 121 CONTINUE -C -C TEST SUBROUTINES SINQI,SINQF AND SINQB -C - CF = 0.25D0/FN - DO 122 I=1,N - Y(I) = XH(I) - 122 CONTINUE - DT = PI/(FN+FN) - DO 124 I=1,N - X(I) = 0.0D0 - ARG = DT*REAL(I,RK) - DO 123 K=1,N - X(I) = X(I)+Y(K)*SIN(REAL(K+K-1,RK)*ARG) - 123 CONTINUE - X(I) = 4.0D0*X(I) - 124 CONTINUE - CALL DSINQI (N,W) - CALL DSINQB (N,Y,W) - SINQBT = 0.0D0 - DO 125 I=1,N - SINQBT = DMAX1(SINQBT,ABS(Y(I)-X(I))) - X(I) = XH(I) - 125 CONTINUE - SINQBT = CF*SINQBT - DO 127 I=1,N - ARG = REAL(I+I-1,RK)*DT - Y(I) = 0.5D0*REAL((-1)**(I+1),RK)*X(N) - DO 126 K=1,NM1 - Y(I) = Y(I)+X(K)*SIN(REAL(K,RK)*ARG) - 126 CONTINUE - Y(I) = Y(I)+Y(I) - 127 CONTINUE - CALL DSINQF (N,X,W) - SINQFT = 0.0D0 - DO 128 I=1,N - SINQFT = DMAX1(SINQFT,ABS(X(I)-Y(I))) - Y(I) = XH(I) - X(I) = XH(I) - 128 CONTINUE - CALL DSINQF (N,Y,W) - CALL DSINQB (N,Y,W) - SINQFB = 0.0D0 - DO 129 I=1,N - SINQFB = DMAX1(SINQFB,ABS(CF*Y(I)-X(I))) - 129 CONTINUE -C -C TEST SUBROUTINES COSQI,COSQF AND COSQB -C - DO 130 I=1,N - Y(I) = XH(I) - 130 CONTINUE - DO 132 I=1,N - X(I) = 0.0D0 - ARG = REAL(I-1,RK)*DT - DO 131 K=1,N - X(I) = X(I)+Y(K)*COS(REAL(K+K-1,RK)*ARG) - 131 CONTINUE - X(I) = 4.0D0*X(I) - 132 CONTINUE - CALL DCOSQI (N,W) - CALL DCOSQB (N,Y,W) - COSQBT = 0.0D0 - DO 133 I=1,N - COSQBT = DMAX1(COSQBT,ABS(X(I)-Y(I))) - X(I) = XH(I) - 133 CONTINUE - COSQBT = CF*COSQBT - DO 135 I=1,N - Y(I) = 0.5D0*X(1) - ARG = REAL(I+I-1,RK)*DT - DO 134 K=2,N - Y(I) = Y(I)+X(K)*COS(REAL(K-1,RK)*ARG) - 134 CONTINUE - Y(I) = Y(I)+Y(I) - 135 CONTINUE - CALL DCOSQF (N,X,W) - COSQFT = 0.0D0 - DO 136 I=1,N - COSQFT = DMAX1(COSQFT,ABS(Y(I)-X(I))) - X(I) = XH(I) - Y(I) = XH(I) - 136 CONTINUE - COSQFT = CF*COSQFT - CALL DCOSQB (N,X,W) - CALL DCOSQF (N,X,W) - COSQFB = 0.0D0 - DO 137 I=1,N - COSQFB = DMAX1(COSQFB,ABS(CF*X(I)-Y(I))) - 137 CONTINUE -C -C TEST PROGRAMS EZFFTI,EZFFTF,EZFFTB -C - CALL DZFFTI(N,W) - DO 138 I=1,N - X(I) = XH(I) - 138 CONTINUE - TPI = 8.0D0*ATAN(1.0D0) - DT = TPI/REAL(N,RK) - NS2 = (N+1)/2 - CF = 2.0D0/REAL(N,RK) - NS2M = NS2-1 - IF (NS2M .LE. 0) GO TO 141 - DO 140 K=1,NS2M - SUM1 = 0.0D0 - SUM2 = 0.0D0 - ARG = REAL(K,RK)*DT - DO 139 I=1,N - ARG1 = REAL(I-1,RK)*ARG - SUM1 = SUM1+X(I)*COS(ARG1) - SUM2 = SUM2+X(I)*SIN(ARG1) - 139 CONTINUE - A(K) = CF*SUM1 - B(K) = CF*SUM2 - 140 CONTINUE - 141 NM1 = N-1 - SUM1 = 0.0D0 - SUM2 = 0.0D0 - DO 142 I=1,NM1,2 - SUM1 = SUM1+X(I) - SUM2 = SUM2+X(I+1) - 142 CONTINUE - IF (MODN .EQ. 1) SUM1 = SUM1+X(N) - AZERO = 0.5D0*CF*(SUM1+SUM2) - IF (MODN .EQ. 0) A(NS2) = 0.5D0*CF*(SUM1-SUM2) - CALL DZFFTF (N,X,AZEROH,AH,BH,W) - DEZF1 = ABS(AZEROH-AZERO) - IF (MODN .EQ. 0) DEZF1 = DMAX1(DEZF1,ABS(A(NS2)-AH(NS2))) - IF (NS2M .LE. 0) GO TO 144 - DO 143 I=1,NS2M - DEZF1 = DMAX1(DEZF1,ABS(AH(I)-A(I)),ABS(BH(I)-B(I))) - 143 CONTINUE - 144 NS2 = N/2 - IF (MODN .EQ. 0) B(NS2) = 0.0D0 - DO 146 I=1,N - SUM = AZERO - ARG1 = REAL(I-1,RK)*DT - DO 145 K=1,NS2 - ARG2 = REAL(K,RK)*ARG1 - SUM = SUM+A(K)*COS(ARG2)+B(K)*SIN(ARG2) - 145 CONTINUE - X(I) = SUM - 146 CONTINUE - CALL DZFFTB (N,Y,AZERO,A,B,W) - DEZB1 = 0.0D0 - DO 147 I=1,N - DEZB1 = DMAX1(DEZB1,ABS(X(I)-Y(I))) - X(I) = XH(I) - 147 CONTINUE - CALL DZFFTF (N,X,AZERO,A,B,W) - CALL DZFFTB (N,Y,AZERO,A,B,W) - DEZFB = 0.0D0 - DO 148 I=1,N - DEZFB = DMAX1(DEZFB,ABS(X(I)-Y(I))) - 148 CONTINUE -C -C TEST CFFTI,CFFTF,CFFTB -C - DO 149 I=1,N - CX(I) =DCMPLX(COS(SQRT2*REAL(I,RK)),SIN(SQRT2*REAL(I*I,RK))) - 149 CONTINUE - DT = (PI+PI)/FN - DO 151 I=1,N - ARG1 = -REAL(I-1,RK)*DT - CY(I) = (0.0D0,0.0D0) - DO 150 K=1,N - ARG2 = REAL(K-1,RK)*ARG1 - CY(I) = CY(I)+DCMPLX(COS(ARG2),SIN(ARG2))*CX(K) - 150 CONTINUE - 151 CONTINUE - CALL ZFFTI (N,W) - CALL ZFFTF (N,CX,W) - DCFFTF = 0.0D0 - DO 152 I=1,N - DCFFTF = DMAX1(DCFFTF,ABS(CX(I)-CY(I))) - CX(I) = CX(I)/FN - 152 CONTINUE - DCFFTF = DCFFTF/FN - DO 154 I=1,N - ARG1 = REAL(I-1,RK)*DT - CY(I) = (0.0D0,0.0D0) - DO 153 K=1,N - ARG2 = REAL(K-1,RK)*ARG1 - CY(I) = CY(I)+DCMPLX(COS(ARG2),SIN(ARG2))*CX(K) - 153 CONTINUE - 154 CONTINUE - CALL ZFFTB (N,CX,W) - DCFFTB = 0.0D0 - DO 155 I=1,N - DCFFTB = DMAX1(DCFFTB,ABS(CX(I)-CY(I))) - CX(I) = CY(I) - 155 CONTINUE - CF = 1.0D0/FN - CALL ZFFTF (N,CX,W) - CALL ZFFTB (N,CX,W) - DCFB = 0.0D0 - DO 156 I=1,N - DCFB = DMAX1(DCFB,ABS(CF*CX(I)-CY(I))) - 156 CONTINUE - WRITE (6,1001) N,RFTF,RFTB,RFTFB,SINTT,SINTFB,COSTT,COSTFB, - 1 SINQFT,SINQBT,SINQFB,COSQFT,COSQBT,COSQFB,DEZF1, - 2 DEZB1,DEZFB,DCFFTF,DCFFTB,DCFB - 157 CONTINUE -C -C -C - 1001 FORMAT (2H0N,I5,8H RFFTF ,E10.3,8H RFFTB ,E10.3,8H RFFTFB , - 1 E10.3,8H SINT ,E10.3,8H SINTFB ,E10.3,8H COST ,E10.3/ - 2 7X,8H COSTFB ,E10.3,8H SINQF ,E10.3,8H SINQB ,E10.3, - 3 8H SINQFB ,E10.3,8H COSQF ,E10.3,8H COSQB ,E10.3/7X, - 4 8H COSQFB ,E10.3,8H DEZF ,E10.3,8H DEZB ,E10.3, - 5 8H DEZFB ,E10.3,8H CFFTF ,E10.3,8H CFFTB ,E10.3/ - 6 7X,8H CFFTFB ,E10.3) -C - END