When using r2c and c2r FFTW in Fortran, are the fo

Question

Blow is a main file

PROGRAM SPHEROID

USE nrtype
USE SUB_INFO
INCLUDE "/usr/local/include/fftw3.f"

INTEGER(I8B) :: plan_forward, plan_backward
INTEGER(I4B) :: i, t, int_N

REAL(DP) :: cth_i, sth_i, real_i, perturbation
REAL(DP) :: PolarEffect, dummy, x1, x2, x3

REAL(DP), DIMENSION(4096) :: dummy1, dummy2, gam, th, ph
REAL(DP), DIMENSION(4096) ::  k1, k2, k3, k4, l1, l2, l3, l4, f_in

COMPLEX(DPC), DIMENSION(2049) :: output1, output2, f_out

CHARACTER(1024) :: baseOutputFilename
CHARACTER(1024) :: outputFile, format_string

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

int_N = 4096

! File Open Section

format_string = '(I5.5)'

! Write the coodinates at t = 0

do i = 1, N
    real_i = real(i)
    gam(i) = 2d0*pi/real_N
    perturbation = 0.01d0*dsin(2d0*pi*real_i/real_N)
    ph(i) = 2d0*pi*real_i/real_N + perturbation
    th(i) = pi/3d0 + perturbation
end do

! Initialization Section for FFTW PLANS

call dfftw_plan_dft_r2c_1d(plan_forward, int_N, f_in, f_out, FFTW_ESTIMATE)
call dfftw_plan_dft_c2r_1d(plan_backward, int_N, f_out, f_in, FFTW_ESTIMATE)

! Runge-Kutta 4th Order Method Section

do t = 1, Iter_N

    call integration(th, ph, gam, k1, l1)

    do i = 1, N
        dummy1(i) = th(i) + 0.5d0*dt*k1(i)
    end do
    do i = 1, N
        dummy2(i) = ph(i) + 0.5d0*dt*l1(i)
    end do

    call integration(dummy1, dummy2, gam, k2, l2)

    do i = 1, N
        dummy1(i) = th(i) + 0.5d0*dt*k2(i)
    end do
    do i = 1, N
        dummy2(i) = ph(i) + 0.5d0*dt*l2(i)
    end do

    call integration(dummy1, dummy2, gam, k3, l3)

    do i = 1, N
        dummy1(i) = th(i) + dt*k3(i)
    end do
    do i = 1, N
        dummy2(i) = ph(i) + dt*l3(i)
    end do

    call integration(dummy1, dummy2, gam, k4, l4)

    do i = 1, N

        cth_i = dcos(th(i))
        sth_i = dsin(th(i))
        PolarEffect = (nv-sv)*dsqrt(1d0+a*sth_i**2) + (nv+sv)*cth_i
        PolarEffect = PolarEffect/(sth_i**2)
        th(i) = th(i) + dt*(k1(i) + 2d0*k2(i) + 2d0*k3(i) + k4(i))/6d0
        ph(i) = ph(i) + dt*(l1(i) + 2d0*l2(i) + 2d0*l3(i) + l4(i))/6d0
        ph(i) = ph(i) + dt*0.25d0*PolarEffect/pi

    end do

!! Fourier Filtering Section

    call dfftw_execute_dft_r2c(plan_forward, th, output1)

    do i = 1, N/2+1
       dummy = abs(output1(i))
        if (dummy.lt.threshhold) then
           output1(i) = dcmplx(0.0d0)
        end if
    end do

    call dfftw_execute_dft_c2r(plan_backward, output1, th)

    do i = 1, N
       th(i) = th(i)/real_N
    end do

    call dfftw_execute_dft_r2c(plan_forward, ph, output2)

    do i = 1, N/2+1
       dummy = abs(output2(i))
        if (dummy.lt.threshhold) then
           output2(i) = dcmplx(0.0d0)
        end if
    end do

    call dfftw_execute_dft_c2r(plan_backward, output2, ph)

    do i = 1, N
       ph(i) = ph(i)/real_N
    end do

!! Data Writing Section

    write(baseOutputFilename, format_string) t
    outputFile = "xyz" // baseOutputFilename
    open(unit=7, file=outputFile)
    outputFile = "Fsptrm" // baseOutputFilename
    open(unit=8, file=outputFile)

    do i = 1, N
        x1 = dsin(th(i))*dcos(ph(i))
        x2 = dsin(th(i))*dsin(ph(i))
        x3 = dsqrt(1d0+a)*dcos(th(i))
        write(7,*) x1, x2, x3
    end do

    do i = 1, N/2+1
        write(8,*) abs(output1(i)), abs(output2(i))
    end do

    close(7)
    close(8)

    do i = 1, N/2+1
        output1(i) = dcmplx(0.0d0)
    end do

    do i = 1, N/2+1
        output2(i) = dcmplx(0.0d0)
    end do

end do

! Destroying Process for FFTW PLANS

call dfftw_destroy_plan(plan_forward)
call dfftw_destroy_plan(plan_backward)

END PROGRAM

Below is a subroutine file for integration

! We implemented Shelly's spectrally accurate convergence method

SUBROUTINE integration(in1,in2,in3,out1,out2)

USE nrtype
USE SUB_INFO

INTEGER(I4B) :: i, j
REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2
REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2
REAL(DP) :: ui, uj, part1, part2, gj, cph, sph
REAL(DP) :: denom, numer, temp

do i = 1, N
    out1(i) = 0d0
end do
do i = 1, N
    out2(i) = 0d0
end do

do i = 1, N

    th_i = in1(i)
    ph_i = in2(i)
    ui = dcos(th_i)
    part1 = dsqrt(1d0+a)/(dsqrt(-a)*ui+dsqrt(1d0+a-a*ui*ui))
    part1 = part1**(dsqrt(-a))
    part2 = (dsqrt(1d0+a-a*ui*ui)+ui)/(dsqrt(1d0+a-a*ui*ui)-ui)
    part2 = dsqrt(part2)

    gi = dsqrt(1d0-ui*ui)*part1*part2

    do j = 1, N

        if (mod(i+j,2).eq.1) then

            th_j = in1(j)
            ph_j = in2(j)
            gam_j = in3(j)

            uj = dcos(th_j)
            part1 = dsqrt(1d0+a)/(dsqrt(-a)*uj+dsqrt(1d0+a-a*uj*uj))
            part1 = part1**(dsqrt(-a))
            part2 = (dsqrt(1d0+a-a*uj*uj)+uj)/(dsqrt(1d0+a-a*uj*uj)-uj)
            part2 = dsqrt(part2)

            gj = dsqrt(1d0-ui*ui)*part1*part2

            cph = dcos(ph_i-ph_j)
            sph = dsin(ph_i-ph_j)

            numer = dsqrt(1d0-uj*uj)*sph
            denom = (gj/gi*(1d0-ui*ui) + gi/gj*(1d0-uj*uj))*0.5d0
            denom = denom - dsqrt((1d0-ui*ui)*(1d0-uj*uj))*cph
            denom = denom + krasny_delta
            v1 = -0.25d0*gam_j*numer/denom/pi

            temp = dsqrt(1d0+(1d0-ui*ui)*a)
            numer = -(gj/gi)*(temp+ui)
            numer = numer + (gi/gj)*((1d0-uj*uj)/(1d0-ui*ui))*(temp-ui)
            numer = numer + 2d0*ui*dsqrt((1d0-uj*uj)/(1d0-ui*ui))*cph
            numer = 0.5d0*numer
            v2 = -0.25d0*gam_j*numer/denom/pi

            out1(i) = out1(i) + 2d0*v1
            out2(i) = out2(i) + 2d0*v2

        end if

    end do

end do

END

Below is a module file

module nrtype
Implicit none
!integer
INTEGER, PARAMETER :: I8B = SELECTED_INT_KIND(20)
INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
!real
INTEGER, PARAMETER :: SP = KIND(1.0)
INTEGER, PARAMETER :: DP = KIND(1.0D0)
!complex
INTEGER, PARAMETER :: SPC = KIND((1.0,1.0))
INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0))
!defualt logical
INTEGER, PARAMETER :: LGT = KIND(.true.)
!mathematical constants
REAL(DP), PARAMETER :: pi = 3.141592653589793238462643383279502884197_dp
!derived data type s for sparse matrices,single and double precision
!User-Defined Constants
INTEGER(I4B), PARAMETER :: N = 4096, Iter_N = 20000
REAL(DP), PARAMETER :: real_N = 4096d0
REAL(DP), PARAMETER :: a = -0.1d0, dt = 0.001d0, krasny_delta = 0.01d0
REAL(DP), PARAMETER :: nv = 0d0, sv = 0d0, threshhold = 0.00000000001d0
!N : The Number of Point Vortices, Iter_N * dt = Total time, dt : Time Step
!krasny_delta : Smoothing Parameter introduced by R.Krasny
!nv : Northern Vortex Strength, sv : Southern Vortex Strength
!a : The Eccentricity in the direction of z , threshhold : Filtering Threshhold
end module nrtype

Below is a subroutine info file

MODULE SUB_INFO

INTERFACE
  SUBROUTINE integration(in1,in2,in3,out1,out2)
  USE nrtype
  INTEGER(I4B) :: i, j
  REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2
  REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2
  REAL(DP) :: ui, uj, part1, part2, gj, cph, sph
  REAL(DP) :: denom, numer, temp
  END SUBROUTINE
END INTERFACE

END MODULE

I compiled them using the below command

gfortran -o p0 -fbounds-check nrtype.f90 spheroid_sub_info.f90 spheroid_sub_integration.f90 spheroid_main.f90 -lfftw3 -lm -Wall -pedantic -pg

nohup ./p0 &

Note that 2049 = 4096 / 2 + 1

When making plan_backward, isn't it correct that we use 2049 instead of 4096 since the dimension of output is 2049?

But when I do that, it blows up. (Blowing up means NAN error)

If I use 4096 in making plan_backward, Everything is fine except that some Fourier coefficients are abnormally big which should not happen.

Please help me use FFTW in Fortran correctly. This issue has discouraged me for a long time.

Answer 1

One issue may be that dfftw_execute_dft_c2r can destroy the content of the input array, as described in this page. The key excerpt is

FFTW_PRESERVE_INPUT specifies that an out-of-place transform must not change its input array. This is ordinarily the default, except for c2r and hc2r (i.e. complex-to-real) transforms for which FFTW_DESTROY_INPUTis the default...

We can verify this, for example, by modifying the sample code by @VladimirF such that it saves data_out to data_save right after the first FFT(r2c) call, and then calculating their difference after the second FFT (c2r) call. So, in the case of OP's code, it seems safer to save output1 and output2 to different arrays before entering the second FFT (c2r).

Answer 2

First, although you claim your example is minimal, it is still pretty large, I have no time to study it.

But I updated my gist code https://gist.github.com/LadaF/73eb430682ef527eea9972ceb96116c5 to show also the backward transform and to answer the title question about the transform dimensions.

The logical size of the transform is the size of the real array (Real-data DFT Array Format) but the complex part is smaller due to inherent symmetries.

But when you make first r2c transform from real array of size n to complex array of size n/2+1. and then an opposite transform back, the real array should be again of size n.

This is my minimal example from the gist:

module FFTW3
  use, intrinsic :: iso_c_binding
  include "fftw3.f03"
end module

    use FFTW3
    implicit none

    integer, parameter :: n = 100

    real(c_double), allocatable :: data_in(:)
    complex(c_double_complex), allocatable :: data_out(:)
    type(c_ptr) :: planf, planb

    allocate(data_in(n))
    allocate(data_out(n/2+1))

    call random_number(data_in)

    planf = fftw_plan_dft_r2c_1d(size(data_in), data_in, data_out, FFTW_ESTIMATE+FFTW_UNALIGNED)
    planb = fftw_plan_dft_c2r_1d(size(data_in), data_out, data_in, FFTW_ESTIMATE+FFTW_UNALIGNED)

    print *, "real input:", real(data_in)

    call fftw_execute_dft_r2c(planf, data_in, data_out)

    print *, "result real part:", real(data_out)

    print *, "result imaginary part:", aimag(data_out)

    call fftw_execute_dft_c2r(planb, data_out, data_in)

    print *, "real output:", real(data_in)/n

    call fftw_destroy_plan(planf)
    call fftw_destroy_plan(planb)
end

Note that I am using the modern Fortran interface. I don't like using the old one.