FFTW 3.3.8: Reversing array dimensions

A minor annoyance in calling FFTW from Fortran is that FFTW’s array cannam@167: dimensions are defined in the C convention (row-major order), while cannam@167: Fortran’s array dimensions are the opposite convention (column-major cannam@167: order). See Multi-dimensional Array Format. This is just a cannam@167: bookkeeping difference, with no effect on performance. The only cannam@167: consequence of this is that, whenever you create an FFTW plan for a cannam@167: multi-dimensional transform, you must always reverse the cannam@167: ordering of the dimensions. cannam@167:

For example, consider the three-dimensional (L × M × N cannam@167: ) arrays: cannam@167:

cannam@167:

  complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out
cannam@167:

To plan a DFT for these arrays using fftw_plan_dft_3d, you could do: cannam@167:

cannam@167:

  plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
cannam@167:

That is, from FFTW’s perspective this is a N × M × L cannam@167: array. cannam@167: No data transposition need occur, as this is only cannam@167: notation. Similarly, to use the more generic routine cannam@167: fftw_plan_dft with the same arrays, you could do: cannam@167:

cannam@167:

  integer(C_INT), dimension(3) :: n = [N,M,L]
cannam@167:   plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
cannam@167:

Note, by the way, that this is different from the legacy Fortran cannam@167: interface (see Fortran-interface routines), which automatically cannam@167: reverses the order of the array dimension for you. Here, you are cannam@167: calling the C interface directly, so there is no “translation” layer. cannam@167:

An important thing to keep in mind is the implication of this for cannam@167: multidimensional real-to-complex transforms (see Multi-Dimensional DFTs of Real Data). In C, a multidimensional real-to-complex DFT cannam@167: chops the last dimension roughly in half (N × M × L cannam@167: real input cannam@167: goes to N × M × L/2+1 cannam@167: complex output). In Fortran, because cannam@167: the array dimension notation is reversed, the first dimension of cannam@167: the complex data is chopped roughly in half. For example consider the cannam@167: ‘r2c’ transform of L × M × N cannam@167: real input in Fortran: cannam@167:

cannam@167:

  type(C_PTR) :: plan
cannam@167:   real(C_DOUBLE), dimension(L,M,N) :: in
cannam@167:   complex(C_DOUBLE_COMPLEX), dimension(L/2+1,M,N) :: out
cannam@167:   plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
cannam@167:   ...
cannam@167:   call fftw_execute_dft_r2c(plan, in, out)
cannam@167:

Alternatively, for an in-place r2c transform, as described in the C cannam@167: documentation we must pad the first dimension of the cannam@167: real input with an extra two entries (which are ignored by FFTW) so as cannam@167: to leave enough space for the complex output. The input is cannam@167: allocated as a 2[L/2+1] × M × N cannam@167: array, even though only cannam@167: L × M × N cannam@167: of it is actually used. In this example, we will cannam@167: allocate the array as a pointer type, using ‘fftw_alloc’ to cannam@167: ensure aligned memory for maximum performance (see Allocating aligned memory in Fortran); this also makes it easy to reference the cannam@167: same memory as both a real array and a complex array. cannam@167:

cannam@167:

  real(C_DOUBLE), pointer :: in(:,:,:)
cannam@167:   complex(C_DOUBLE_COMPLEX), pointer :: out(:,:,:)
cannam@167:   type(C_PTR) :: plan, data
cannam@167:   data = fftw_alloc_complex(int((L/2+1) * M * N, C_SIZE_T))
cannam@167:   call c_f_pointer(data, in, [2*(L/2+1),M,N])
cannam@167:   call c_f_pointer(data, out, [L/2+1,M,N])
cannam@167:   plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
cannam@167:   ...
cannam@167:   call fftw_execute_dft_r2c(plan, in, out)
cannam@167:   ...
cannam@167:   call fftw_destroy_plan(plan)
cannam@167:   call fftw_free(data)
cannam@167:

7.2 Reversing array dimensions