CODEWITHSUNDEEP
fortranfftfortran90subroutinefftw
Ziegl
Ziegl

Reputation: 147

Passing an FFTW plan to a Fortran subroutine

I have the following, seemingly simple issue, but I haven't been able to find any information on it:

In fortran90, is it possible to pass an FFTW plan to a subroutine, so as to allow reusing a plan in different "places" of the code? The naive approach of simply passing it as an INTEGER*8 variable yields mayhem, suggesting the code needs more than the plain integer plan to execute the Fourier Transform.

Upvotes: 2

Views: 1403

Answers (2)

max
max

Reputation: 564

To this end, I usually write a Fortran module in which I put plans and allocatable arrays that will contain the in and out of FFTW3. Then, I can simply USE the module wherever I want.

Here is an example for 3dimensional FFT of real data:

MODULE modfft  

    IMPLICIT NONE  

    TYPE :: fftw3Needs  
        INTEGER(i4b) :: plan_forward, plan_backward ! plans for forward and inverse fft
        REAL(dp), ALLOCATABLE, DIMENSION (:,:,:) :: in_forward, out_backward ! input of FFT and output (result) of inverse FFT  
        COMPLEX(dp), ALLOCATABLE, DIMENSION (:,:,:) :: out_forward, in_backward ! output (result) of FFT and input of inverse FFT          
    END TYPE  
    TYPE (fftw3Needs) :: fftw3  

CONTAINS

    SUBROUTINE init
        CALL init_fftw3_plans! Where you initialize your plans with, e.g., dfftw_plan_dft_r2c_3d.
    END SUBROUTINE init

    SUBROUTINE deallocate_everything_fft
        CALL dfftw_destroy_plan ( fftw3%plan_forward ) ! destroy FFTW3 plans 
        CALL dfftw_destroy_plan ( fftw3%plan_backward ) ! destroy FFTW3 plans 
    END SUBROUTINE deallocate_everything_fft


END MODULE modfft

Note that this uses the old Fortran interface.

Upvotes: 1

Vladimir F Героям слава
Vladimir F Героям слава

Reputation: 60078

There are two different Fortran interfaces to FFTW. I assume FFTW3 as version 2 is deprecated.

I strongly recommend the modern Fortran interface that uses Fortran 2003 interoperability with C. The plane is then of type type(c_ptr) which is defined in the intrinsic module iso_c_binding.

In the legacy interface, the plan is indeed integer*8 on the most common 64-bit platforms.

There are many examples how to use the plan in the manual. See http://www.fftw.org/doc/Fortran-Examples.html#Fortran-Examples for the legacy interface and http://fftw.org/doc/Overview-of-Fortran-interface.html#Overview-of-Fortran-interface for the modern interface.

Whatever interface you choose, you can pass the plan variable as a dummy argument to subroutines, functions and wherever you need it to use or store.

If you encounter problems executing the plan that worked before, make sure the working arrays it refers to exist and are properly allocated. Local unsaved arrays of subroutines are not safe in this regard.

Be cautious of dummy argument arrays. They may be valid only during execution of the subroutine which created the plan because they were passed as a copy. The target argument does not help 100% because it may still be valid only inside the defining procedure.

--Important edit--

There is another problem explained in http://www.fftw.org/doc/FFTW-Execution-in-Fortran.html#FFTW-Execution-in-Fortran

Try to use

call dfftw_execute_dft(plan, in, out)

(in the legacy interface, there is an equivalent in the modern one). It may improve also the possible problem I explained two paragraphs before, because new pointers will be passed and will be OK even in the event of a copy of the original array.

Upvotes: 2

Related Questions