Real to Half-complex 3d Cluster FFT mkl implementation

[giorgos17]

Greetings,

I would like to ask if there is any information available about plans to implement real to half-complex Cluster FFT routines in future versions of mkl. If I am not mistaken, such routines are currently unavailable.

My implementation of the real to half-complex cluster FFT based on the complex to complex cluster FFT mkl routine spends an unacceptable amount of time in communicating complex numbers between positive and negative wavevectors in order to make a conjugate-symmetric complex input for the mkl routine.

Is there any better way to do this than the obvious one I am currently using?

[Dmitry_B_Intel]

Hi giorgos17,

Could you specify what do you mean by 'half-complex 3d'? This term may have different meaning than what you intended to say. For example, in FFTW this term refers to peculiar layout of the data in the frequency domain of 1D real-to-complex transform. For one 1D half-complex result is laid out as all real parts followed by all imaginary parts in the reverse order (R0,R1,R2,...I3,I2,I1). How do you want it be extended to 3D?

Thanks
Dima

[giorgos17]

Hi Dima,

I use the following half-complex format:

for odd i only:

r(i,j,k) = Re(c(floor(i/2)+1,j,k))

r(i+1,j,k)=Im(c(floor(i/2)+1,j,k))

where c is the complex array and r is the corresponding half-complex array.

i runs from 1 to n, so that only positive frequencies are stored in the half-complex

array, and their symmetric conjugates in the negative frequencies are dropped.

Giorgos

[Vladimir_Petrov__Int]

Hi Giorgos,

Actually MKL does have the functionality you are asking for. What it lacks are examples illustrating this type of transforms via Dfti-like family of functions. Such examples will be added in one of our future releases.

Meanwhile you can get a notion of how to do it by looking at the following C example

examples\fftw2x_cdft\source\wrappers_r3d.c

on how to do it in C via our FFTW MPI wrappers.

Best regards,

-Vladimir

[giorgos17]

Thank you Vladimir, is there a simple way I can use these routines form Fortran 90?

Best regards,

Giorgos

[Vladimir_Petrov__Int]

Giorgos,

If you prefer to use the FFTW MPI interface then use the following routines provided by MKL wrappers:

rfftw3d_f77_mpi_create_plan

rfftwnd_f77_mpi_local_sizes

rfftwnd_f77_mpi

rfftwnd_f77_mpi_destroy_plan
.

If you choose the Dfti-like interface you just build the module from the file include/mkl_cdft.f90 and use the functions:

DftiCreateDescriptorDM

DftiSetValueDM/DftiGetValueDM

DftiCommitDescriptorDM

DftiComputeForwardDM/DftiComputeBackwardDM

DftiFreeDescriptorDM

.

Best regards,

-Vladimir

[giorgos17]

I prefer the Dfti-like interface, but it is impossible to call DftiComputeForwardDM/DftiComputeBackwardDM

with a real array as argument, there is no matching interface.

[Dmitry_B_Intel]

You are right, CFFT doesn't provide Compute functions that would accept real array. This will be fixed.

At this point, however, tricking the compute function into accepting real arrays as complex shouldwork. Unfortunately, Fortran doesn't providelegal means to storage associate dynamically allocated real array withcomplex array (let me know if you know how to do this).So I'd suggest you tryingto bypass type checking.

Thanks
Dima

[giorgos17]

No, as far as i know there is no way to storage-associate an allocatable array, but using an external procedure

to bypass type checking should work.

Do I have to set some particular value via DftiSetValueDM, or should I just use defaults?

Thanks

Giorgos

[Vladimir_Petrov__Int]

Giorgos,

You will definitely have to call DftiSetValueDM(h, DFTI_PLACEMENT, DFTI_NOT_INPLACE) if you want the transform to be done out-of-place.

Best regards,

-Vladimir

[giorgos17]

Thanks, I will try your suggestion.