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). 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)