Multi-dimensional FFTs in Fortran using allocatable arrays

[chambes]

The MKL DFTI uses 1-D arrays for the computation and suggests using the EQUIVALENCE statement to transform a multi-dimensional array into a 1-D array to pass to the DFTI functions.

How can one do a multi-dimensional FFT using MKL_DFTI on an allocatable array since an allocatable array cannot be used with the EQUIVALENCE statement?

[Dmitry_B_Intel]

Hi chambes,

You are right, EQUIVALNCE cannot be used to pass mD-arrays to DftiCompute functions if the arrays are allocatable or dummy arguments to the function/subroutine. A workaround to this problem is to pass the 1st column as1D subarray. AssumingX is declared as

real :: x(:,:,:)

(it can be a dummy argument or an allocatable array), the followingusewill work in many cases:

status = DftiComputeForward(hand, x(:,1,1))

Dima

[chambes]

That compiled without error, but the resulting FT was not correct. It looked like half of the 3D matrix was unchanged, while the other half had a version of the FT'ed data. I am running a rather large 3D FFT (256x256x256). It works fine if I allocate a 1D array and assign it using the RESHAPE command, but I would like to avoid increasing my memory use like that. Your method assumes the data is contiguous across the whole array - is that always correct?

Thanks for your help!

[Dmitry_B_Intel · Accepted solution]

Chambes,

That half of the 3D matrix remained unchanged may indicatemis-typedarray (real vs complex) and a version of MKL older than MKL 10.1 Update 1. Could you check what type of transform you do and what MKL version do you use?

I see that RESHAPE can do the job, and I agree it better be avoided for memory saving reason. Fortunately, the MKL interface to FFTs allows to describe nearly anypossible layout of the data, not necessarily contiguous. The layout should be properly set via DFTI_INPUT_STRIDES/DFTI_OUTPUT_STRIDES.

For the workaround of a previous postto apply, the reference passed to DftiCompute should point to the actual data and the descriptor should properly describe the layout of the data.

Simple example (256x256x256 complex-to-complex transform, in-place):

complex, allocatable :: X(:,:,:)
allocate( X( 300, 400, 500) ) ! reservemore space than needed, for illustration
status = DftiCreateDescriptor(hand, DFTI_SINGLE, DFTI_COMPLEX, 3, (/256,256,256/))
status = DftiSetValue(hand, DFTI_INPUT_STRIDES, (/0,1,300,300*400/))
status = DftiCommitDescriptor(hand)
! Input is X(i,j,k) where i,j,k are all in the range 1...256
status = DftiComputeForward(hand, X(:,1,1))
! Nowoutput is X(i,j,k) where i,j,kin the range 1...256

Dima

[chambes]

It works after setting DFTI_INPUT_STRIDES. Thanks!