Add support for different real kinds

Makefile

@@ -3,21 +3,24 @@
 LIB = dfftpack
 
 FC = gfortran
-FFLAGS = -O2
+FFLAGS = -O2 -fPIC
 
 export LIB
 export FC
 export FFLAGS
 
-.PHONY: all clean test
+.PHONY: build clean test
 
-all:
-	$(MAKE) -f Makefile --directory=src
-	$(MAKE) -f Makefile --directory=test
+build:
+	$(MAKE) -f Makefile $@ --directory=src
 
-test:
-	$(MAKE) -f Makefile --directory=test
+test: build
+	$(MAKE) -f Makefile $@ --directory=test
+
+bench: build
+	$(MAKE) -f Makefile $@ --directory=example
 
 clean:
-	$(MAKE) -f Makefile clean --directory=src
-	$(MAKE) -f Makefile clean --directory=test
+	$(MAKE) -f Makefile $@ --directory=src
+	$(MAKE) -f Makefile $@ --directory=test
+	$(MAKE) -f Makefile $@ --directory=example

README.md

@@ -31,6 +31,7 @@ fftpack = { git="https://github.com/fortran-lang/fftpack.git" }
 Alternatively, you can build using provided `Makefile`:
 ```bash
 make
+make test
 ```
 
 ## Links

example/Makefile

@@ -0,0 +1,18 @@
+SRCF90 = \
+	bench1.f90
+
+OBJ = $(SRCF90:%.f90=%.o) 
+
+build: bench1.x
+
+bench: build
+	./bench1.x
+
+bench1.x: $(OBJ)
+	$(FC) $(FFLAGS) $(OBJ) -L../src -l$(LIB) -I../src -o $@
+
+%.o: %.f90
+	$(FC) $(FFLAGS) -I../src -c $<
+
+clean:
+	rm -f -r *.o *.x

src/Makefile

@@ -1,4 +1,4 @@
-SRCF = \
+SRCF90 = \
 	zfftb.f90\
 	cfftb1.f90\
 	zfftf.f90\
@@ -47,9 +47,7 @@ SRCF = \
 	dsinqi.f90\
 	dsint.f90\
 	sint1.f90\
-	dsinti.f90
-
-SRCF90 = \
+	dsinti.f90\
 	fftpack.f90\
 	fftpack_fft.f90\
 	fftpack_ifft.f90\
@@ -59,21 +57,28 @@ SRCF90 = \
 	fftpack_ifftshift.f90\
 	fftpack_qct.f90\
 	fftpack_iqct.f90\
-	fftpack_dct.f90\
-	rk.f90
+	fftpack_dct.f90
+
+SRCFPP = \
+	rk.F90
 
-OBJF := $(SRCF:.f90=.o)
-OBJF90 := $(SRCF90:.f90=.o)
+OBJ = $(SRCF90:.f90=.o)
+OBJ += $(SRCFPP:.F90=.o)
 
-lib$(LIB).a: $(OBJF) $(OBJF90)
-	ar -rcs lib$(LIB).a $(OBJF) $(OBJF90)
+build: lib$(LIB).a lib$(LIB).so
 
-shared: $(OBJ)
+lib$(LIB).a: $(OBJ)
+	ar -rcs lib$(LIB).a $(OBJ)
+
+lib$(LIB).so: $(OBJ)
 	$(FC) -shared -o lib$(LIB).so $(OBJ)
 
 clean:
 	rm -f -r *.o *.a *.so *.mod *.smod
 
+%.o: %.F90
+	$(FC) $(FFLAGS) -c $<
+
 %.o: %.f90
 	$(FC) $(FFLAGS) -c $<
 

src/rk.F90

@@ -0,0 +1,16 @@
+!> fftpack kind
+module fftpack_kind
+    implicit none
+
+    !> fftpack real kind
+#if defined(fftpack_sp)
+    integer, parameter :: rk = selected_real_kind(6)
+#elif defined(fftpack_xdp)
+    integer, parameter :: rk = selected_real_kind(18)
+#elif defined(fftpack_qp)
+    integer, parameter :: rk = selected_real_kind(33)
+#else
+    integer, parameter :: rk = selected_real_kind(15)
+#endif
+
+end module fftpack_kind

test/Makefile

@@ -1,35 +1,39 @@
 FETCH = curl -L
 
-SRC = \
+SRCF90 = \
 	test_fftpack_fft.f90 \
 	test_fftpack_rfft.f90 \
-	test_fftpack_qct.f90 \
-	test_fftpack_dct.f90 \
 	test_fftpack_utils.f90 \
-	test_fftpack.f90 \
+	test_fftpack.f90
+
+SRCFPP = \
+	test_fftpack_qct.F90 \
+	test_fftpack_dct.F90 \
 	testdrive.F90
-
-OBJ = $(SRC:.f90=.o)
-OBJ := $(OBJ:.F90=.o)
 
-all: tstfft \
-	test_fftpack
+OBJ = $(SRCF90:%.f90=%.o) 
+OBJ += $(SRCFPP:%.F90=%.o)
+
+build: tstfft.x \
+	test_fftpack.x
+
+test: build
+	./tstfft.x
+	./test_fftpack.x
 
 # Orginal test
-tstfft: tstfft.f
-	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
-	time ./tstfft.x
+tstfft.x: tstfft.f
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o $@
 
 # `fftpack` fft routines
-test_fftpack: $(OBJ)
-	$(FC) $(FFLAGS) $(OBJ) -L../src -l$(LIB) -I../src -o [email protected]
-	./test_fftpack.x
+test_fftpack.x: $(OBJ)
+	$(FC) $(FFLAGS) $(OBJ) -L../src -l$(LIB) -I../src -o $@
 
 testdrive.F90:
 	$(FETCH) https://github.com/fortran-lang/test-drive/raw/v0.4.0/src/testdrive.F90 > $@
 
 %.o: %.F90
-	$(FC) $(FFLAGS) -c $<
+	$(FC) $(FFLAGS) -I../src -c $<
 
 %.o: %.f90
 	$(FC) $(FFLAGS) -I../src -c $<

test/test_fftpack_dct.f90 → test/test_fftpack_dct.F90

@@ -7,6 +7,12 @@ module test_fftpack_dct
 
     public :: collect_dct
 
+#if defined(fftpack_sp)
+    real(kind=rk) :: eps = 1.0e-5_rk
+#else
+    real(kind=rk) :: eps = 1.0e-10_rk
+#endif
+
 contains
 
     !> Collect all exported unit tests
@@ -26,15 +32,16 @@ subroutine test_classic_dct(error)
         type(error_type), allocatable, intent(out) :: error
         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(error, all(x == [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk]), "`dcosti` failed.")
+        call check(error, sum(abs(x - [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk])) < eps, &
+                   "`dcosti` failed.")
         if (allocated(error)) return
 
         call dcost(4, x, w)
-        call check(error, all(x/(2.0_rk*(4.0_rk - 1.0_rk)) == [real(kind=rk) :: 1, 2, 3, 4]), "`dcost` failed.")
+        call check(error, sum(abs(x/(2.0_rk*(4.0_rk - 1.0_rk)) - &
+                                  [real(kind=rk) :: 1, 2, 3, 4])) < eps, "`dcost` failed.")
 
     end subroutine test_classic_dct
 
@@ -46,23 +53,25 @@ subroutine test_modernized_dct(error)
         if (allocated(error)) return
         call check(error, all(dct(x, 3) == dct(x)), "`dct(x, 3)` failed.")
         if (allocated(error)) return
-        call check(error, all(dct(x, 4) == [real(kind=rk) :: -3, -3.0000000000000036_rk, 15, 33]), "`dct(x, 4)` failed.")
+        call check(error, sum(abs(dct(x, 4) - [real(kind=rk) :: -3, -3.0000000000000036_rk, 15, 33])) &
+                   < eps, "`dct(x, 4)` failed.")
 
     end subroutine test_modernized_dct
 
     subroutine test_modernized_idct(error)
         type(error_type), allocatable, intent(out) :: error
-        real(kind=rk) :: eps = 1.0e-10_rk
         real(kind=rk) :: x(4) = [1, 2, 3, 4]
 
-        call check(error, all(idct(dct(x))/(2.0_rk*(4.0_rk - 1.0_rk)) == [real(kind=rk) :: 1, 2, 3, 4]), &
+        call check(error, sum(abs(idct(dct(x))/(2.0_rk*(4.0_rk - 1.0_rk)) - &
+                                  [real(kind=rk) :: 1, 2, 3, 4])) < eps, &
                    "`idct(dct(x))/(2.0_rk*(4.0_rk-1.0_rk))` failed.")
         if (allocated(error)) return
         call check(error, all(idct(dct(x), 2)/(2.0_rk*(2.0_rk - 1.0_rk)) == [real(kind=rk) :: 5.5, 9.5]), &
                    "`idct(dct(x), 2)/(2.0_rk*(2.0_rk-1.0_rk))` failed.")
         if (allocated(error)) return
-        call check(error, 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]), &
+        call check(error, sum(abs(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])) < eps, &
                    "`idct(dct(x, 2), 4)/(2.0_rk*(4.0_rk-1.0_rk))` failed.")
 
     end subroutine test_modernized_idct

test/test_fftpack_qct.f90 → test/test_fftpack_qct.F90

@@ -7,6 +7,12 @@ module test_fftpack_qct
 
     public :: collect_qct
 
+#if defined(fftpack_sp)
+    real(kind=rk) :: eps = 1.0e-5_rk
+#else
+    real(kind=rk) :: eps = 1.0e-10_rk
+#endif
+
 contains
 
     !> Collect all exported unit tests
@@ -26,7 +32,6 @@ subroutine test_classic_qct(error)
         type(error_type), allocatable, intent(out) :: error
         real(kind=rk) :: w(3*4 + 15)
         real(kind=rk) :: x(4) = [1, 2, 3, 4]
-        real(kind=rk) :: eps = 1.0e-10_rk
 
         call dcosqi(4, w)
         call dcosqf(4, x, w)
@@ -42,7 +47,6 @@ end subroutine test_classic_qct
 
     subroutine test_modernized_qct(error)
         type(error_type), allocatable, intent(out) :: error
-        real(kind=rk) :: eps = 1.0e-10_rk
         real(kind=rk) :: x(3) = [9, -9, 3]
 
         call check(error, sum(abs(qct(x, 2) - [-3.7279220613578570_rk, 21.727922061357859_rk])) < eps, &
@@ -59,7 +63,6 @@ end subroutine test_modernized_qct
 
     subroutine test_modernized_iqct(error)
         type(error_type), allocatable, intent(out) :: error
-        real(kind=rk) :: eps = 1.0e-10_rk
         real(kind=rk) :: x(4) = [1, 2, 3, 4]
 
         call check(error, sum(abs(iqct(qct(x))/(4.0_rk*4.0_rk) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, &