Next: FFTW Fortran type reference, Previous: Overview of Fortran interface, Up: Calling FFTW from Modern Fortran [Contents][Index]

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). See 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 × M × 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 × M × 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 (see 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 (see Multi-Dimensional DFTs of Real Data). In C, a multidimensional real-to-complex DFT
chops the last dimension roughly in half (N × M × L
real input
goes to N × M × 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 × M × 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] × M × N
array, even though only
L × M × N
of it is actually used. In this example, we will
allocate the array as a pointer type, using ‘`fftw_alloc`’ to
ensure aligned memory for maximum performance (see 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)