From 8b59c6181c758bd133ae54725931cbe2f25d2696 Mon Sep 17 00:00:00 2001 From: "Harris M. Snyder" Date: Sat, 11 Sep 2021 23:09:32 -0400 Subject: [PATCH 1/3] the kind used by fftpack can now be chosen at build time by modifying the defintion of 'rk' in 'src/rk.f' --- src/cfftb1.f | 3 +- src/cfftf1.f | 3 +- src/cffti1.f | 3 +- src/cosqb1.f | 3 +- src/cosqf1.f | 3 +- src/dcosqb.f | 3 +- src/dcosqf.f | 3 +- src/dcosqi.f | 3 +- src/dcost.f | 3 +- src/dcosti.f | 3 +- src/dfftb.f | 3 +- src/dfftf.f | 3 +- src/dffti.f | 3 +- src/dsinqb.f | 3 +- src/dsinqf.f | 3 +- src/dsinqi.f | 3 +- src/dsint.f | 3 +- src/dsinti.f | 3 +- src/dzfftb.f | 3 +- src/dzfftf.f | 3 +- src/dzffti.f | 3 +- src/ezfft1.f | 3 +- src/fftpack.f90 | 155 ++++++++++++++++---------------- src/fftpack_fft.f90 | 12 +-- src/fftpack_fftshift.f90 | 20 ++--- src/fftpack_ifft.f90 | 12 +-- src/fftpack_ifftshift.f90 | 20 ++--- src/fftpack_iqct.f90 | 12 +-- src/fftpack_irfft.f90 | 12 +-- src/fftpack_qct.f90 | 12 +-- src/fftpack_rfft.f90 | 12 +-- src/passb.f | 3 +- src/passb2.f | 3 +- src/passb3.f | 3 +- src/passb4.f | 3 +- src/passb5.f | 3 +- src/passf.f | 3 +- src/passf2.f | 3 +- src/passf3.f | 3 +- src/passf4.f | 3 +- src/passf5.f | 3 +- src/radb2.f | 3 +- src/radb3.f | 3 +- src/radb4.f | 3 +- src/radb5.f | 3 +- src/radbg.f | 3 +- src/radf2.f | 3 +- src/radf3.f | 3 +- src/radf4.f | 3 +- src/radf5.f | 3 +- src/radfg.f | 3 +- src/rfftb1.f | 3 +- src/rfftf1.f | 3 +- src/rffti1.f | 3 +- src/rk.f | 4 + src/sint1.f | 3 +- src/zfftb.f | 3 +- src/zfftf.f | 3 +- src/zffti.f | 3 +- test/test_fftpack_dcosq.f90 | 15 ++-- test/test_fftpack_dfft.f90 | 11 +-- test/test_fftpack_dzfft.f90 | 17 ++-- test/test_fftpack_fft.f90 | 11 +-- test/test_fftpack_fftshift.f90 | 30 ++++--- test/test_fftpack_ifft.f90 | 15 ++-- test/test_fftpack_ifftshift.f90 | 30 ++++--- test/test_fftpack_iqct.f90 | 25 +++--- test/test_fftpack_irfft.f90 | 15 ++-- test/test_fftpack_qct.f90 | 13 +-- test/test_fftpack_rfft.f90 | 11 +-- test/test_fftpack_zfft.f90 | 14 +-- 71 files changed, 347 insertions(+), 278 deletions(-) create mode 100644 src/rk.f diff --git a/src/cfftb1.f b/src/cfftb1.f index 9faebe2..bd2a2c4 100644 --- a/src/cfftb1.f +++ b/src/cfftb1.f @@ -1,5 +1,6 @@ SUBROUTINE CFFTB1 (N,C,CH,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) NF = IFAC(2) NA = 0 diff --git a/src/cfftf1.f b/src/cfftf1.f index 2eacd2a..ac847b6 100644 --- a/src/cfftf1.f +++ b/src/cfftf1.f @@ -1,5 +1,6 @@ SUBROUTINE CFFTF1 (N,C,CH,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) NF = IFAC(2) NA = 0 diff --git a/src/cffti1.f b/src/cffti1.f index e1f45e9..923cdfe 100644 --- a/src/cffti1.f +++ b/src/cffti1.f @@ -1,5 +1,6 @@ SUBROUTINE CFFTI1 (N,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/ NL = N diff --git a/src/cosqb1.f b/src/cosqb1.f index 5c8b6ec..1515d9e 100644 --- a/src/cosqb1.f +++ b/src/cosqb1.f @@ -1,5 +1,6 @@ SUBROUTINE COSQB1 (N,X,W,XH) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(1) ,W(1) ,XH(1) NS2 = (N+1)/2 NP2 = N+2 diff --git a/src/cosqf1.f b/src/cosqf1.f index 87a0612..409cb95 100644 --- a/src/cosqf1.f +++ b/src/cosqf1.f @@ -1,5 +1,6 @@ SUBROUTINE COSQF1 (N,X,W,XH) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(1) ,W(1) ,XH(1) NS2 = (N+1)/2 NP2 = N+2 diff --git a/src/dcosqb.f b/src/dcosqb.f index 341d8c3..e054c6c 100644 --- a/src/dcosqb.f +++ b/src/dcosqb.f @@ -1,5 +1,6 @@ SUBROUTINE DCOSQB (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(*) ,WSAVE(*) DATA TSQRT2 /2.82842712474619009760D0/ IF (N-2) 101,102,103 diff --git a/src/dcosqf.f b/src/dcosqf.f index ca585bb..911657e 100644 --- a/src/dcosqf.f +++ b/src/dcosqf.f @@ -1,5 +1,6 @@ SUBROUTINE DCOSQF (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(*) ,WSAVE(*) DATA SQRT2 /1.41421356237309504880D0/ IF (N-2) 102,101,103 diff --git a/src/dcosqi.f b/src/dcosqi.f index 215e5c1..9c09127 100644 --- a/src/dcosqi.f +++ b/src/dcosqi.f @@ -1,5 +1,6 @@ SUBROUTINE DCOSQI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) DATA PIH /1.57079632679489661923D0/ DT = PIH/FLOAT(N) diff --git a/src/dcost.f b/src/dcost.f index fbf90f2..77544ab 100644 --- a/src/dcost.f +++ b/src/dcost.f @@ -1,5 +1,6 @@ SUBROUTINE DCOST (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(*) ,WSAVE(*) NM1 = N-1 NP1 = N+1 diff --git a/src/dcosti.f b/src/dcosti.f index 4c6d7bb..4993339 100644 --- a/src/dcosti.f +++ b/src/dcosti.f @@ -1,5 +1,6 @@ SUBROUTINE DCOSTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) DATA PI /3.14159265358979323846D0/ IF (N .LE. 3) RETURN diff --git a/src/dfftb.f b/src/dfftb.f index e4adfca..dabef1a 100644 --- a/src/dfftb.f +++ b/src/dfftb.f @@ -1,5 +1,6 @@ SUBROUTINE DFFTB (N,R,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION R(1) ,WSAVE(1) IF (N .EQ. 1) RETURN CALL RFFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1)) diff --git a/src/dfftf.f b/src/dfftf.f index 7511a66..49dae32 100644 --- a/src/dfftf.f +++ b/src/dfftf.f @@ -1,5 +1,6 @@ SUBROUTINE DFFTF (N,R,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION R(1) ,WSAVE(1) IF (N .EQ. 1) RETURN CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1)) diff --git a/src/dffti.f b/src/dffti.f index d70250c..1166966 100644 --- a/src/dffti.f +++ b/src/dffti.f @@ -1,5 +1,6 @@ SUBROUTINE DFFTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) IF (N .EQ. 1) RETURN CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1)) diff --git a/src/dsinqb.f b/src/dsinqb.f index d8d7835..ae11841 100644 --- a/src/dsinqb.f +++ b/src/dsinqb.f @@ -1,5 +1,6 @@ SUBROUTINE DSINQB (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(1) ,WSAVE(1) IF (N .GT. 1) GO TO 101 X(1) = 4.0D0*X(1) diff --git a/src/dsinqf.f b/src/dsinqf.f index 1b257ea..7ff4de7 100644 --- a/src/dsinqf.f +++ b/src/dsinqf.f @@ -1,5 +1,6 @@ SUBROUTINE DSINQF (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(1) ,WSAVE(1) IF (N .EQ. 1) RETURN NS2 = N/2 diff --git a/src/dsinqi.f b/src/dsinqi.f index c4897c2..0aa3beb 100644 --- a/src/dsinqi.f +++ b/src/dsinqi.f @@ -1,5 +1,6 @@ SUBROUTINE DSINQI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) CALL DCOSQI (N,WSAVE) RETURN diff --git a/src/dsint.f b/src/dsint.f index feae4e0..12313f1 100644 --- a/src/dsint.f +++ b/src/dsint.f @@ -1,5 +1,6 @@ SUBROUTINE DSINT (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION X(1) ,WSAVE(1) NP1 = N+1 IW1 = N/2+1 diff --git a/src/dsinti.f b/src/dsinti.f index e655c0a..fbb283c 100644 --- a/src/dsinti.f +++ b/src/dsinti.f @@ -1,5 +1,6 @@ SUBROUTINE DSINTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) DATA PI /3.14159265358979323846D0/ IF (N .LE. 1) RETURN diff --git a/src/dzfftb.f b/src/dzfftb.f index 07ced9b..045f59b 100644 --- a/src/dzfftb.f +++ b/src/dzfftb.f @@ -1,5 +1,6 @@ SUBROUTINE DZFFTB (N,R,AZERO,A,B,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION R(*) ,A(*) ,B(*) ,WSAVE(*) IF (N-2) 101,102,103 101 R(1) = AZERO diff --git a/src/dzfftf.f b/src/dzfftf.f index f4c86de..13b4f6b 100644 --- a/src/dzfftf.f +++ b/src/dzfftf.f @@ -1,8 +1,9 @@ SUBROUTINE DZFFTF (N,R,AZERO,A,B,WSAVE) + USE fftpack_kind C C VERSION 3 JUNE 1979 C - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION R(*) ,A(*) ,B(*) ,WSAVE(*) IF (N-2) 101,102,103 101 AZERO = R(1) diff --git a/src/dzffti.f b/src/dzffti.f index 8985149..c0a6f69 100644 --- a/src/dzffti.f +++ b/src/dzffti.f @@ -1,5 +1,6 @@ SUBROUTINE DZFFTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) IF (N .EQ. 1) RETURN CALL EZFFT1 (N,WSAVE(2*N+1),WSAVE(3*N+1)) diff --git a/src/ezfft1.f b/src/ezfft1.f index 230944a..2600c89 100644 --- a/src/ezfft1.f +++ b/src/ezfft1.f @@ -1,5 +1,6 @@ SUBROUTINE EZFFT1 (N,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ 1 ,TPI/6.28318530717958647692D0/ diff --git a/src/fftpack.f90 b/src/fftpack.f90 index 6a3b9f7..ae54319 100644 --- a/src/fftpack.f90 +++ b/src/fftpack.f90 @@ -1,10 +1,9 @@ module fftpack + use fftpack_kind implicit none private - integer, parameter :: dp = kind(1.0d0) - public :: dp public :: zffti, zfftf, zfftb public :: fft, ifft public :: fftshift, ifftshift @@ -24,9 +23,9 @@ module fftpack !> Initialize `zfftf` and `zfftb`. !> ([Specification](../page/specs/fftpack.html#zffti)) pure subroutine zffti(n, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(out) :: wsave(*) + real(kind=rk), intent(out) :: wsave(*) end subroutine zffti !> Version: experimental @@ -34,10 +33,10 @@ end subroutine zffti !> Forward transform of a double complex periodic sequence. !> ([Specification](../page/specs/fftpack.html#zfftf)) pure subroutine zfftf(n, c, wsave) - import dp + import rk integer, intent(in) :: n - complex(kind=dp), intent(inout) :: c(*) - real(kind=dp), intent(in) :: wsave(*) + complex(kind=rk), intent(inout) :: c(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine zfftf !> Version: experimental @@ -45,10 +44,10 @@ end subroutine zfftf !> Unnormalized inverse of `zfftf`. !> ([Specification](../page/specs/fftpack.html#zfftb)) pure subroutine zfftb(n, c, wsave) - import dp + import rk integer, intent(in) :: n - complex(kind=dp), intent(inout) :: c(*) - real(kind=dp), intent(in) :: wsave(*) + complex(kind=rk), intent(inout) :: c(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine zfftb !> Version: experimental @@ -56,9 +55,9 @@ end subroutine zfftb !> Initialize `dfftf` and `dfftb`. !> ([Specification](../page/specs/fftpack.html#dffti)) pure subroutine dffti(n, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(out) :: wsave(*) + real(kind=rk), intent(out) :: wsave(*) end subroutine dffti !> Version: experimental @@ -66,10 +65,10 @@ end subroutine dffti !> Forward transform of a double real periodic sequence. !> ([Specification](../page/specs/fftpack.html#dfftf)) pure subroutine dfftf(n, r, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(inout) :: r(*) - real(kind=dp), intent(in) :: wsave(*) + real(kind=rk), intent(inout) :: r(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine dfftf !> Version: experimental @@ -77,10 +76,10 @@ end subroutine dfftf !> Unnormalized inverse of `dfftf`. !> ([Specification](../page/specs/fftpack.html#dfftb)) pure subroutine dfftb(n, r, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(inout) :: r(*) - real(kind=dp), intent(in) :: wsave(*) + real(kind=rk), intent(inout) :: r(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine dfftb !> Version: experimental @@ -88,9 +87,9 @@ end subroutine dfftb !> Initialize `dzfftf` and `dzfftb`. !> ([Specification](../page/specs/fftpack.html#dzffti)) pure subroutine dzffti(n, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(out) :: wsave(*) + real(kind=rk), intent(out) :: wsave(*) end subroutine dzffti !> Version: experimental @@ -98,12 +97,12 @@ end subroutine dzffti !> Simplified forward transform of a double real periodic sequence. !> ([Specification](../page/specs/fftpack.html#dzfftf)) pure subroutine dzfftf(n, r, azero, a, b, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(in) :: r(*) - real(kind=dp), intent(out) :: azero - real(kind=dp), intent(out) :: a(*), b(*) - real(kind=dp), intent(in) :: wsave(*) + real(kind=rk), intent(in) :: r(*) + real(kind=rk), intent(out) :: azero + real(kind=rk), intent(out) :: a(*), b(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine dzfftf !> Version: experimental @@ -111,12 +110,12 @@ end subroutine dzfftf !> Unnormalized inverse of `dzfftf`. !> ([Specification](../page/specs/fftpack.html#dzfftb)) pure subroutine dzfftb(n, r, azero, a, b, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(out) :: r(*) - real(kind=dp), intent(in) :: azero - real(kind=dp), intent(in) :: a(*), b(*) - real(kind=dp), intent(in) :: wsave(*) + real(kind=rk), intent(out) :: r(*) + real(kind=rk), intent(in) :: azero + real(kind=rk), intent(in) :: a(*), b(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine dzfftb !> Version: experimental @@ -124,9 +123,9 @@ end subroutine dzfftb !> Initialize `dcosqf` and `dcosqb`. !> ([Specification](../page/specs/fftpack.html#dcosqi)) pure subroutine dcosqi(n, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(out) :: wsave(*) + real(kind=rk), intent(out) :: wsave(*) end subroutine dcosqi !> Version: experimental @@ -134,10 +133,10 @@ end subroutine dcosqi !> Forward transform of quarter wave data. !> ([Specification](../page/specs/fftpack.html#dcosqf)) pure subroutine dcosqf(n, x, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(inout) :: x(*) - real(kind=dp), intent(in) :: wsave(*) + real(kind=rk), intent(inout) :: x(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine dcosqf !> Version: experimental @@ -145,10 +144,10 @@ end subroutine dcosqf !> Unnormalized inverse of `dcosqf`. !> ([Specification](../page/specs/fftpack.html#dcosqb)) pure subroutine dcosqb(n, x, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(inout) :: x(*) - real(kind=dp), intent(in) :: wsave(*) + real(kind=rk), intent(inout) :: x(*) + real(kind=rk), intent(in) :: wsave(*) end subroutine dcosqb end interface @@ -158,11 +157,11 @@ end subroutine dcosqb !> Forward transform of a double complex periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#fft)) interface fft - pure module function fft_dp(x, n) result(result) - complex(kind=dp), intent(in) :: x(:) + pure module function fft_rk(x, n) result(result) + complex(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - complex(kind=dp), allocatable :: result(:) - end function fft_dp + complex(kind=rk), allocatable :: result(:) + end function fft_rk end interface fft !> Version: experimental @@ -170,11 +169,11 @@ end function fft_dp !> Backward transform of a double complex periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#ifft)) interface ifft - pure module function ifft_dp(x, n) result(result) - complex(kind=dp), intent(in) :: x(:) + pure module function ifft_rk(x, n) result(result) + complex(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - complex(kind=dp), allocatable :: result(:) - end function ifft_dp + complex(kind=rk), allocatable :: result(:) + end function ifft_rk end interface ifft !> Version: experimental @@ -182,11 +181,11 @@ end function ifft_dp !> Forward transform of a double real periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#rfft)) interface rfft - pure module function rfft_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function rfft_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) - end function rfft_dp + real(kind=rk), allocatable :: result(:) + end function rfft_rk end interface rfft !> Version: experimental @@ -194,11 +193,11 @@ end function rfft_dp !> Backward transform of a double real periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#irfft)) interface irfft - pure module function irfft_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function irfft_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) - end function irfft_dp + real(kind=rk), allocatable :: result(:) + end function irfft_rk end interface irfft !> Version: experimental @@ -206,11 +205,11 @@ end function irfft_dp !> Forward transform of quarter wave data. !> ([Specifiction](../page/specs/fftpack.html#qct)) interface qct - pure module function qct_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function qct_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) - end function qct_dp + real(kind=rk), allocatable :: result(:) + end function qct_rk end interface qct !> Version: experimental @@ -218,11 +217,11 @@ end function qct_dp !> Backward transform of quarter wave data. !> ([Specifiction](../page/specs/fftpack.html#iqct)) interface iqct - pure module function iqct_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function iqct_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) - end function iqct_dp + real(kind=rk), allocatable :: result(:) + end function iqct_rk end interface iqct !> Version: experimental @@ -230,14 +229,14 @@ end function iqct_dp !> Shifts zero-frequency component to center of spectrum. !> ([Specifiction](../page/specs/fftpack.html#fftshift)) interface fftshift - pure module function fftshift_cdp(x) result(result) - complex(kind=dp), intent(in) :: x(:) - complex(kind=dp), allocatable :: result(:) - end function fftshift_cdp - pure module function fftshift_rdp(x) result(result) - real(kind=dp), intent(in) :: x(:) - real(kind=dp), allocatable :: result(:) - end function fftshift_rdp + pure module function fftshift_crk(x) result(result) + complex(kind=rk), intent(in) :: x(:) + complex(kind=rk), allocatable :: result(:) + end function fftshift_crk + pure module function fftshift_rrk(x) result(result) + real(kind=rk), intent(in) :: x(:) + real(kind=rk), allocatable :: result(:) + end function fftshift_rrk end interface fftshift !> Version: experimental @@ -245,14 +244,14 @@ end function fftshift_rdp !> Shifts zero-frequency component to beginning of spectrum. !> ([Specifiction](../page/specs/fftpack.html#ifftshift)) interface ifftshift - pure module function ifftshift_cdp(x) result(result) - complex(kind=dp), intent(in) :: x(:) - complex(kind=dp), allocatable :: result(:) - end function ifftshift_cdp - pure module function ifftshift_rdp(x) result(result) - real(kind=dp), intent(in) :: x(:) - real(kind=dp), allocatable :: result(:) - end function ifftshift_rdp + pure module function ifftshift_crk(x) result(result) + complex(kind=rk), intent(in) :: x(:) + complex(kind=rk), allocatable :: result(:) + end function ifftshift_crk + pure module function ifftshift_rrk(x) result(result) + real(kind=rk), intent(in) :: x(:) + real(kind=rk), allocatable :: result(:) + end function ifftshift_rrk end interface ifftshift end module fftpack diff --git a/src/fftpack_fft.f90 b/src/fftpack_fft.f90 index b6e362f..bf21eb5 100644 --- a/src/fftpack_fft.f90 +++ b/src/fftpack_fft.f90 @@ -3,20 +3,20 @@ contains !> Forward transform of a double complex periodic sequence. - pure module function fft_dp(x, n) result(result) - complex(kind=dp), intent(in) :: x(:) + pure module function fft_rk(x, n) result(result) + complex(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - complex(kind=dp), allocatable :: result(:) + complex(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i - real(kind=dp), allocatable :: wsave(:) + real(kind=rk), allocatable :: wsave(:) if (present(n)) then lenseq = n if (lenseq <= size(x)) then result = x(:lenseq) else if (lenseq > size(x)) then - result = [x, ((0.0_dp, 0.0_dp), i=1, lenseq - size(x))] + result = [x, ((0.0_rk, 0.0_rk), i=1, lenseq - size(x))] end if else lenseq = size(x) @@ -31,6 +31,6 @@ pure module function fft_dp(x, n) result(result) !> Forward transformation call zfftf(lenseq, result, wsave) - end function fft_dp + end function fft_rk end submodule fftpack_fft diff --git a/src/fftpack_fftshift.f90 b/src/fftpack_fftshift.f90 index 8cb937b..a22dfc6 100644 --- a/src/fftpack_fftshift.f90 +++ b/src/fftpack_fftshift.f90 @@ -3,21 +3,21 @@ contains !> Shifts zero-frequency component to center of spectrum for `complex` type. - pure module function fftshift_cdp(x) result(result) - complex(kind=dp), intent(in) :: x(:) - complex(kind=dp), allocatable :: result(:) + pure module function fftshift_crk(x) result(result) + complex(kind=rk), intent(in) :: x(:) + complex(kind=rk), allocatable :: result(:) - result = cshift(x, shift=-floor(0.5_dp*size(x))) + result = cshift(x, shift=-floor(0.5_rk*size(x))) - end function fftshift_cdp + end function fftshift_crk !> Shifts zero-frequency component to center of spectrum for `real` type. - pure module function fftshift_rdp(x) result(result) - real(kind=dp), intent(in) :: x(:) - real(kind=dp), allocatable :: result(:) + pure module function fftshift_rrk(x) result(result) + real(kind=rk), intent(in) :: x(:) + real(kind=rk), allocatable :: result(:) - result = cshift(x, shift=-floor(0.5_dp*size(x))) + result = cshift(x, shift=-floor(0.5_rk*size(x))) - end function fftshift_rdp + end function fftshift_rrk end submodule fftpack_fftshift diff --git a/src/fftpack_ifft.f90 b/src/fftpack_ifft.f90 index 9e6dea0..40068d4 100644 --- a/src/fftpack_ifft.f90 +++ b/src/fftpack_ifft.f90 @@ -3,20 +3,20 @@ contains !> Backward transform of a double complex periodic sequence. - pure module function ifft_dp(x, n) result(result) - complex(kind=dp), intent(in) :: x(:) + pure module function ifft_rk(x, n) result(result) + complex(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - complex(kind=dp), allocatable :: result(:) + complex(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i - real(kind=dp), allocatable :: wsave(:) + real(kind=rk), allocatable :: wsave(:) if (present(n)) then lenseq = n if (lenseq <= size(x)) then result = x(:lenseq) else if (lenseq > size(x)) then - result = [x, ((0.0_dp, 0.0_dp), i=1, lenseq - size(x))] + result = [x, ((0.0_rk, 0.0_rk), i=1, lenseq - size(x))] end if else lenseq = size(x) @@ -31,6 +31,6 @@ pure module function ifft_dp(x, n) result(result) !> Backward transformation call zfftb(lenseq, result, wsave) - end function ifft_dp + end function ifft_rk end submodule fftpack_ifft diff --git a/src/fftpack_ifftshift.f90 b/src/fftpack_ifftshift.f90 index da8132a..49830a7 100644 --- a/src/fftpack_ifftshift.f90 +++ b/src/fftpack_ifftshift.f90 @@ -3,21 +3,21 @@ contains !> Shifts zero-frequency component to beginning of spectrum for `complex` type. - pure module function ifftshift_cdp(x) result(result) - complex(kind=dp), intent(in) :: x(:) - complex(kind=dp), allocatable :: result(:) + pure module function ifftshift_crk(x) result(result) + complex(kind=rk), intent(in) :: x(:) + complex(kind=rk), allocatable :: result(:) - result = cshift(x, shift=-ceiling(0.5_dp*size(x))) + result = cshift(x, shift=-ceiling(0.5_rk*size(x))) - end function ifftshift_cdp + end function ifftshift_crk !> Shifts zero-frequency component to beginning of spectrum for `real` type. - pure module function ifftshift_rdp(x) result(result) - real(kind=dp), intent(in) :: x(:) - real(kind=dp), allocatable :: result(:) + pure module function ifftshift_rrk(x) result(result) + real(kind=rk), intent(in) :: x(:) + real(kind=rk), allocatable :: result(:) - result = cshift(x, shift=-ceiling(0.5_dp*size(x))) + result = cshift(x, shift=-ceiling(0.5_rk*size(x))) - end function ifftshift_rdp + end function ifftshift_rrk end submodule fftpack_ifftshift diff --git a/src/fftpack_iqct.f90 b/src/fftpack_iqct.f90 index f116fef..1ae2a20 100644 --- a/src/fftpack_iqct.f90 +++ b/src/fftpack_iqct.f90 @@ -3,20 +3,20 @@ contains !> Backward transform of quarter wave data. - pure module function iqct_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function iqct_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) + real(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i - real(kind=dp), allocatable :: wsave(:) + real(kind=rk), allocatable :: wsave(:) if (present(n)) then lenseq = n if (lenseq <= size(x)) then result = x(:lenseq) else if (lenseq > size(x)) then - result = [x, (0.0_dp, i=1, lenseq - size(x))] + result = [x, (0.0_rk, i=1, lenseq - size(x))] end if else lenseq = size(x) @@ -31,6 +31,6 @@ pure module function iqct_dp(x, n) result(result) !> Backward transformation call dcosqb(lenseq, result, wsave) - end function iqct_dp + end function iqct_rk end submodule fftpack_iqct diff --git a/src/fftpack_irfft.f90 b/src/fftpack_irfft.f90 index eeb6d34..97565cb 100644 --- a/src/fftpack_irfft.f90 +++ b/src/fftpack_irfft.f90 @@ -3,20 +3,20 @@ contains !> Backward transform of a double real periodic sequence. - pure module function irfft_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function irfft_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) + real(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i - real(kind=dp), allocatable :: wsave(:) + real(kind=rk), allocatable :: wsave(:) if (present(n)) then lenseq = n if (lenseq <= size(x)) then result = x(:lenseq) else if (lenseq > size(x)) then - result = [x, (0.0_dp, i=1, lenseq - size(x))] + result = [x, (0.0_rk, i=1, lenseq - size(x))] end if else lenseq = size(x) @@ -31,6 +31,6 @@ pure module function irfft_dp(x, n) result(result) !> Backward transformation call dfftb(lenseq, result, wsave) - end function irfft_dp + end function irfft_rk end submodule fftpack_irfft diff --git a/src/fftpack_qct.f90 b/src/fftpack_qct.f90 index 582918b..ceb2f6b 100644 --- a/src/fftpack_qct.f90 +++ b/src/fftpack_qct.f90 @@ -3,20 +3,20 @@ contains !> Forward transform of quarter wave data. - pure module function qct_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function qct_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) + real(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i - real(kind=dp), allocatable :: wsave(:) + real(kind=rk), allocatable :: wsave(:) if (present(n)) then lenseq = n if (lenseq <= size(x)) then result = x(:lenseq) else if (lenseq > size(x)) then - result = [x, (0.0_dp, i=1, lenseq - size(x))] + result = [x, (0.0_rk, i=1, lenseq - size(x))] end if else lenseq = size(x) @@ -31,6 +31,6 @@ pure module function qct_dp(x, n) result(result) !> Forward transformation call dcosqf(lenseq, result, wsave) - end function qct_dp + end function qct_rk end submodule fftpack_qct diff --git a/src/fftpack_rfft.f90 b/src/fftpack_rfft.f90 index 6cbcdac..a0c6260 100644 --- a/src/fftpack_rfft.f90 +++ b/src/fftpack_rfft.f90 @@ -3,20 +3,20 @@ contains !> Forward transform of a double real periodic sequence. - pure module function rfft_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function rfft_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) + real(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i - real(kind=dp), allocatable :: wsave(:) + real(kind=rk), allocatable :: wsave(:) if (present(n)) then lenseq = n if (lenseq <= size(x)) then result = x(:lenseq) else if (lenseq > size(x)) then - result = [x, (0.0_dp, i=1, lenseq - size(x))] + result = [x, (0.0_rk, i=1, lenseq - size(x))] end if else lenseq = size(x) @@ -31,6 +31,6 @@ pure module function rfft_dp(x, n) result(result) !> Forward transformation call dfftf(lenseq, result, wsave) - end function rfft_dp + end function rfft_rk end submodule fftpack_rfft diff --git a/src/passb.f b/src/passb.f index 43dba87..f347d5c 100644 --- a/src/passb.f +++ b/src/passb.f @@ -1,5 +1,6 @@ SUBROUTINE PASSB (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,WA(1) ,C2(IDL1,IP), 2 CH2(IDL1,IP) diff --git a/src/passb2.f b/src/passb2.f index 22f6279..7fcdc45 100644 --- a/src/passb2.f +++ b/src/passb2.f @@ -1,5 +1,6 @@ SUBROUTINE PASSB2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1 WA1(1) IF (IDO .GT. 2) GO TO 102 diff --git a/src/passb3.f b/src/passb3.f index ee8ef10..613b528 100644 --- a/src/passb3.f +++ b/src/passb3.f @@ -1,5 +1,6 @@ SUBROUTINE PASSB3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1 WA1(1) ,WA2(1) C *** TAUI IS SQRT(3)/2 *** diff --git a/src/passb4.f b/src/passb4.f index fd9fbab..06c7b13 100644 --- a/src/passb4.f +++ b/src/passb4.f @@ -1,5 +1,6 @@ SUBROUTINE PASSB4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1 WA1(1) ,WA2(1) ,WA3(1) IF (IDO .NE. 2) GO TO 102 diff --git a/src/passb5.f b/src/passb5.f index f8a6f53..48c6e69 100644 --- a/src/passb5.f +++ b/src/passb5.f @@ -1,5 +1,6 @@ SUBROUTINE PASSB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) C *** TR11=COS(2*PI/5), TI11=SIN(2*PI/5) diff --git a/src/passf.f b/src/passf.f index cb97d12..1689851 100644 --- a/src/passf.f +++ b/src/passf.f @@ -1,5 +1,6 @@ SUBROUTINE PASSF (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,WA(1) ,C2(IDL1,IP), 2 CH2(IDL1,IP) diff --git a/src/passf2.f b/src/passf2.f index 765b69f..cddf984 100644 --- a/src/passf2.f +++ b/src/passf2.f @@ -1,5 +1,6 @@ SUBROUTINE PASSF2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1 WA1(1) IF (IDO .GT. 2) GO TO 102 diff --git a/src/passf3.f b/src/passf3.f index d3547ca..a64d207 100644 --- a/src/passf3.f +++ b/src/passf3.f @@ -1,5 +1,6 @@ SUBROUTINE PASSF3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1 WA1(1) ,WA2(1) C *** TAUI IS -SQRT(3)/2 *** diff --git a/src/passf4.f b/src/passf4.f index d14ad61..09daafe 100644 --- a/src/passf4.f +++ b/src/passf4.f @@ -1,5 +1,6 @@ SUBROUTINE PASSF4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1 WA1(1) ,WA2(1) ,WA3(1) IF (IDO .NE. 2) GO TO 102 diff --git a/src/passf5.f b/src/passf5.f index 9e56a3b..6e3da88 100644 --- a/src/passf5.f +++ b/src/passf5.f @@ -1,5 +1,6 @@ SUBROUTINE PASSF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) C *** TR11=COS(2*PI/5), TI11=-SIN(2*PI/5) diff --git a/src/radb2.f b/src/radb2.f index e8f18a1..8a05aed 100644 --- a/src/radb2.f +++ b/src/radb2.f @@ -1,5 +1,6 @@ SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1 WA1(1) DO 101 K=1,L1 diff --git a/src/radb3.f b/src/radb3.f index 6f2b689..b5722e3 100644 --- a/src/radb3.f +++ b/src/radb3.f @@ -1,5 +1,6 @@ SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1 WA1(1) ,WA2(1) C *** TAUI IS SQRT(3)/2 *** diff --git a/src/radb4.f b/src/radb4.f index 6917fb5..72ff6a4 100644 --- a/src/radb4.f +++ b/src/radb4.f @@ -1,5 +1,6 @@ SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1 WA1(1) ,WA2(1) ,WA3(1) DATA SQRT2 /1.41421356237309504880D0/ diff --git a/src/radb5.f b/src/radb5.f index e1a8664..3b90b46 100644 --- a/src/radb5.f +++ b/src/radb5.f @@ -1,5 +1,6 @@ SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) C *** TR11=COS(2*PI/5), TI11=SIN(2*PI/5) diff --git a/src/radbg.f b/src/radbg.f index deaed58..47472d6 100644 --- a/src/radbg.f +++ b/src/radbg.f @@ -1,5 +1,6 @@ SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(1) diff --git a/src/radf2.f b/src/radf2.f index 831682f..184da2d 100644 --- a/src/radf2.f +++ b/src/radf2.f @@ -1,5 +1,6 @@ SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) , 1 WA1(1) DO 101 K=1,L1 diff --git a/src/radf3.f b/src/radf3.f index f6f402f..078b6b2 100644 --- a/src/radf3.f +++ b/src/radf3.f @@ -1,5 +1,6 @@ SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) , 1 WA1(1) ,WA2(1) C *** TAUI IS -SQRT(3)/2 *** diff --git a/src/radf4.f b/src/radf4.f index 1a33060..c7305f1 100644 --- a/src/radf4.f +++ b/src/radf4.f @@ -1,5 +1,6 @@ SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) , 1 WA1(1) ,WA2(1) ,WA3(1) DATA HSQT2 /0.70710678118654752440D0/ diff --git a/src/radf5.f b/src/radf5.f index e45521f..bd79c6d 100644 --- a/src/radf5.f +++ b/src/radf5.f @@ -1,5 +1,6 @@ SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) , 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) DATA TR11,TI11,TR12,TI12 /0.3090169943749474241D0, diff --git a/src/radfg.f b/src/radfg.f index f0ad326..a7d91c6 100644 --- a/src/radfg.f +++ b/src/radfg.f @@ -1,5 +1,6 @@ SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(1) diff --git a/src/rfftb1.f b/src/rfftb1.f index be45445..36f122c 100644 --- a/src/rfftb1.f +++ b/src/rfftb1.f @@ -1,5 +1,6 @@ SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) NF = IFAC(2) NA = 0 diff --git a/src/rfftf1.f b/src/rfftf1.f index 46c3acd..bbe6520 100644 --- a/src/rfftf1.f +++ b/src/rfftf1.f @@ -1,5 +1,6 @@ SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) NF = IFAC(2) NA = 1 diff --git a/src/rffti1.f b/src/rffti1.f index 0fb3b87..5dd1f41 100644 --- a/src/rffti1.f +++ b/src/rffti1.f @@ -1,5 +1,6 @@ SUBROUTINE RFFTI1 (N,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ NL = N diff --git a/src/rk.f b/src/rk.f new file mode 100644 index 0000000..24ab713 --- /dev/null +++ b/src/rk.f @@ -0,0 +1,4 @@ + module fftpack_kind + implicit none + integer,parameter :: rk = kind(1.0d0) + end module diff --git a/src/sint1.f b/src/sint1.f index e6e328c..3619df6 100644 --- a/src/sint1.f +++ b/src/sint1.f @@ -1,5 +1,6 @@ SUBROUTINE SINT1(N,WAR,WAS,XH,X,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WAR(*),WAS(*),X(*),XH(*),IFAC(*) DATA SQRT3 /1.73205080756887729352D0/ DO 100 I=1,N diff --git a/src/zfftb.f b/src/zfftb.f index e4f5204..a2af5ee 100644 --- a/src/zfftb.f +++ b/src/zfftb.f @@ -1,5 +1,6 @@ SUBROUTINE ZFFTB (N,C,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION C(1) ,WSAVE(1) IF (N .EQ. 1) RETURN IW1 = N+N+1 diff --git a/src/zfftf.f b/src/zfftf.f index d16e7e4..78cb869 100644 --- a/src/zfftf.f +++ b/src/zfftf.f @@ -1,5 +1,6 @@ SUBROUTINE ZFFTF (N,C,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION C(1) ,WSAVE(1) IF (N .EQ. 1) RETURN IW1 = N+N+1 diff --git a/src/zffti.f b/src/zffti.f index df69b15..ee64c84 100644 --- a/src/zffti.f +++ b/src/zffti.f @@ -1,5 +1,6 @@ SUBROUTINE ZFFTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + USE fftpack_kind + IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) IF (N .EQ. 1) RETURN IW1 = N+N+1 diff --git a/test/test_fftpack_dcosq.f90 b/test/test_fftpack_dcosq.f90 index b5735f5..7169be6 100644 --- a/test/test_fftpack_dcosq.f90 +++ b/test/test_fftpack_dcosq.f90 @@ -12,17 +12,18 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_dcosq_real - use fftpack, only: dcosqi, dcosqf, dcosqb, dp - real(kind=dp) :: w(3*4 + 15) - real(kind=dp) :: x(4) = [1, 2, 3, 4] - real(kind=dp) :: eps = 1.0e-10_dp + use fftpack, only: dcosqi, dcosqf, dcosqb + use fftpack_kind + real(kind=rk) :: w(3*4 + 15) + real(kind=rk) :: x(4) = [1, 2, 3, 4] + real(kind=rk) :: eps = 1.0e-10_rk call dcosqi(4, w) call dcosqf(4, x, w) - call check(sum(abs(x - [11.999626276085150_dp, -9.1029432177492193_dp, & - 2.6176618435106480_dp, -1.5143449018465791_dp])) < eps, msg="`dcosqf` failed.") + call check(sum(abs(x - [11.999626276085150_rk, -9.1029432177492193_rk, & + 2.6176618435106480_rk, -1.5143449018465791_rk])) < eps, msg="`dcosqf` failed.") call dcosqb(4, x, w) !! - call check(sum(abs(x/(4.0_dp*4.0_dp) - [real(kind=dp) :: 1, 2, 3, 4])) < eps, msg="`dcosqb` failed.") + call check(sum(abs(x/(4.0_rk*4.0_rk) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, msg="`dcosqb` failed.") end subroutine test_fftpack_dcosq_real diff --git a/test/test_fftpack_dfft.f90 b/test/test_fftpack_dfft.f90 index 1b69479..bbad26a 100644 --- a/test/test_fftpack_dfft.f90 +++ b/test/test_fftpack_dfft.f90 @@ -12,20 +12,21 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_dfft() - use fftpack, only: dffti, dfftf, dfftb, dp + use fftpack, only: dffti, dfftf, dfftb + use fftpack_kind - real(kind=dp) :: x(4) - real(kind=dp) :: w(31) + real(kind=rk) :: x(4) + real(kind=rk) :: w(31) x = [1, 2, 3, 4] call dffti(4, w) call dfftf(4, x, w) - call check(all(x == [real(kind=dp) :: 10, -2, 2, -2]), & + call check(all(x == [real(kind=rk) :: 10, -2, 2, -2]), & msg="`dfftf` failed.") call dfftb(4, x, w) - call check(all(x/4.0_dp == [real(kind=dp) :: 1, 2, 3, 4]), & + call check(all(x/4.0_rk == [real(kind=rk) :: 1, 2, 3, 4]), & msg="`dfftb` failed.") end subroutine test_fftpack_dfft diff --git a/test/test_fftpack_dzfft.f90 b/test/test_fftpack_dzfft.f90 index 0464feb..31f3ebc 100644 --- a/test/test_fftpack_dzfft.f90 +++ b/test/test_fftpack_dzfft.f90 @@ -12,21 +12,22 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_dzfft - use fftpack, only: dzffti, dzfftf, dzfftb, dp + use fftpack, only: dzffti, dzfftf, dzfftb + use fftpack_kind - real(kind=dp) :: x(4) = [1, 2, 3, 4] - real(kind=dp) :: w(3*4 + 15) - real(kind=dp) :: azero, a(4/2), b(4/2) + real(kind=rk) :: x(4) = [1, 2, 3, 4] + real(kind=rk) :: w(3*4 + 15) + real(kind=rk) :: azero, a(4/2), b(4/2) call dzffti(4, w) call dzfftf(4, x, azero, a, b, w) - call check(azero == 2.5_dp, msg="azero == 2.5_dp failed.") - call check(all(a == [-1.0_dp, -0.5_dp]), msg="all(a == [-1.0, -0.5]) failed.") - call check(all(b == [-1.0_dp, 0.0_dp]), msg="all(b == [-1.0, 0.0]) failed.") + call check(azero == 2.5_rk, msg="azero == 2.5_rk failed.") + call check(all(a == [-1.0_rk, -0.5_rk]), msg="all(a == [-1.0, -0.5]) failed.") + call check(all(b == [-1.0_rk, 0.0_rk]), msg="all(b == [-1.0, 0.0]) failed.") x = 0 call dzfftb(4, x, azero, a, b, w) - call check(all(x == [real(kind=dp) :: 1, 2, 3, 4]), msg="all(x = [real(kind=dp) :: 1, 2, 3, 4]) failed.") + call check(all(x == [real(kind=rk) :: 1, 2, 3, 4]), msg="all(x = [real(kind=rk) :: 1, 2, 3, 4]) failed.") end subroutine test_fftpack_dzfft diff --git a/test/test_fftpack_fft.f90 b/test/test_fftpack_fft.f90 index 0e2387d..1770065 100644 --- a/test/test_fftpack_fft.f90 +++ b/test/test_fftpack_fft.f90 @@ -12,16 +12,17 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_fft - use fftpack, only: fft, dp - real(kind=dp) :: eps = 1.0e-10_dp + use fftpack, only: fft + use fftpack_kind + real(kind=rk) :: eps = 1.0e-10_rk - complex(kind=dp) :: x(3) = [1.0_dp, 2.0_dp, 3.0_dp] + complex(kind=rk) :: x(3) = [1.0_rk, 2.0_rk, 3.0_rk] - call check(sum(abs(fft(x, 2) - [(3.0_dp, 0.0_dp), (-1.0_dp, 0.0_dp)])) < eps, & + call check(sum(abs(fft(x, 2) - [(3.0_rk, 0.0_rk), (-1.0_rk, 0.0_rk)])) < eps, & msg="`fft(x, 2)` failed.") call check(sum(abs(fft(x, 3) - fft(x))) < eps, & msg="`fft(x, 3)` failed.") - call check(sum(abs(fft(x, 4) - [(6.0_dp, 0.0_dp), (-2.0_dp, -2.0_dp), (2.0_dp, 0.0_dp), (-2.0_dp, 2.0_dp)])) < eps, & + call check(sum(abs(fft(x, 4) - [(6.0_rk, 0.0_rk), (-2.0_rk, -2.0_rk), (2.0_rk, 0.0_rk), (-2.0_rk, 2.0_rk)])) < eps, & msg="`fft(x, 4)` failed.") end subroutine test_fftpack_fft diff --git a/test/test_fftpack_fftshift.f90 b/test/test_fftpack_fftshift.f90 index 32de4bf..2b90b3c 100644 --- a/test/test_fftpack_fftshift.f90 +++ b/test/test_fftpack_fftshift.f90 @@ -13,28 +13,30 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_fftshift_complex - use fftpack, only: fftshift, dp + use fftpack, only: fftshift + use fftpack_kind - complex(kind=dp) :: xeven(4) = [1, 2, 3, 4] - complex(kind=dp) :: xodd(5) = [1, 2, 3, 4, 5] + complex(kind=rk) :: xeven(4) = [1, 2, 3, 4] + complex(kind=rk) :: xodd(5) = [1, 2, 3, 4, 5] - call check(all(fftshift(xeven) == [complex(kind=dp) :: 3, 4, 1, 2]), & - msg="all(fftshift(xeven) == [complex(kind=dp) :: 3, 4, 1, 2]) failed.") - call check(all(fftshift(xodd) == [complex(kind=dp) :: 4, 5, 1, 2, 3]), & - msg="all(fftshift(xodd) == [complex(kind=dp) :: 4, 5, 1, 2, 3]) failed.") + call check(all(fftshift(xeven) == [complex(kind=rk) :: 3, 4, 1, 2]), & + msg="all(fftshift(xeven) == [complex(kind=rk) :: 3, 4, 1, 2]) failed.") + call check(all(fftshift(xodd) == [complex(kind=rk) :: 4, 5, 1, 2, 3]), & + msg="all(fftshift(xodd) == [complex(kind=rk) :: 4, 5, 1, 2, 3]) failed.") end subroutine test_fftpack_fftshift_complex subroutine test_fftpack_fftshift_real - use fftpack, only: fftshift, dp + use fftpack, only: fftshift + use fftpack_kind - real(kind=dp) :: xeven(4) = [1, 2, 3, 4] - real(kind=dp) :: xodd(5) = [1, 2, 3, 4, 5] + real(kind=rk) :: xeven(4) = [1, 2, 3, 4] + real(kind=rk) :: xodd(5) = [1, 2, 3, 4, 5] - call check(all(fftshift(xeven) == [real(kind=dp) :: 3, 4, 1, 2]), & - msg="all(fftshift(xeven) == [real(kind=dp) :: 3, 4, 1, 2]) failed.") - call check(all(fftshift(xodd) == [real(kind=dp) :: 4, 5, 1, 2, 3]), & - msg="all(fftshift(xodd) == [real(kind=dp) :: 4, 5, 1, 2, 3]) failed.") + call check(all(fftshift(xeven) == [real(kind=rk) :: 3, 4, 1, 2]), & + msg="all(fftshift(xeven) == [real(kind=rk) :: 3, 4, 1, 2]) failed.") + call check(all(fftshift(xodd) == [real(kind=rk) :: 4, 5, 1, 2, 3]), & + msg="all(fftshift(xodd) == [real(kind=rk) :: 4, 5, 1, 2, 3]) failed.") end subroutine test_fftpack_fftshift_real diff --git a/test/test_fftpack_ifft.f90 b/test/test_fftpack_ifft.f90 index b052b95..d7c7f73 100644 --- a/test/test_fftpack_ifft.f90 +++ b/test/test_fftpack_ifft.f90 @@ -12,16 +12,17 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_ifft - use fftpack, only: fft, ifft, dp - real(kind=dp) :: eps = 1.0e-10_dp + use fftpack, only: fft, ifft + use fftpack_kind + real(kind=rk) :: eps = 1.0e-10_rk - complex(kind=dp) :: x(4) = [1, 2, 3, 4] + complex(kind=rk) :: x(4) = [1, 2, 3, 4] - call check(sum(abs(ifft(fft(x))/4.0_dp - [complex(kind=dp) :: 1, 2, 3, 4])) < eps, & - msg="`ifft(fft(x))/4.0_dp` failed.") - call check(sum(abs(ifft(fft(x), 2) - [complex(kind=dp) ::(8, 2), (12, -2)])) < eps, & + call check(sum(abs(ifft(fft(x))/4.0_rk - [complex(kind=rk) :: 1, 2, 3, 4])) < eps, & + msg="`ifft(fft(x))/4.0_rk` failed.") + call check(sum(abs(ifft(fft(x), 2) - [complex(kind=rk) ::(8, 2), (12, -2)])) < eps, & msg="`ifft(fft(x), 2)` failed.") - call check(sum(abs(ifft(fft(x, 2), 4) - [complex(kind=dp) ::(2, 0), (3, -1), (4, 0), (3, 1)])) < eps, & + call check(sum(abs(ifft(fft(x, 2), 4) - [complex(kind=rk) ::(2, 0), (3, -1), (4, 0), (3, 1)])) < eps, & msg="`ifft(fft(x, 2), 4)` failed.") end subroutine test_fftpack_ifft diff --git a/test/test_fftpack_ifftshift.f90 b/test/test_fftpack_ifftshift.f90 index 640fd92..88f56ed 100644 --- a/test/test_fftpack_ifftshift.f90 +++ b/test/test_fftpack_ifftshift.f90 @@ -13,30 +13,32 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_ifftshift_complex - use fftpack, only: ifftshift, dp + use fftpack, only: ifftshift + use fftpack_kind integer :: i - complex(kind=dp) :: xeven(4) = [3, 4, 1, 2] - complex(kind=dp) :: xodd(5) = [4, 5, 1, 2, 3] + complex(kind=rk) :: xeven(4) = [3, 4, 1, 2] + complex(kind=rk) :: xodd(5) = [4, 5, 1, 2, 3] - call check(all(ifftshift(xeven) == [complex(kind=dp) ::(i, i=1, 4)]), & - msg="all(ifftshift(xeven) == [complex(kind=dp) ::(i, i=1, 4)]) failed.") - call check(all(ifftshift(xodd) == [complex(kind=dp) ::(i, i=1, 5)]), & - msg="all(ifftshift(xodd) == [complex(kind=dp) ::(i, i=1, 5)]) failed.") + call check(all(ifftshift(xeven) == [complex(kind=rk) ::(i, i=1, 4)]), & + msg="all(ifftshift(xeven) == [complex(kind=rk) ::(i, i=1, 4)]) failed.") + call check(all(ifftshift(xodd) == [complex(kind=rk) ::(i, i=1, 5)]), & + msg="all(ifftshift(xodd) == [complex(kind=rk) ::(i, i=1, 5)]) failed.") end subroutine test_fftpack_ifftshift_complex subroutine test_fftpack_ifftshift_real - use fftpack, only: ifftshift, dp + use fftpack, only: ifftshift + use fftpack_kind integer :: i - real(kind=dp) :: xeven(4) = [3, 4, 1, 2] - real(kind=dp) :: xodd(5) = [4, 5, 1, 2, 3] + real(kind=rk) :: xeven(4) = [3, 4, 1, 2] + real(kind=rk) :: xodd(5) = [4, 5, 1, 2, 3] - call check(all(ifftshift(xeven) == [real(kind=dp) ::(i, i=1, 4)]), & - msg="all(ifftshift(xeven) == [real(kind=dp) ::(i, i=1, 4)]) failed.") - call check(all(ifftshift(xodd) == [real(kind=dp) ::(i, i=1, 5)]), & - msg="all(ifftshift(xodd) == [real(kind=dp) ::(i, i=1, 5)]) failed.") + call check(all(ifftshift(xeven) == [real(kind=rk) ::(i, i=1, 4)]), & + msg="all(ifftshift(xeven) == [real(kind=rk) ::(i, i=1, 4)]) failed.") + call check(all(ifftshift(xodd) == [real(kind=rk) ::(i, i=1, 5)]), & + msg="all(ifftshift(xodd) == [real(kind=rk) ::(i, i=1, 5)]) failed.") end subroutine test_fftpack_ifftshift_real diff --git a/test/test_fftpack_iqct.f90 b/test/test_fftpack_iqct.f90 index 7f5d97c..62f535b 100644 --- a/test/test_fftpack_iqct.f90 +++ b/test/test_fftpack_iqct.f90 @@ -12,18 +12,19 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_iqct - use fftpack, only: qct, iqct, dp - real(kind=dp) :: eps = 1.0e-10_dp - - real(kind=dp) :: x(4) = [1, 2, 3, 4] - - call check(sum(abs(iqct(qct(x))/(4.0_dp*4.0_dp) - [real(kind=dp) :: 1, 2, 3, 4])) < eps, & - msg="`iqct(qct(x)/(4.0_dp*4.0_dp)` failed.") - call check(sum(abs(iqct(qct(x), 2)/(4.0_dp*2.0_dp) - [1.4483415291679655_dp, 7.4608849947753271_dp])) < eps, & - msg="`iqct(qct(x), 2)/(4.0_dp*2.0_dp)` failed.") - call check(sum(abs(iqct(qct(x, 2), 4)/(4.0_dp*4.0_dp) - [0.5_dp, 0.70932417358418376_dp, 1.0_dp, & - 0.78858050747473762_dp])) < eps, & - msg="`iqct(qct(x, 2), 4)/(4.0_dp*4.0_dp)` failed.") + use fftpack, only: qct, iqct + use fftpack_kind + real(kind=rk) :: eps = 1.0e-10_rk + + real(kind=rk) :: x(4) = [1, 2, 3, 4] + + call check(sum(abs(iqct(qct(x))/(4.0_rk*4.0_rk) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & + msg="`iqct(qct(x)/(4.0_rk*4.0_rk)` failed.") + call check(sum(abs(iqct(qct(x), 2)/(4.0_rk*2.0_rk) - [1.4483415291679655_rk, 7.4608849947753271_rk])) < eps, & + msg="`iqct(qct(x), 2)/(4.0_rk*2.0_rk)` failed.") + call check(sum(abs(iqct(qct(x, 2), 4)/(4.0_rk*4.0_rk) - [0.5_rk, 0.70932417358418376_rk, 1.0_rk, & + 0.78858050747473762_rk])) < eps, & + msg="`iqct(qct(x, 2), 4)/(4.0_rk*4.0_rk)` failed.") end subroutine test_fftpack_iqct diff --git a/test/test_fftpack_irfft.f90 b/test/test_fftpack_irfft.f90 index 7ac0112..0c7e970 100644 --- a/test/test_fftpack_irfft.f90 +++ b/test/test_fftpack_irfft.f90 @@ -12,16 +12,17 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_irfft - use fftpack, only: rfft, irfft, dp - real(kind=dp) :: eps = 1.0e-10_dp + use fftpack, only: rfft, irfft + use fftpack_kind + real(kind=rk) :: eps = 1.0e-10_rk - real(kind=dp) :: x(4) = [1, 2, 3, 4] + real(kind=rk) :: x(4) = [1, 2, 3, 4] - call check(sum(abs(irfft(rfft(x))/4.0_dp - [real(kind=dp) :: 1, 2, 3, 4])) < eps, & - msg="`irfft(rfft(x))/4.0_dp` failed.") - call check(sum(abs(irfft(rfft(x), 2) - [real(kind=dp) :: 8, 12])) < eps, & + call check(sum(abs(irfft(rfft(x))/4.0_rk - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & + msg="`irfft(rfft(x))/4.0_rk` failed.") + call check(sum(abs(irfft(rfft(x), 2) - [real(kind=rk) :: 8, 12])) < eps, & msg="`irfft(rfft(x), 2)` failed.") - call check(sum(abs(irfft(rfft(x, 2), 4) - [real(kind=dp) :: 1, 3, 5, 3])) < eps, & + call check(sum(abs(irfft(rfft(x, 2), 4) - [real(kind=rk) :: 1, 3, 5, 3])) < eps, & msg="`irfft(rfft(x, 2), 4)` failed.") end subroutine test_fftpack_irfft diff --git a/test/test_fftpack_qct.f90 b/test/test_fftpack_qct.f90 index af92545..850d989 100644 --- a/test/test_fftpack_qct.f90 +++ b/test/test_fftpack_qct.f90 @@ -12,17 +12,18 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_qct - use fftpack, only: qct, dp - real(kind=dp) :: eps = 1.0e-10_dp + use fftpack, only: qct + use fftpack_kind + real(kind=rk) :: eps = 1.0e-10_rk - real(kind=dp) :: x(3) = [9, -9, 3] + real(kind=rk) :: x(3) = [9, -9, 3] - call check(sum(abs(qct(x, 2) - [-3.7279220613578570_dp, 21.727922061357859_dp])) < eps, & + call check(sum(abs(qct(x, 2) - [-3.7279220613578570_rk, 21.727922061357859_rk])) < eps, & msg="`qct(x, 2)` failed.") call check(sum(abs(qct(x, 3) - qct(x))) < eps, & msg="`qct(x,3)` failed.") - call check(sum(abs(qct(x, 4) - [-3.3871908980838743_dp, -2.1309424696909023_dp, & - 11.645661095452331_dp, 29.872472272322447_dp])) < eps, & + call check(sum(abs(qct(x, 4) - [-3.3871908980838743_rk, -2.1309424696909023_rk, & + 11.645661095452331_rk, 29.872472272322447_rk])) < eps, & msg="`qct(x, 4)` failed.") end subroutine test_fftpack_qct diff --git a/test/test_fftpack_rfft.f90 b/test/test_fftpack_rfft.f90 index c9079b9..cf37192 100644 --- a/test/test_fftpack_rfft.f90 +++ b/test/test_fftpack_rfft.f90 @@ -12,16 +12,17 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_rfft - use fftpack, only: rfft, dp - real(kind=dp) :: eps = 1.0e-10_dp + use fftpack, only: rfft + use fftpack_kind + real(kind=rk) :: eps = 1.0e-10_rk - real(kind=dp) :: x(3) = [9, -9, 3] + real(kind=rk) :: x(3) = [9, -9, 3] - call check(sum(abs(rfft(x, 2) - [real(kind=dp) :: 0, 18])) < eps, & + call check(sum(abs(rfft(x, 2) - [real(kind=rk) :: 0, 18])) < eps, & msg="`rfft(x, 2)` failed.") call check(sum(abs(rfft(x, 3) - rfft(x))) < eps, & msg="`rfft(x, 3)` failed.") - call check(sum(abs(rfft(x, 4) - [real(kind=dp) :: 3, 6, 9, 21])) < eps, & + call check(sum(abs(rfft(x, 4) - [real(kind=rk) :: 3, 6, 9, 21])) < eps, & msg="`rfft(x, 4)` failed.") end subroutine test_fftpack_rfft diff --git a/test/test_fftpack_zfft.f90 b/test/test_fftpack_zfft.f90 index 7e30646..b45ab3b 100644 --- a/test/test_fftpack_zfft.f90 +++ b/test/test_fftpack_zfft.f90 @@ -12,18 +12,22 @@ subroutine check(condition, msg) end subroutine check subroutine test_fftpack_zfft() - use fftpack, only: zffti, zfftf, zfftb, dp + use fftpack_kind - complex(kind=dp) :: x(4) = [1, 2, 3, 4] - real(kind=dp) :: w(31) + use fftpack_kind + use fftpack, only: zffti, zfftf, zfftb + use fftpack_kind + + complex(kind=rk) :: x(4) = [1, 2, 3, 4] + real(kind=rk) :: w(31) call zffti(4, w) call zfftf(4, x, w) - call check(all(x == [complex(kind=dp) ::(10, 0), (-2, 2), (-2, 0), (-2, -2)]), & + call check(all(x == [complex(kind=rk) ::(10, 0), (-2, 2), (-2, 0), (-2, -2)]), & msg="`zfftf` failed.") call zfftb(4, x, w) - call check(all(x/4.0_dp == [complex(kind=dp) ::(1, 0), (2, 0), (3, 0), (4, 0)]), & + call check(all(x/4.0_rk == [complex(kind=rk) ::(1, 0), (2, 0), (3, 0), (4, 0)]), & msg="`zfftb` failed.") end subroutine test_fftpack_zfft From 7cc1dbe51000be19d61f8aa4e3fdacf5b9e4eaaa Mon Sep 17 00:00:00 2001 From: "Harris M. Snyder" Date: Sat, 11 Sep 2021 23:14:32 -0400 Subject: [PATCH 2/3] a few cleanup fixes to previous commit --- src/fftpack.f90 | 14 +++++++------- src/fftpack_fft.f90 | 2 +- src/fftpack_ifft.f90 | 2 +- src/fftpack_irfft.f90 | 2 +- src/fftpack_rfft.f90 | 2 +- test/tstfft.f | 5 +++-- 6 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/fftpack.f90 b/src/fftpack.f90 index ae54319..5a27c95 100644 --- a/src/fftpack.f90 +++ b/src/fftpack.f90 @@ -30,7 +30,7 @@ end subroutine zffti !> Version: experimental !> - !> Forward transform of a double complex periodic sequence. + !> Forward transform of a complex periodic sequence. !> ([Specification](../page/specs/fftpack.html#zfftf)) pure subroutine zfftf(n, c, wsave) import rk @@ -62,7 +62,7 @@ end subroutine dffti !> Version: experimental !> - !> Forward transform of a double real periodic sequence. + !> Forward transform of a real periodic sequence. !> ([Specification](../page/specs/fftpack.html#dfftf)) pure subroutine dfftf(n, r, wsave) import rk @@ -94,7 +94,7 @@ end subroutine dzffti !> Version: experimental !> - !> Simplified forward transform of a double real periodic sequence. + !> Simplified forward transform of a real periodic sequence. !> ([Specification](../page/specs/fftpack.html#dzfftf)) pure subroutine dzfftf(n, r, azero, a, b, wsave) import rk @@ -154,7 +154,7 @@ end subroutine dcosqb !> Version: experimental !> - !> Forward transform of a double complex periodic sequence. + !> Forward transform of a complex periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#fft)) interface fft pure module function fft_rk(x, n) result(result) @@ -166,7 +166,7 @@ end function fft_rk !> Version: experimental !> - !> Backward transform of a double complex periodic sequence. + !> Backward transform of a complex periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#ifft)) interface ifft pure module function ifft_rk(x, n) result(result) @@ -178,7 +178,7 @@ end function ifft_rk !> Version: experimental !> - !> Forward transform of a double real periodic sequence. + !> Forward transform of a real periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#rfft)) interface rfft pure module function rfft_rk(x, n) result(result) @@ -190,7 +190,7 @@ end function rfft_rk !> Version: experimental !> - !> Backward transform of a double real periodic sequence. + !> Backward transform of a real periodic sequence. !> ([Specifiction](../page/specs/fftpack.html#irfft)) interface irfft pure module function irfft_rk(x, n) result(result) diff --git a/src/fftpack_fft.f90 b/src/fftpack_fft.f90 index bf21eb5..8f02ee0 100644 --- a/src/fftpack_fft.f90 +++ b/src/fftpack_fft.f90 @@ -2,7 +2,7 @@ contains - !> Forward transform of a double complex periodic sequence. + !> Forward transform of a complex periodic sequence. pure module function fft_rk(x, n) result(result) complex(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n diff --git a/src/fftpack_ifft.f90 b/src/fftpack_ifft.f90 index 40068d4..680e64b 100644 --- a/src/fftpack_ifft.f90 +++ b/src/fftpack_ifft.f90 @@ -2,7 +2,7 @@ contains - !> Backward transform of a double complex periodic sequence. + !> Backward transform of a complex periodic sequence. pure module function ifft_rk(x, n) result(result) complex(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n diff --git a/src/fftpack_irfft.f90 b/src/fftpack_irfft.f90 index 97565cb..13cdb05 100644 --- a/src/fftpack_irfft.f90 +++ b/src/fftpack_irfft.f90 @@ -2,7 +2,7 @@ contains - !> Backward transform of a double real periodic sequence. + !> Backward transform of a real periodic sequence. pure module function irfft_rk(x, n) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n diff --git a/src/fftpack_rfft.f90 b/src/fftpack_rfft.f90 index a0c6260..2d10767 100644 --- a/src/fftpack_rfft.f90 +++ b/src/fftpack_rfft.f90 @@ -2,7 +2,7 @@ contains - !> Forward transform of a double real periodic sequence. + !> Forward transform of a real periodic sequence. pure module function rfft_rk(x, n) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n diff --git a/test/tstfft.f b/test/tstfft.f index 5111c69..24ced38 100644 --- a/test/tstfft.f +++ b/test/tstfft.f @@ -53,11 +53,12 @@ PROGRAM TSTFFT C - IMPLICIT DOUBLE PRECISION (A-H,O-Z) + 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) - DOUBLE COMPLEX CX ,CY + 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 From 1c23efa5296ded27c7dbac880d023ad245631bf8 Mon Sep 17 00:00:00 2001 From: "Harris M. Snyder" Date: Sun, 12 Sep 2021 22:35:02 -0400 Subject: [PATCH 3/3] Changed FLOAT() calls to REAL(,RK) --- src/cffti1.f | 4 ++-- src/dcosqi.f | 2 +- src/dcosti.f | 2 +- src/dsinti.f | 2 +- src/dzfftf.f | 2 +- src/ezfft1.f | 4 ++-- src/radbg.f | 2 +- src/radfg.f | 2 +- src/rffti1.f | 4 ++-- test/tstfft.f | 64 +++++++++++++++++++++++++-------------------------- 10 files changed, 44 insertions(+), 44 deletions(-) diff --git a/src/cffti1.f b/src/cffti1.f index 923cdfe..f962cd3 100644 --- a/src/cffti1.f +++ b/src/cffti1.f @@ -28,7 +28,7 @@ SUBROUTINE CFFTI1 (N,WA,IFAC) IFAC(1) = N IFAC(2) = NF TPI = 6.28318530717958647692D0 - ARGH = TPI/FLOAT(N) + ARGH = TPI/REAL(N,RK) I = 2 L1 = 1 DO 110 K1=1,NF @@ -44,7 +44,7 @@ SUBROUTINE CFFTI1 (N,WA,IFAC) WA(I) = 0.0D0 LD = LD+L1 FI = 0.0D0 - ARGLD = FLOAT(LD)*ARGH + ARGLD = REAL(LD,RK)*ARGH DO 108 II=4,IDOT,2 I = I+2 FI = FI+1.D0 diff --git a/src/dcosqi.f b/src/dcosqi.f index 9c09127..b4cc80d 100644 --- a/src/dcosqi.f +++ b/src/dcosqi.f @@ -3,7 +3,7 @@ SUBROUTINE DCOSQI (N,WSAVE) IMPLICIT REAL(RK) (A-H,O-Z) DIMENSION WSAVE(1) DATA PIH /1.57079632679489661923D0/ - DT = PIH/FLOAT(N) + DT = PIH/REAL(N,RK) FK = 0.0D0 DO 101 K=1,N FK = FK+1.0D0 diff --git a/src/dcosti.f b/src/dcosti.f index 4993339..1a28918 100644 --- a/src/dcosti.f +++ b/src/dcosti.f @@ -7,7 +7,7 @@ SUBROUTINE DCOSTI (N,WSAVE) NM1 = N-1 NP1 = N+1 NS2 = N/2 - DT = PI/FLOAT(NM1) + DT = PI/REAL(NM1,RK) FK = 0.0D0 DO 101 K=2,NS2 KC = NP1-K diff --git a/src/dsinti.f b/src/dsinti.f index fbb283c..aaae25b 100644 --- a/src/dsinti.f +++ b/src/dsinti.f @@ -6,7 +6,7 @@ SUBROUTINE DSINTI (N,WSAVE) IF (N .LE. 1) RETURN NS2 = N/2 NP1 = N+1 - DT = PI/FLOAT(NP1) + DT = PI/REAL(NP1,RK) DO 101 K=1,NS2 WSAVE(K) = 2.0D0*SIN(K*DT) 101 CONTINUE diff --git a/src/dzfftf.f b/src/dzfftf.f index 13b4f6b..2781e61 100644 --- a/src/dzfftf.f +++ b/src/dzfftf.f @@ -15,7 +15,7 @@ SUBROUTINE DZFFTF (N,R,AZERO,A,B,WSAVE) WSAVE(I) = R(I) 104 CONTINUE CALL DFFTF (N,WSAVE,WSAVE(N+1)) - CF = 2.0D0/FLOAT(N) + CF = 2.0D0/REAL(N,RK) CFM = -CF AZERO = 0.5D0*CF*WSAVE(1) NS2 = (N+1)/2 diff --git a/src/ezfft1.f b/src/ezfft1.f index 2600c89..8c21abe 100644 --- a/src/ezfft1.f +++ b/src/ezfft1.f @@ -28,7 +28,7 @@ SUBROUTINE EZFFT1 (N,WA,IFAC) 107 IF (NL .NE. 1) GO TO 104 IFAC(1) = N IFAC(2) = NF - ARGH = TPI/FLOAT(N) + ARGH = TPI/REAL(N,RK) IS = 0 NFM1 = NF-1 L1 = 1 @@ -38,7 +38,7 @@ SUBROUTINE EZFFT1 (N,WA,IFAC) L2 = L1*IP IDO = N/L2 IPM = IP-1 - ARG1 = FLOAT(L1)*ARGH + ARG1 = REAL(L1,RK)*ARGH CH1 = 1.0D0 SH1 = 0.0D0 DCH1 = COS(ARG1) diff --git a/src/radbg.f b/src/radbg.f index 47472d6..c38e501 100644 --- a/src/radbg.f +++ b/src/radbg.f @@ -5,7 +5,7 @@ SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(1) DATA TPI/6.28318530717958647692D0/ - ARG = TPI/FLOAT(IP) + ARG = TPI/REAL(IP,RK) DCP = COS(ARG) DSP = SIN(ARG) IDP2 = IDO+2 diff --git a/src/radfg.f b/src/radfg.f index a7d91c6..eac2006 100644 --- a/src/radfg.f +++ b/src/radfg.f @@ -5,7 +5,7 @@ SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(1) DATA TPI/6.28318530717958647692D0/ - ARG = TPI/FLOAT(IP) + ARG = TPI/REAL(IP,RK) DCP = COS(ARG) DSP = SIN(ARG) IPPH = (IP+1)/2 diff --git a/src/rffti1.f b/src/rffti1.f index 5dd1f41..7d4a8e0 100644 --- a/src/rffti1.f +++ b/src/rffti1.f @@ -28,7 +28,7 @@ SUBROUTINE RFFTI1 (N,WA,IFAC) IFAC(1) = N IFAC(2) = NF TPI = 6.28318530717958647692D0 - ARGH = TPI/FLOAT(N) + ARGH = TPI/REAL(N,RK) IS = 0 NFM1 = NF-1 L1 = 1 @@ -42,7 +42,7 @@ SUBROUTINE RFFTI1 (N,WA,IFAC) DO 109 J=1,IPM LD = LD+L1 I = IS - ARGLD = FLOAT(LD)*ARGH + ARGLD = REAL(LD,RK)*ARGH FI = 0.0D0 DO 108 II=3,IDO,2 I = I+2 diff --git a/test/tstfft.f b/test/tstfft.f index 24ced38..36a9a3f 100644 --- a/test/tstfft.f +++ b/test/tstfft.f @@ -65,12 +65,12 @@ PROGRAM TSTFFT DO 157 NZ=1,NNS N = ND(NZ) MODN = MOD(N,2) - FN = FLOAT(N) + FN = REAL(N,RK) TFN = FN+FN NP1 = N+1 NM1 = N-1 DO 101 J=1,NP1 - X(J) = SIN(FLOAT(J)*SQRT2) + X(J) = SIN(REAL(J,RK)*SQRT2) Y(J) = X(J) XH(J) = X(J) 101 CONTINUE @@ -85,9 +85,9 @@ PROGRAM TSTFFT DO 103 K=2,NS2 SUM1 = 0.0D0 SUM2 = 0.0D0 - ARG = FLOAT(K-1)*DT + ARG = REAL(K-1,RK)*DT DO 102 I=1,N - ARG1 = FLOAT(I-1)*ARG + ARG1 = REAL(I-1,RK)*ARG SUM1 = SUM1+X(I)*COS(ARG1) SUM2 = SUM2+X(I)*SIN(ARG1) 102 CONTINUE @@ -112,13 +112,13 @@ PROGRAM TSTFFT RFTF = RFTF/FN DO 109 I=1,N SUM = 0.5D0*X(1) - ARG = FLOAT(I-1)*DT + ARG = REAL(I-1,RK)*DT IF (NS2 .LT. 2) GO TO 108 DO 107 K=2,NS2 - ARG1 = FLOAT(K-1)*ARG + 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*FLOAT((-1)**(I-1))*X(N) + 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) @@ -144,9 +144,9 @@ PROGRAM TSTFFT 112 CONTINUE DO 114 I=1,NM1 Y(I) = 0.0D0 - ARG1 = FLOAT(I)*DT + ARG1 = REAL(I,RK)*DT DO 113 K=1,NM1 - Y(I) = Y(I)+X(K)*SIN(FLOAT(K)*ARG1) + Y(I) = Y(I)+X(K)*SIN(REAL(K,RK)*ARG1) 113 CONTINUE Y(I) = Y(I)+Y(I) 114 CONTINUE @@ -173,10 +173,10 @@ PROGRAM TSTFFT X(I) = XH(I) 117 CONTINUE DO 119 I=1,NP1 - Y(I) = 0.5D0*(X(1)+FLOAT((-1)**(I+1))*X(N+1)) - ARG = FLOAT(I-1)*DT + 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(FLOAT(K-1)*ARG) + Y(I) = Y(I)+X(K)*COS(REAL(K-1,RK)*ARG) 118 CONTINUE Y(I) = Y(I)+Y(I) 119 CONTINUE @@ -205,9 +205,9 @@ PROGRAM TSTFFT DT = PI/(FN+FN) DO 124 I=1,N X(I) = 0.0D0 - ARG = DT*FLOAT(I) + ARG = DT*REAL(I,RK) DO 123 K=1,N - X(I) = X(I)+Y(K)*SIN(FLOAT(K+K-1)*ARG) + X(I) = X(I)+Y(K)*SIN(REAL(K+K-1,RK)*ARG) 123 CONTINUE X(I) = 4.0D0*X(I) 124 CONTINUE @@ -220,10 +220,10 @@ PROGRAM TSTFFT 125 CONTINUE SINQBT = CF*SINQBT DO 127 I=1,N - ARG = FLOAT(I+I-1)*DT - Y(I) = 0.5D0*FLOAT((-1)**(I+1))*X(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(FLOAT(K)*ARG) + Y(I) = Y(I)+X(K)*SIN(REAL(K,RK)*ARG) 126 CONTINUE Y(I) = Y(I)+Y(I) 127 CONTINUE @@ -248,9 +248,9 @@ PROGRAM TSTFFT 130 CONTINUE DO 132 I=1,N X(I) = 0.0D0 - ARG = FLOAT(I-1)*DT + ARG = REAL(I-1,RK)*DT DO 131 K=1,N - X(I) = X(I)+Y(K)*COS(FLOAT(K+K-1)*ARG) + X(I) = X(I)+Y(K)*COS(REAL(K+K-1,RK)*ARG) 131 CONTINUE X(I) = 4.0D0*X(I) 132 CONTINUE @@ -264,9 +264,9 @@ PROGRAM TSTFFT COSQBT = CF*COSQBT DO 135 I=1,N Y(I) = 0.5D0*X(1) - ARG = FLOAT(I+I-1)*DT + ARG = REAL(I+I-1,RK)*DT DO 134 K=2,N - Y(I) = Y(I)+X(K)*COS(FLOAT(K-1)*ARG) + Y(I) = Y(I)+X(K)*COS(REAL(K-1,RK)*ARG) 134 CONTINUE Y(I) = Y(I)+Y(I) 135 CONTINUE @@ -292,17 +292,17 @@ PROGRAM TSTFFT X(I) = XH(I) 138 CONTINUE TPI = 8.0D0*ATAN(1.0D0) - DT = TPI/FLOAT(N) + DT = TPI/REAL(N,RK) NS2 = (N+1)/2 - CF = 2.0D0/FLOAT(N) + 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 = FLOAT(K)*DT + ARG = REAL(K,RK)*DT DO 139 I=1,N - ARG1 = FLOAT(I-1)*ARG + ARG1 = REAL(I-1,RK)*ARG SUM1 = SUM1+X(I)*COS(ARG1) SUM2 = SUM2+X(I)*SIN(ARG1) 139 CONTINUE @@ -330,9 +330,9 @@ PROGRAM TSTFFT IF (MODN .EQ. 0) B(NS2) = 0.0D0 DO 146 I=1,N SUM = AZERO - ARG1 = FLOAT(I-1)*DT + ARG1 = REAL(I-1,RK)*DT DO 145 K=1,NS2 - ARG2 = FLOAT(K)*ARG1 + ARG2 = REAL(K,RK)*ARG1 SUM = SUM+A(K)*COS(ARG2)+B(K)*SIN(ARG2) 145 CONTINUE X(I) = SUM @@ -353,14 +353,14 @@ PROGRAM TSTFFT C TEST CFFTI,CFFTF,CFFTB C DO 149 I=1,N - CX(I) = DCMPLX(COS(SQRT2*FLOAT(I)),SIN(SQRT2*FLOAT(I*I))) + 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 = -FLOAT(I-1)*DT + ARG1 = -REAL(I-1,RK)*DT CY(I) = (0.0D0,0.0D0) DO 150 K=1,N - ARG2 = FLOAT(K-1)*ARG1 + ARG2 = REAL(K-1,RK)*ARG1 CY(I) = CY(I)+DCMPLX(COS(ARG2),SIN(ARG2))*CX(K) 150 CONTINUE 151 CONTINUE @@ -373,10 +373,10 @@ PROGRAM TSTFFT 152 CONTINUE DCFFTF = DCFFTF/FN DO 154 I=1,N - ARG1 = FLOAT(I-1)*DT + ARG1 = REAL(I-1,RK)*DT CY(I) = (0.0D0,0.0D0) DO 153 K=1,N - ARG2 = FLOAT(K-1)*ARG1 + ARG2 = REAL(K-1,RK)*ARG1 CY(I) = CY(I)+DCMPLX(COS(ARG2),SIN(ARG2))*CX(K) 153 CONTINUE 154 CONTINUE