Merge pull request #18 from jacobwilliams/refactor

Refactor

Author: certik

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" <test_name, see `fpm.toml` or list>
+fpm build
+fpm test --list
+fpm test <test_name, see `fpm.toml` or list>
 ```
 To use `fftpack` within your `fpm` project, add the following to your `fpm.toml` file:
 ```toml

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

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
+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

src/cfftb1.f

@@ -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