Info Catalog fftw3: Overview of Fortran interface fftw3: Calling FFTW from Modern Fortran fftw3: FFTW Fortran type reference

fftw3: Reversing array dimensions

 
 7.2 Reversing array dimensions
 ==============================
 
 A minor annoyance in calling FFTW from Fortran is that FFTW's array
 dimensions are defined in the C convention (row-major order), while
 Fortran's array dimensions are the opposite convention (column-major
 order).  Multi-dimensional Array Format.  This is just a
 bookkeeping difference, with no effect on performance.  The only
 consequence of this is that, whenever you create an FFTW plan for a
 multi-dimensional transform, you must always _reverse the ordering of
 the dimensions_.
 
    For example, consider the three-dimensional (L x M x N ) arrays:
 
        complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out
 
    To plan a DFT for these arrays using 'fftw_plan_dft_3d', you could
 do:
 
        plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
 
    That is, from FFTW's perspective this is a N x M x L array.  _No data
 transposition need occur_, as this is _only notation_.  Similarly, to
 use the more generic routine 'fftw_plan_dft' with the same arrays, you
 could do:
 
        integer(C_INT), dimension(3) :: n = [N,M,L]
        plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
 
    Note, by the way, that this is different from the legacy Fortran
 interface (Fortran-interface routines), which automatically
 reverses the order of the array dimension for you.  Here, you are
 calling the C interface directly, so there is no "translation" layer.
 
    An important thing to keep in mind is the implication of this for
 multidimensional real-to-complex transforms (Multi-Dimensional
 DFTs of Real Data).  In C, a multidimensional real-to-complex DFT
 chops the last dimension roughly in half (N x M x L real input goes to N
 x M x L/2+1 complex output).  In Fortran, because the array dimension
 notation is reversed, the _first_ dimension of the complex data is
 chopped roughly in half.  For example consider the 'r2c' transform of L
 x M x N real input in Fortran:
 
        type(C_PTR) :: plan
        real(C_DOUBLE), dimension(L,M,N) :: in
        complex(C_DOUBLE_COMPLEX), dimension(L/2+1,M,N) :: out
        plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
        ...
        call fftw_execute_dft_r2c(plan, in, out)
 
    Alternatively, for an in-place r2c transform, as described in the C
 documentation we must _pad_ the _first_ dimension of the real input with
 an extra two entries (which are ignored by FFTW) so as to leave enough
 space for the complex output.  The input is _allocated_ as a 2[L/2+1] x
 M x N array, even though only L x M x N of it is actually used.  In this
 example, we will allocate the array as a pointer type, using
 Allocating aligned memory in Fortran::); this also makes it easy to
 reference the same memory as both a real array and a complex array.
 
        real(C_DOUBLE), pointer :: in(:,:,:)
        complex(C_DOUBLE_COMPLEX), pointer :: out(:,:,:)
        type(C_PTR) :: plan, data
        data = fftw_alloc_complex(int((L/2+1) * M * N, C_SIZE_T))
        call c_f_pointer(data, in, [2*(L/2+1),M,N])
        call c_f_pointer(data, out, [L/2+1,M,N])
        plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
        ...
        call fftw_execute_dft_r2c(plan, in, out)
        ...
        call fftw_destroy_plan(plan)
        call fftw_free(data)
 

Info Catalog fftw3: Overview of Fortran interface fftw3: Calling FFTW from Modern Fortran fftw3: FFTW Fortran type reference