add (i)fft2 interfaces.

add `stdlib/forlab` dependencies for example/test checks.

Author: zoziha

doc/specs/fftpack.md

@@ -281,6 +281,105 @@ program demo_ifft
 end program demo_ifft
 ```
 
+### `fft2`
+
+#### Description
+
+Computes the rank-2 forward complex discrete fourier transform (the fourier analysis).   
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Syntax
+
+`result = [[fftpack(module):fft2(interface)]](x)`
+
+#### Argument
+
+`x`: Shall be a `complex` and rank-2 array.
+This argument is `intent(in)`. 
+
+#### Return value
+
+Returns a `complex` and rank-2 array, the Discrete Fourier Transform (DFT) of `x`.
+
+#### Notes
+
+Within numerical accuracy, `x == ifft2(fft2(x))/(size(x, 1)*size(x, 2))`.
+
+#### Example
+
+```fortran
+program demo_fftpack_fft2
+
+    use fftpack, only: fft2
+    use forlab_io, only: disp
+
+    complex(kind=8), allocatable :: x(:, :)
+    x = reshape(spread([1, 2, 3, 4], 1, 4), [4, 4])
+    call disp(fft2(x), "fft2(x): ")
+    
+    !! (40.00,0.000)            (-8.000,8.000)            (-8.000,0.000)           (-8.000,-8.000)
+    !! (0.000,0.000)             (0.000,0.000)             (0.000,0.000)             (0.000,0.000)
+    !! (0.000,0.000)             (0.000,0.000)             (0.000,0.000)             (0.000,0.000)
+    !! (0.000,0.000)             (0.000,0.000)             (0.000,0.000)             (0.000,0.000)
+
+end program demo_fftpack_fft2
+
+```
+
+### `ifft2`
+
+#### Description
+
+Unnormalized inverse of `fft2`.
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Syntax
+
+`result = [[fftpack(module):ifft(interface)]](x)`
+
+#### Argument
+
+`x`: Shall be a `complex` and rank-2 array.
+This argument is `intent(in)`. 
+
+#### Return value
+
+Returns a `complex` and rank-2 array, the unnormalized inverse Discrete Fourier Transform (DFT) of `x`.
+
+#### Example
+
+```fortran
+program demo_fftpack_ifft2
+
+    use fftpack, only: fft2, ifft2
+    use forlab_io, only: disp
+
+    complex(kind=8), allocatable :: x(:, :)
+    x = reshape(spread([1, 2, 3, 4], 1, 4), [4, 4])
+    call disp(ifft2(fft2(x))/(size(x, 1)*size(x, 2)), "ifft2(fft2(x))/(size(x, 1)*size(x, 2)): ")
+
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+
+end program demo_fftpack_ifft2
+```
+
 ## Fourier transform of double real periodic sequences
 ### `dffti`
 

example/demo_fftpack_fft2.f90

@@ -0,0 +1,15 @@
+program demo_fftpack_fft2
+
+    use fftpack, only: fft2
+    use forlab_io, only: disp
+
+    complex(kind=8), allocatable :: x(:, :)
+    x = reshape(spread([1, 2, 3, 4], 1, 4), [4, 4])
+    call disp(fft2(x), "fft2(x): ")
+    
+    !! (40.00,0.000)            (-8.000,8.000)            (-8.000,0.000)           (-8.000,-8.000)
+    !! (0.000,0.000)             (0.000,0.000)             (0.000,0.000)             (0.000,0.000)
+    !! (0.000,0.000)             (0.000,0.000)             (0.000,0.000)             (0.000,0.000)
+    !! (0.000,0.000)             (0.000,0.000)             (0.000,0.000)             (0.000,0.000)
+
+end program demo_fftpack_fft2

example/demo_fftpack_ifft2.f90

@@ -0,0 +1,15 @@
+program demo_fftpack_ifft2
+
+    use fftpack, only: fft2, ifft2
+    use forlab_io, only: disp
+
+    complex(kind=8), allocatable :: x(:, :)
+    x = reshape(spread([1, 2, 3, 4], 1, 4), [4, 4])
+    call disp(ifft2(fft2(x))/(size(x, 1)*size(x, 2)), "ifft2(fft2(x))/(size(x, 1)*size(x, 2)): ")
+
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+    !! (1.000,0.000)             (2.000,0.000)             (3.000,0.000)             (4.000,0.000)
+
+end program demo_fftpack_ifft2

fpm.toml

@@ -15,6 +15,21 @@ auto-executables = false
 auto-tests = false
 auto-examples = false
 
+# examples
+[[example]]
+name = "fftpack_fft2"
+source-dir = "example"
+main = "demo_fftpack_fft2.f90"
+[[example]]
+name = "fftpack_ifft2"
+source-dir = "example"
+main = "demo_fftpack_ifft2.f90"
+
+[dev-dependencies]
+stdlib = { git="https://github.com/LKedward/stdlib-fpm.git" }
+forlab = { git="https://github.com/fortran-fans/forlab.git" }
+
+
 # Original test
 [[test]]
 name = "tstfft"
@@ -32,6 +47,11 @@ name = "fftpack_fft"
 source-dir = "test"
 main = "test_fftpack_fft.f90"
 
+[[test]]
+name = "fftpack_fft2"
+source-dir = "test"
+main = "test_fftpack_fft2.f90"
+
 [[test]]
 name = "fftpack_ifft"
 source-dir = "test"

src/fftpack.f90

@@ -5,7 +5,7 @@ module fftpack
     private
 
     public :: zffti, zfftf, zfftb
-    public :: fft, ifft
+    public :: fft, ifft, fft2, ifft2
     public :: fftshift, ifftshift
 
     public :: dffti, dfftf, dfftb
@@ -176,6 +176,28 @@ pure module function ifft_dp(x, n) result(result)
         end function ifft_dp
     end interface ifft
 
+    !> Version: experimental
+    !>
+    !> Forward transform of a rank-2 and double complex periodic sequence.
+    !> ([Specifiction](../page/specs/fftpack.html#fft2))
+    interface fft2
+        pure module function fft2_dp(x) result(result)
+            complex(kind=dp), intent(in) :: x(:, :)
+            complex(kind=dp), allocatable :: result(:, :)
+        end function fft2_dp
+    end interface fft2
+
+    !> Version: experimental
+    !>
+    !> Backward transform of a rank-2 and double complex periodic sequence.
+    !> ([Specifiction](../page/specs/fftpack.html#ifft2))
+    interface ifft2
+        pure module function ifft2_dp(x) result(result)
+            complex(kind=dp), intent(in) :: x(:, :)
+            complex(kind=dp), allocatable :: result(:, :)
+        end function ifft2_dp
+    end interface ifft2
+
     !> Version: experimental
     !>
     !> Forward transform of a double real periodic sequence.

src/fftpack_fft2.f90

@@ -0,0 +1,51 @@
+submodule(fftpack) fftpack_fft2
+
+contains
+
+    !> Forward transform of a rank-2 double complex periodic sequence.
+    !> `fft2` can achieve parallel acceleration through `OpenMP`.
+    pure module function fft2_dp(x) result(result)
+        complex(kind=dp), intent(in) :: x(:, :)
+        complex(kind=dp), allocatable :: result(:, :)
+
+        integer :: lendim(rank(x))
+        integer :: lenseq, lensav, i
+        real(kind=dp), allocatable :: wsave(:)
+
+        lendim = shape(x)
+        result = x
+
+        !! FFT is performed on each dimension of the array in turn, 
+        !!  and the scale of the array remains the same after each FFT.
+
+        !* First diemnsion:
+
+        !> Initialize FFT
+        lenseq = lendim(1)
+        lensav = 4*lenseq + 15
+        allocate(wsave(lensav))
+        call zffti(lenseq, wsave)
+
+        !> Forward transformation
+        do i = 1, lendim(2)
+            call zfftf(lenseq, result(:, i), wsave)
+        end do
+
+        deallocate(wsave)
+
+        !* Second dimension:
+
+        !> Initialize FFT
+        lenseq = lendim(2)
+        lensav = 4*lenseq + 15
+        allocate(wsave(lensav))
+        call zffti(lenseq, wsave)
+
+        !> Forward transformation
+        do i = 1, lendim(1)
+            call zfftf(lenseq, result(i, :), wsave)
+        end do
+
+    end function fft2_dp
+
+end submodule fftpack_fft2

src/fftpack_ifft2.f90

@@ -0,0 +1,51 @@
+submodule(fftpack) fftpack_ifft2
+
+contains
+
+    !> Backward transform of a rank-2 double complex periodic sequence.
+    !> `ifft2` can achieve parallel acceleration through `OpenMP`.
+    pure module function ifft2_dp(x) result(result)
+        complex(kind=dp), intent(in) :: x(:, :)
+        complex(kind=dp), allocatable :: result(:, :)
+
+        integer :: lendim(rank(x))
+        integer :: lenseq, lensav, i
+        real(kind=dp), allocatable :: wsave(:)
+
+        lendim = shape(x)
+        result = x
+        
+        !! FFT is performed on each dimension of the array in turn, 
+        !!  and the scale of the array remains the same after each FFT.
+
+        !* First diemnsion:
+
+        !> Initialize FFT
+        lenseq = lendim(1)
+        lensav = 4*lenseq + 15
+        allocate(wsave(lensav))
+        call zffti(lenseq, wsave)
+
+        !> Backward transformation
+        do i = 1, lendim(2)
+            call zfftb(lenseq, result(:, i), wsave)
+        end do
+
+        deallocate(wsave)
+
+        !* Second dimension:
+
+        !> Initialize FFT
+        lenseq = lendim(2)
+        lensav = 4*lenseq + 15
+        allocate(wsave(lensav))
+        call zffti(lenseq, wsave)
+
+        !> Backward transformation
+        do i = 1, lendim(1)
+            call zfftb(lenseq, result(i, :), wsave)
+        end do
+
+    end function ifft2_dp
+
+end submodule fftpack_ifft2

test/test_fftpack_fft2.f90

@@ -0,0 +1,12 @@
+program test_fftpack_fft2
+    use fftpack, only: fft2, ifft2
+    use forlab_math, only: is_close
+    use stdlib_error, only: check
+
+    complex(kind=8), allocatable :: x(:, :)
+    x = reshape(spread([1, 2, 3, 4], 1, 4), [4, 4])
+    
+    call check(all(is_close(ifft2(fft2(x))/16.0_8, x)), msg="`ifft2(fft2(x))` failed.")
+    print *, "All tests in `test_fftpack_fft2` passed."
+
+end program test_fftpack_fft2