Merge pull request #51 from fortran-lang/complex-fft-example

Add complex-fft example; fix comment & rename real-fft example

.github/workflows/fpm.yml

@@ -8,13 +8,13 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-latest, macos-11]
-        gcc_v: [10] # Version of GFortran we want to use.
+        os: [ubuntu-latest, macos-12]
+        gcc_v: [11] # Version of GFortran we want to use.
         include:
         - os: ubuntu-latest
           os-arch: linux-x86_64
 
-        - os: macos-11
+        - os: macos-12
           os-arch: macos-x86_64
 
     env:
@@ -25,24 +25,35 @@ jobs:
     - name: Checkout code
       uses: actions/checkout@v1
 
-    - name: Install GFortran macOS
-      if: contains(matrix.os, 'macos')
-      run: |
-        ln -s /usr/local/bin/gfortran-${GCC_V} /usr/local/bin/gfortran
-        which gfortran-${GCC_V}
-        which gfortran
-
     - name: Install GFortran Linux
       if: contains(matrix.os, 'ubuntu')
       run: |
         sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-${GCC_V} 100 \
         --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${GCC_V} \
         --slave /usr/bin/gcov gcov /usr/bin/gcov-${GCC_V}
-    
+
+    # Backport gfortran shared libraries to version 9 folder. This is necessary because the macOS release of fpm 
+    # 0.10.0 used for bootstrapping has these paths hardcoded in the executable.
+    # See https://github.com/fortran-lang/fpm/pull/1061
+    - name: MacOS patch libgfortran
+      if: contains(matrix.os, 'macos')
+      run: |
+        ln -s /usr/local/bin/gfortran-${GCC_V} /usr/local/bin/gfortran
+        which gfortran-${GCC_V}
+        which gfortran
+        mkdir /usr/local/opt/gcc@10
+        mkdir /usr/local/opt/gcc@10/lib
+        mkdir /usr/local/opt/gcc@10/lib/gcc
+        mkdir /usr/local/opt/gcc@10/lib/gcc/10
+        mkdir /usr/local/lib/gcc/10
+        ln -fs /usr/local/opt/gcc@${GCC_V}/lib/gcc/${GCC_V}/libquadmath.0.dylib /usr/local/opt/gcc@10/lib/gcc/10/libquadmath.0.dylib 
+        ln -fs /usr/local/opt/gcc@${GCC_V}/lib/gcc/${GCC_V}/libgfortran.5.dylib /usr/local/opt/gcc@10/lib/gcc/10/libgfortran.5.dylib
+        ln -fs /usr/local/lib/gcc/${GCC_V}/libgcc_s.1.dylib /usr/local/lib/gcc/10/libgcc_s.1.dylib
+
     - name: Install fpm
-      uses: fortran-lang/setup-fpm@v3
+      uses: fortran-lang/setup-fpm@v5
       with:
-        fpm-version: 'v0.8.2'
+        fpm-version: 'v0.9.0'
 
     - name: Build fftpack
       run: |

example/CMakeLists.txt

@@ -9,3 +9,9 @@ target_link_libraries(bench3 fftpack)
 
 add_executable(rfft_example)
 target_link_libraries(rfft_example fftpack)
+
+add_executable(complex_transforms)
+target_link_libraries(complex_transforms fftpack)
+
+add_executable(real_transforms)
+target_link_libraries(real_transforms fftpack)

example/complex_transforms.f90

@@ -0,0 +1,46 @@
+program complex_transforms
+  !! This program invokes fftpack's fft function to compute the forward transform of a complex function
+  !! and the inverse transform of the result.  An assertion verifies the expected result of the forward 
+  !! transform according to the element layout described at [1].  A second assertion checks that the
+  !! inverse transform recovers the original function.
+  !!
+  !! [1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.fft.html#scipy.fftpack.fft
+  use fftpack, only: fft, ifft
+  implicit none
+  integer j, k
+  integer, parameter :: N = 8
+  double precision, parameter :: two_pi = 2.D0*acos(-1.D0), tolerance = 1.0D-06, f_avg = 3.D0, zero=0.D0
+  double precision, parameter :: x(0:N-1) = [(two_pi*dble(j)/dble(N), j=0,N-1)]
+  integer, parameter :: rk = kind(two_pi)
+  complex(rk), parameter :: f(0:N-1) = f_avg + cos(x)
+    !! sample f(x) = 3 + cos(x) uniformly on [0,2*pi)
+    !!             = 3 + (exp(i*x) - exp(-i*x))/2  
+    !! which yields the Fourier coefficients
+    !!                {  3, k = 0
+    !!       f_hat =  {  1/2, k = 1
+    !!                {  1/2, k = -1
+    !!                {  0, otherwise
+  complex(rk), dimension(0:N-1) :: f_round_trip, fft_f
+  character(len=*), parameter :: real_format = "(a,*(g10.4,:,1x))" !! space-separated values
+  character(len=*), parameter :: complex_format= "(a,*('(',g10.4,',',g10.4,')',:,1x)))" !! space-separated complex values
+
+  call assert(mod(N,2)==0, "the algorithm below requires even N")
+
+  fft_f(:) = fft(f)/dble(N)
+  f_round_trip(:) = ifft(fft_f)
+
+  print complex_format, "f = ", f
+  print complex_format, "fft_f = ", fft_f
+  print complex_format, "f_round_trip = ",f_round_trip
+
+  !call assert(all(abs(f_round_trip - f) < tolerance), "inverse of forward FFT must yield the original function")
+
+contains
+
+  pure subroutine assert(assertion, description)
+    logical, intent(in) :: assertion
+    character(len=*), intent(in) :: description
+    if (.not. assertion) error stop description
+  end subroutine
+
+end program

example/real-forward-transform.f90 → example/real_transforms.f90

@@ -1,8 +1,8 @@
-program forward_transform_of_real_function
+program real_transforms
   !! This program invokes fftpack's rrft function to compute the forward transform of a real function
   !! and constructs the resulting complex Fourier coefficients by (re)organizing and normalizing the
   !! rfft result according to array element layout described at [1].  The program also demonstrates
-  !! the inverse transform of the raw rrft result to recover the original function.
+  !! the inverse transform of the normalized rrft result to recover the original function.
   !!
   !! [1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.rfft.html#scipy.fftpack.rfft
   use fftpack, only: rfft, irfft