Combine dct and qct interfaces

src/CMakeLists.txt

@@ -31,9 +31,7 @@ set(FFTPACK_SOURCES
     ${dir}/fftpack_fftshift.f90
     ${dir}/fftpack_ifft.f90
     ${dir}/fftpack_ifftshift.f90
-    ${dir}/fftpack_iqct.f90
     ${dir}/fftpack_irfft.f90
-    ${dir}/fftpack_qct.f90
     ${dir}/fftpack_rfft.f90
     ${dir}/fftpack_utils.f90
     ${dir}/passb.f90

src/Makefile

@@ -57,8 +57,6 @@ SRCF90 = \
 	fftpack_irfft.f90\
 	fftpack_fftshift.f90\
 	fftpack_ifftshift.f90\
-	fftpack_qct.f90\
-	fftpack_iqct.f90\
 	fftpack_dct.f90\
 	rk.f90\
 	fftpack_utils.f90
@@ -82,8 +80,6 @@ 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

src/fftpack.f90

@@ -15,8 +15,6 @@ module fftpack
     public :: dzffti, dzfftf, dzfftb
 
     public :: dcosqi, dcosqf, dcosqb
-    public :: qct, iqct
-
     public :: dcosti, dcost
     public :: dct, idct
 
@@ -246,46 +244,28 @@ end function irfft_rk
 
     !> Version: experimental
     !>
-    !> Forward transform of quarter wave data.
-    !> ([Specifiction](../page/specs/fftpack.html#qct))
-    interface qct
-        pure module function qct_rk(x, n) result(result)
-            real(kind=rk), intent(in) :: x(:)
-            integer, intent(in), optional :: n
-            real(kind=rk), allocatable :: result(:)
-        end function qct_rk
-    end interface qct
-
-    !> Version: experimental
-    !>
-    !> Backward transform of quarter wave data.
-    !> ([Specifiction](../page/specs/fftpack.html#iqct))
-    interface iqct
-        pure module function iqct_rk(x, n) result(result)
-            real(kind=rk), intent(in) :: x(:)
-            integer, intent(in), optional :: n
-            real(kind=rk), allocatable :: result(:)
-        end function iqct_rk
-    end interface iqct
-
-    !> Version: experimental
-    !>
-    !> Discrete fourier cosine (forward) transform of an even sequence.
+    !> Dsicrete cosine transforms.
     !> ([Specification](../page/specs/fftpack.html#dct))
     interface dct
-        pure module function dct_rk(x, n) result(result)
+        pure module function dct_rk(x, n, t) result(result)
             real(kind=rk), intent(in) :: x(:)
             integer, intent(in), optional :: n
+            integer, intent(in), optional :: t
             real(kind=rk), allocatable :: result(:)
         end function dct_rk
     end interface dct
 
     !> Version: experimental
     !>
-    !> Discrete fourier cosine (backward) transform of an even sequence.
+    !> Inverse discrete cosine transforms.
     !> ([Specification](../page/specs/fftpack.html#idct))
     interface idct
-        module procedure :: dct_rk
+    pure module function idct_rk(x, n, t) result(result)
+        real(kind=rk), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        integer, intent(in), optional :: t
+        real(kind=rk), allocatable :: result(:)
+        end function idct_rk
     end interface idct
 
     !> Version: experimental

src/fftpack_dct.f90

@@ -2,10 +2,11 @@
 
 contains
 
-    !> Discrete fourier cosine transform of an even sequence.
-    pure module function dct_rk(x, n) result(result)
+    !> Discrete cosine transforms of types 1, 2, 3.
+    pure module function dct_rk(x, n, t) result(result)
         real(kind=rk), intent(in) :: x(:)
         integer, intent(in), optional :: n
+        integer, intent(in), optional :: t
         real(kind=rk), allocatable :: result(:)
 
         integer :: lenseq, lensav, i
@@ -23,14 +24,80 @@ pure module function dct_rk(x, n) result(result)
             result = x
         end if
 
-        !> Initialize FFT
-        lensav = 3*lenseq + 15
-        allocate (wsave(lensav))
-        call dcosti(lenseq, wsave)
-
-        !> Discrete fourier cosine transformation
-        call dcost(lenseq, result, wsave)
+        ! Default to DCT-2
+        if (.not.present(t)) then
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosqi(lenseq, wsave)
+            call dcosqb(lenseq, result, wsave)
+            return
+        end if
 
+        if (t == 1) then  ! DCT-1
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosti(lenseq, wsave)
+            call dcost(lenseq, result, wsave)
+        else if (t == 2) then  ! DCT-2
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosqi(lenseq, wsave)
+            call dcosqb(lenseq, result, wsave)
+        else if (t == 3) then  ! DCT-3
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosqi(lenseq, wsave)
+            call dcosqf(lenseq, result, wsave)
+        end if
     end function dct_rk
 
+    !> Inverse discrete cosine transforms of types 1, 2, 3.
+    pure module function idct_rk(x, n, t) result(result)
+        real(kind=rk), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        integer, intent(in), optional :: t
+        real(kind=rk), allocatable :: result(:)
+
+        integer :: lenseq, lensav, i
+        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_rk, i=1, lenseq - size(x))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        ! Default to t=2; inverse DCT-2 is DCT-3
+        if (.not.present(t)) then
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosqi(lenseq, wsave)
+            call dcosqf(lenseq, result, wsave)
+            return
+        end if
+
+        if (t == 1) then  ! inverse DCT-1 is DCT-1
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosti(lenseq, wsave)
+            call dcost(lenseq, result, wsave)
+        else if (t == 2) then  ! inverse DCT-2 is DCT-3
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosqi(lenseq, wsave)
+            call dcosqf(lenseq, result, wsave)
+        else if (t == 3) then  ! inverse DCT-3 is DCT-2
+            lensav = 3*lenseq + 15
+            allocate (wsave(lensav))
+            call dcosqi(lenseq, wsave)
+            call dcosqb(lenseq, result, wsave)
+        end if
+    end function idct_rk
+
 end submodule fftpack_dct

test/CMakeLists.txt

@@ -14,7 +14,6 @@ endmacro()
 set(FFTPACK_TEST_SOURCES
     test_fftpack_dct.f90
     test_fftpack_fft.f90
-    test_fftpack_qct.f90
     test_fftpack_rfft.f90
     test_fftpack_utils.f90
     test_fftpack.f90

test/Makefile

@@ -3,12 +3,11 @@ FETCH = curl -L
 SRC = \
 	test_fftpack_fft.f90 \
 	test_fftpack_rfft.f90 \
-	test_fftpack_qct.f90 \
 	test_fftpack_dct.f90 \
 	test_fftpack_utils.f90 \
 	test_fftpack.f90 \
 	testdrive.F90
-	
+
 OBJ = $(SRC:.f90=.o)
 OBJ := $(OBJ:.F90=.o)
 
@@ -24,10 +23,10 @@ tstfft: tstfft.f
 test_fftpack: $(OBJ)
 	$(FC) $(FFLAGS) $(OBJ) -L../src -l$(LIB) -I../src -o $@.x
 	./test_fftpack.x
-	
+
 testdrive.F90:
 	$(FETCH) https://github.com/fortran-lang/test-drive/raw/v0.4.0/src/testdrive.F90 > $@
-	
+
 %.o: %.F90
 	$(FC) $(FFLAGS) -c $<
 
@@ -36,14 +35,12 @@ testdrive.F90:
 
 test_fftpack.o: test_fftpack_fft.o \
 	test_fftpack_rfft.o \
-	test_fftpack_qct.o \
 	test_fftpack_dct.o \
 	test_fftpack_utils.o \
 	testdrive.o
 
 test_fftpack_fft.o: testdrive.o
 test_fftpack_rfft.o: testdrive.o
-test_fftpack_qct.o: testdrive.o
 test_fftpack_dct.o: testdrive.o
 test_fftpack_utils.o: testdrive.o
 

test/test_fftpack.f90

@@ -3,7 +3,6 @@ program test_fftpack
     use testdrive, only: run_testsuite, new_testsuite, testsuite_type
     use test_fftpack_fft, only: collect_fft
     use test_fftpack_rfft, only: collect_rfft
-    use test_fftpack_qct, only: collect_qct
     use test_fftpack_dct, only: collect_dct
     use test_fftpack_utils, only: collect_utils
     implicit none
@@ -16,7 +15,6 @@ program test_fftpack
     testsuites = [ &
                  new_testsuite("fft", collect_fft), &
                  new_testsuite("rfft", collect_rfft), &
-                 new_testsuite("qct", collect_qct), &
                  new_testsuite("dct", collect_dct), &
                  new_testsuite("utils", collect_utils) &
                  ]