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