In order to obtain maximum performance in FFTW, you should store your
data in arrays that have been specially aligned in memory (see SIMD alignment and fftw_malloc). Enforcing alignment also permits you to
safely use the new-array execute functions (see New-array Execute Functions) to apply a given plan to more than one pair of in/out
arrays. Unfortunately, standard Fortran arrays do not provide
any alignment guarantees. The only way to allocate aligned
memory in standard Fortran is to allocate it with an external C
function, like the
fftw_alloc_complex functions. Fortunately, Fortran 2003 provides
a simple way to associate such allocated memory with a standard Fortran
array pointer that you can then use normally.
We therefore recommend allocating all your input/output arrays using the following technique:
arr, to your array of the desired type and dimensions. For example,
real(C_DOUBLE), pointer :: a(:,:)for a 2d real array, or
complex(C_DOUBLE_COMPLEX), pointer :: a(:,:,:)for a 3d complex array.
integer(C_SIZE_T). You can either declare a variable of this type, e.g.
integer(C_SIZE_T) :: sz, to store the number of elements to allocate, or you can use the
int(..., C_SIZE_T)intrinsic function. e.g. set
sz = L * M * Nor use
int(L * M * N, C_SIZE_T)for an L × M × N array.
type(C_PTR) :: pto hold the return value from FFTW's allocation routine. Set
p = fftw_alloc_real(sz)for a real array, or
p = fftw_alloc_complex(sz)for a complex array.
arrwith the allocated memory
pusing the standard
call c_f_pointer(p, arr, [...dimensions...]), where
[...dimensions...])are an array of the dimensions of the array (in the usual Fortran order). e.g.
call c_f_pointer(p, arr, [L,M,N])for an L × M × N array. (Alternatively, you can omit the dimensions argument if you specified the shape explicitly when declaring
arr.) You can now use
arras a usual multidimensional array.
For example, here is how we would allocate an L × M 2d real array:
real(C_DOUBLE), pointer :: arr(:,:) type(C_PTR) :: p p = fftw_alloc_real(int(L * M, C_SIZE_T)) call c_f_pointer(p, arr, [L,M]) ...use arr and arr(i,j) as usual... call fftw_free(p)
and here is an L × M × N 3d complex array:
complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:) type(C_PTR) :: p p = fftw_alloc_complex(int(L * M * N, C_SIZE_T)) call c_f_pointer(p, arr, [L,M,N]) ...use arr and arr(i,j,k) as usual... call fftw_free(p)
See Reversing array dimensions for an example allocating a single array and associating both real and complex array pointers with it, for in-place real-to-complex transforms.