Add (i)dst/(i)qst interfaces

fpm.toml

@@ -72,6 +72,21 @@ name = "fftpack_iqct"
 source-dir = "test"
 main = "test_fftpack_iqct.f90"
 
+[[test]]
+name = "fftpack_dsinq"
+source-dir = "test"
+main = "test_fftpack_dsinq.f90"
+
+[[test]]
+name = "fftpack_qst"
+source-dir = "test"
+main = "test_fftpack_qst.f90"
+
+[[test]]
+name = "fftpack_iqst"
+source-dir = "test"
+main = "test_fftpack_iqst.f90"
+
 [[test]]
 name = "fftpack_dcost"
 source-dir = "test"
@@ -82,6 +97,16 @@ name = "fftpack_dct"
 source-dir = "test"
 main = "test_fftpack_dct.f90"
 
+[[test]]
+name = "fftpack_dsint"
+source-dir = "test"
+main = "test_fftpack_dsint.f90"
+
+[[test]]
+name = "fftpack_dst"
+source-dir = "test"
+main = "test_fftpack_dst.f90"
+
 # `fftpack` utility routines
 [[test]]
 name = "fftpack_fftshift"

src/Makefile

@@ -59,7 +59,10 @@ SRCF90 = \
 	fftpack_ifftshift.f90\
 	fftpack_qct.f90\
 	fftpack_iqct.f90\
-	fftpack_dct.f90
+	fftpack_qst.f90\
+	fftpack_iqst.f90\
+	fftpack_dct.f90\
+	fftpack_dst.f90
 
 OBJF := $(SRCF:.f=.o)
 OBJF90 := $(SRCF90:.f90=.o)
@@ -82,6 +85,9 @@ fftpack_rfft.o: fftpack.o
 fftpack_irfft.o: fftpack.o
 fftpack_qct.o: fftpack.o
 fftpack_iqct.o: fftpack.o
+fftpack_qst.o: fftpack.o
+fftpack_iqst.o: fftpack.o
 fftpack_dct.o: fftpack.o
+fftpack_dst.o: fftpack.o
 fftpack_fftshift.o: fftpack.o
 fftpack_ifftshift.o: fftpack.o

src/fftpack.f90

@@ -16,9 +16,15 @@ module fftpack
     public :: dcosqi, dcosqf, dcosqb
     public :: qct, iqct
 
+    public :: dsinqi, dsinqf, dsinqb
+    public :: qst, iqst
+
     public :: dcosti, dcost
     public :: dct, idct
 
+    public :: dsinti, dsint
+    public :: dst, idst
+
     interface
 
         !> Version: experimental
@@ -153,6 +159,38 @@ pure subroutine dcosqb(n, x, wsave)
             real(kind=dp), intent(in) :: wsave(*)
         end subroutine dcosqb
 
+        !> Version: experimental
+        !>
+        !> Initialize `dsinqf` and `dsinqb`.
+        !> ([Specification](../page/specs/fftpack.html#dsinqi))
+        pure subroutine dsinqi(n, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(out) :: wsave(*)
+        end subroutine dsinqi
+
+        !> Version: experimental
+        !>
+        !> Forward transform of quarter wave data.
+        !> ([Specification](../page/specs/fftpack.html#dsinqf))
+        pure subroutine dsinqf(n, x, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(inout) :: x(*)
+            real(kind=dp), intent(in) :: wsave(*)
+        end subroutine dsinqf
+
+        !> Version: experimental
+        !>
+        !> Unnormalized inverse of `dsinqf`.
+        !> ([Specification](../page/specs/fftpack.html#dsinqb))
+        pure subroutine dsinqb(n, x, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(inout) :: x(*)
+            real(kind=dp), intent(in) :: wsave(*)
+        end subroutine dsinqb
+
         !> Version: experimental
         !>
         !> Initialize `dcost`. ([Specification](../page/specs/fftpack.html#dcosti))
@@ -173,6 +211,26 @@ pure subroutine dcost(n, x, wsave)
             real(kind=dp), intent(in) :: wsave(*)
         end subroutine dcost
 
+        !> Version: experimental
+        !>
+        !> Initialize `dsint`. ([Specification](../page/specs/fftpack.html#dsinti))
+        pure subroutine dsinti(n, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(out) :: wsave(*)
+        end subroutine dsinti
+
+        !> Version: experimental
+        !>
+        !> Discrete fourier sine transform of an odd sequence.
+        !> ([Specification](../page/specs/fftpack.html#dsint))
+        pure subroutine dsint(n, x, wsave)
+            import dp
+            integer, intent(in) :: n
+            real(kind=dp), intent(inout) :: x(*)
+            real(kind=dp), intent(in) :: wsave(*)
+        end subroutine dsint
+
     end interface
 
     !> Version: experimental
@@ -225,7 +283,7 @@ end function irfft_dp
 
     !> Version: experimental
     !>
-    !> Forward transform of quarter wave data.
+    !> Forward cosine-transform of quarter wave data.
     !> ([Specifiction](../page/specs/fftpack.html#qct))
     interface qct
         pure module function qct_dp(x, n) result(result)
@@ -237,7 +295,7 @@ end function qct_dp
 
     !> Version: experimental
     !>
-    !> Backward transform of quarter wave data.
+    !> Backward cosine-transform of quarter wave data.
     !> ([Specifiction](../page/specs/fftpack.html#iqct))
     interface iqct
         pure module function iqct_dp(x, n) result(result)
@@ -247,6 +305,30 @@ pure module function iqct_dp(x, n) result(result)
         end function iqct_dp
     end interface iqct
 
+    !> Version: experimental
+    !>
+    !> Forward sine-transform of quarter wave data.
+    !> ([Specifiction](../page/specs/fftpack.html#qst))
+    interface qst
+        pure module function qst_dp(x, n) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            real(kind=dp), allocatable :: result(:)
+        end function qst_dp
+    end interface qst
+
+    !> Version: experimental
+    !>
+    !> Backward sine-transform of quarter wave data.
+    !> ([Specifiction](../page/specs/fftpack.html#iqst))
+    interface iqst
+        pure module function iqst_dp(x, n) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            real(kind=dp), allocatable :: result(:)
+        end function iqst_dp
+    end interface iqst
+
     !> Version: experimental
     !>
     !> Discrete fourier cosine (forward) transform of an even sequence.
@@ -267,6 +349,26 @@ end function dct_dp
         module procedure :: dct_dp
     end interface idct
 
+    !> Version: experimental
+    !>
+    !> Discrete fourier sine (forward) transform of an odd sequence.
+    !> ([Specification](../page/specs/fftpack.html#dst))
+    interface dst
+        pure module function dst_dp(x, n) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            real(kind=dp), allocatable :: result(:)
+        end function dst_dp
+    end interface dst
+
+    !> Version: experimental
+    !>
+    !> Discrete fourier sine (backward) transform of an even sequence.
+    !> ([Specification](../page/specs/fftpack.html#idst))
+    interface idst
+        module procedure :: dst_dp
+    end interface idst
+
     !> Version: experimental
     !>
     !> Shifts zero-frequency component to center of spectrum.

src/fftpack_dst.f90

@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_dst
+
+contains
+
+    !> Discrete fourier sine transform of an odd sequence.
+    pure module function dst_dp(x, n) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        real(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, i
+        real(kind=dp), 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))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = int(2.5_dp*lenseq + 15)
+        allocate (wsave(lensav))
+        call dsinti(lenseq, wsave)
+
+        !> Discrete fourier sine transformation
+        call dsint(lenseq, result, wsave)
+
+    end function dst_dp
+
+end submodule fftpack_dst

src/fftpack_iqct.f90

@@ -2,7 +2,7 @@
 
 contains
 
-    !> Backward transform of quarter wave data.
+    !> Backward cosine-transform of quarter wave data.
     pure module function iqct_dp(x, n) result(result)
         real(kind=dp), intent(in) :: x(:)
         integer, intent(in), optional :: n

src/fftpack_iqst.f90

@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_iqst
+
+contains
+
+    !> Backward sine-transform of quarter wave data.
+    pure module function iqst_dp(x, n) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        real(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, i
+        real(kind=dp), 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))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = 3*lenseq + 15
+        allocate (wsave(lensav))
+        call dsinqi(lenseq, wsave)
+
+        !> Backward transformation
+        call dsinqb(lenseq, result, wsave)
+
+    end function iqst_dp
+
+end submodule fftpack_iqst

src/fftpack_qct.f90

@@ -2,7 +2,7 @@
 
 contains
 
-    !> Forward transform of quarter wave data.
+    !> Forward cosine-transform of quarter wave data.
     pure module function qct_dp(x, n) result(result)
         real(kind=dp), intent(in) :: x(:)
         integer, intent(in), optional :: n

src/fftpack_qst.f90

@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_qst
+
+contains
+
+    !> Forward sine-transform of quarter wave data.
+    pure module function qst_dp(x, n) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        real(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, i
+        real(kind=dp), 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))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = 3*lenseq + 15
+        allocate (wsave(lensav))
+        call dsinqi(lenseq, wsave)
+
+        !> Forward transformation
+        call dsinqf(lenseq, result, wsave)
+
+    end function qst_dp
+
+end submodule fftpack_qst

test/Makefile

@@ -9,8 +9,13 @@ all: tstfft \
 	fftpack_dcosq \
 	fftpack_qct \
 	fftpack_iqct \
+	fftpack_dsinq \
+	fftpack_qst \
+	fftpack_iqst \
 	fftpack_dcost \
-	fftpack_dct
+	fftpack_dct \
+	fftpack_dsint \
+	fftpack_dst
 
 # Orginal test
 tstfft: tstfft.o
@@ -50,6 +55,18 @@ fftpack_iqct: test_fftpack_iqct.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
 	./fftpack_iqct.x
 
+fftpack_dsinq: test_fftpack_dsinq.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
+	./fftpack_dsinq.x
+
+fftpack_qst: test_fftpack_qst.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
+	./fftpack_qst.x
+
+fftpack_iqst: test_fftpack_iqst.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
+	./fftpack_iqst.x
+
 fftpack_dcost: test_fftpack_dcost.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
 	./fftpack_dcost.x
@@ -58,6 +75,14 @@ fftpack_dct: test_fftpack_dct.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
 	./fftpack_dct.x
 
+fftpack_dsint: test_fftpack_dsint.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
+	./fftpack_dsint.x
+
+fftpack_dst: test_fftpack_dst.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
+	./fftpack_dst.x
+
 # `fftpack` utility routines
 fftpack_fftshift: test_fftpack_fftshift.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]