DFT complex conjugate even data structure

[Thomas_D_1]

I'm trying to do a comparison between MKLs DFT for a complex --> complex forward DFT vs a real --> real forward DFT. In both cases my starting array is entirely real so the real --> real transform should be able to take advantage of complex conjugate even symmetries. The complex transform is giving me the expected results. However, for the real transform the output seems completely jumbled. It's possible I'm not applying the strides correctly but the code I'm using is practically copied directly from the MKL examples so that seems unlikely. From the documentation the formatting for the data in the output array is suppose to be CCE. For all the other formats the documentation provides an explicit (and pictorial) description of how the data is laid out. But I can't seem to find anything describing how the data for the CCE format is distributed. Is there a reference that covers this?

[Evgueni_P_Intel]

Hi Thomas D.,

If the 1D CCE format is used, then only the first N/2+1 complex elements of the transformed real sequence are stored in the output array.

For the multi-dimensional CCE format, please read the following thread on this forum.

https://software.intel.com/en-us/forums/topic/288202

Evgueni.

[Thomas_D_1]

Thanks for the response Evgueni,

If I understand the post you linked correctly, that post is showing how to go from a half sized complex array that has complex conjugate even symmetry to the full complex array. My question if slightly different. I'm trying to understand how the real and imaginary parts of the CCE format are distributed in a strictly real output array. My thought originally was that they would alternate R(i,j,..) I(i,j,...) along the first index, giving an effective complex array of size [N1/2,N2,N3,..], however from my testing that doesn't appear to be the case (unless of course I've made a silly mistake with my input and output stides).

[Evgueni_P_Intel]

By default, Intel MKL stores the real and imaginary parts of a complex number next to each other in the output array, the real part being closer to address 0.

Here is an example of the default strides for a 3D real-to-complex transform.

MKL_LONG N[] = {16, 32, 64};
// default strides for out-of-place transforms
MKL_LONG default_oop_is[] = {0, N[1]*N[2]      , N[2]    , 1};
MKL_LONG default_oop_os[] = {0, N[1]*(N[2]/2+1), N[2]/2+1, 1};
// default strides for in-place transforms
MKL_LONG default_ip_is[] = {0, N[1]*(N[2]/2+1)*2, (N[2]/2+1)*2, 1};
MKL_LONG default_ip_os[] = {0, N[1]*(N[2]/2+1)  , N[2]/2+1    , 1};

By default, the size of the output array is N[0]*N[1]*(N[2]/2+1) complex numbers.
 

[Thomas_D_1]

Thanks again Evgueni,

It turns out my problem came from not setting DFTI_OUTPUT_DISTANCE correctly. The transform I was using was a 3+1 D transform, and while I had set the input and output strides correctly as you described in your previous post, I was only setting the DFTI_INPUT_DISTANCE. Since this was a complex conjugate even, real to real, in place transform the meant setting

! system size = [N1,N2,N3,M]
real :: array(2*(N1/2+1),N2,N3,M)
! create descriptor
status = DftiCreateDescriptor(hand, DFTI_SINGLE, DFTI_REAL, 3, [N1,N2,N3])
! set transform to use complex-conjugate-even symmetry
statsu = DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)
! set strides
status = DftiSetValue(hand, DFTI_INPUT_STRIDES , [0, 1, 2*(N1/2+1), 2*(N1/2+1)*N2])
status = DftiSetValue(hand, DFTI_OUTPUT_STRIDES, [0, 1,   (N1/2+1),   (N1/2+1)*N2])
! set number of transforms
status = DftiSetValue(hand, DFTI_NUMBER_OF_TRANSFORMS, M)
! set distance for doing multiple transforms
status = DftiSetValue(hand, DFTI_INPUT_DISTANCE , 2*(N1/2+1)*N2*N3)
status = DftiSetValue(hand, DFTI_OUTPUT_DISTANCE,   (N1/2+1)*N2*N3) ! this was the line I was missing...