diff --git a/README.md b/README.md index ec1486b..91845e2 100644 --- a/README.md +++ b/README.md @@ -12,12 +12,12 @@ cd fftpack ``` ### Build with [fortran-lang/fpm](https://github.com/fortran-lang/fpm) -Fortran Package Manager (fpm) is a great package manager and build system for Fortran. +Fortran Package Manager (fpm) is a package manager and build system for Fortran. You can build using provided `fpm.toml`: ```bash -fpm build --flag "-O2" -fpm test --flag "-O2" --list -fpm test --flag "-O2" +fpm build +fpm test --list +fpm test ``` To use `fftpack` within your `fpm` project, add the following to your `fpm.toml` file: ```toml diff --git a/example/bench1.f90 b/example/bench1.f90 index af9c7f9..c46496d 100644 --- a/example/bench1.f90 +++ b/example/bench1.f90 @@ -1,10 +1,10 @@ program bench1 use fftpack, only: zffti, zfftf, zfftb +use fftpack_kind, only: rk implicit none -integer, parameter :: dp = kind(0.d0) -complex(dp), allocatable :: z(:) -real(dp), allocatable :: w(:), x(:) -real(dp) :: err, time_init, time_forward, time_backward, t1, t2 +complex(rk), allocatable :: z(:) +real(rk), allocatable :: w(:), x(:) +real(rk) :: err, time_init, time_forward, time_backward, t1, t2 integer :: N N = 1024*1014*16 @@ -32,7 +32,7 @@ program bench1 time_backward = t2-t1 print *, "Done" -err = maxval(abs(x-real(z/N,dp))) +err = maxval(abs(x-real(z/N,rk))) print * print *, "Error: ", err print *, "Init time: ", time_init diff --git a/src/Makefile b/src/Makefile index 2a9f61f..7b868a6 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,53 +1,53 @@ SRCF = \ - zfftb.f\ - cfftb1.f\ - zfftf.f\ - cfftf1.f\ - zffti.f\ - cffti1.f\ - dcosqb.f\ - cosqb1.f\ - dcosqf.f\ - cosqf1.f\ - dcosqi.f\ - dcost.f\ - dcosti.f\ - ezfft1.f\ - dzfftb.f\ - dzfftf.f\ - dzffti.f\ - passb.f\ - passb2.f\ - passb3.f\ - passb4.f\ - passb5.f\ - passf.f\ - passf2.f\ - passf3.f\ - passf4.f\ - passf5.f\ - radb2.f\ - radb3.f\ - radb4.f\ - radb5.f\ - radbg.f\ - radf2.f\ - radf3.f\ - radf4.f\ - radf5.f\ - radfg.f\ - dfftb.f\ - rfftb1.f\ - dfftf.f\ - rfftf1.f\ - dffti.f\ - rffti1.f\ - dsinqb.f\ - dsinqf.f\ - dsinqi.f\ - dsint.f\ - sint1.f\ - dsinti.f + zfftb.f90\ + cfftb1.f90\ + zfftf.f90\ + cfftf1.f90\ + zffti.f90\ + cffti1.f90\ + dcosqb.f90\ + cosqb1.f90\ + dcosqf.f90\ + cosqf1.f90\ + dcosqi.f90\ + dcost.f90\ + dcosti.f90\ + ezfft1.f90\ + dzfftb.f90\ + dzfftf.f90\ + dzffti.f90\ + passb.f90\ + passb2.f90\ + passb3.f90\ + passb4.f90\ + passb5.f90\ + passf.f90\ + passf2.f90\ + passf3.f90\ + passf4.f90\ + passf5.f90\ + radb2.f90\ + radb3.f90\ + radb4.f90\ + radb5.f90\ + radbg.f90\ + radf2.f90\ + radf3.f90\ + radf4.f90\ + radf5.f90\ + radfg.f90\ + dfftb.f90\ + rfftb1.f90\ + dfftf.f90\ + rfftf1.f90\ + dffti.f90\ + rffti1.f90\ + dsinqb.f90\ + dsinqf.f90\ + dsinqi.f90\ + dsint.f90\ + sint1.f90\ + dsinti.f90 SRCF90 = \ fftpack.f90\ @@ -59,9 +59,10 @@ SRCF90 = \ fftpack_ifftshift.f90\ fftpack_qct.f90\ fftpack_iqct.f90\ - fftpack_dct.f90 + fftpack_dct.f90\ + rk.f90 -OBJF := $(SRCF:.f=.o) +OBJF := $(SRCF:.f90=.o) OBJF90 := $(SRCF90:.f90=.o) lib$(LIB).a: $(OBJF) $(OBJF90) @@ -76,12 +77,62 @@ clean: %.o: %.f90 $(FC) $(FFLAGS) -c $< -fftpack_fft.o: fftpack.o -fftpack_ifft.o: fftpack.o -fftpack_rfft.o: fftpack.o -fftpack_irfft.o: fftpack.o -fftpack_qct.o: fftpack.o -fftpack_iqct.o: fftpack.o -fftpack_dct.o: fftpack.o -fftpack_fftshift.o: fftpack.o -fftpack_ifftshift.o: fftpack.o \ No newline at end of file +fftpack_fft.o: fftpack.o rk.o +fftpack_ifft.o: fftpack.o rk.o +fftpack_rfft.o: fftpack.o rk.o +fftpack_irfft.o: fftpack.o rk.o +fftpack_qct.o: fftpack.o rk.o +fftpack_iqct.o: fftpack.o rk.o +fftpack_dct.o: fftpack.o rk.o +fftpack_fftshift.o: fftpack.o rk.o +fftpack_ifftshift.o: fftpack.o rk.o + +zfftb.f90: rk.o +cfftb1.f90: rk.o +zfftf.f90: rk.o +cfftf1.f90: rk.o +zffti.f90: rk.o +cffti1.f90: rk.o +dcosqb.f90: rk.o +cosqb1.f90: rk.o +dcosqf.f90: rk.o +cosqf1.f90: rk.o +dcosqi.f90: rk.o +dcost.f90: rk.o +dcosti.f90: rk.o +ezfft1.f90: rk.o +dzfftb.f90: rk.o +dzfftf.f90: rk.o +dzffti.f90: rk.o +passb.f90: rk.o +passb2.f90: rk.o +passb3.f90: rk.o +passb4.f90: rk.o +passb5.f90: rk.o +passf.f90: rk.o +passf2.f90: rk.o +passf3.f90: rk.o +passf4.f90: rk.o +passf5.f90: rk.o +radb2.f90: rk.o +radb3.f90: rk.o +radb4.f90: rk.o +radb5.f90: rk.o +radbg.f90: rk.o +radf2.f90: rk.o +radf3.f90: rk.o +radf4.f90: rk.o +radf5.f90: rk.o +radfg.f90: rk.o +dfftb.f90: rk.o +rfftb1.f90: rk.o +dfftf.f90: rk.o +rfftf1.f90: rk.o +dffti.f90: rk.o +rffti1.f90: rk.o +dsinqb.f90: rk.o +dsinqf.f90: rk.o +dsinqi.f90: rk.o +dsint.f90: rk.o +sint1.f90: rk.o +dsinti.f90: rk.o diff --git a/src/cfftb1.f b/src/cfftb1.f deleted file mode 100644 index 9faebe2..0000000 --- a/src/cfftb1.f +++ /dev/null @@ -1,62 +0,0 @@ - SUBROUTINE CFFTB1 (N,C,CH,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) - NF = IFAC(2) - NA = 0 - L1 = 1 - IW = 1 - DO 116 K1=1,NF - IP = IFAC(K1+2) - L2 = IP*L1 - IDO = N/L2 - IDOT = IDO+IDO - IDL1 = IDOT*L1 - IF (IP .NE. 4) GO TO 103 - IX2 = IW+IDOT - IX3 = IX2+IDOT - IF (NA .NE. 0) GO TO 101 - CALL PASSB4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) - GO TO 102 - 101 CALL PASSB4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) - 102 NA = 1-NA - GO TO 115 - 103 IF (IP .NE. 2) GO TO 106 - IF (NA .NE. 0) GO TO 104 - CALL PASSB2 (IDOT,L1,C,CH,WA(IW)) - GO TO 105 - 104 CALL PASSB2 (IDOT,L1,CH,C,WA(IW)) - 105 NA = 1-NA - GO TO 115 - 106 IF (IP .NE. 3) GO TO 109 - IX2 = IW+IDOT - IF (NA .NE. 0) GO TO 107 - CALL PASSB3 (IDOT,L1,C,CH,WA(IW),WA(IX2)) - GO TO 108 - 107 CALL PASSB3 (IDOT,L1,CH,C,WA(IW),WA(IX2)) - 108 NA = 1-NA - GO TO 115 - 109 IF (IP .NE. 5) GO TO 112 - IX2 = IW+IDOT - IX3 = IX2+IDOT - IX4 = IX3+IDOT - IF (NA .NE. 0) GO TO 110 - CALL PASSB5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) - GO TO 111 - 110 CALL PASSB5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) - 111 NA = 1-NA - GO TO 115 - 112 IF (NA .NE. 0) GO TO 113 - CALL PASSB (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) - GO TO 114 - 113 CALL PASSB (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) - 114 IF (NAC .NE. 0) NA = 1-NA - 115 L1 = L2 - IW = IW+(IP-1)*IDOT - 116 CONTINUE - IF (NA .EQ. 0) RETURN - N2 = N+N - DO 117 I=1,N2 - C(I) = CH(I) - 117 CONTINUE - RETURN - END diff --git a/src/cfftb1.f90 b/src/cfftb1.f90 new file mode 100644 index 0000000..43f5b54 --- /dev/null +++ b/src/cfftb1.f90 @@ -0,0 +1,68 @@ + subroutine cfftb1(n,c,Ch,Wa,Ifac) + use fftpack_kind + implicit none + real(rk) :: c , Ch , Wa + integer :: i , idl1 , ido , idot , Ifac , ip , iw , ix2 , ix3 , ix4, & + k1 , l1 , l2 , n , n2 , na , nac , nf + dimension Ch(*) , c(*) , Wa(*) , Ifac(*) + nf = Ifac(2) + na = 0 + l1 = 1 + iw = 1 + do k1 = 1 , nf + ip = Ifac(k1+2) + l2 = ip*l1 + ido = n/l2 + idot = ido + ido + idl1 = idot*l1 + if ( ip==4 ) then + ix2 = iw + idot + ix3 = ix2 + idot + if ( na/=0 ) then + call passb4(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3)) + else + call passb4(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3)) + endif + na = 1 - na + elseif ( ip==2 ) then + if ( na/=0 ) then + call passb2(idot,l1,Ch,c,Wa(iw)) + else + call passb2(idot,l1,c,Ch,Wa(iw)) + endif + na = 1 - na + elseif ( ip==3 ) then + ix2 = iw + idot + if ( na/=0 ) then + call passb3(idot,l1,Ch,c,Wa(iw),Wa(ix2)) + else + call passb3(idot,l1,c,Ch,Wa(iw),Wa(ix2)) + endif + na = 1 - na + elseif ( ip/=5 ) then + if ( na/=0 ) then + call passb(nac,idot,ip,l1,idl1,Ch,Ch,Ch,c,c,Wa(iw)) + else + call passb(nac,idot,ip,l1,idl1,c,c,c,Ch,Ch,Wa(iw)) + endif + if ( nac/=0 ) na = 1 - na + else + ix2 = iw + idot + ix3 = ix2 + idot + ix4 = ix3 + idot + if ( na/=0 ) then + call passb5(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) + else + call passb5(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) + endif + na = 1 - na + endif + l1 = l2 + iw = iw + (ip-1)*idot + enddo + if ( na==0 ) return + n2 = n + n + do i = 1 , n2 + c(i) = Ch(i) + enddo + end subroutine cfftb1 \ No newline at end of file diff --git a/src/cfftf1.f b/src/cfftf1.f deleted file mode 100644 index 2eacd2a..0000000 --- a/src/cfftf1.f +++ /dev/null @@ -1,62 +0,0 @@ - SUBROUTINE CFFTF1 (N,C,CH,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) - NF = IFAC(2) - NA = 0 - L1 = 1 - IW = 1 - DO 116 K1=1,NF - IP = IFAC(K1+2) - L2 = IP*L1 - IDO = N/L2 - IDOT = IDO+IDO - IDL1 = IDOT*L1 - IF (IP .NE. 4) GO TO 103 - IX2 = IW+IDOT - IX3 = IX2+IDOT - IF (NA .NE. 0) GO TO 101 - CALL PASSF4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) - GO TO 102 - 101 CALL PASSF4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) - 102 NA = 1-NA - GO TO 115 - 103 IF (IP .NE. 2) GO TO 106 - IF (NA .NE. 0) GO TO 104 - CALL PASSF2 (IDOT,L1,C,CH,WA(IW)) - GO TO 105 - 104 CALL PASSF2 (IDOT,L1,CH,C,WA(IW)) - 105 NA = 1-NA - GO TO 115 - 106 IF (IP .NE. 3) GO TO 109 - IX2 = IW+IDOT - IF (NA .NE. 0) GO TO 107 - CALL PASSF3 (IDOT,L1,C,CH,WA(IW),WA(IX2)) - GO TO 108 - 107 CALL PASSF3 (IDOT,L1,CH,C,WA(IW),WA(IX2)) - 108 NA = 1-NA - GO TO 115 - 109 IF (IP .NE. 5) GO TO 112 - IX2 = IW+IDOT - IX3 = IX2+IDOT - IX4 = IX3+IDOT - IF (NA .NE. 0) GO TO 110 - CALL PASSF5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) - GO TO 111 - 110 CALL PASSF5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) - 111 NA = 1-NA - GO TO 115 - 112 IF (NA .NE. 0) GO TO 113 - CALL PASSF (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) - GO TO 114 - 113 CALL PASSF (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) - 114 IF (NAC .NE. 0) NA = 1-NA - 115 L1 = L2 - IW = IW+(IP-1)*IDOT - 116 CONTINUE - IF (NA .EQ. 0) RETURN - N2 = N+N - DO 117 I=1,N2 - C(I) = CH(I) - 117 CONTINUE - RETURN - END diff --git a/src/cfftf1.f90 b/src/cfftf1.f90 new file mode 100644 index 0000000..0139f39 --- /dev/null +++ b/src/cfftf1.f90 @@ -0,0 +1,68 @@ + subroutine cfftf1(n,c,Ch,Wa,Ifac) + use fftpack_kind + implicit none + real(rk) :: c , Ch , Wa + integer :: i , idl1 , ido , idot , Ifac , ip , iw , ix2 , ix3 , ix4, & + k1 , l1 , l2 , n , n2 , na , nac , nf + dimension Ch(*) , c(*) , Wa(*) , Ifac(*) + nf = Ifac(2) + na = 0 + l1 = 1 + iw = 1 + do k1 = 1 , nf + ip = Ifac(k1+2) + l2 = ip*l1 + ido = n/l2 + idot = ido + ido + idl1 = idot*l1 + if ( ip==4 ) then + ix2 = iw + idot + ix3 = ix2 + idot + if ( na/=0 ) then + call passf4(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3)) + else + call passf4(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3)) + endif + na = 1 - na + elseif ( ip==2 ) then + if ( na/=0 ) then + call passf2(idot,l1,Ch,c,Wa(iw)) + else + call passf2(idot,l1,c,Ch,Wa(iw)) + endif + na = 1 - na + elseif ( ip==3 ) then + ix2 = iw + idot + if ( na/=0 ) then + call passf3(idot,l1,Ch,c,Wa(iw),Wa(ix2)) + else + call passf3(idot,l1,c,Ch,Wa(iw),Wa(ix2)) + endif + na = 1 - na + elseif ( ip/=5 ) then + if ( na/=0 ) then + call passf(nac,idot,ip,l1,idl1,Ch,Ch,Ch,c,c,Wa(iw)) + else + call passf(nac,idot,ip,l1,idl1,c,c,c,Ch,Ch,Wa(iw)) + endif + if ( nac/=0 ) na = 1 - na + else + ix2 = iw + idot + ix3 = ix2 + idot + ix4 = ix3 + idot + if ( na/=0 ) then + call passf5(idot,l1,Ch,c,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) + else + call passf5(idot,l1,c,Ch,Wa(iw),Wa(ix2),Wa(ix3),Wa(ix4)) + endif + na = 1 - na + endif + l1 = l2 + iw = iw + (ip-1)*idot + enddo + if ( na==0 ) return + n2 = n + n + do i = 1 , n2 + c(i) = Ch(i) + enddo + end subroutine cfftf1 \ No newline at end of file diff --git a/src/cffti1.f b/src/cffti1.f deleted file mode 100644 index e1f45e9..0000000 --- a/src/cffti1.f +++ /dev/null @@ -1,61 +0,0 @@ - SUBROUTINE CFFTI1 (N,WA,IFAC) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) - DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/ - NL = N - NF = 0 - J = 0 - 101 J = J+1 - IF (J-4) 102,102,103 - 102 NTRY = NTRYH(J) - GO TO 104 - 103 NTRY = NTRY+2 - 104 NQ = NL/NTRY - NR = NL-NTRY*NQ - IF (NR) 101,105,101 - 105 NF = NF+1 - IFAC(NF+2) = NTRY - NL = NQ - IF (NTRY .NE. 2) GO TO 107 - IF (NF .EQ. 1) GO TO 107 - DO 106 I=2,NF - IB = NF-I+2 - IFAC(IB+2) = IFAC(IB+1) - 106 CONTINUE - IFAC(3) = 2 - 107 IF (NL .NE. 1) GO TO 104 - IFAC(1) = N - IFAC(2) = NF - TPI = 6.28318530717958647692D0 - ARGH = TPI/FLOAT(N) - I = 2 - L1 = 1 - DO 110 K1=1,NF - IP = IFAC(K1+2) - LD = 0 - L2 = L1*IP - IDO = N/L2 - IDOT = IDO+IDO+2 - IPM = IP-1 - DO 109 J=1,IPM - I1 = I - WA(I-1) = 1.0D0 - WA(I) = 0.0D0 - LD = LD+L1 - FI = 0.0D0 - ARGLD = FLOAT(LD)*ARGH - DO 108 II=4,IDOT,2 - I = I+2 - FI = FI+1.D0 - ARG = FI*ARGLD - WA(I-1) = COS(ARG) - WA(I) = SIN(ARG) - 108 CONTINUE - IF (IP .LE. 5) GO TO 109 - WA(I1-1) = WA(I-1) - WA(I1) = WA(I) - 109 CONTINUE - L1 = L2 - 110 CONTINUE - RETURN - END diff --git a/src/cffti1.f90 b/src/cffti1.f90 new file mode 100644 index 0000000..1335697 --- /dev/null +++ b/src/cffti1.f90 @@ -0,0 +1,68 @@ + subroutine cffti1(n,Wa,Ifac) + use fftpack_kind + implicit none + real(rk) :: arg , argh , argld , fi , Wa + integer :: i , i1 , ib , ido , idot , Ifac , ii , ip , ipm , j , k1, & + l1 , l2 , ld , n , nf , nl , nq , nr , ntry + dimension Wa(*) , Ifac(*) + integer,dimension(4),parameter :: ntryh = [3 , 4 , 2 , 5] + real(rk),parameter :: tpi = 2.0_rk * acos(-1.0_rk) ! 2 * pi + nl = n + nf = 0 + j = 0 + 100 j = j + 1 + if ( j<=4 ) then + ntry = ntryh(j) + else + ntry = ntry + 2 + endif + 200 nq = nl/ntry + nr = nl - ntry*nq + if ( nr/=0 ) goto 100 + nf = nf + 1 + Ifac(nf+2) = ntry + nl = nq + if ( ntry==2 ) then + if ( nf/=1 ) then + do i = 2 , nf + ib = nf - i + 2 + Ifac(ib+2) = Ifac(ib+1) + enddo + Ifac(3) = 2 + endif + endif + if ( nl/=1 ) goto 200 + Ifac(1) = n + Ifac(2) = nf + argh = tpi/real(n, rk) + i = 2 + l1 = 1 + do k1 = 1 , nf + ip = Ifac(k1+2) + ld = 0 + l2 = l1*ip + ido = n/l2 + idot = ido + ido + 2 + ipm = ip - 1 + do j = 1 , ipm + i1 = i + Wa(i-1) = 1.0_rk + Wa(i) = 0.0_rk + ld = ld + l1 + fi = 0.0_rk + argld = real(ld, rk)*argh + do ii = 4 , idot , 2 + i = i + 2 + fi = fi + 1.0_rk + arg = fi*argld + Wa(i-1) = cos(arg) + Wa(i) = sin(arg) + enddo + if ( ip>5 ) then + Wa(i1-1) = Wa(i-1) + Wa(i1) = Wa(i) + endif + enddo + l1 = l2 + enddo + end subroutine cffti1 \ No newline at end of file diff --git a/src/cosqb1.f b/src/cosqb1.f deleted file mode 100644 index 5c8b6ec..0000000 --- a/src/cosqb1.f +++ /dev/null @@ -1,28 +0,0 @@ - SUBROUTINE COSQB1 (N,X,W,XH) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(1) ,W(1) ,XH(1) - NS2 = (N+1)/2 - NP2 = N+2 - DO 101 I=3,N,2 - XIM1 = X(I-1)+X(I) - X(I) = X(I)-X(I-1) - X(I-1) = XIM1 - 101 CONTINUE - X(1) = X(1)+X(1) - MODN = MOD(N,2) - IF (MODN .EQ. 0) X(N) = X(N)+X(N) - CALL DFFTB (N,X,XH) - DO 102 K=2,NS2 - KC = NP2-K - XH(K) = W(K-1)*X(KC)+W(KC-1)*X(K) - XH(KC) = W(K-1)*X(K)-W(KC-1)*X(KC) - 102 CONTINUE - IF (MODN .EQ. 0) X(NS2+1) = W(NS2)*(X(NS2+1)+X(NS2+1)) - DO 103 K=2,NS2 - KC = NP2-K - X(K) = XH(K)+XH(KC) - X(KC) = XH(K)-XH(KC) - 103 CONTINUE - X(1) = X(1)+X(1) - RETURN - END diff --git a/src/cosqb1.f90 b/src/cosqb1.f90 new file mode 100644 index 0000000..11d7f4a --- /dev/null +++ b/src/cosqb1.f90 @@ -0,0 +1,30 @@ + subroutine cosqb1(n,x,w,Xh) + use fftpack_kind + implicit none + integer :: i , k , kc , modn , n , np2 , ns2 + real(rk) :: w , x , Xh , xim1 + dimension x(*) , w(*) , Xh(*) + ns2 = (n+1)/2 + np2 = n + 2 + do i = 3 , n , 2 + xim1 = x(i-1) + x(i) + x(i) = x(i) - x(i-1) + x(i-1) = xim1 + enddo + x(1) = x(1) + x(1) + modn = mod(n,2) + if ( modn==0 ) x(n) = x(n) + x(n) + call dfftb(n,x,Xh) + do k = 2 , ns2 + kc = np2 - k + Xh(k) = w(k-1)*x(kc) + w(kc-1)*x(k) + Xh(kc) = w(k-1)*x(k) - w(kc-1)*x(kc) + enddo + if ( modn==0 ) x(ns2+1) = w(ns2)*(x(ns2+1)+x(ns2+1)) + do k = 2 , ns2 + kc = np2 - k + x(k) = Xh(k) + Xh(kc) + x(kc) = Xh(k) - Xh(kc) + enddo + x(1) = x(1) + x(1) + end subroutine cosqb1 \ No newline at end of file diff --git a/src/cosqf1.f b/src/cosqf1.f deleted file mode 100644 index 87a0612..0000000 --- a/src/cosqf1.f +++ /dev/null @@ -1,26 +0,0 @@ - SUBROUTINE COSQF1 (N,X,W,XH) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(1) ,W(1) ,XH(1) - NS2 = (N+1)/2 - NP2 = N+2 - DO 101 K=2,NS2 - KC = NP2-K - XH(K) = X(K)+X(KC) - XH(KC) = X(K)-X(KC) - 101 CONTINUE - MODN = MOD(N,2) - IF (MODN .EQ. 0) XH(NS2+1) = X(NS2+1)+X(NS2+1) - DO 102 K=2,NS2 - KC = NP2-K - X(K) = W(K-1)*XH(KC)+W(KC-1)*XH(K) - X(KC) = W(K-1)*XH(K)-W(KC-1)*XH(KC) - 102 CONTINUE - IF (MODN .EQ. 0) X(NS2+1) = W(NS2)*XH(NS2+1) - CALL DFFTF (N,X,XH) - DO 103 I=3,N,2 - XIM1 = X(I-1)-X(I) - X(I) = X(I-1)+X(I) - X(I-1) = XIM1 - 103 CONTINUE - RETURN - END diff --git a/src/cosqf1.f90 b/src/cosqf1.f90 new file mode 100644 index 0000000..8ef39f4 --- /dev/null +++ b/src/cosqf1.f90 @@ -0,0 +1,28 @@ + subroutine cosqf1(n,x,w,Xh) + use fftpack_kind + implicit none + integer :: i , k , kc , modn , n , np2 , ns2 + real(rk) :: w , x , Xh , xim1 + dimension x(*) , w(*) , Xh(*) + ns2 = (n+1)/2 + np2 = n + 2 + do k = 2 , ns2 + kc = np2 - k + Xh(k) = x(k) + x(kc) + Xh(kc) = x(k) - x(kc) + enddo + modn = mod(n,2) + if ( modn==0 ) Xh(ns2+1) = x(ns2+1) + x(ns2+1) + do k = 2 , ns2 + kc = np2 - k + x(k) = w(k-1)*Xh(kc) + w(kc-1)*Xh(k) + x(kc) = w(k-1)*Xh(k) - w(kc-1)*Xh(kc) + enddo + if ( modn==0 ) x(ns2+1) = w(ns2)*Xh(ns2+1) + call dfftf(n,x,Xh) + do i = 3 , n , 2 + xim1 = x(i-1) - x(i) + x(i) = x(i-1) + x(i) + x(i-1) = xim1 + enddo + end subroutine cosqf1 \ No newline at end of file diff --git a/src/dcosqb.f b/src/dcosqb.f deleted file mode 100644 index 341d8c3..0000000 --- a/src/dcosqb.f +++ /dev/null @@ -1,14 +0,0 @@ - SUBROUTINE DCOSQB (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(*) ,WSAVE(*) - DATA TSQRT2 /2.82842712474619009760D0/ - IF (N-2) 101,102,103 - 101 X(1) = 4.0D0*X(1) - RETURN - 102 X1 = 4.0D0*(X(1)+X(2)) - X(2) = TSQRT2*(X(1)-X(2)) - X(1) = X1 - RETURN - 103 CALL COSQB1 (N,X,WSAVE,WSAVE(N+1)) - RETURN - END diff --git a/src/dcosqb.f90 b/src/dcosqb.f90 new file mode 100644 index 0000000..1ab8de7 --- /dev/null +++ b/src/dcosqb.f90 @@ -0,0 +1,19 @@ + subroutine dcosqb(n,x,Wsave) + use fftpack_kind + implicit none + integer :: n + real(rk) :: Wsave , x , x1 + dimension x(*) , Wsave(*) + real(rk),parameter :: tsqrt2 = 2.0_rk * sqrt(2.0_rk) + if ( n<2 ) then + x(1) = 4.0_rk*x(1) + return + elseif ( n==2 ) then + x1 = 4.0_rk*(x(1)+x(2)) + x(2) = tsqrt2*(x(1)-x(2)) + x(1) = x1 + return + else + call cosqb1(n,x,Wsave,Wsave(n+1)) + endif + end subroutine dcosqb \ No newline at end of file diff --git a/src/dcosqf.f b/src/dcosqf.f deleted file mode 100644 index ca585bb..0000000 --- a/src/dcosqf.f +++ /dev/null @@ -1,12 +0,0 @@ - SUBROUTINE DCOSQF (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(*) ,WSAVE(*) - DATA SQRT2 /1.41421356237309504880D0/ - IF (N-2) 102,101,103 - 101 TSQX = SQRT2*X(2) - X(2) = X(1)-TSQX - X(1) = X(1)+TSQX - 102 RETURN - 103 CALL COSQF1 (N,X,WSAVE,WSAVE(N+1)) - RETURN - END diff --git a/src/dcosqf.f90 b/src/dcosqf.f90 new file mode 100644 index 0000000..27f4ceb --- /dev/null +++ b/src/dcosqf.f90 @@ -0,0 +1,17 @@ + subroutine dcosqf(n,x,Wsave) + use fftpack_kind + implicit none + integer :: n + real(rk) :: tsqx , Wsave , x + dimension x(*) , Wsave(*) + real(rk),parameter :: sqrt2 = sqrt(2.0_rk) + if ( n<2 ) then + return + elseif ( n==2 ) then + tsqx = sqrt2*x(2) + x(2) = x(1) - tsqx + x(1) = x(1) + tsqx + else + call cosqf1(n,x,Wsave,Wsave(n+1)) + endif + end subroutine dcosqf \ No newline at end of file diff --git a/src/dcosqi.f b/src/dcosqi.f deleted file mode 100644 index 215e5c1..0000000 --- a/src/dcosqi.f +++ /dev/null @@ -1,13 +0,0 @@ - SUBROUTINE DCOSQI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION WSAVE(1) - DATA PIH /1.57079632679489661923D0/ - DT = PIH/FLOAT(N) - FK = 0.0D0 - DO 101 K=1,N - FK = FK+1.0D0 - WSAVE(K) = COS(FK*DT) - 101 CONTINUE - CALL DFFTI (N,WSAVE(N+1)) - RETURN - END diff --git a/src/dcosqi.f90 b/src/dcosqi.f90 new file mode 100644 index 0000000..1faf8d0 --- /dev/null +++ b/src/dcosqi.f90 @@ -0,0 +1,15 @@ + subroutine dcosqi(n,Wsave) + use fftpack_kind + implicit none + real(rk) :: dt , fk , Wsave + integer :: k , n + dimension Wsave(*) + real(rk),parameter :: pih = acos(-1.0_rk) / 2.0_rk ! pi / 2 + dt = pih/real(n, rk) + fk = 0.0_rk + do k = 1 , n + fk = fk + 1.0_rk + Wsave(k) = cos(fk*dt) + enddo + call dffti(n,Wsave(n+1)) + end subroutine dcosqi \ No newline at end of file diff --git a/src/dcost.f b/src/dcost.f deleted file mode 100644 index fbf90f2..0000000 --- a/src/dcost.f +++ /dev/null @@ -1,43 +0,0 @@ - SUBROUTINE DCOST (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(*) ,WSAVE(*) - NM1 = N-1 - NP1 = N+1 - NS2 = N/2 - IF (N-2) 106,101,102 - 101 X1H = X(1)+X(2) - X(2) = X(1)-X(2) - X(1) = X1H - RETURN - 102 IF (N .GT. 3) GO TO 103 - X1P3 = X(1)+X(3) - TX2 = X(2)+X(2) - X(2) = X(1)-X(3) - X(1) = X1P3+TX2 - X(3) = X1P3-TX2 - RETURN - 103 C1 = X(1)-X(N) - X(1) = X(1)+X(N) - DO 104 K=2,NS2 - KC = NP1-K - T1 = X(K)+X(KC) - T2 = X(K)-X(KC) - C1 = C1+WSAVE(KC)*T2 - T2 = WSAVE(K)*T2 - X(K) = T1-T2 - X(KC) = T1+T2 - 104 CONTINUE - MODN = MOD(N,2) - IF (MODN .NE. 0) X(NS2+1) = X(NS2+1)+X(NS2+1) - CALL DFFTF (NM1,X,WSAVE(N+1)) - XIM2 = X(2) - X(2) = C1 - DO 105 I=4,N,2 - XI = X(I) - X(I) = X(I-2)-X(I-1) - X(I-1) = XIM2 - XIM2 = XI - 105 CONTINUE - IF (MODN .NE. 0) X(N) = XIM2 - 106 RETURN - END diff --git a/src/dcost.f90 b/src/dcost.f90 new file mode 100644 index 0000000..fdc431d --- /dev/null +++ b/src/dcost.f90 @@ -0,0 +1,48 @@ + subroutine dcost(n,x,Wsave) + use fftpack_kind + implicit none + real(rk) :: c1 , t1 , t2 , tx2 , Wsave , x , x1h , x1p3 , & + xi , xim2 + integer :: i , k , kc , modn , n , nm1 , np1 , ns2 + dimension x(*) , Wsave(*) + nm1 = n - 1 + np1 = n + 1 + ns2 = n/2 + if ( n<2 ) return + if ( n==2 ) then + x1h = x(1) + x(2) + x(2) = x(1) - x(2) + x(1) = x1h + return + elseif ( n>3 ) then + c1 = x(1) - x(n) + x(1) = x(1) + x(n) + do k = 2 , ns2 + kc = np1 - k + t1 = x(k) + x(kc) + t2 = x(k) - x(kc) + c1 = c1 + Wsave(kc)*t2 + t2 = Wsave(k)*t2 + x(k) = t1 - t2 + x(kc) = t1 + t2 + enddo + modn = mod(n,2) + if ( modn/=0 ) x(ns2+1) = x(ns2+1) + x(ns2+1) + call dfftf(nm1,x,Wsave(n+1)) + xim2 = x(2) + x(2) = c1 + do i = 4 , n , 2 + xi = x(i) + x(i) = x(i-2) - x(i-1) + x(i-1) = xim2 + xim2 = xi + enddo + if ( modn/=0 ) x(n) = xim2 + return + endif + x1p3 = x(1) + x(3) + tx2 = x(2) + x(2) + x(2) = x(1) - x(3) + x(1) = x1p3 + tx2 + x(3) = x1p3 - tx2 + end subroutine dcost \ No newline at end of file diff --git a/src/dcosti.f b/src/dcosti.f deleted file mode 100644 index 4c6d7bb..0000000 --- a/src/dcosti.f +++ /dev/null @@ -1,19 +0,0 @@ - SUBROUTINE DCOSTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION WSAVE(1) - DATA PI /3.14159265358979323846D0/ - IF (N .LE. 3) RETURN - NM1 = N-1 - NP1 = N+1 - NS2 = N/2 - DT = PI/FLOAT(NM1) - FK = 0.0D0 - DO 101 K=2,NS2 - KC = NP1-K - FK = FK+1.0D0 - WSAVE(K) = 2.0D0*SIN(FK*DT) - WSAVE(KC) = 2.0D0*COS(FK*DT) - 101 CONTINUE - CALL DFFTI (NM1,WSAVE(N+1)) - RETURN - END diff --git a/src/dcosti.f90 b/src/dcosti.f90 new file mode 100644 index 0000000..ed4fc1d --- /dev/null +++ b/src/dcosti.f90 @@ -0,0 +1,21 @@ + subroutine dcosti(n,Wsave) + use fftpack_kind + implicit none + real(rk) :: dt , fk , Wsave + integer :: k , kc , n , nm1 , np1 , ns2 + dimension Wsave(*) + real(rk),parameter :: pi = acos(-1.0_rk) + if ( n<=3 ) return + nm1 = n - 1 + np1 = n + 1 + ns2 = n/2 + dt = pi/real(nm1, rk) + fk = 0.0_rk + do k = 2 , ns2 + kc = np1 - k + fk = fk + 1.0_rk + Wsave(k) = 2.0_rk*sin(fk*dt) + Wsave(kc) = 2.0_rk*cos(fk*dt) + enddo + call dffti(nm1,Wsave(n+1)) + end subroutine dcosti \ No newline at end of file diff --git a/src/dfftb.f b/src/dfftb.f deleted file mode 100644 index e4adfca..0000000 --- a/src/dfftb.f +++ /dev/null @@ -1,7 +0,0 @@ - SUBROUTINE DFFTB (N,R,WSAVE) - IMPLICIT DOUBLE PRECISION (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)) - RETURN - END diff --git a/src/dfftb.f90 b/src/dfftb.f90 new file mode 100644 index 0000000..045edc1 --- /dev/null +++ b/src/dfftb.f90 @@ -0,0 +1,9 @@ + subroutine dfftb(n,r,Wsave) + use fftpack_kind + implicit none + integer :: n + real(rk) :: r , Wsave + dimension r(*) , Wsave(*) + if ( n==1 ) return + call rfftb1(n,r,Wsave,Wsave(n+1),Wsave(2*n+1)) + end subroutine dfftb \ No newline at end of file diff --git a/src/dfftf.f b/src/dfftf.f deleted file mode 100644 index 7511a66..0000000 --- a/src/dfftf.f +++ /dev/null @@ -1,7 +0,0 @@ - SUBROUTINE DFFTF (N,R,WSAVE) - IMPLICIT DOUBLE PRECISION (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)) - RETURN - END diff --git a/src/dfftf.f90 b/src/dfftf.f90 new file mode 100644 index 0000000..d23437e --- /dev/null +++ b/src/dfftf.f90 @@ -0,0 +1,9 @@ + subroutine dfftf(n,r,Wsave) + use fftpack_kind + implicit none + integer :: n + real(rk) :: r , Wsave + dimension r(*) , Wsave(*) + if ( n==1 ) return + call rfftf1(n,r,Wsave,Wsave(n+1),Wsave(2*n+1)) + end subroutine dfftf \ No newline at end of file diff --git a/src/dffti.f b/src/dffti.f deleted file mode 100644 index d70250c..0000000 --- a/src/dffti.f +++ /dev/null @@ -1,7 +0,0 @@ - SUBROUTINE DFFTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION WSAVE(1) - IF (N .EQ. 1) RETURN - CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1)) - RETURN - END diff --git a/src/dffti.f90 b/src/dffti.f90 new file mode 100644 index 0000000..d9f07db --- /dev/null +++ b/src/dffti.f90 @@ -0,0 +1,9 @@ + subroutine dffti(n,Wsave) + use fftpack_kind + implicit none + integer :: n + real(rk) :: Wsave + dimension Wsave(*) + if ( n==1 ) return + call rffti1(n,Wsave(n+1),Wsave(2*n+1)) + end subroutine dffti \ No newline at end of file diff --git a/src/dsinqb.f b/src/dsinqb.f deleted file mode 100644 index d8d7835..0000000 --- a/src/dsinqb.f +++ /dev/null @@ -1,19 +0,0 @@ - SUBROUTINE DSINQB (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(1) ,WSAVE(1) - IF (N .GT. 1) GO TO 101 - X(1) = 4.0D0*X(1) - RETURN - 101 NS2 = N/2 - DO 102 K=2,N,2 - X(K) = -X(K) - 102 CONTINUE - CALL DCOSQB (N,X,WSAVE) - DO 103 K=1,NS2 - KC = N-K - XHOLD = X(K) - X(K) = X(KC+1) - X(KC+1) = XHOLD - 103 CONTINUE - RETURN - END diff --git a/src/dsinqb.f90 b/src/dsinqb.f90 new file mode 100644 index 0000000..1832e6b --- /dev/null +++ b/src/dsinqb.f90 @@ -0,0 +1,23 @@ + subroutine dsinqb(n,x,Wsave) + use fftpack_kind + implicit none + integer :: k , kc , n , ns2 + real(rk) :: Wsave , x , xhold + dimension x(*) , Wsave(*) + if ( n>1 ) then + ns2 = n/2 + do k = 2 , n , 2 + x(k) = -x(k) + enddo + call dcosqb(n,x,Wsave) + do k = 1 , ns2 + kc = n - k + xhold = x(k) + x(k) = x(kc+1) + x(kc+1) = xhold + enddo + return + endif + x(1) = 4.0_rk*x(1) + return + end subroutine dsinqb \ No newline at end of file diff --git a/src/dsinqf.f b/src/dsinqf.f deleted file mode 100644 index 1b257ea..0000000 --- a/src/dsinqf.f +++ /dev/null @@ -1,17 +0,0 @@ - SUBROUTINE DSINQF (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(1) ,WSAVE(1) - IF (N .EQ. 1) RETURN - NS2 = N/2 - DO 101 K=1,NS2 - KC = N-K - XHOLD = X(K) - X(K) = X(KC+1) - X(KC+1) = XHOLD - 101 CONTINUE - CALL DCOSQF (N,X,WSAVE) - DO 102 K=2,N,2 - X(K) = -X(K) - 102 CONTINUE - RETURN - END diff --git a/src/dsinqf.f90 b/src/dsinqf.f90 new file mode 100644 index 0000000..66e7312 --- /dev/null +++ b/src/dsinqf.f90 @@ -0,0 +1,19 @@ + subroutine dsinqf(n,x,Wsave) + use fftpack_kind + implicit none + integer :: k , kc , n , ns2 + real(rk) :: Wsave , x , xhold + dimension x(*) , Wsave(*) + if ( n==1 ) return + ns2 = n/2 + do k = 1 , ns2 + kc = n - k + xhold = x(k) + x(k) = x(kc+1) + x(kc+1) = xhold + enddo + call dcosqf(n,x,Wsave) + do k = 2 , n , 2 + x(k) = -x(k) + enddo + end subroutine dsinqf \ No newline at end of file diff --git a/src/dsinqi.f b/src/dsinqi.f deleted file mode 100644 index c4897c2..0000000 --- a/src/dsinqi.f +++ /dev/null @@ -1,6 +0,0 @@ - SUBROUTINE DSINQI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION WSAVE(1) - CALL DCOSQI (N,WSAVE) - RETURN - END diff --git a/src/dsinqi.f90 b/src/dsinqi.f90 new file mode 100644 index 0000000..d95e2db --- /dev/null +++ b/src/dsinqi.f90 @@ -0,0 +1,8 @@ + subroutine dsinqi(n,Wsave) + use fftpack_kind + implicit none + integer :: n + real(rk) :: Wsave + dimension Wsave(*) + call dcosqi(n,Wsave) + end subroutine dsinqi \ No newline at end of file diff --git a/src/dsint.f b/src/dsint.f deleted file mode 100644 index feae4e0..0000000 --- a/src/dsint.f +++ /dev/null @@ -1,10 +0,0 @@ - SUBROUTINE DSINT (N,X,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION X(1) ,WSAVE(1) - NP1 = N+1 - IW1 = N/2+1 - IW2 = IW1+NP1 - IW3 = IW2+NP1 - CALL SINT1(N,X,WSAVE,WSAVE(IW1),WSAVE(IW2),WSAVE(IW3)) - RETURN - END diff --git a/src/dsint.f90 b/src/dsint.f90 new file mode 100644 index 0000000..8c97f2f --- /dev/null +++ b/src/dsint.f90 @@ -0,0 +1,12 @@ + subroutine dsint(n,x,Wsave) + use fftpack_kind + implicit none + integer :: iw1 , iw2 , iw3 , n , np1 + real(rk) :: Wsave , x + dimension x(*) , Wsave(*) + np1 = n + 1 + iw1 = n/2 + 1 + iw2 = iw1 + np1 + iw3 = iw2 + np1 + call sint1(n,x,Wsave,Wsave(iw1),Wsave(iw2),Wsave(iw3)) + end subroutine dsint \ No newline at end of file diff --git a/src/dsinti.f b/src/dsinti.f deleted file mode 100644 index e655c0a..0000000 --- a/src/dsinti.f +++ /dev/null @@ -1,14 +0,0 @@ - SUBROUTINE DSINTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION WSAVE(1) - DATA PI /3.14159265358979323846D0/ - IF (N .LE. 1) RETURN - NS2 = N/2 - NP1 = N+1 - DT = PI/FLOAT(NP1) - DO 101 K=1,NS2 - WSAVE(K) = 2.0D0*SIN(K*DT) - 101 CONTINUE - CALL DFFTI (NP1,WSAVE(NS2+1)) - RETURN - END diff --git a/src/dsinti.f90 b/src/dsinti.f90 new file mode 100644 index 0000000..2069076 --- /dev/null +++ b/src/dsinti.f90 @@ -0,0 +1,16 @@ + subroutine dsinti(n,Wsave) + use fftpack_kind + implicit none + real(rk) :: dt , Wsave + integer :: k , n , np1 , ns2 + dimension Wsave(*) + real(rk),parameter :: pi = acos(-1.0_rk) + if ( n<=1 ) return + ns2 = n/2 + np1 = n + 1 + dt = pi/real(np1, rk) + do k = 1 , ns2 + Wsave(k) = 2.0_rk*sin(k*dt) + enddo + call dffti(np1,Wsave(ns2+1)) + end subroutine dsinti \ No newline at end of file diff --git a/src/dzfftb.f b/src/dzfftb.f deleted file mode 100644 index 07ced9b..0000000 --- a/src/dzfftb.f +++ /dev/null @@ -1,19 +0,0 @@ - SUBROUTINE DZFFTB (N,R,AZERO,A,B,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION R(*) ,A(*) ,B(*) ,WSAVE(*) - IF (N-2) 101,102,103 - 101 R(1) = AZERO - RETURN - 102 R(1) = AZERO+A(1) - R(2) = AZERO-A(1) - RETURN - 103 NS2 = (N-1)/2 - DO 104 I=1,NS2 - R(2*I) = 0.5D0*A(I) - R(2*I+1) = -0.5D0*B(I) - 104 CONTINUE - R(1) = AZERO - IF (MOD(N,2) .EQ. 0) R(N) = A(NS2+1) - CALL DFFTB (N,R,WSAVE(N+1)) - RETURN - END diff --git a/src/dzfftb.f90 b/src/dzfftb.f90 new file mode 100644 index 0000000..e02dcc1 --- /dev/null +++ b/src/dzfftb.f90 @@ -0,0 +1,24 @@ + subroutine dzfftb(n,r,Azero,a,b,Wsave) + use fftpack_kind + implicit none + real(rk) :: a , Azero , b , r , Wsave + integer :: i , n , ns2 + dimension r(*) , a(*) , b(*) , Wsave(*) + if ( n<2 ) then + r(1) = Azero + return + elseif ( n==2 ) then + r(1) = Azero + a(1) + r(2) = Azero - a(1) + return + else + ns2 = (n-1)/2 + do i = 1 , ns2 + r(2*i) = 0.5_rk*a(i) + r(2*i+1) = -0.5_rk*b(i) + enddo + r(1) = Azero + if ( mod(n,2)==0 ) r(n) = a(ns2+1) + call dfftb(n,r,Wsave(n+1)) + endif + end subroutine dzfftb \ No newline at end of file diff --git a/src/dzfftf.f b/src/dzfftf.f deleted file mode 100644 index f4c86de..0000000 --- a/src/dzfftf.f +++ /dev/null @@ -1,30 +0,0 @@ - SUBROUTINE DZFFTF (N,R,AZERO,A,B,WSAVE) -C -C VERSION 3 JUNE 1979 -C - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION R(*) ,A(*) ,B(*) ,WSAVE(*) - IF (N-2) 101,102,103 - 101 AZERO = R(1) - RETURN - 102 AZERO = 0.5D0*(R(1)+R(2)) - A(1) = 0.5D0*(R(1)-R(2)) - RETURN - 103 DO 104 I=1,N - WSAVE(I) = R(I) - 104 CONTINUE - CALL DFFTF (N,WSAVE,WSAVE(N+1)) - CF = 2.0D0/FLOAT(N) - CFM = -CF - AZERO = 0.5D0*CF*WSAVE(1) - NS2 = (N+1)/2 - NS2M = NS2-1 - DO 105 I=1,NS2M - A(I) = CF*WSAVE(2*I) - B(I) = CFM*WSAVE(2*I+1) - 105 CONTINUE - IF (MOD(N,2) .EQ. 1) RETURN - A(NS2) = 0.5D0*CF*WSAVE(N) - B(NS2) = 0.0D0 - RETURN - END diff --git a/src/dzfftf.f90 b/src/dzfftf.f90 new file mode 100644 index 0000000..d122b9b --- /dev/null +++ b/src/dzfftf.f90 @@ -0,0 +1,35 @@ + subroutine dzfftf(n,r,Azero,a,b,Wsave) +! +! VERSION 3 JUNE 1979 +! + use fftpack_kind + implicit none + real(rk) :: a , Azero , b , cf , cfm , r , Wsave + integer :: i , n , ns2 , ns2m + dimension r(*) , a(*) , b(*) , Wsave(*) + if ( n<2 ) then + Azero = r(1) + return + elseif ( n==2 ) then + Azero = 0.5_rk*(r(1)+r(2)) + a(1) = 0.5_rk*(r(1)-r(2)) + return + else + do i = 1 , n + Wsave(i) = r(i) + enddo + call dfftf(n,Wsave,Wsave(n+1)) + cf = 2.0_rk/real(n, rk) + cfm = -cf + Azero = 0.5_rk*cf*Wsave(1) + ns2 = (n+1)/2 + ns2m = ns2 - 1 + do i = 1 , ns2m + a(i) = cf*Wsave(2*i) + b(i) = cfm*Wsave(2*i+1) + enddo + if ( mod(n,2)==1 ) return + a(ns2) = 0.5_rk*cf*Wsave(n) + b(ns2) = 0.0_rk + endif + end subroutine dzfftf \ No newline at end of file diff --git a/src/dzffti.f b/src/dzffti.f deleted file mode 100644 index 8985149..0000000 --- a/src/dzffti.f +++ /dev/null @@ -1,7 +0,0 @@ - SUBROUTINE DZFFTI (N,WSAVE) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION WSAVE(1) - IF (N .EQ. 1) RETURN - CALL EZFFT1 (N,WSAVE(2*N+1),WSAVE(3*N+1)) - RETURN - END diff --git a/src/dzffti.f90 b/src/dzffti.f90 new file mode 100644 index 0000000..8b51821 --- /dev/null +++ b/src/dzffti.f90 @@ -0,0 +1,9 @@ + subroutine dzffti(n,Wsave) + use fftpack_kind + implicit none + integer :: n + real(rk) :: Wsave + dimension Wsave(*) + if ( n==1 ) return + call ezfft1(n,Wsave(2*n+1),Wsave(3*n+1)) + end subroutine dzffti \ No newline at end of file diff --git a/src/ezfft1.f b/src/ezfft1.f deleted file mode 100644 index 230944a..0000000 --- a/src/ezfft1.f +++ /dev/null @@ -1,63 +0,0 @@ - SUBROUTINE EZFFT1 (N,WA,IFAC) - IMPLICIT DOUBLE PRECISION (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/ - NL = N - NF = 0 - J = 0 - 101 J = J+1 - IF (J-4) 102,102,103 - 102 NTRY = NTRYH(J) - GO TO 104 - 103 NTRY = NTRY+2 - 104 NQ = NL/NTRY - NR = NL-NTRY*NQ - IF (NR) 101,105,101 - 105 NF = NF+1 - IFAC(NF+2) = NTRY - NL = NQ - IF (NTRY .NE. 2) GO TO 107 - IF (NF .EQ. 1) GO TO 107 - DO 106 I=2,NF - IB = NF-I+2 - IFAC(IB+2) = IFAC(IB+1) - 106 CONTINUE - IFAC(3) = 2 - 107 IF (NL .NE. 1) GO TO 104 - IFAC(1) = N - IFAC(2) = NF - ARGH = TPI/FLOAT(N) - IS = 0 - NFM1 = NF-1 - L1 = 1 - IF (NFM1 .EQ. 0) RETURN - DO 111 K1=1,NFM1 - IP = IFAC(K1+2) - L2 = L1*IP - IDO = N/L2 - IPM = IP-1 - ARG1 = FLOAT(L1)*ARGH - CH1 = 1.0D0 - SH1 = 0.0D0 - DCH1 = COS(ARG1) - DSH1 = SIN(ARG1) - DO 110 J=1,IPM - CH1H = DCH1*CH1-DSH1*SH1 - SH1 = DCH1*SH1+DSH1*CH1 - CH1 = CH1H - I = IS+2 - WA(I-1) = CH1 - WA(I) = SH1 - IF (IDO .LT. 5) GO TO 109 - DO 108 II=5,IDO,2 - I = I+2 - WA(I-1) = CH1*WA(I-3)-SH1*WA(I-2) - WA(I) = CH1*WA(I-2)+SH1*WA(I-3) - 108 CONTINUE - 109 IS = IS+IDO - 110 CONTINUE - L1 = L2 - 111 CONTINUE - RETURN - END diff --git a/src/ezfft1.f90 b/src/ezfft1.f90 new file mode 100644 index 0000000..45e1c7d --- /dev/null +++ b/src/ezfft1.f90 @@ -0,0 +1,71 @@ + subroutine ezfft1(n,Wa,Ifac) + use fftpack_kind + implicit none + real(rk) :: arg1 , argh , ch1 , ch1h , dch1 , dsh1 , sh1 , & + Wa + integer :: i , ib , ido , Ifac , ii , ip , ipm , is , j , k1 , l1 , & + l2 , n , nf , nfm1 , nl , nq , nr , ntry + dimension Wa(*) , Ifac(*) + integer,dimension(4),parameter :: ntryh = [4 , 2 , 3 , 5] + real(rk),parameter :: tpi = 2.0_rk * acos(-1.0_rk) ! 2 * pi + nl = n + nf = 0 + j = 0 + 100 j = j + 1 + if ( j<=4 ) then + ntry = ntryh(j) + else + ntry = ntry + 2 + endif + 200 nq = nl/ntry + nr = nl - ntry*nq + if ( nr/=0 ) goto 100 + nf = nf + 1 + Ifac(nf+2) = ntry + nl = nq + if ( ntry==2 ) then + if ( nf/=1 ) then + do i = 2 , nf + ib = nf - i + 2 + Ifac(ib+2) = Ifac(ib+1) + enddo + Ifac(3) = 2 + endif + endif + if ( nl/=1 ) goto 200 + Ifac(1) = n + Ifac(2) = nf + argh = tpi/real(n, rk) + is = 0 + nfm1 = nf - 1 + l1 = 1 + if ( nfm1==0 ) return + do k1 = 1 , nfm1 + ip = Ifac(k1+2) + l2 = l1*ip + ido = n/l2 + ipm = ip - 1 + arg1 = real(l1, rk)*argh + ch1 = 1.0_rk + sh1 = 0.0_rk + dch1 = cos(arg1) + dsh1 = sin(arg1) + do j = 1 , ipm + ch1h = dch1*ch1 - dsh1*sh1 + sh1 = dch1*sh1 + dsh1*ch1 + ch1 = ch1h + i = is + 2 + Wa(i-1) = ch1 + Wa(i) = sh1 + if ( ido>=5 ) then + do ii = 5 , ido , 2 + i = i + 2 + Wa(i-1) = ch1*Wa(i-3) - sh1*Wa(i-2) + Wa(i) = ch1*Wa(i-2) + sh1*Wa(i-3) + enddo + endif + is = is + ido + enddo + l1 = l2 + enddo + end subroutine ezfft1 \ No newline at end of file diff --git a/src/fftpack.f90 b/src/fftpack.f90 index 6d07caf..b6b9659 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 @@ -19,6 +18,8 @@ module fftpack public :: dcosti, dcost public :: dct, idct + + public :: rk interface @@ -27,20 +28,20 @@ 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 !> - !> 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 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 @@ -48,10 +49,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 @@ -59,20 +60,20 @@ 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 !> - !> 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 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 @@ -80,10 +81,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 @@ -91,22 +92,22 @@ 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 !> - !> 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 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 @@ -114,12 +115,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 @@ -127,9 +128,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 @@ -137,10 +138,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 @@ -148,19 +149,19 @@ 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 !> Version: experimental !> !> Initialize `dcost`. ([Specification](../page/specs/fftpack.html#dcosti)) pure subroutine dcosti(n, wsave) - import dp + import rk integer, intent(in) :: n - real(kind=dp), intent(out) :: wsave(*) + real(kind=rk), intent(out) :: wsave(*) end subroutine dcosti !> Version: experimental @@ -168,60 +169,60 @@ end subroutine dcosti !> Discrete fourier cosine transform of an even sequence. !> ([Specification](../page/specs/fftpack.html#dcost)) pure subroutine dcost(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 dcost end interface !> 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_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 !> - !> 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_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 !> - !> 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_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 !> - !> 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_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 @@ -229,11 +230,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 @@ -241,11 +242,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 @@ -253,11 +254,11 @@ end function iqct_dp !> Discrete fourier cosine (forward) transform of an even sequence. !> ([Specification](../page/specs/fftpack.html#dct)) interface dct - pure module function dct_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function dct_rk(x, n) result(result) + real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - real(kind=dp), allocatable :: result(:) - end function dct_dp + real(kind=rk), allocatable :: result(:) + end function dct_rk end interface dct !> Version: experimental @@ -265,7 +266,7 @@ end function dct_dp !> Discrete fourier cosine (backward) transform of an even sequence. !> ([Specification](../page/specs/fftpack.html#idct)) interface idct - module procedure :: dct_dp + module procedure :: dct_rk end interface idct !> Version: experimental @@ -273,14 +274,14 @@ end function dct_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 @@ -288,14 +289,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_dct.f90 b/src/fftpack_dct.f90 index 85f3539..88a46d9 100644 --- a/src/fftpack_dct.f90 +++ b/src/fftpack_dct.f90 @@ -3,20 +3,20 @@ contains !> Discrete fourier cosine transform of an even sequence. - pure module function dct_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + pure module function dct_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 dct_dp(x, n) result(result) !> Discrete fourier cosine transformation call dcost(lenseq, result, wsave) - end function dct_dp + end function dct_rk end submodule fftpack_dct diff --git a/src/fftpack_fft.f90 b/src/fftpack_fft.f90 index b6e362f..8f02ee0 100644 --- a/src/fftpack_fft.f90 +++ b/src/fftpack_fft.f90 @@ -2,21 +2,21 @@ contains - !> Forward transform of a double complex periodic sequence. - pure module function fft_dp(x, n) result(result) - complex(kind=dp), intent(in) :: x(:) + !> 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 - 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..680e64b 100644 --- a/src/fftpack_ifft.f90 +++ b/src/fftpack_ifft.f90 @@ -2,21 +2,21 @@ contains - !> Backward transform of a double complex periodic sequence. - pure module function ifft_dp(x, n) result(result) - complex(kind=dp), intent(in) :: x(:) + !> 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 - 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..13cdb05 100644 --- a/src/fftpack_irfft.f90 +++ b/src/fftpack_irfft.f90 @@ -2,21 +2,21 @@ contains - !> Backward transform of a double real periodic sequence. - pure module function irfft_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + !> 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 - 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..2d10767 100644 --- a/src/fftpack_rfft.f90 +++ b/src/fftpack_rfft.f90 @@ -2,21 +2,21 @@ contains - !> Forward transform of a double real periodic sequence. - pure module function rfft_dp(x, n) result(result) - real(kind=dp), intent(in) :: x(:) + !> 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 - 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 deleted file mode 100644 index 43dba87..0000000 --- a/src/passb.f +++ /dev/null @@ -1,117 +0,0 @@ - SUBROUTINE PASSB (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (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) - IDOT = IDO/2 - NT = IP*IDL1 - IPP2 = IP+2 - IPPH = (IP+1)/2 - IDP = IP*IDO -C - IF (IDO .LT. L1) GO TO 106 - DO 103 J=2,IPPH - JC = IPP2-J - DO 102 K=1,L1 - DO 101 I=1,IDO - CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) - CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) - 101 CONTINUE - 102 CONTINUE - 103 CONTINUE - DO 105 K=1,L1 - DO 104 I=1,IDO - CH(I,K,1) = CC(I,1,K) - 104 CONTINUE - 105 CONTINUE - GO TO 112 - 106 DO 109 J=2,IPPH - JC = IPP2-J - DO 108 I=1,IDO - DO 107 K=1,L1 - CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) - CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) - 107 CONTINUE - 108 CONTINUE - 109 CONTINUE - DO 111 I=1,IDO - DO 110 K=1,L1 - CH(I,K,1) = CC(I,1,K) - 110 CONTINUE - 111 CONTINUE - 112 IDL = 2-IDO - INC = 0 - DO 116 L=2,IPPH - LC = IPP2-L - IDL = IDL+IDO - DO 113 IK=1,IDL1 - C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2) - C2(IK,LC) = WA(IDL)*CH2(IK,IP) - 113 CONTINUE - IDLJ = IDL - INC = INC+IDO - DO 115 J=3,IPPH - JC = IPP2-J - IDLJ = IDLJ+INC - IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP - WAR = WA(IDLJ-1) - WAI = WA(IDLJ) - DO 114 IK=1,IDL1 - C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J) - C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC) - 114 CONTINUE - 115 CONTINUE - 116 CONTINUE - DO 118 J=2,IPPH - DO 117 IK=1,IDL1 - CH2(IK,1) = CH2(IK,1)+CH2(IK,J) - 117 CONTINUE - 118 CONTINUE - DO 120 J=2,IPPH - JC = IPP2-J - DO 119 IK=2,IDL1,2 - CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC) - CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC) - CH2(IK,J) = C2(IK,J)+C2(IK-1,JC) - CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC) - 119 CONTINUE - 120 CONTINUE - NAC = 1 - IF (IDO .EQ. 2) RETURN - NAC = 0 - DO 121 IK=1,IDL1 - C2(IK,1) = CH2(IK,1) - 121 CONTINUE - DO 123 J=2,IP - DO 122 K=1,L1 - C1(1,K,J) = CH(1,K,J) - C1(2,K,J) = CH(2,K,J) - 122 CONTINUE - 123 CONTINUE - IF (IDOT .GT. L1) GO TO 127 - IDIJ = 0 - DO 126 J=2,IP - IDIJ = IDIJ+2 - DO 125 I=4,IDO,2 - IDIJ = IDIJ+2 - DO 124 K=1,L1 - C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) - C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) - 124 CONTINUE - 125 CONTINUE - 126 CONTINUE - RETURN - 127 IDJ = 2-IDO - DO 130 J=2,IP - IDJ = IDJ+IDO - DO 129 K=1,L1 - IDIJ = IDJ - DO 128 I=4,IDO,2 - IDIJ = IDIJ+2 - C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) - C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) - 128 CONTINUE - 129 CONTINUE - 130 CONTINUE - RETURN - END diff --git a/src/passb.f90 b/src/passb.f90 new file mode 100644 index 0000000..12c3928 --- /dev/null +++ b/src/passb.f90 @@ -0,0 +1,125 @@ + subroutine passb(Nac,Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) + use fftpack_kind + implicit none + real(rk) :: c1 , c2 , Cc , Ch , Ch2 , Wa , wai , war + integer :: i , idij , idj , idl , Idl1 , idlj , Ido , idot , idp , & + ik , inc , Ip , ipp2 , ipph , j , jc , k , l , l1 , lc + integer :: Nac , nt + dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , Wa(*) , & + c2(Idl1,Ip) , Ch2(Idl1,Ip) + idot = Ido/2 + nt = Ip*Idl1 + ipp2 = Ip + 2 + ipph = (Ip+1)/2 + idp = Ip*Ido +! + if ( Idoidp ) idlj = idlj - idp + war = Wa(idlj-1) + wai = Wa(idlj) + do ik = 1 , Idl1 + c2(ik,l) = c2(ik,l) + war*Ch2(ik,j) + c2(ik,lc) = c2(ik,lc) + wai*Ch2(ik,jc) + enddo + enddo + enddo + do j = 2 , ipph + do ik = 1 , Idl1 + Ch2(ik,1) = Ch2(ik,1) + Ch2(ik,j) + enddo + enddo + do j = 2 , ipph + jc = ipp2 - j + do ik = 2 , Idl1 , 2 + Ch2(ik-1,j) = c2(ik-1,j) - c2(ik,jc) + Ch2(ik-1,jc) = c2(ik-1,j) + c2(ik,jc) + Ch2(ik,j) = c2(ik,j) + c2(ik-1,jc) + Ch2(ik,jc) = c2(ik,j) - c2(ik-1,jc) + enddo + enddo + Nac = 1 + if ( Ido==2 ) return + Nac = 0 + do ik = 1 , Idl1 + c2(ik,1) = Ch2(ik,1) + enddo + do j = 2 , Ip + do k = 1 , l1 + c1(1,k,j) = Ch(1,k,j) + c1(2,k,j) = Ch(2,k,j) + enddo + enddo + if ( idot>l1 ) then + idj = 2 - Ido + do j = 2 , Ip + idj = idj + Ido + do k = 1 , l1 + idij = idj + do i = 4 , Ido , 2 + idij = idij + 2 + c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij) & + *Ch(i,k,j) + c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij) & + *Ch(i-1,k,j) + enddo + enddo + enddo + return + endif + idij = 0 + do j = 2 , Ip + idij = idij + 2 + do i = 4 , Ido , 2 + idij = idij + 2 + do k = 1 , l1 + c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij)*Ch(i,k,j) + c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij)*Ch(i-1,k,j) + enddo + enddo + enddo + return + end subroutine passb \ No newline at end of file diff --git a/src/passb2.f b/src/passb2.f deleted file mode 100644 index 22f6279..0000000 --- a/src/passb2.f +++ /dev/null @@ -1,24 +0,0 @@ - SUBROUTINE PASSB2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , - 1 WA1(1) - IF (IDO .GT. 2) GO TO 102 - DO 101 K=1,L1 - CH(1,K,1) = CC(1,1,K)+CC(1,2,K) - CH(1,K,2) = CC(1,1,K)-CC(1,2,K) - CH(2,K,1) = CC(2,1,K)+CC(2,2,K) - CH(2,K,2) = CC(2,1,K)-CC(2,2,K) - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) - TR2 = CC(I-1,1,K)-CC(I-1,2,K) - CH(I,K,1) = CC(I,1,K)+CC(I,2,K) - TI2 = CC(I,1,K)-CC(I,2,K) - CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2 - CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passb2.f90 b/src/passb2.f90 new file mode 100644 index 0000000..f75bf41 --- /dev/null +++ b/src/passb2.f90 @@ -0,0 +1,26 @@ + subroutine passb2(Ido,l1,Cc,Ch,Wa1) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 + integer :: i , Ido , k , l1 + dimension Cc(Ido,2,l1) , Ch(Ido,l1,2) , Wa1(*) + if ( Ido>2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + Ch(i-1,k,1) = Cc(i-1,1,k) + Cc(i-1,2,k) + tr2 = Cc(i-1,1,k) - Cc(i-1,2,k) + Ch(i,k,1) = Cc(i,1,k) + Cc(i,2,k) + ti2 = Cc(i,1,k) - Cc(i,2,k) + Ch(i,k,2) = Wa1(i-1)*ti2 + Wa1(i)*tr2 + Ch(i-1,k,2) = Wa1(i-1)*tr2 - Wa1(i)*ti2 + enddo + enddo + else + do k = 1 , l1 + Ch(1,k,1) = Cc(1,1,k) + Cc(1,2,k) + Ch(1,k,2) = Cc(1,1,k) - Cc(1,2,k) + Ch(2,k,1) = Cc(2,1,k) + Cc(2,2,k) + Ch(2,k,2) = Cc(2,1,k) - Cc(2,2,k) + enddo + end if + end subroutine passb2 \ No newline at end of file diff --git a/src/passb3.f b/src/passb3.f deleted file mode 100644 index ee8ef10..0000000 --- a/src/passb3.f +++ /dev/null @@ -1,44 +0,0 @@ - SUBROUTINE PASSB3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , - 1 WA1(1) ,WA2(1) -C *** TAUI IS SQRT(3)/2 *** - DATA TAUR,TAUI /-0.5D0,0.86602540378443864676D0/ - IF (IDO .NE. 2) GO TO 102 - DO 101 K=1,L1 - TR2 = CC(1,2,K)+CC(1,3,K) - CR2 = CC(1,1,K)+TAUR*TR2 - CH(1,K,1) = CC(1,1,K)+TR2 - TI2 = CC(2,2,K)+CC(2,3,K) - CI2 = CC(2,1,K)+TAUR*TI2 - CH(2,K,1) = CC(2,1,K)+TI2 - CR3 = TAUI*(CC(1,2,K)-CC(1,3,K)) - CI3 = TAUI*(CC(2,2,K)-CC(2,3,K)) - CH(1,K,2) = CR2-CI3 - CH(1,K,3) = CR2+CI3 - CH(2,K,2) = CI2+CR3 - CH(2,K,3) = CI2-CR3 - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - TR2 = CC(I-1,2,K)+CC(I-1,3,K) - CR2 = CC(I-1,1,K)+TAUR*TR2 - CH(I-1,K,1) = CC(I-1,1,K)+TR2 - TI2 = CC(I,2,K)+CC(I,3,K) - CI2 = CC(I,1,K)+TAUR*TI2 - CH(I,K,1) = CC(I,1,K)+TI2 - CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) - CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) - DR2 = CR2-CI3 - DR3 = CR2+CI3 - DI2 = CI2+CR3 - DI3 = CI2-CR3 - CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 - CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 - CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 - CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passb3.f90 b/src/passb3.f90 new file mode 100644 index 0000000..0161f73 --- /dev/null +++ b/src/passb3.f90 @@ -0,0 +1,47 @@ + subroutine passb3(Ido,l1,Cc,Ch,Wa1,Wa2) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , cr2 , cr3 , di2 , di3 , & + dr2 , dr3 , ti2 , tr2 , Wa1 , Wa2 + integer :: i , Ido , k , l1 + dimension Cc(Ido,3,l1) , Ch(Ido,l1,3) , Wa1(*) , Wa2(*) + real(rk),parameter :: taur = -0.5_rk + real(rk),parameter :: taui = sqrt(3.0_rk) / 2.0_rk + if ( Ido/=2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + tr2 = Cc(i-1,2,k) + Cc(i-1,3,k) + cr2 = Cc(i-1,1,k) + taur*tr2 + Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + ti2 = Cc(i,2,k) + Cc(i,3,k) + ci2 = Cc(i,1,k) + taur*ti2 + Ch(i,k,1) = Cc(i,1,k) + ti2 + cr3 = taui*(Cc(i-1,2,k)-Cc(i-1,3,k)) + ci3 = taui*(Cc(i,2,k)-Cc(i,3,k)) + dr2 = cr2 - ci3 + dr3 = cr2 + ci3 + di2 = ci2 + cr3 + di3 = ci2 - cr3 + Ch(i,k,2) = Wa1(i-1)*di2 + Wa1(i)*dr2 + Ch(i-1,k,2) = Wa1(i-1)*dr2 - Wa1(i)*di2 + Ch(i,k,3) = Wa2(i-1)*di3 + Wa2(i)*dr3 + Ch(i-1,k,3) = Wa2(i-1)*dr3 - Wa2(i)*di3 + enddo + enddo + else + do k = 1 , l1 + tr2 = Cc(1,2,k) + Cc(1,3,k) + cr2 = Cc(1,1,k) + taur*tr2 + Ch(1,k,1) = Cc(1,1,k) + tr2 + ti2 = Cc(2,2,k) + Cc(2,3,k) + ci2 = Cc(2,1,k) + taur*ti2 + Ch(2,k,1) = Cc(2,1,k) + ti2 + cr3 = taui*(Cc(1,2,k)-Cc(1,3,k)) + ci3 = taui*(Cc(2,2,k)-Cc(2,3,k)) + Ch(1,k,2) = cr2 - ci3 + Ch(1,k,3) = cr2 + ci3 + Ch(2,k,2) = ci2 + cr3 + Ch(2,k,3) = ci2 - cr3 + enddo + end if + end subroutine passb3 \ No newline at end of file diff --git a/src/passb4.f b/src/passb4.f deleted file mode 100644 index fd9fbab..0000000 --- a/src/passb4.f +++ /dev/null @@ -1,52 +0,0 @@ - SUBROUTINE PASSB4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (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 - DO 101 K=1,L1 - TI1 = CC(2,1,K)-CC(2,3,K) - TI2 = CC(2,1,K)+CC(2,3,K) - TR4 = CC(2,4,K)-CC(2,2,K) - TI3 = CC(2,2,K)+CC(2,4,K) - TR1 = CC(1,1,K)-CC(1,3,K) - TR2 = CC(1,1,K)+CC(1,3,K) - TI4 = CC(1,2,K)-CC(1,4,K) - TR3 = CC(1,2,K)+CC(1,4,K) - CH(1,K,1) = TR2+TR3 - CH(1,K,3) = TR2-TR3 - CH(2,K,1) = TI2+TI3 - CH(2,K,3) = TI2-TI3 - CH(1,K,2) = TR1+TR4 - CH(1,K,4) = TR1-TR4 - CH(2,K,2) = TI1+TI4 - CH(2,K,4) = TI1-TI4 - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - TI1 = CC(I,1,K)-CC(I,3,K) - TI2 = CC(I,1,K)+CC(I,3,K) - TI3 = CC(I,2,K)+CC(I,4,K) - TR4 = CC(I,4,K)-CC(I,2,K) - TR1 = CC(I-1,1,K)-CC(I-1,3,K) - TR2 = CC(I-1,1,K)+CC(I-1,3,K) - TI4 = CC(I-1,2,K)-CC(I-1,4,K) - TR3 = CC(I-1,2,K)+CC(I-1,4,K) - CH(I-1,K,1) = TR2+TR3 - CR3 = TR2-TR3 - CH(I,K,1) = TI2+TI3 - CI3 = TI2-TI3 - CR2 = TR1+TR4 - CR4 = TR1-TR4 - CI2 = TI1+TI4 - CI4 = TI1-TI4 - CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2 - CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2 - CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3 - CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3 - CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4 - CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passb4.f90 b/src/passb4.f90 new file mode 100644 index 0000000..0c78a1d --- /dev/null +++ b/src/passb4.f90 @@ -0,0 +1,56 @@ + subroutine passb4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & + & ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3 , tr4 , & + & Wa1 , Wa2 , Wa3 + integer :: i , Ido , k , l1 + dimension Cc(Ido,4,l1) , Ch(Ido,l1,4) , Wa1(*) , Wa2(*) , Wa3(*) + if ( Ido/=2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + ti1 = Cc(i,1,k) - Cc(i,3,k) + ti2 = Cc(i,1,k) + Cc(i,3,k) + ti3 = Cc(i,2,k) + Cc(i,4,k) + tr4 = Cc(i,4,k) - Cc(i,2,k) + tr1 = Cc(i-1,1,k) - Cc(i-1,3,k) + tr2 = Cc(i-1,1,k) + Cc(i-1,3,k) + ti4 = Cc(i-1,2,k) - Cc(i-1,4,k) + tr3 = Cc(i-1,2,k) + Cc(i-1,4,k) + Ch(i-1,k,1) = tr2 + tr3 + cr3 = tr2 - tr3 + Ch(i,k,1) = ti2 + ti3 + ci3 = ti2 - ti3 + cr2 = tr1 + tr4 + cr4 = tr1 - tr4 + ci2 = ti1 + ti4 + ci4 = ti1 - ti4 + Ch(i-1,k,2) = Wa1(i-1)*cr2 - Wa1(i)*ci2 + Ch(i,k,2) = Wa1(i-1)*ci2 + Wa1(i)*cr2 + Ch(i-1,k,3) = Wa2(i-1)*cr3 - Wa2(i)*ci3 + Ch(i,k,3) = Wa2(i-1)*ci3 + Wa2(i)*cr3 + Ch(i-1,k,4) = Wa3(i-1)*cr4 - Wa3(i)*ci4 + Ch(i,k,4) = Wa3(i-1)*ci4 + Wa3(i)*cr4 + enddo + enddo + else + do k = 1 , l1 + ti1 = Cc(2,1,k) - Cc(2,3,k) + ti2 = Cc(2,1,k) + Cc(2,3,k) + tr4 = Cc(2,4,k) - Cc(2,2,k) + ti3 = Cc(2,2,k) + Cc(2,4,k) + tr1 = Cc(1,1,k) - Cc(1,3,k) + tr2 = Cc(1,1,k) + Cc(1,3,k) + ti4 = Cc(1,2,k) - Cc(1,4,k) + tr3 = Cc(1,2,k) + Cc(1,4,k) + Ch(1,k,1) = tr2 + tr3 + Ch(1,k,3) = tr2 - tr3 + Ch(2,k,1) = ti2 + ti3 + Ch(2,k,3) = ti2 - ti3 + Ch(1,k,2) = tr1 + tr4 + Ch(1,k,4) = tr1 - tr4 + Ch(2,k,2) = ti1 + ti4 + Ch(2,k,4) = ti1 - ti4 + enddo + end if + end subroutine passb4 \ No newline at end of file diff --git a/src/passb5.f b/src/passb5.f deleted file mode 100644 index f8a6f53..0000000 --- a/src/passb5.f +++ /dev/null @@ -1,79 +0,0 @@ - SUBROUTINE PASSB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (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) -C *** TR12=COS(4*PI/5), TI12=SIN(4*PI/5) - DATA TR11,TI11,TR12,TI12 /0.3090169943749474241D0, - + 0.95105651629515357212D0, - + -0.8090169943749474241D0,0.58778525229247312917D0/ - IF (IDO .NE. 2) GO TO 102 - DO 101 K=1,L1 - TI5 = CC(2,2,K)-CC(2,5,K) - TI2 = CC(2,2,K)+CC(2,5,K) - TI4 = CC(2,3,K)-CC(2,4,K) - TI3 = CC(2,3,K)+CC(2,4,K) - TR5 = CC(1,2,K)-CC(1,5,K) - TR2 = CC(1,2,K)+CC(1,5,K) - TR4 = CC(1,3,K)-CC(1,4,K) - TR3 = CC(1,3,K)+CC(1,4,K) - CH(1,K,1) = CC(1,1,K)+TR2+TR3 - CH(2,K,1) = CC(2,1,K)+TI2+TI3 - CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 - CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3 - CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 - CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3 - CR5 = TI11*TR5+TI12*TR4 - CI5 = TI11*TI5+TI12*TI4 - CR4 = TI12*TR5-TI11*TR4 - CI4 = TI12*TI5-TI11*TI4 - CH(1,K,2) = CR2-CI5 - CH(1,K,5) = CR2+CI5 - CH(2,K,2) = CI2+CR5 - CH(2,K,3) = CI3+CR4 - CH(1,K,3) = CR3-CI4 - CH(1,K,4) = CR3+CI4 - CH(2,K,4) = CI3-CR4 - CH(2,K,5) = CI2-CR5 - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - TI5 = CC(I,2,K)-CC(I,5,K) - TI2 = CC(I,2,K)+CC(I,5,K) - TI4 = CC(I,3,K)-CC(I,4,K) - TI3 = CC(I,3,K)+CC(I,4,K) - TR5 = CC(I-1,2,K)-CC(I-1,5,K) - TR2 = CC(I-1,2,K)+CC(I-1,5,K) - TR4 = CC(I-1,3,K)-CC(I-1,4,K) - TR3 = CC(I-1,3,K)+CC(I-1,4,K) - CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 - CH(I,K,1) = CC(I,1,K)+TI2+TI3 - CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 - CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 - CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 - CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 - CR5 = TI11*TR5+TI12*TR4 - CI5 = TI11*TI5+TI12*TI4 - CR4 = TI12*TR5-TI11*TR4 - CI4 = TI12*TI5-TI11*TI4 - DR3 = CR3-CI4 - DR4 = CR3+CI4 - DI3 = CI3+CR4 - DI4 = CI3-CR4 - DR5 = CR2+CI5 - DR2 = CR2-CI5 - DI5 = CI2-CR5 - DI2 = CI2+CR5 - CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 - CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 - CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 - CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 - CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4 - CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4 - CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5 - CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passb5.f90 b/src/passb5.f90 new file mode 100644 index 0000000..acdfbfe --- /dev/null +++ b/src/passb5.f90 @@ -0,0 +1,86 @@ + subroutine passb5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & + cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & + dr4 , dr5 + real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & + tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 + integer :: i , Ido , k , l1 + dimension Cc(Ido,5,l1) , Ch(Ido,l1,5) , Wa1(*) , Wa2(*) , Wa3(*), & + Wa4(*) + real(rk),parameter :: pi = acos(-1.0_rk) + real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti11 = sin(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti12 = sin(4.0_rk * pi / 5.0_rk) + if ( Ido/=2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + ti5 = Cc(i,2,k) - Cc(i,5,k) + ti2 = Cc(i,2,k) + Cc(i,5,k) + ti4 = Cc(i,3,k) - Cc(i,4,k) + ti3 = Cc(i,3,k) + Cc(i,4,k) + tr5 = Cc(i-1,2,k) - Cc(i-1,5,k) + tr2 = Cc(i-1,2,k) + Cc(i-1,5,k) + tr4 = Cc(i-1,3,k) - Cc(i-1,4,k) + tr3 = Cc(i-1,3,k) + Cc(i-1,4,k) + Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + tr3 + Ch(i,k,1) = Cc(i,1,k) + ti2 + ti3 + cr2 = Cc(i-1,1,k) + tr11*tr2 + tr12*tr3 + ci2 = Cc(i,1,k) + tr11*ti2 + tr12*ti3 + cr3 = Cc(i-1,1,k) + tr12*tr2 + tr11*tr3 + ci3 = Cc(i,1,k) + tr12*ti2 + tr11*ti3 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + dr3 = cr3 - ci4 + dr4 = cr3 + ci4 + di3 = ci3 + cr4 + di4 = ci3 - cr4 + dr5 = cr2 + ci5 + dr2 = cr2 - ci5 + di5 = ci2 - cr5 + di2 = ci2 + cr5 + Ch(i-1,k,2) = Wa1(i-1)*dr2 - Wa1(i)*di2 + Ch(i,k,2) = Wa1(i-1)*di2 + Wa1(i)*dr2 + Ch(i-1,k,3) = Wa2(i-1)*dr3 - Wa2(i)*di3 + Ch(i,k,3) = Wa2(i-1)*di3 + Wa2(i)*dr3 + Ch(i-1,k,4) = Wa3(i-1)*dr4 - Wa3(i)*di4 + Ch(i,k,4) = Wa3(i-1)*di4 + Wa3(i)*dr4 + Ch(i-1,k,5) = Wa4(i-1)*dr5 - Wa4(i)*di5 + Ch(i,k,5) = Wa4(i-1)*di5 + Wa4(i)*dr5 + enddo + enddo + else + do k = 1 , l1 + ti5 = Cc(2,2,k) - Cc(2,5,k) + ti2 = Cc(2,2,k) + Cc(2,5,k) + ti4 = Cc(2,3,k) - Cc(2,4,k) + ti3 = Cc(2,3,k) + Cc(2,4,k) + tr5 = Cc(1,2,k) - Cc(1,5,k) + tr2 = Cc(1,2,k) + Cc(1,5,k) + tr4 = Cc(1,3,k) - Cc(1,4,k) + tr3 = Cc(1,3,k) + Cc(1,4,k) + Ch(1,k,1) = Cc(1,1,k) + tr2 + tr3 + Ch(2,k,1) = Cc(2,1,k) + ti2 + ti3 + cr2 = Cc(1,1,k) + tr11*tr2 + tr12*tr3 + ci2 = Cc(2,1,k) + tr11*ti2 + tr12*ti3 + cr3 = Cc(1,1,k) + tr12*tr2 + tr11*tr3 + ci3 = Cc(2,1,k) + tr12*ti2 + tr11*ti3 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + Ch(1,k,2) = cr2 - ci5 + Ch(1,k,5) = cr2 + ci5 + Ch(2,k,2) = ci2 + cr5 + Ch(2,k,3) = ci3 + cr4 + Ch(1,k,3) = cr3 - ci4 + Ch(1,k,4) = cr3 + ci4 + Ch(2,k,4) = ci3 - cr4 + Ch(2,k,5) = ci2 - cr5 + enddo + end if + end subroutine passb5 \ No newline at end of file diff --git a/src/passf.f b/src/passf.f deleted file mode 100644 index cb97d12..0000000 --- a/src/passf.f +++ /dev/null @@ -1,117 +0,0 @@ - SUBROUTINE PASSF (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (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) - IDOT = IDO/2 - NT = IP*IDL1 - IPP2 = IP+2 - IPPH = (IP+1)/2 - IDP = IP*IDO -C - IF (IDO .LT. L1) GO TO 106 - DO 103 J=2,IPPH - JC = IPP2-J - DO 102 K=1,L1 - DO 101 I=1,IDO - CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) - CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) - 101 CONTINUE - 102 CONTINUE - 103 CONTINUE - DO 105 K=1,L1 - DO 104 I=1,IDO - CH(I,K,1) = CC(I,1,K) - 104 CONTINUE - 105 CONTINUE - GO TO 112 - 106 DO 109 J=2,IPPH - JC = IPP2-J - DO 108 I=1,IDO - DO 107 K=1,L1 - CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) - CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) - 107 CONTINUE - 108 CONTINUE - 109 CONTINUE - DO 111 I=1,IDO - DO 110 K=1,L1 - CH(I,K,1) = CC(I,1,K) - 110 CONTINUE - 111 CONTINUE - 112 IDL = 2-IDO - INC = 0 - DO 116 L=2,IPPH - LC = IPP2-L - IDL = IDL+IDO - DO 113 IK=1,IDL1 - C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2) - C2(IK,LC) = -WA(IDL)*CH2(IK,IP) - 113 CONTINUE - IDLJ = IDL - INC = INC+IDO - DO 115 J=3,IPPH - JC = IPP2-J - IDLJ = IDLJ+INC - IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP - WAR = WA(IDLJ-1) - WAI = WA(IDLJ) - DO 114 IK=1,IDL1 - C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J) - C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC) - 114 CONTINUE - 115 CONTINUE - 116 CONTINUE - DO 118 J=2,IPPH - DO 117 IK=1,IDL1 - CH2(IK,1) = CH2(IK,1)+CH2(IK,J) - 117 CONTINUE - 118 CONTINUE - DO 120 J=2,IPPH - JC = IPP2-J - DO 119 IK=2,IDL1,2 - CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC) - CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC) - CH2(IK,J) = C2(IK,J)+C2(IK-1,JC) - CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC) - 119 CONTINUE - 120 CONTINUE - NAC = 1 - IF (IDO .EQ. 2) RETURN - NAC = 0 - DO 121 IK=1,IDL1 - C2(IK,1) = CH2(IK,1) - 121 CONTINUE - DO 123 J=2,IP - DO 122 K=1,L1 - C1(1,K,J) = CH(1,K,J) - C1(2,K,J) = CH(2,K,J) - 122 CONTINUE - 123 CONTINUE - IF (IDOT .GT. L1) GO TO 127 - IDIJ = 0 - DO 126 J=2,IP - IDIJ = IDIJ+2 - DO 125 I=4,IDO,2 - IDIJ = IDIJ+2 - DO 124 K=1,L1 - C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J) - C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J) - 124 CONTINUE - 125 CONTINUE - 126 CONTINUE - RETURN - 127 IDJ = 2-IDO - DO 130 J=2,IP - IDJ = IDJ+IDO - DO 129 K=1,L1 - IDIJ = IDJ - DO 128 I=4,IDO,2 - IDIJ = IDIJ+2 - C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J) - C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J) - 128 CONTINUE - 129 CONTINUE - 130 CONTINUE - RETURN - END diff --git a/src/passf.f90 b/src/passf.f90 new file mode 100644 index 0000000..ebf9278 --- /dev/null +++ b/src/passf.f90 @@ -0,0 +1,124 @@ + subroutine passf(Nac,Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) + use fftpack_kind + implicit none + real(rk) :: c1 , c2 , Cc , Ch , Ch2 , Wa , wai , war + integer :: i , idij , idj , idl , Idl1 , idlj , Ido , idot , idp , & + ik , inc , Ip , ipp2 , ipph , j , jc , k , l , l1 , lc + integer :: Nac , nt + dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , Wa(*) , & + & c2(Idl1,Ip) , Ch2(Idl1,Ip) + idot = Ido/2 + nt = Ip*Idl1 + ipp2 = Ip + 2 + ipph = (Ip+1)/2 + idp = Ip*Ido +! + if ( Idoidp ) idlj = idlj - idp + war = Wa(idlj-1) + wai = Wa(idlj) + do ik = 1 , Idl1 + c2(ik,l) = c2(ik,l) + war*Ch2(ik,j) + c2(ik,lc) = c2(ik,lc) - wai*Ch2(ik,jc) + enddo + enddo + enddo + do j = 2 , ipph + do ik = 1 , Idl1 + Ch2(ik,1) = Ch2(ik,1) + Ch2(ik,j) + enddo + enddo + do j = 2 , ipph + jc = ipp2 - j + do ik = 2 , Idl1 , 2 + Ch2(ik-1,j) = c2(ik-1,j) - c2(ik,jc) + Ch2(ik-1,jc) = c2(ik-1,j) + c2(ik,jc) + Ch2(ik,j) = c2(ik,j) + c2(ik-1,jc) + Ch2(ik,jc) = c2(ik,j) - c2(ik-1,jc) + enddo + enddo + Nac = 1 + if ( Ido==2 ) return + Nac = 0 + do ik = 1 , Idl1 + c2(ik,1) = Ch2(ik,1) + enddo + do j = 2 , Ip + do k = 1 , l1 + c1(1,k,j) = Ch(1,k,j) + c1(2,k,j) = Ch(2,k,j) + enddo + enddo + if ( idot>l1 ) then + idj = 2 - Ido + do j = 2 , Ip + idj = idj + Ido + do k = 1 , l1 + idij = idj + do i = 4 , Ido , 2 + idij = idij + 2 + c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) + Wa(idij) & + *Ch(i,k,j) + c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) - Wa(idij) & + *Ch(i-1,k,j) + enddo + enddo + enddo + else + idij = 0 + do j = 2 , Ip + idij = idij + 2 + do i = 4 , Ido , 2 + idij = idij + 2 + do k = 1 , l1 + c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) + Wa(idij)*Ch(i,k,j) + c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) - Wa(idij)*Ch(i-1,k,j) + enddo + enddo + enddo + end if + end subroutine passf \ No newline at end of file diff --git a/src/passf2.f b/src/passf2.f deleted file mode 100644 index 765b69f..0000000 --- a/src/passf2.f +++ /dev/null @@ -1,24 +0,0 @@ - SUBROUTINE PASSF2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , - 1 WA1(1) - IF (IDO .GT. 2) GO TO 102 - DO 101 K=1,L1 - CH(1,K,1) = CC(1,1,K)+CC(1,2,K) - CH(1,K,2) = CC(1,1,K)-CC(1,2,K) - CH(2,K,1) = CC(2,1,K)+CC(2,2,K) - CH(2,K,2) = CC(2,1,K)-CC(2,2,K) - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) - TR2 = CC(I-1,1,K)-CC(I-1,2,K) - CH(I,K,1) = CC(I,1,K)+CC(I,2,K) - TI2 = CC(I,1,K)-CC(I,2,K) - CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2 - CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passf2.f90 b/src/passf2.f90 new file mode 100644 index 0000000..4c9f17c --- /dev/null +++ b/src/passf2.f90 @@ -0,0 +1,26 @@ + subroutine passf2(Ido,l1,Cc,Ch,Wa1) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 + integer :: i , Ido , k , l1 + dimension Cc(Ido,2,l1) , Ch(Ido,l1,2) , Wa1(*) + if ( Ido>2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + Ch(i-1,k,1) = Cc(i-1,1,k) + Cc(i-1,2,k) + tr2 = Cc(i-1,1,k) - Cc(i-1,2,k) + Ch(i,k,1) = Cc(i,1,k) + Cc(i,2,k) + ti2 = Cc(i,1,k) - Cc(i,2,k) + Ch(i,k,2) = Wa1(i-1)*ti2 - Wa1(i)*tr2 + Ch(i-1,k,2) = Wa1(i-1)*tr2 + Wa1(i)*ti2 + enddo + enddo + else + do k = 1 , l1 + Ch(1,k,1) = Cc(1,1,k) + Cc(1,2,k) + Ch(1,k,2) = Cc(1,1,k) - Cc(1,2,k) + Ch(2,k,1) = Cc(2,1,k) + Cc(2,2,k) + Ch(2,k,2) = Cc(2,1,k) - Cc(2,2,k) + enddo + end if + end subroutine passf2 \ No newline at end of file diff --git a/src/passf3.f b/src/passf3.f deleted file mode 100644 index d3547ca..0000000 --- a/src/passf3.f +++ /dev/null @@ -1,44 +0,0 @@ - SUBROUTINE PASSF3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , - 1 WA1(1) ,WA2(1) -C *** TAUI IS -SQRT(3)/2 *** - DATA TAUR,TAUI /-0.5D0,-0.86602540378443864676D0/ - IF (IDO .NE. 2) GO TO 102 - DO 101 K=1,L1 - TR2 = CC(1,2,K)+CC(1,3,K) - CR2 = CC(1,1,K)+TAUR*TR2 - CH(1,K,1) = CC(1,1,K)+TR2 - TI2 = CC(2,2,K)+CC(2,3,K) - CI2 = CC(2,1,K)+TAUR*TI2 - CH(2,K,1) = CC(2,1,K)+TI2 - CR3 = TAUI*(CC(1,2,K)-CC(1,3,K)) - CI3 = TAUI*(CC(2,2,K)-CC(2,3,K)) - CH(1,K,2) = CR2-CI3 - CH(1,K,3) = CR2+CI3 - CH(2,K,2) = CI2+CR3 - CH(2,K,3) = CI2-CR3 - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - TR2 = CC(I-1,2,K)+CC(I-1,3,K) - CR2 = CC(I-1,1,K)+TAUR*TR2 - CH(I-1,K,1) = CC(I-1,1,K)+TR2 - TI2 = CC(I,2,K)+CC(I,3,K) - CI2 = CC(I,1,K)+TAUR*TI2 - CH(I,K,1) = CC(I,1,K)+TI2 - CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) - CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) - DR2 = CR2-CI3 - DR3 = CR2+CI3 - DI2 = CI2+CR3 - DI3 = CI2-CR3 - CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 - CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 - CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 - CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passf3.f90 b/src/passf3.f90 new file mode 100644 index 0000000..165b286 --- /dev/null +++ b/src/passf3.f90 @@ -0,0 +1,47 @@ + subroutine passf3(Ido,l1,Cc,Ch,Wa1,Wa2) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , cr2 , cr3 , di2 , di3 , & + & dr2 , dr3 , ti2 , tr2 , Wa1 , Wa2 + integer :: i , Ido , k , l1 + dimension Cc(Ido,3,l1) , Ch(Ido,l1,3) , Wa1(*) , Wa2(*) + real(rk),parameter :: taur = -0.5_rk + real(rk),parameter :: taui = -sqrt(3.0_rk) / 2.0_rk + if ( Ido/=2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + tr2 = Cc(i-1,2,k) + Cc(i-1,3,k) + cr2 = Cc(i-1,1,k) + taur*tr2 + Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + ti2 = Cc(i,2,k) + Cc(i,3,k) + ci2 = Cc(i,1,k) + taur*ti2 + Ch(i,k,1) = Cc(i,1,k) + ti2 + cr3 = taui*(Cc(i-1,2,k)-Cc(i-1,3,k)) + ci3 = taui*(Cc(i,2,k)-Cc(i,3,k)) + dr2 = cr2 - ci3 + dr3 = cr2 + ci3 + di2 = ci2 + cr3 + di3 = ci2 - cr3 + Ch(i,k,2) = Wa1(i-1)*di2 - Wa1(i)*dr2 + Ch(i-1,k,2) = Wa1(i-1)*dr2 + Wa1(i)*di2 + Ch(i,k,3) = Wa2(i-1)*di3 - Wa2(i)*dr3 + Ch(i-1,k,3) = Wa2(i-1)*dr3 + Wa2(i)*di3 + enddo + enddo + else + do k = 1 , l1 + tr2 = Cc(1,2,k) + Cc(1,3,k) + cr2 = Cc(1,1,k) + taur*tr2 + Ch(1,k,1) = Cc(1,1,k) + tr2 + ti2 = Cc(2,2,k) + Cc(2,3,k) + ci2 = Cc(2,1,k) + taur*ti2 + Ch(2,k,1) = Cc(2,1,k) + ti2 + cr3 = taui*(Cc(1,2,k)-Cc(1,3,k)) + ci3 = taui*(Cc(2,2,k)-Cc(2,3,k)) + Ch(1,k,2) = cr2 - ci3 + Ch(1,k,3) = cr2 + ci3 + Ch(2,k,2) = ci2 + cr3 + Ch(2,k,3) = ci2 - cr3 + enddo + end if + end subroutine passf3 \ No newline at end of file diff --git a/src/passf4.f b/src/passf4.f deleted file mode 100644 index d14ad61..0000000 --- a/src/passf4.f +++ /dev/null @@ -1,52 +0,0 @@ - SUBROUTINE PASSF4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (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 - DO 101 K=1,L1 - TI1 = CC(2,1,K)-CC(2,3,K) - TI2 = CC(2,1,K)+CC(2,3,K) - TR4 = CC(2,2,K)-CC(2,4,K) - TI3 = CC(2,2,K)+CC(2,4,K) - TR1 = CC(1,1,K)-CC(1,3,K) - TR2 = CC(1,1,K)+CC(1,3,K) - TI4 = CC(1,4,K)-CC(1,2,K) - TR3 = CC(1,2,K)+CC(1,4,K) - CH(1,K,1) = TR2+TR3 - CH(1,K,3) = TR2-TR3 - CH(2,K,1) = TI2+TI3 - CH(2,K,3) = TI2-TI3 - CH(1,K,2) = TR1+TR4 - CH(1,K,4) = TR1-TR4 - CH(2,K,2) = TI1+TI4 - CH(2,K,4) = TI1-TI4 - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - TI1 = CC(I,1,K)-CC(I,3,K) - TI2 = CC(I,1,K)+CC(I,3,K) - TI3 = CC(I,2,K)+CC(I,4,K) - TR4 = CC(I,2,K)-CC(I,4,K) - TR1 = CC(I-1,1,K)-CC(I-1,3,K) - TR2 = CC(I-1,1,K)+CC(I-1,3,K) - TI4 = CC(I-1,4,K)-CC(I-1,2,K) - TR3 = CC(I-1,2,K)+CC(I-1,4,K) - CH(I-1,K,1) = TR2+TR3 - CR3 = TR2-TR3 - CH(I,K,1) = TI2+TI3 - CI3 = TI2-TI3 - CR2 = TR1+TR4 - CR4 = TR1-TR4 - CI2 = TI1+TI4 - CI4 = TI1-TI4 - CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2 - CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2 - CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3 - CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3 - CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4 - CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passf4.f90 b/src/passf4.f90 new file mode 100644 index 0000000..110dea8 --- /dev/null +++ b/src/passf4.f90 @@ -0,0 +1,56 @@ + subroutine passf4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & + & ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3 , tr4 , & + & Wa1 , Wa2 , Wa3 + integer :: i , Ido , k , l1 + dimension Cc(Ido,4,l1) , Ch(Ido,l1,4) , Wa1(*) , Wa2(*) , Wa3(*) + if ( Ido/=2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + ti1 = Cc(i,1,k) - Cc(i,3,k) + ti2 = Cc(i,1,k) + Cc(i,3,k) + ti3 = Cc(i,2,k) + Cc(i,4,k) + tr4 = Cc(i,2,k) - Cc(i,4,k) + tr1 = Cc(i-1,1,k) - Cc(i-1,3,k) + tr2 = Cc(i-1,1,k) + Cc(i-1,3,k) + ti4 = Cc(i-1,4,k) - Cc(i-1,2,k) + tr3 = Cc(i-1,2,k) + Cc(i-1,4,k) + Ch(i-1,k,1) = tr2 + tr3 + cr3 = tr2 - tr3 + Ch(i,k,1) = ti2 + ti3 + ci3 = ti2 - ti3 + cr2 = tr1 + tr4 + cr4 = tr1 - tr4 + ci2 = ti1 + ti4 + ci4 = ti1 - ti4 + Ch(i-1,k,2) = Wa1(i-1)*cr2 + Wa1(i)*ci2 + Ch(i,k,2) = Wa1(i-1)*ci2 - Wa1(i)*cr2 + Ch(i-1,k,3) = Wa2(i-1)*cr3 + Wa2(i)*ci3 + Ch(i,k,3) = Wa2(i-1)*ci3 - Wa2(i)*cr3 + Ch(i-1,k,4) = Wa3(i-1)*cr4 + Wa3(i)*ci4 + Ch(i,k,4) = Wa3(i-1)*ci4 - Wa3(i)*cr4 + enddo + enddo + else + do k = 1 , l1 + ti1 = Cc(2,1,k) - Cc(2,3,k) + ti2 = Cc(2,1,k) + Cc(2,3,k) + tr4 = Cc(2,2,k) - Cc(2,4,k) + ti3 = Cc(2,2,k) + Cc(2,4,k) + tr1 = Cc(1,1,k) - Cc(1,3,k) + tr2 = Cc(1,1,k) + Cc(1,3,k) + ti4 = Cc(1,4,k) - Cc(1,2,k) + tr3 = Cc(1,2,k) + Cc(1,4,k) + Ch(1,k,1) = tr2 + tr3 + Ch(1,k,3) = tr2 - tr3 + Ch(2,k,1) = ti2 + ti3 + Ch(2,k,3) = ti2 - ti3 + Ch(1,k,2) = tr1 + tr4 + Ch(1,k,4) = tr1 - tr4 + Ch(2,k,2) = ti1 + ti4 + Ch(2,k,4) = ti1 - ti4 + enddo + end if + end subroutine passf4 \ No newline at end of file diff --git a/src/passf5.f b/src/passf5.f deleted file mode 100644 index 9e56a3b..0000000 --- a/src/passf5.f +++ /dev/null @@ -1,79 +0,0 @@ - SUBROUTINE PASSF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (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) -C *** TR12=-COS(4*PI/5), TI12=-SIN(4*PI/5) - DATA TR11,TI11,TR12,TI12 /0.3090169943749474241D0, - + -0.95105651629515357212D0, - 1 -0.8090169943749474241D0, -0.58778525229247312917D0/ - IF (IDO .NE. 2) GO TO 102 - DO 101 K=1,L1 - TI5 = CC(2,2,K)-CC(2,5,K) - TI2 = CC(2,2,K)+CC(2,5,K) - TI4 = CC(2,3,K)-CC(2,4,K) - TI3 = CC(2,3,K)+CC(2,4,K) - TR5 = CC(1,2,K)-CC(1,5,K) - TR2 = CC(1,2,K)+CC(1,5,K) - TR4 = CC(1,3,K)-CC(1,4,K) - TR3 = CC(1,3,K)+CC(1,4,K) - CH(1,K,1) = CC(1,1,K)+TR2+TR3 - CH(2,K,1) = CC(2,1,K)+TI2+TI3 - CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 - CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3 - CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 - CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3 - CR5 = TI11*TR5+TI12*TR4 - CI5 = TI11*TI5+TI12*TI4 - CR4 = TI12*TR5-TI11*TR4 - CI4 = TI12*TI5-TI11*TI4 - CH(1,K,2) = CR2-CI5 - CH(1,K,5) = CR2+CI5 - CH(2,K,2) = CI2+CR5 - CH(2,K,3) = CI3+CR4 - CH(1,K,3) = CR3-CI4 - CH(1,K,4) = CR3+CI4 - CH(2,K,4) = CI3-CR4 - CH(2,K,5) = CI2-CR5 - 101 CONTINUE - RETURN - 102 DO 104 K=1,L1 - DO 103 I=2,IDO,2 - TI5 = CC(I,2,K)-CC(I,5,K) - TI2 = CC(I,2,K)+CC(I,5,K) - TI4 = CC(I,3,K)-CC(I,4,K) - TI3 = CC(I,3,K)+CC(I,4,K) - TR5 = CC(I-1,2,K)-CC(I-1,5,K) - TR2 = CC(I-1,2,K)+CC(I-1,5,K) - TR4 = CC(I-1,3,K)-CC(I-1,4,K) - TR3 = CC(I-1,3,K)+CC(I-1,4,K) - CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 - CH(I,K,1) = CC(I,1,K)+TI2+TI3 - CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 - CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 - CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 - CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 - CR5 = TI11*TR5+TI12*TR4 - CI5 = TI11*TI5+TI12*TI4 - CR4 = TI12*TR5-TI11*TR4 - CI4 = TI12*TI5-TI11*TI4 - DR3 = CR3-CI4 - DR4 = CR3+CI4 - DI3 = CI3+CR4 - DI4 = CI3-CR4 - DR5 = CR2+CI5 - DR2 = CR2-CI5 - DI5 = CI2-CR5 - DI2 = CI2+CR5 - CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 - CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 - CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 - CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 - CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4 - CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4 - CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5 - CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5 - 103 CONTINUE - 104 CONTINUE - RETURN - END diff --git a/src/passf5.f90 b/src/passf5.f90 new file mode 100644 index 0000000..6c29b44 --- /dev/null +++ b/src/passf5.f90 @@ -0,0 +1,86 @@ + subroutine passf5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & + cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & + dr4 , dr5 + real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & + tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 + integer :: i , Ido , k , l1 + dimension Cc(Ido,5,l1) , Ch(Ido,l1,5) , Wa1(*) , Wa2(*) , Wa3(*), & + Wa4(*) + real(rk),parameter :: pi = acos(-1.0_rk) + real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti11 = -sin(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti12 = -sin(4.0_rk * pi / 5.0_rk) + if ( Ido/=2 ) then + do k = 1 , l1 + do i = 2 , Ido , 2 + ti5 = Cc(i,2,k) - Cc(i,5,k) + ti2 = Cc(i,2,k) + Cc(i,5,k) + ti4 = Cc(i,3,k) - Cc(i,4,k) + ti3 = Cc(i,3,k) + Cc(i,4,k) + tr5 = Cc(i-1,2,k) - Cc(i-1,5,k) + tr2 = Cc(i-1,2,k) + Cc(i-1,5,k) + tr4 = Cc(i-1,3,k) - Cc(i-1,4,k) + tr3 = Cc(i-1,3,k) + Cc(i-1,4,k) + Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + tr3 + Ch(i,k,1) = Cc(i,1,k) + ti2 + ti3 + cr2 = Cc(i-1,1,k) + tr11*tr2 + tr12*tr3 + ci2 = Cc(i,1,k) + tr11*ti2 + tr12*ti3 + cr3 = Cc(i-1,1,k) + tr12*tr2 + tr11*tr3 + ci3 = Cc(i,1,k) + tr12*ti2 + tr11*ti3 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + dr3 = cr3 - ci4 + dr4 = cr3 + ci4 + di3 = ci3 + cr4 + di4 = ci3 - cr4 + dr5 = cr2 + ci5 + dr2 = cr2 - ci5 + di5 = ci2 - cr5 + di2 = ci2 + cr5 + Ch(i-1,k,2) = Wa1(i-1)*dr2 + Wa1(i)*di2 + Ch(i,k,2) = Wa1(i-1)*di2 - Wa1(i)*dr2 + Ch(i-1,k,3) = Wa2(i-1)*dr3 + Wa2(i)*di3 + Ch(i,k,3) = Wa2(i-1)*di3 - Wa2(i)*dr3 + Ch(i-1,k,4) = Wa3(i-1)*dr4 + Wa3(i)*di4 + Ch(i,k,4) = Wa3(i-1)*di4 - Wa3(i)*dr4 + Ch(i-1,k,5) = Wa4(i-1)*dr5 + Wa4(i)*di5 + Ch(i,k,5) = Wa4(i-1)*di5 - Wa4(i)*dr5 + enddo + enddo + else + do k = 1 , l1 + ti5 = Cc(2,2,k) - Cc(2,5,k) + ti2 = Cc(2,2,k) + Cc(2,5,k) + ti4 = Cc(2,3,k) - Cc(2,4,k) + ti3 = Cc(2,3,k) + Cc(2,4,k) + tr5 = Cc(1,2,k) - Cc(1,5,k) + tr2 = Cc(1,2,k) + Cc(1,5,k) + tr4 = Cc(1,3,k) - Cc(1,4,k) + tr3 = Cc(1,3,k) + Cc(1,4,k) + Ch(1,k,1) = Cc(1,1,k) + tr2 + tr3 + Ch(2,k,1) = Cc(2,1,k) + ti2 + ti3 + cr2 = Cc(1,1,k) + tr11*tr2 + tr12*tr3 + ci2 = Cc(2,1,k) + tr11*ti2 + tr12*ti3 + cr3 = Cc(1,1,k) + tr12*tr2 + tr11*tr3 + ci3 = Cc(2,1,k) + tr12*ti2 + tr11*ti3 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + Ch(1,k,2) = cr2 - ci5 + Ch(1,k,5) = cr2 + ci5 + Ch(2,k,2) = ci2 + cr5 + Ch(2,k,3) = ci3 + cr4 + Ch(1,k,3) = cr3 - ci4 + Ch(1,k,4) = cr3 + ci4 + Ch(2,k,4) = ci3 - cr4 + Ch(2,k,5) = ci2 - cr5 + enddo + end if + end subroutine passf5 \ No newline at end of file diff --git a/src/radb2.f b/src/radb2.f deleted file mode 100644 index e8f18a1..0000000 --- a/src/radb2.f +++ /dev/null @@ -1,28 +0,0 @@ - SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , - 1 WA1(1) - DO 101 K=1,L1 - CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K) - CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K) - 101 CONTINUE - IF (IDO-2) 107,105,102 - 102 IDP2 = IDO+2 - DO 104 K=1,L1 - DO 103 I=3,IDO,2 - IC = IDP2-I - CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K) - TR2 = CC(I-1,1,K)-CC(IC-1,2,K) - CH(I,K,1) = CC(I,1,K)-CC(IC,2,K) - TI2 = CC(I,1,K)+CC(IC,2,K) - CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2 - CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2 - 103 CONTINUE - 104 CONTINUE - IF (MOD(IDO,2) .EQ. 1) RETURN - 105 DO 106 K=1,L1 - CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K) - CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K)) - 106 CONTINUE - 107 RETURN - END diff --git a/src/radb2.f90 b/src/radb2.f90 new file mode 100644 index 0000000..b37708a --- /dev/null +++ b/src/radb2.f90 @@ -0,0 +1,31 @@ + subroutine radb2(Ido,l1,Cc,Ch,Wa1) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Cc(Ido,2,l1) , Ch(Ido,l1,2) , Wa1(*) + do k = 1 , l1 + Ch(1,k,1) = Cc(1,1,k) + Cc(Ido,2,k) + Ch(1,k,2) = Cc(1,1,k) - Cc(Ido,2,k) + enddo + if ( Ido<2 ) return + if ( Ido/=2 ) then + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + Ch(i-1,k,1) = Cc(i-1,1,k) + Cc(ic-1,2,k) + tr2 = Cc(i-1,1,k) - Cc(ic-1,2,k) + Ch(i,k,1) = Cc(i,1,k) - Cc(ic,2,k) + ti2 = Cc(i,1,k) + Cc(ic,2,k) + Ch(i-1,k,2) = Wa1(i-2)*tr2 - Wa1(i-1)*ti2 + Ch(i,k,2) = Wa1(i-2)*ti2 + Wa1(i-1)*tr2 + enddo + enddo + if ( mod(Ido,2)==1 ) return + endif + do k = 1 , l1 + Ch(Ido,k,1) = Cc(Ido,1,k) + Cc(Ido,1,k) + Ch(Ido,k,2) = -(Cc(1,2,k)+Cc(1,2,k)) + enddo + end subroutine radb2 \ No newline at end of file diff --git a/src/radb3.f b/src/radb3.f deleted file mode 100644 index 6f2b689..0000000 --- a/src/radb3.f +++ /dev/null @@ -1,39 +0,0 @@ - SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , - 1 WA1(1) ,WA2(1) -C *** TAUI IS SQRT(3)/2 *** - DATA TAUR,TAUI /-0.5D0,0.86602540378443864676D0/ - DO 101 K=1,L1 - TR2 = CC(IDO,2,K)+CC(IDO,2,K) - CR2 = CC(1,1,K)+TAUR*TR2 - CH(1,K,1) = CC(1,1,K)+TR2 - CI3 = TAUI*(CC(1,3,K)+CC(1,3,K)) - CH(1,K,2) = CR2-CI3 - CH(1,K,3) = CR2+CI3 - 101 CONTINUE - IF (IDO .EQ. 1) RETURN - IDP2 = IDO+2 - DO 103 K=1,L1 - DO 102 I=3,IDO,2 - IC = IDP2-I - TR2 = CC(I-1,3,K)+CC(IC-1,2,K) - CR2 = CC(I-1,1,K)+TAUR*TR2 - CH(I-1,K,1) = CC(I-1,1,K)+TR2 - TI2 = CC(I,3,K)-CC(IC,2,K) - CI2 = CC(I,1,K)+TAUR*TI2 - CH(I,K,1) = CC(I,1,K)+TI2 - CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K)) - CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K)) - DR2 = CR2-CI3 - DR3 = CR2+CI3 - DI2 = CI2+CR3 - DI3 = CI2-CR3 - CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2 - CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2 - CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3 - CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3 - 102 CONTINUE - 103 CONTINUE - RETURN - END diff --git a/src/radb3.f90 b/src/radb3.f90 new file mode 100644 index 0000000..f2c3e23 --- /dev/null +++ b/src/radb3.f90 @@ -0,0 +1,41 @@ + subroutine radb3(Ido,l1,Cc,Ch,Wa1,Wa2) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , cr2 , cr3 , di2 , di3 , & + dr2 , dr3 , ti2 , tr2 , Wa1 , Wa2 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Cc(Ido,3,l1) , Ch(Ido,l1,3) , Wa1(*) , Wa2(*) + real(rk),parameter :: taur = - 0.5_rk + real(rk),parameter :: taui = sqrt(3.0_rk) / 2.0_rk + do k = 1 , l1 + tr2 = Cc(Ido,2,k) + Cc(Ido,2,k) + cr2 = Cc(1,1,k) + taur*tr2 + Ch(1,k,1) = Cc(1,1,k) + tr2 + ci3 = taui*(Cc(1,3,k)+Cc(1,3,k)) + Ch(1,k,2) = cr2 - ci3 + Ch(1,k,3) = cr2 + ci3 + enddo + if ( Ido==1 ) return + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + tr2 = Cc(i-1,3,k) + Cc(ic-1,2,k) + cr2 = Cc(i-1,1,k) + taur*tr2 + Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + ti2 = Cc(i,3,k) - Cc(ic,2,k) + ci2 = Cc(i,1,k) + taur*ti2 + Ch(i,k,1) = Cc(i,1,k) + ti2 + cr3 = taui*(Cc(i-1,3,k)-Cc(ic-1,2,k)) + ci3 = taui*(Cc(i,3,k)+Cc(ic,2,k)) + dr2 = cr2 - ci3 + dr3 = cr2 + ci3 + di2 = ci2 + cr3 + di3 = ci2 - cr3 + Ch(i-1,k,2) = Wa1(i-2)*dr2 - Wa1(i-1)*di2 + Ch(i,k,2) = Wa1(i-2)*di2 + Wa1(i-1)*dr2 + Ch(i-1,k,3) = Wa2(i-2)*dr3 - Wa2(i-1)*di3 + Ch(i,k,3) = Wa2(i-2)*di3 + Wa2(i-1)*dr3 + enddo + enddo + end subroutine radb3 \ No newline at end of file diff --git a/src/radb4.f b/src/radb4.f deleted file mode 100644 index 6917fb5..0000000 --- a/src/radb4.f +++ /dev/null @@ -1,58 +0,0 @@ - SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , - 1 WA1(1) ,WA2(1) ,WA3(1) - DATA SQRT2 /1.41421356237309504880D0/ - DO 101 K=1,L1 - TR1 = CC(1,1,K)-CC(IDO,4,K) - TR2 = CC(1,1,K)+CC(IDO,4,K) - TR3 = CC(IDO,2,K)+CC(IDO,2,K) - TR4 = CC(1,3,K)+CC(1,3,K) - CH(1,K,1) = TR2+TR3 - CH(1,K,2) = TR1-TR4 - CH(1,K,3) = TR2-TR3 - CH(1,K,4) = TR1+TR4 - 101 CONTINUE - IF (IDO-2) 107,105,102 - 102 IDP2 = IDO+2 - DO 104 K=1,L1 - DO 103 I=3,IDO,2 - IC = IDP2-I - TI1 = CC(I,1,K)+CC(IC,4,K) - TI2 = CC(I,1,K)-CC(IC,4,K) - TI3 = CC(I,3,K)-CC(IC,2,K) - TR4 = CC(I,3,K)+CC(IC,2,K) - TR1 = CC(I-1,1,K)-CC(IC-1,4,K) - TR2 = CC(I-1,1,K)+CC(IC-1,4,K) - TI4 = CC(I-1,3,K)-CC(IC-1,2,K) - TR3 = CC(I-1,3,K)+CC(IC-1,2,K) - CH(I-1,K,1) = TR2+TR3 - CR3 = TR2-TR3 - CH(I,K,1) = TI2+TI3 - CI3 = TI2-TI3 - CR2 = TR1-TR4 - CR4 = TR1+TR4 - CI2 = TI1+TI4 - CI4 = TI1-TI4 - CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2 - CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2 - CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3 - CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3 - CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4 - CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4 - 103 CONTINUE - 104 CONTINUE - IF (MOD(IDO,2) .EQ. 1) RETURN - 105 CONTINUE - DO 106 K=1,L1 - TI1 = CC(1,2,K)+CC(1,4,K) - TI2 = CC(1,4,K)-CC(1,2,K) - TR1 = CC(IDO,1,K)-CC(IDO,3,K) - TR2 = CC(IDO,1,K)+CC(IDO,3,K) - CH(IDO,K,1) = TR2+TR2 - CH(IDO,K,2) = SQRT2*(TR1-TI1) - CH(IDO,K,3) = TI2+TI2 - CH(IDO,K,4) = -SQRT2*(TR1+TI1) - 106 CONTINUE - 107 RETURN - END diff --git a/src/radb4.f90 b/src/radb4.f90 new file mode 100644 index 0000000..e6fadf5 --- /dev/null +++ b/src/radb4.f90 @@ -0,0 +1,62 @@ + subroutine radb4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & + ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3, & + tr4 , Wa1 , Wa2 , Wa3 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Cc(Ido,4,l1) , Ch(Ido,l1,4) , Wa1(*) , Wa2(*) , Wa3(*) + real(rk),parameter :: sqrt2 = sqrt(2.0_rk) + do k = 1 , l1 + tr1 = Cc(1,1,k) - Cc(Ido,4,k) + tr2 = Cc(1,1,k) + Cc(Ido,4,k) + tr3 = Cc(Ido,2,k) + Cc(Ido,2,k) + tr4 = Cc(1,3,k) + Cc(1,3,k) + Ch(1,k,1) = tr2 + tr3 + Ch(1,k,2) = tr1 - tr4 + Ch(1,k,3) = tr2 - tr3 + Ch(1,k,4) = tr1 + tr4 + enddo + if ( Ido<2 ) return + if ( Ido/=2 ) then + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + ti1 = Cc(i,1,k) + Cc(ic,4,k) + ti2 = Cc(i,1,k) - Cc(ic,4,k) + ti3 = Cc(i,3,k) - Cc(ic,2,k) + tr4 = Cc(i,3,k) + Cc(ic,2,k) + tr1 = Cc(i-1,1,k) - Cc(ic-1,4,k) + tr2 = Cc(i-1,1,k) + Cc(ic-1,4,k) + ti4 = Cc(i-1,3,k) - Cc(ic-1,2,k) + tr3 = Cc(i-1,3,k) + Cc(ic-1,2,k) + Ch(i-1,k,1) = tr2 + tr3 + cr3 = tr2 - tr3 + Ch(i,k,1) = ti2 + ti3 + ci3 = ti2 - ti3 + cr2 = tr1 - tr4 + cr4 = tr1 + tr4 + ci2 = ti1 + ti4 + ci4 = ti1 - ti4 + Ch(i-1,k,2) = Wa1(i-2)*cr2 - Wa1(i-1)*ci2 + Ch(i,k,2) = Wa1(i-2)*ci2 + Wa1(i-1)*cr2 + Ch(i-1,k,3) = Wa2(i-2)*cr3 - Wa2(i-1)*ci3 + Ch(i,k,3) = Wa2(i-2)*ci3 + Wa2(i-1)*cr3 + Ch(i-1,k,4) = Wa3(i-2)*cr4 - Wa3(i-1)*ci4 + Ch(i,k,4) = Wa3(i-2)*ci4 + Wa3(i-1)*cr4 + enddo + enddo + if ( mod(Ido,2)==1 ) return + endif + do k = 1 , l1 + ti1 = Cc(1,2,k) + Cc(1,4,k) + ti2 = Cc(1,4,k) - Cc(1,2,k) + tr1 = Cc(Ido,1,k) - Cc(Ido,3,k) + tr2 = Cc(Ido,1,k) + Cc(Ido,3,k) + Ch(Ido,k,1) = tr2 + tr2 + Ch(Ido,k,2) = sqrt2*(tr1-ti1) + Ch(Ido,k,3) = ti2 + ti2 + Ch(Ido,k,4) = -sqrt2*(tr1+ti1) + enddo + end subroutine radb4 \ No newline at end of file diff --git a/src/radb5.f b/src/radb5.f deleted file mode 100644 index e1a8664..0000000 --- a/src/radb5.f +++ /dev/null @@ -1,67 +0,0 @@ - SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (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) -C *** TR12=COS(4*PI/5), TI12=SIN(4*PI/5) - DATA TR11,TI11,TR12,TI12 /0.3090169943749474241D0, - + 0.95105651629515357212D0, - + -0.8090169943749474241D0,0.58778525229247312917D0/ - DO 101 K=1,L1 - TI5 = CC(1,3,K)+CC(1,3,K) - TI4 = CC(1,5,K)+CC(1,5,K) - TR2 = CC(IDO,2,K)+CC(IDO,2,K) - TR3 = CC(IDO,4,K)+CC(IDO,4,K) - CH(1,K,1) = CC(1,1,K)+TR2+TR3 - CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 - CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 - CI5 = TI11*TI5+TI12*TI4 - CI4 = TI12*TI5-TI11*TI4 - CH(1,K,2) = CR2-CI5 - CH(1,K,3) = CR3-CI4 - CH(1,K,4) = CR3+CI4 - CH(1,K,5) = CR2+CI5 - 101 CONTINUE - IF (IDO .EQ. 1) RETURN - IDP2 = IDO+2 - DO 103 K=1,L1 - DO 102 I=3,IDO,2 - IC = IDP2-I - TI5 = CC(I,3,K)+CC(IC,2,K) - TI2 = CC(I,3,K)-CC(IC,2,K) - TI4 = CC(I,5,K)+CC(IC,4,K) - TI3 = CC(I,5,K)-CC(IC,4,K) - TR5 = CC(I-1,3,K)-CC(IC-1,2,K) - TR2 = CC(I-1,3,K)+CC(IC-1,2,K) - TR4 = CC(I-1,5,K)-CC(IC-1,4,K) - TR3 = CC(I-1,5,K)+CC(IC-1,4,K) - CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 - CH(I,K,1) = CC(I,1,K)+TI2+TI3 - CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 - CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 - CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 - CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 - CR5 = TI11*TR5+TI12*TR4 - CI5 = TI11*TI5+TI12*TI4 - CR4 = TI12*TR5-TI11*TR4 - CI4 = TI12*TI5-TI11*TI4 - DR3 = CR3-CI4 - DR4 = CR3+CI4 - DI3 = CI3+CR4 - DI4 = CI3-CR4 - DR5 = CR2+CI5 - DR2 = CR2-CI5 - DI5 = CI2-CR5 - DI2 = CI2+CR5 - CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2 - CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2 - CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3 - CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3 - CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4 - CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4 - CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5 - CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5 - 102 CONTINUE - 103 CONTINUE - RETURN - END diff --git a/src/radb5.f90 b/src/radb5.f90 new file mode 100644 index 0000000..e90a4d8 --- /dev/null +++ b/src/radb5.f90 @@ -0,0 +1,73 @@ + subroutine radb5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & + cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & + dr4 , dr5 + real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & + tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Cc(Ido,5,l1) , Ch(Ido,l1,5) , Wa1(*) , Wa2(*) , Wa3(*), & + Wa4(*) + real(rk),parameter :: pi = acos(-1.0_rk) + real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti11 = sin(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti12 = sin(4.0_rk * pi / 5.0_rk) + do k = 1 , l1 + ti5 = Cc(1,3,k) + Cc(1,3,k) + ti4 = Cc(1,5,k) + Cc(1,5,k) + tr2 = Cc(Ido,2,k) + Cc(Ido,2,k) + tr3 = Cc(Ido,4,k) + Cc(Ido,4,k) + Ch(1,k,1) = Cc(1,1,k) + tr2 + tr3 + cr2 = Cc(1,1,k) + tr11*tr2 + tr12*tr3 + cr3 = Cc(1,1,k) + tr12*tr2 + tr11*tr3 + ci5 = ti11*ti5 + ti12*ti4 + ci4 = ti12*ti5 - ti11*ti4 + Ch(1,k,2) = cr2 - ci5 + Ch(1,k,3) = cr3 - ci4 + Ch(1,k,4) = cr3 + ci4 + Ch(1,k,5) = cr2 + ci5 + enddo + if ( Ido==1 ) return + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + ti5 = Cc(i,3,k) + Cc(ic,2,k) + ti2 = Cc(i,3,k) - Cc(ic,2,k) + ti4 = Cc(i,5,k) + Cc(ic,4,k) + ti3 = Cc(i,5,k) - Cc(ic,4,k) + tr5 = Cc(i-1,3,k) - Cc(ic-1,2,k) + tr2 = Cc(i-1,3,k) + Cc(ic-1,2,k) + tr4 = Cc(i-1,5,k) - Cc(ic-1,4,k) + tr3 = Cc(i-1,5,k) + Cc(ic-1,4,k) + Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + tr3 + Ch(i,k,1) = Cc(i,1,k) + ti2 + ti3 + cr2 = Cc(i-1,1,k) + tr11*tr2 + tr12*tr3 + ci2 = Cc(i,1,k) + tr11*ti2 + tr12*ti3 + cr3 = Cc(i-1,1,k) + tr12*tr2 + tr11*tr3 + ci3 = Cc(i,1,k) + tr12*ti2 + tr11*ti3 + cr5 = ti11*tr5 + ti12*tr4 + ci5 = ti11*ti5 + ti12*ti4 + cr4 = ti12*tr5 - ti11*tr4 + ci4 = ti12*ti5 - ti11*ti4 + dr3 = cr3 - ci4 + dr4 = cr3 + ci4 + di3 = ci3 + cr4 + di4 = ci3 - cr4 + dr5 = cr2 + ci5 + dr2 = cr2 - ci5 + di5 = ci2 - cr5 + di2 = ci2 + cr5 + Ch(i-1,k,2) = Wa1(i-2)*dr2 - Wa1(i-1)*di2 + Ch(i,k,2) = Wa1(i-2)*di2 + Wa1(i-1)*dr2 + Ch(i-1,k,3) = Wa2(i-2)*dr3 - Wa2(i-1)*di3 + Ch(i,k,3) = Wa2(i-2)*di3 + Wa2(i-1)*dr3 + Ch(i-1,k,4) = Wa3(i-2)*dr4 - Wa3(i-1)*di4 + Ch(i,k,4) = Wa3(i-2)*di4 + Wa3(i-1)*dr4 + Ch(i-1,k,5) = Wa4(i-2)*dr5 - Wa4(i-1)*di5 + Ch(i,k,5) = Wa4(i-2)*di5 + Wa4(i-1)*dr5 + enddo + enddo + end subroutine radb5 \ No newline at end of file diff --git a/src/radbg.f b/src/radbg.f deleted file mode 100644 index deaed58..0000000 --- a/src/radbg.f +++ /dev/null @@ -1,161 +0,0 @@ - SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (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) - DATA TPI/6.28318530717958647692D0/ - ARG = TPI/FLOAT(IP) - DCP = COS(ARG) - DSP = SIN(ARG) - IDP2 = IDO+2 - NBD = (IDO-1)/2 - IPP2 = IP+2 - IPPH = (IP+1)/2 - IF (IDO .LT. L1) GO TO 103 - DO 102 K=1,L1 - DO 101 I=1,IDO - CH(I,K,1) = CC(I,1,K) - 101 CONTINUE - 102 CONTINUE - GO TO 106 - 103 DO 105 I=1,IDO - DO 104 K=1,L1 - CH(I,K,1) = CC(I,1,K) - 104 CONTINUE - 105 CONTINUE - 106 DO 108 J=2,IPPH - JC = IPP2-J - J2 = J+J - DO 107 K=1,L1 - CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K) - CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K) - 107 CONTINUE - 108 CONTINUE - IF (IDO .EQ. 1) GO TO 116 - IF (NBD .LT. L1) GO TO 112 - DO 111 J=2,IPPH - JC = IPP2-J - DO 110 K=1,L1 - DO 109 I=3,IDO,2 - IC = IDP2-I - CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K) - CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K) - CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K) - CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K) - 109 CONTINUE - 110 CONTINUE - 111 CONTINUE - GO TO 116 - 112 DO 115 J=2,IPPH - JC = IPP2-J - DO 114 I=3,IDO,2 - IC = IDP2-I - DO 113 K=1,L1 - CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K) - CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K) - CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K) - CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K) - 113 CONTINUE - 114 CONTINUE - 115 CONTINUE - 116 AR1 = 1.0D0 - AI1 = 0.0D0 - DO 120 L=2,IPPH - LC = IPP2-L - AR1H = DCP*AR1-DSP*AI1 - AI1 = DCP*AI1+DSP*AR1 - AR1 = AR1H - DO 117 IK=1,IDL1 - C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2) - C2(IK,LC) = AI1*CH2(IK,IP) - 117 CONTINUE - DC2 = AR1 - DS2 = AI1 - AR2 = AR1 - AI2 = AI1 - DO 119 J=3,IPPH - JC = IPP2-J - AR2H = DC2*AR2-DS2*AI2 - AI2 = DC2*AI2+DS2*AR2 - AR2 = AR2H - DO 118 IK=1,IDL1 - C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J) - C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC) - 118 CONTINUE - 119 CONTINUE - 120 CONTINUE - DO 122 J=2,IPPH - DO 121 IK=1,IDL1 - CH2(IK,1) = CH2(IK,1)+CH2(IK,J) - 121 CONTINUE - 122 CONTINUE - DO 124 J=2,IPPH - JC = IPP2-J - DO 123 K=1,L1 - CH(1,K,J) = C1(1,K,J)-C1(1,K,JC) - CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC) - 123 CONTINUE - 124 CONTINUE - IF (IDO .EQ. 1) GO TO 132 - IF (NBD .LT. L1) GO TO 128 - DO 127 J=2,IPPH - JC = IPP2-J - DO 126 K=1,L1 - DO 125 I=3,IDO,2 - CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC) - CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC) - CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC) - CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC) - 125 CONTINUE - 126 CONTINUE - 127 CONTINUE - GO TO 132 - 128 DO 131 J=2,IPPH - JC = IPP2-J - DO 130 I=3,IDO,2 - DO 129 K=1,L1 - CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC) - CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC) - CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC) - CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC) - 129 CONTINUE - 130 CONTINUE - 131 CONTINUE - 132 CONTINUE - IF (IDO .EQ. 1) RETURN - DO 133 IK=1,IDL1 - C2(IK,1) = CH2(IK,1) - 133 CONTINUE - DO 135 J=2,IP - DO 134 K=1,L1 - C1(1,K,J) = CH(1,K,J) - 134 CONTINUE - 135 CONTINUE - IF (NBD .GT. L1) GO TO 139 - IS = -IDO - DO 138 J=2,IP - IS = IS+IDO - IDIJ = IS - DO 137 I=3,IDO,2 - IDIJ = IDIJ+2 - DO 136 K=1,L1 - C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) - C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) - 136 CONTINUE - 137 CONTINUE - 138 CONTINUE - GO TO 143 - 139 IS = -IDO - DO 142 J=2,IP - IS = IS+IDO - DO 141 K=1,L1 - IDIJ = IS - DO 140 I=3,IDO,2 - IDIJ = IDIJ+2 - C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) - C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) - 140 CONTINUE - 141 CONTINUE - 142 CONTINUE - 143 RETURN - END diff --git a/src/radbg.f90 b/src/radbg.f90 new file mode 100644 index 0000000..55a6af4 --- /dev/null +++ b/src/radbg.f90 @@ -0,0 +1,174 @@ + subroutine radbg(Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) + use fftpack_kind + implicit none + real(rk) :: ai1 , ai2 , ar1 , ar1h , ar2 , ar2h , arg , c1 , & + c2 , Cc , Ch , Ch2 , dc2 , dcp , ds2 , dsp , & + Wa + integer :: i , ic , idij , Idl1 , Ido , idp2 , ik , Ip , ipp2 , & + ipph , is , j , j2 , jc , k , l , l1 , lc , nbd + dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , & + c2(Idl1,Ip) , Ch2(Idl1,Ip) , Wa(*) + real(rk),parameter :: tpi = 2*acos(-1.0_rk) ! 2 * pi + arg = tpi/real(Ip, rk) + dcp = cos(arg) + dsp = sin(arg) + idp2 = Ido + 2 + nbd = (Ido-1)/2 + ipp2 = Ip + 2 + ipph = (Ip+1)/2 + if ( Idol1 ) then + is = -Ido + do j = 2 , Ip + is = is + Ido + do k = 1 , l1 + idij = is + do i = 3 , Ido , 2 + idij = idij + 2 + c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij) & + *Ch(i,k,j) + c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij) & + *Ch(i-1,k,j) + enddo + enddo + enddo + else + is = -Ido + do j = 2 , Ip + is = is + Ido + idij = is + do i = 3 , Ido , 2 + idij = idij + 2 + do k = 1 , l1 + c1(i-1,k,j) = Wa(idij-1)*Ch(i-1,k,j) - Wa(idij) & + *Ch(i,k,j) + c1(i,k,j) = Wa(idij-1)*Ch(i,k,j) + Wa(idij) & + *Ch(i-1,k,j) + enddo + enddo + enddo + endif + end subroutine radbg \ No newline at end of file diff --git a/src/radf2.f b/src/radf2.f deleted file mode 100644 index 831682f..0000000 --- a/src/radf2.f +++ /dev/null @@ -1,28 +0,0 @@ - SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) , - 1 WA1(1) - DO 101 K=1,L1 - CH(1,1,K) = CC(1,K,1)+CC(1,K,2) - CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2) - 101 CONTINUE - IF (IDO-2) 107,105,102 - 102 IDP2 = IDO+2 - DO 104 K=1,L1 - DO 103 I=3,IDO,2 - IC = IDP2-I - TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) - TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) - CH(I,1,K) = CC(I,K,1)+TI2 - CH(IC,2,K) = TI2-CC(I,K,1) - CH(I-1,1,K) = CC(I-1,K,1)+TR2 - CH(IC-1,2,K) = CC(I-1,K,1)-TR2 - 103 CONTINUE - 104 CONTINUE - IF (MOD(IDO,2) .EQ. 1) RETURN - 105 DO 106 K=1,L1 - CH(1,2,K) = -CC(IDO,K,2) - CH(IDO,1,K) = CC(IDO,K,1) - 106 CONTINUE - 107 RETURN - END diff --git a/src/radf2.f90 b/src/radf2.f90 new file mode 100644 index 0000000..51a6fc4 --- /dev/null +++ b/src/radf2.f90 @@ -0,0 +1,31 @@ + subroutine radf2(Ido,l1,Cc,Ch,Wa1) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ti2 , tr2 , Wa1 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Ch(Ido,2,l1) , Cc(Ido,l1,2) , Wa1(*) + do k = 1 , l1 + Ch(1,1,k) = Cc(1,k,1) + Cc(1,k,2) + Ch(Ido,2,k) = Cc(1,k,1) - Cc(1,k,2) + enddo + if ( Ido<2 ) return + if ( Ido/=2 ) then + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + tr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) + ti2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) + Ch(i,1,k) = Cc(i,k,1) + ti2 + Ch(ic,2,k) = ti2 - Cc(i,k,1) + Ch(i-1,1,k) = Cc(i-1,k,1) + tr2 + Ch(ic-1,2,k) = Cc(i-1,k,1) - tr2 + enddo + enddo + if ( mod(Ido,2)==1 ) return + endif + do k = 1 , l1 + Ch(1,2,k) = -Cc(Ido,k,2) + Ch(Ido,1,k) = Cc(Ido,k,1) + enddo + end subroutine radf2 \ No newline at end of file diff --git a/src/radf3.f b/src/radf3.f deleted file mode 100644 index f6f402f..0000000 --- a/src/radf3.f +++ /dev/null @@ -1,37 +0,0 @@ - SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) , - 1 WA1(1) ,WA2(1) -C *** TAUI IS -SQRT(3)/2 *** - DATA TAUR,TAUI /-0.5D0,0.86602540378443864676D0/ - DO 101 K=1,L1 - CR2 = CC(1,K,2)+CC(1,K,3) - CH(1,1,K) = CC(1,K,1)+CR2 - CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2)) - CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2 - 101 CONTINUE - IF (IDO .EQ. 1) RETURN - IDP2 = IDO+2 - DO 103 K=1,L1 - DO 102 I=3,IDO,2 - IC = IDP2-I - DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) - DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) - DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) - DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) - CR2 = DR2+DR3 - CI2 = DI2+DI3 - CH(I-1,1,K) = CC(I-1,K,1)+CR2 - CH(I,1,K) = CC(I,K,1)+CI2 - TR2 = CC(I-1,K,1)+TAUR*CR2 - TI2 = CC(I,K,1)+TAUR*CI2 - TR3 = TAUI*(DI2-DI3) - TI3 = TAUI*(DR3-DR2) - CH(I-1,3,K) = TR2+TR3 - CH(IC-1,2,K) = TR2-TR3 - CH(I,3,K) = TI2+TI3 - CH(IC,2,K) = TI3-TI2 - 102 CONTINUE - 103 CONTINUE - RETURN - END diff --git a/src/radf3.f90 b/src/radf3.f90 new file mode 100644 index 0000000..472489a --- /dev/null +++ b/src/radf3.f90 @@ -0,0 +1,40 @@ + subroutine radf3(Ido,l1,Cc,Ch,Wa1,Wa2) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , cr2 , di2 , di3 , dr2 , dr3 , & + ti2 , ti3 , tr2 , tr3 , Wa1 , Wa2 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Ch(Ido,3,l1) , Cc(Ido,l1,3) , Wa1(*) , Wa2(*) + real(rk),parameter :: taur = -0.5_rk + ! note: original comment said this was -SQRT(3)/2 but value was 0.86602540378443864676d0 + real(rk),parameter :: taui = sqrt(3.0_rk) / 2.0_rk + do k = 1 , l1 + cr2 = Cc(1,k,2) + Cc(1,k,3) + Ch(1,1,k) = Cc(1,k,1) + cr2 + Ch(1,3,k) = taui*(Cc(1,k,3)-Cc(1,k,2)) + Ch(Ido,2,k) = Cc(1,k,1) + taur*cr2 + enddo + if ( Ido==1 ) return + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + dr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) + di2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) + dr3 = Wa2(i-2)*Cc(i-1,k,3) + Wa2(i-1)*Cc(i,k,3) + di3 = Wa2(i-2)*Cc(i,k,3) - Wa2(i-1)*Cc(i-1,k,3) + cr2 = dr2 + dr3 + ci2 = di2 + di3 + Ch(i-1,1,k) = Cc(i-1,k,1) + cr2 + Ch(i,1,k) = Cc(i,k,1) + ci2 + tr2 = Cc(i-1,k,1) + taur*cr2 + ti2 = Cc(i,k,1) + taur*ci2 + tr3 = taui*(di2-di3) + ti3 = taui*(dr3-dr2) + Ch(i-1,3,k) = tr2 + tr3 + Ch(ic-1,2,k) = tr2 - tr3 + Ch(i,3,k) = ti2 + ti3 + Ch(ic,2,k) = ti3 - ti2 + enddo + enddo + end subroutine radf3 \ No newline at end of file diff --git a/src/radf4.f b/src/radf4.f deleted file mode 100644 index 1a33060..0000000 --- a/src/radf4.f +++ /dev/null @@ -1,54 +0,0 @@ - SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3) - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) , - 1 WA1(1) ,WA2(1) ,WA3(1) - DATA HSQT2 /0.70710678118654752440D0/ - DO 101 K=1,L1 - TR1 = CC(1,K,2)+CC(1,K,4) - TR2 = CC(1,K,1)+CC(1,K,3) - CH(1,1,K) = TR1+TR2 - CH(IDO,4,K) = TR2-TR1 - CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3) - CH(1,3,K) = CC(1,K,4)-CC(1,K,2) - 101 CONTINUE - IF (IDO-2) 107,105,102 - 102 IDP2 = IDO+2 - DO 104 K=1,L1 - DO 103 I=3,IDO,2 - IC = IDP2-I - CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) - CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) - CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) - CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) - CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) - CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) - TR1 = CR2+CR4 - TR4 = CR4-CR2 - TI1 = CI2+CI4 - TI4 = CI2-CI4 - TI2 = CC(I,K,1)+CI3 - TI3 = CC(I,K,1)-CI3 - TR2 = CC(I-1,K,1)+CR3 - TR3 = CC(I-1,K,1)-CR3 - CH(I-1,1,K) = TR1+TR2 - CH(IC-1,4,K) = TR2-TR1 - CH(I,1,K) = TI1+TI2 - CH(IC,4,K) = TI1-TI2 - CH(I-1,3,K) = TI4+TR3 - CH(IC-1,2,K) = TR3-TI4 - CH(I,3,K) = TR4+TI3 - CH(IC,2,K) = TR4-TI3 - 103 CONTINUE - 104 CONTINUE - IF (MOD(IDO,2) .EQ. 1) RETURN - 105 CONTINUE - DO 106 K=1,L1 - TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4)) - TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4)) - CH(IDO,1,K) = TR1+CC(IDO,K,1) - CH(IDO,3,K) = CC(IDO,K,1)-TR1 - CH(1,2,K) = TI1-CC(IDO,K,3) - CH(1,4,K) = TI1+CC(IDO,K,3) - 106 CONTINUE - 107 RETURN - END diff --git a/src/radf4.f90 b/src/radf4.f90 new file mode 100644 index 0000000..c0b62bb --- /dev/null +++ b/src/radf4.f90 @@ -0,0 +1,58 @@ + subroutine radf4(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , cr2 , cr3 , cr4 , & + ti1 , ti2 , ti3 , ti4 , tr1 , tr2 , tr3, & + tr4 , Wa1 , Wa2 , Wa3 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Cc(Ido,l1,4) , Ch(Ido,4,l1) , Wa1(*) , Wa2(*) , Wa3(*) + real(rk),parameter :: hsqt2 = sqrt(2.0_rk) / 2.0_rk + do k = 1 , l1 + tr1 = Cc(1,k,2) + Cc(1,k,4) + tr2 = Cc(1,k,1) + Cc(1,k,3) + Ch(1,1,k) = tr1 + tr2 + Ch(Ido,4,k) = tr2 - tr1 + Ch(Ido,2,k) = Cc(1,k,1) - Cc(1,k,3) + Ch(1,3,k) = Cc(1,k,4) - Cc(1,k,2) + enddo + if ( Ido<2 ) return + if ( Ido/=2 ) then + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + cr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) + ci2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) + cr3 = Wa2(i-2)*Cc(i-1,k,3) + Wa2(i-1)*Cc(i,k,3) + ci3 = Wa2(i-2)*Cc(i,k,3) - Wa2(i-1)*Cc(i-1,k,3) + cr4 = Wa3(i-2)*Cc(i-1,k,4) + Wa3(i-1)*Cc(i,k,4) + ci4 = Wa3(i-2)*Cc(i,k,4) - Wa3(i-1)*Cc(i-1,k,4) + tr1 = cr2 + cr4 + tr4 = cr4 - cr2 + ti1 = ci2 + ci4 + ti4 = ci2 - ci4 + ti2 = Cc(i,k,1) + ci3 + ti3 = Cc(i,k,1) - ci3 + tr2 = Cc(i-1,k,1) + cr3 + tr3 = Cc(i-1,k,1) - cr3 + Ch(i-1,1,k) = tr1 + tr2 + Ch(ic-1,4,k) = tr2 - tr1 + Ch(i,1,k) = ti1 + ti2 + Ch(ic,4,k) = ti1 - ti2 + Ch(i-1,3,k) = ti4 + tr3 + Ch(ic-1,2,k) = tr3 - ti4 + Ch(i,3,k) = tr4 + ti3 + Ch(ic,2,k) = tr4 - ti3 + enddo + enddo + if ( mod(Ido,2)==1 ) return + endif + do k = 1 , l1 + ti1 = -hsqt2*(Cc(Ido,k,2)+Cc(Ido,k,4)) + tr1 = hsqt2*(Cc(Ido,k,2)-Cc(Ido,k,4)) + Ch(Ido,1,k) = tr1 + Cc(Ido,k,1) + Ch(Ido,3,k) = Cc(Ido,k,1) - tr1 + Ch(1,2,k) = ti1 - Cc(Ido,k,3) + Ch(1,4,k) = ti1 + Cc(Ido,k,3) + enddo + end subroutine radf4 \ No newline at end of file diff --git a/src/radf5.f b/src/radf5.f deleted file mode 100644 index e45521f..0000000 --- a/src/radf5.f +++ /dev/null @@ -1,61 +0,0 @@ - SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) - IMPLICIT DOUBLE PRECISION (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, - + 0.95105651629515357212D0, - 1 -0.8090169943749474241D0, 0.58778525229247312917D0/ - DO 101 K=1,L1 - CR2 = CC(1,K,5)+CC(1,K,2) - CI5 = CC(1,K,5)-CC(1,K,2) - CR3 = CC(1,K,4)+CC(1,K,3) - CI4 = CC(1,K,4)-CC(1,K,3) - CH(1,1,K) = CC(1,K,1)+CR2+CR3 - CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3 - CH(1,3,K) = TI11*CI5+TI12*CI4 - CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3 - CH(1,5,K) = TI12*CI5-TI11*CI4 - 101 CONTINUE - IF (IDO .EQ. 1) RETURN - IDP2 = IDO+2 - DO 103 K=1,L1 - DO 102 I=3,IDO,2 - IC = IDP2-I - DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) - DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) - DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) - DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) - DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) - DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) - DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5) - DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5) - CR2 = DR2+DR5 - CI5 = DR5-DR2 - CR5 = DI2-DI5 - CI2 = DI2+DI5 - CR3 = DR3+DR4 - CI4 = DR4-DR3 - CR4 = DI3-DI4 - CI3 = DI3+DI4 - CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3 - CH(I,1,K) = CC(I,K,1)+CI2+CI3 - TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3 - TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3 - TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3 - TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3 - TR5 = TI11*CR5+TI12*CR4 - TI5 = TI11*CI5+TI12*CI4 - TR4 = TI12*CR5-TI11*CR4 - TI4 = TI12*CI5-TI11*CI4 - CH(I-1,3,K) = TR2+TR5 - CH(IC-1,2,K) = TR2-TR5 - CH(I,3,K) = TI2+TI5 - CH(IC,2,K) = TI5-TI2 - CH(I-1,5,K) = TR3+TR4 - CH(IC-1,4,K) = TR3-TR4 - CH(I,5,K) = TI3+TI4 - CH(IC,4,K) = TI4-TI3 - 102 CONTINUE - 103 CONTINUE - RETURN - END diff --git a/src/radf5.f90 b/src/radf5.f90 new file mode 100644 index 0000000..5fbcb11 --- /dev/null +++ b/src/radf5.f90 @@ -0,0 +1,69 @@ + subroutine radf5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4) + use fftpack_kind + implicit none + real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , & + cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , & + dr4 , dr5 + real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, & + tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4 + integer :: i , ic , Ido , idp2 , k , l1 + dimension Cc(Ido,l1,5) , Ch(Ido,5,l1) , Wa1(*) , Wa2(*) , Wa3(*), & + Wa4(*) + real(rk),parameter :: pi = acos(-1.0_rk) + real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti11 = sin(2.0_rk * pi / 5.0_rk) + real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk) + real(rk),parameter :: ti12 = sin(4.0_rk * pi / 5.0_rk) + do k = 1 , l1 + cr2 = Cc(1,k,5) + Cc(1,k,2) + ci5 = Cc(1,k,5) - Cc(1,k,2) + cr3 = Cc(1,k,4) + Cc(1,k,3) + ci4 = Cc(1,k,4) - Cc(1,k,3) + Ch(1,1,k) = Cc(1,k,1) + cr2 + cr3 + Ch(Ido,2,k) = Cc(1,k,1) + tr11*cr2 + tr12*cr3 + Ch(1,3,k) = ti11*ci5 + ti12*ci4 + Ch(Ido,4,k) = Cc(1,k,1) + tr12*cr2 + tr11*cr3 + Ch(1,5,k) = ti12*ci5 - ti11*ci4 + enddo + if ( Ido==1 ) return + idp2 = Ido + 2 + do k = 1 , l1 + do i = 3 , Ido , 2 + ic = idp2 - i + dr2 = Wa1(i-2)*Cc(i-1,k,2) + Wa1(i-1)*Cc(i,k,2) + di2 = Wa1(i-2)*Cc(i,k,2) - Wa1(i-1)*Cc(i-1,k,2) + dr3 = Wa2(i-2)*Cc(i-1,k,3) + Wa2(i-1)*Cc(i,k,3) + di3 = Wa2(i-2)*Cc(i,k,3) - Wa2(i-1)*Cc(i-1,k,3) + dr4 = Wa3(i-2)*Cc(i-1,k,4) + Wa3(i-1)*Cc(i,k,4) + di4 = Wa3(i-2)*Cc(i,k,4) - Wa3(i-1)*Cc(i-1,k,4) + dr5 = Wa4(i-2)*Cc(i-1,k,5) + Wa4(i-1)*Cc(i,k,5) + di5 = Wa4(i-2)*Cc(i,k,5) - Wa4(i-1)*Cc(i-1,k,5) + cr2 = dr2 + dr5 + ci5 = dr5 - dr2 + cr5 = di2 - di5 + ci2 = di2 + di5 + cr3 = dr3 + dr4 + ci4 = dr4 - dr3 + cr4 = di3 - di4 + ci3 = di3 + di4 + Ch(i-1,1,k) = Cc(i-1,k,1) + cr2 + cr3 + Ch(i,1,k) = Cc(i,k,1) + ci2 + ci3 + tr2 = Cc(i-1,k,1) + tr11*cr2 + tr12*cr3 + ti2 = Cc(i,k,1) + tr11*ci2 + tr12*ci3 + tr3 = Cc(i-1,k,1) + tr12*cr2 + tr11*cr3 + ti3 = Cc(i,k,1) + tr12*ci2 + tr11*ci3 + tr5 = ti11*cr5 + ti12*cr4 + ti5 = ti11*ci5 + ti12*ci4 + tr4 = ti12*cr5 - ti11*cr4 + ti4 = ti12*ci5 - ti11*ci4 + Ch(i-1,3,k) = tr2 + tr5 + Ch(ic-1,2,k) = tr2 - tr5 + Ch(i,3,k) = ti2 + ti5 + Ch(ic,2,k) = ti5 - ti2 + Ch(i-1,5,k) = tr3 + tr4 + Ch(ic-1,4,k) = tr3 - tr4 + Ch(i,5,k) = ti3 + ti4 + Ch(ic,4,k) = ti4 - ti3 + enddo + enddo + end subroutine radf5 \ No newline at end of file diff --git a/src/radfg.f b/src/radfg.f deleted file mode 100644 index f0ad326..0000000 --- a/src/radfg.f +++ /dev/null @@ -1,167 +0,0 @@ - SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) - IMPLICIT DOUBLE PRECISION (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) - DATA TPI/6.28318530717958647692D0/ - ARG = TPI/FLOAT(IP) - DCP = COS(ARG) - DSP = SIN(ARG) - IPPH = (IP+1)/2 - IPP2 = IP+2 - IDP2 = IDO+2 - NBD = (IDO-1)/2 - IF (IDO .EQ. 1) GO TO 119 - DO 101 IK=1,IDL1 - CH2(IK,1) = C2(IK,1) - 101 CONTINUE - DO 103 J=2,IP - DO 102 K=1,L1 - CH(1,K,J) = C1(1,K,J) - 102 CONTINUE - 103 CONTINUE - IF (NBD .GT. L1) GO TO 107 - IS = -IDO - DO 106 J=2,IP - IS = IS+IDO - IDIJ = IS - DO 105 I=3,IDO,2 - IDIJ = IDIJ+2 - DO 104 K=1,L1 - CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) - CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) - 104 CONTINUE - 105 CONTINUE - 106 CONTINUE - GO TO 111 - 107 IS = -IDO - DO 110 J=2,IP - IS = IS+IDO - DO 109 K=1,L1 - IDIJ = IS - DO 108 I=3,IDO,2 - IDIJ = IDIJ+2 - CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) - CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) - 108 CONTINUE - 109 CONTINUE - 110 CONTINUE - 111 IF (NBD .LT. L1) GO TO 115 - DO 114 J=2,IPPH - JC = IPP2-J - DO 113 K=1,L1 - DO 112 I=3,IDO,2 - C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) - C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) - C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) - C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) - 112 CONTINUE - 113 CONTINUE - 114 CONTINUE - GO TO 121 - 115 DO 118 J=2,IPPH - JC = IPP2-J - DO 117 I=3,IDO,2 - DO 116 K=1,L1 - C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) - C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) - C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) - C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) - 116 CONTINUE - 117 CONTINUE - 118 CONTINUE - GO TO 121 - 119 DO 120 IK=1,IDL1 - C2(IK,1) = CH2(IK,1) - 120 CONTINUE - 121 DO 123 J=2,IPPH - JC = IPP2-J - DO 122 K=1,L1 - C1(1,K,J) = CH(1,K,J)+CH(1,K,JC) - C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J) - 122 CONTINUE - 123 CONTINUE -C - AR1 = 1.0D0 - AI1 = 0.0D0 - DO 127 L=2,IPPH - LC = IPP2-L - AR1H = DCP*AR1-DSP*AI1 - AI1 = DCP*AI1+DSP*AR1 - AR1 = AR1H - DO 124 IK=1,IDL1 - CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2) - CH2(IK,LC) = AI1*C2(IK,IP) - 124 CONTINUE - DC2 = AR1 - DS2 = AI1 - AR2 = AR1 - AI2 = AI1 - DO 126 J=3,IPPH - JC = IPP2-J - AR2H = DC2*AR2-DS2*AI2 - AI2 = DC2*AI2+DS2*AR2 - AR2 = AR2H - DO 125 IK=1,IDL1 - CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J) - CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC) - 125 CONTINUE - 126 CONTINUE - 127 CONTINUE - DO 129 J=2,IPPH - DO 128 IK=1,IDL1 - CH2(IK,1) = CH2(IK,1)+C2(IK,J) - 128 CONTINUE - 129 CONTINUE -C - IF (IDO .LT. L1) GO TO 132 - DO 131 K=1,L1 - DO 130 I=1,IDO - CC(I,1,K) = CH(I,K,1) - 130 CONTINUE - 131 CONTINUE - GO TO 135 - 132 DO 134 I=1,IDO - DO 133 K=1,L1 - CC(I,1,K) = CH(I,K,1) - 133 CONTINUE - 134 CONTINUE - 135 DO 137 J=2,IPPH - JC = IPP2-J - J2 = J+J - DO 136 K=1,L1 - CC(IDO,J2-2,K) = CH(1,K,J) - CC(1,J2-1,K) = CH(1,K,JC) - 136 CONTINUE - 137 CONTINUE - IF (IDO .EQ. 1) RETURN - IF (NBD .LT. L1) GO TO 141 - DO 140 J=2,IPPH - JC = IPP2-J - J2 = J+J - DO 139 K=1,L1 - DO 138 I=3,IDO,2 - IC = IDP2-I - CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) - CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) - CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) - CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) - 138 CONTINUE - 139 CONTINUE - 140 CONTINUE - RETURN - 141 DO 144 J=2,IPPH - JC = IPP2-J - J2 = J+J - DO 143 I=3,IDO,2 - IC = IDP2-I - DO 142 K=1,L1 - CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) - CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) - CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) - CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) - 142 CONTINUE - 143 CONTINUE - 144 CONTINUE - RETURN - END diff --git a/src/radfg.f90 b/src/radfg.f90 new file mode 100644 index 0000000..589a466 --- /dev/null +++ b/src/radfg.f90 @@ -0,0 +1,180 @@ + subroutine radfg(Ido,Ip,l1,Idl1,Cc,c1,c2,Ch,Ch2,Wa) + use fftpack_kind + implicit none + real(rk) :: ai1 , ai2 , ar1 , ar1h , ar2 , ar2h , arg , c1 , & + c2 , Cc , Ch , Ch2 , dc2 , dcp , ds2 , dsp , & + Wa + integer :: i , ic , idij , Idl1 , Ido , idp2 , ik , Ip , ipp2 , & + ipph , is , j , j2 , jc , k , l , l1 , lc , nbd + dimension Ch(Ido,l1,Ip) , Cc(Ido,Ip,l1) , c1(Ido,l1,Ip) , & + c2(Idl1,Ip) , Ch2(Idl1,Ip) , Wa(*) + real(rk),parameter :: tpi = 2.0_rk * acos(-1.0_rk) ! 2 * pi + arg = tpi/real(Ip, rk) + dcp = cos(arg) + dsp = sin(arg) + ipph = (Ip+1)/2 + ipp2 = Ip + 2 + idp2 = Ido + 2 + nbd = (Ido-1)/2 + if ( Ido==1 ) then + do ik = 1 , Idl1 + c2(ik,1) = Ch2(ik,1) + enddo + else + do ik = 1 , Idl1 + Ch2(ik,1) = c2(ik,1) + enddo + do j = 2 , Ip + do k = 1 , l1 + Ch(1,k,j) = c1(1,k,j) + enddo + enddo + if ( nbd>l1 ) then + is = -Ido + do j = 2 , Ip + is = is + Ido + do k = 1 , l1 + idij = is + do i = 3 , Ido , 2 + idij = idij + 2 + Ch(i-1,k,j) = Wa(idij-1)*c1(i-1,k,j) + Wa(idij) & + *c1(i,k,j) + Ch(i,k,j) = Wa(idij-1)*c1(i,k,j) - Wa(idij) & + *c1(i-1,k,j) + enddo + enddo + enddo + else + is = -Ido + do j = 2 , Ip + is = is + Ido + idij = is + do i = 3 , Ido , 2 + idij = idij + 2 + do k = 1 , l1 + Ch(i-1,k,j) = Wa(idij-1)*c1(i-1,k,j) + Wa(idij) & + *c1(i,k,j) + Ch(i,k,j) = Wa(idij-1)*c1(i,k,j) - Wa(idij) & + *c1(i-1,k,j) + enddo + enddo + enddo + endif + if ( nbd real64 - 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_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 dcosti(4, w) call dcost(4, x, w) - call check(all(x == [real(kind=dp) :: 15, -4, 0, -1.0000000000000009_dp]), msg="`dcosti` failed.") + call check(all(x == [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk]), msg="`dcosti` failed.") call dcost(4, x, w) - call check(all(x/(2.0_dp*(4.0_dp - 1.0_dp)) == [real(kind=dp) :: 1, 2, 3, 4]), msg="`dcost` failed.") + call check(all(x/(2.0_rk*(4.0_rk - 1.0_rk)) == [real(kind=rk) :: 1, 2, 3, 4]), msg="`dcost` failed.") end subroutine test_fftpack_dcost_real diff --git a/test/test_fftpack_dct.f90 b/test/test_fftpack_dct.f90 index ef1c6f8..f2a2307 100644 --- a/test/test_fftpack_dct.f90 +++ b/test/test_fftpack_dct.f90 @@ -14,29 +14,29 @@ end subroutine check subroutine test_fftpack_dct use fftpack, only: dct - use iso_fortran_env, only: dp => real64 + use fftpack_kind - real(kind=dp) :: x(3) = [9, -9, 3] + real(kind=rk) :: x(3) = [9, -9, 3] - call check(all(dct(x, 2) == [real(kind=dp) :: 0, 18]), msg="`dct(x, 2)` failed.") + call check(all(dct(x, 2) == [real(kind=rk) :: 0, 18]), msg="`dct(x, 2)` failed.") call check(all(dct(x, 3) == dct(x)), msg="`dct(x, 3)` failed.") - call check(all(dct(x, 4) == [real(kind=dp) :: -3, -3.0000000000000036_dp, 15, 33]), msg="`dct(x, 4)` failed.") + call check(all(dct(x, 4) == [real(kind=rk) :: -3, -3.0000000000000036_rk, 15, 33]), msg="`dct(x, 4)` failed.") end subroutine test_fftpack_dct subroutine test_fftpack_idct use fftpack, only: dct, idct - use iso_fortran_env, only: dp => real64 - real(kind=dp) :: eps = 1.0e-10_dp - real(kind=dp) :: x(4) = [1, 2, 3, 4] - - call check(all(idct(dct(x))/(2.0_dp*(4.0_dp - 1.0_dp)) == [real(kind=dp) :: 1, 2, 3, 4]), & - msg="`idct(dct(x))/(2.0_dp*(4.0_dp-1.0_dp))` failed.") - call check(all(idct(dct(x), 2)/(2.0_dp*(2.0_dp - 1.0_dp)) == [real(kind=dp) :: 5.5, 9.5]), & - msg="`idct(dct(x), 2)/(2.0_dp*(2.0_dp-1.0_dp))` failed.") - call check(all(idct(dct(x, 2), 4)/(2.0_dp*(4.0_dp - 1.0_dp)) == & - [0.16666666666666666_dp, 0.33333333333333331_dp, 0.66666666666666663_dp, 0.83333333333333315_dp]), & - msg="`idct(dct(x, 2), 4)/(2.0_dp*(4.0_dp-1.0_dp))` failed.") + use iso_fortran_env, only: rk => real64 + real(kind=rk) :: eps = 1.0e-10_rk + real(kind=rk) :: x(4) = [1, 2, 3, 4] + + call check(all(idct(dct(x))/(2.0_rk*(4.0_rk - 1.0_rk)) == [real(kind=rk) :: 1, 2, 3, 4]), & + msg="`idct(dct(x))/(2.0_rk*(4.0_rk-1.0_rk))` failed.") + call check(all(idct(dct(x), 2)/(2.0_rk*(2.0_rk - 1.0_rk)) == [real(kind=rk) :: 5.5, 9.5]), & + msg="`idct(dct(x), 2)/(2.0_rk*(2.0_rk-1.0_rk))` failed.") + call check(all(idct(dct(x, 2), 4)/(2.0_rk*(4.0_rk - 1.0_rk)) == & + [0.16666666666666666_rk, 0.33333333333333331_rk, 0.66666666666666663_rk, 0.83333333333333315_rk]), & + msg="`idct(dct(x, 2), 4)/(2.0_rk*(4.0_rk-1.0_rk))` failed.") end subroutine test_fftpack_idct 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 diff --git a/test/tstfft.f b/test/tstfft.f index 5111c69..36a9a3f 100644 --- a/test/tstfft.f +++ b/test/tstfft.f @@ -53,23 +53,24 @@ 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 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 @@ -84,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 @@ -111,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) @@ -143,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 @@ -172,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 @@ -204,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 @@ -219,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 @@ -247,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 @@ -263,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 @@ -291,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 @@ -329,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 @@ -352,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 @@ -372,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