topic Re: FFT examples wrong? in Intel® oneAPI Math Kernel Library

FFT examples wrong?

scrognoid — Thu, 07 May 2009 22:49:59 GMT

Is there not a problem in the Fortran FFT code examples:

! Fortran example.
! 1D complex to complex, and real to conjugate-even
Use MKL_DFTI
Complex :: X_in(32)
Complex :: X_out(32)
Real :: Y_in(32)
Real :: Y_out(34)
type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle, My_Desc2_Handle
Integer :: Status
...put input data into X_in(1),...,X_in(32); Y_in(1),...,Y_in(32)
! Perform a complex to complex transform
Status = DftiCreateDescriptor( My_Desc1_Handle, DFTI_SINGLE,
DFTI_COMPLEX, 1, 32 )
Status = DftiSetValue( My_Desc1_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
Status = DftiCommitDescriptor( My_Desc1_Handle )
Status = DftiComputeForward( My_Desc1_Handle, X_in, X_out )
Status = DftiFreeDescriptor(My_Desc1_Handle)
! result is given by {X_out(1),X_out(2),...,X_out(32)}
! Perform a real to complex conjugate-even transform
Status = DftiCreateDescriptor(My_Desc2_Handle, DFTI_SINGLE,
DFTI_REAL, 1, 32)
Status = DftiSetValue( My_Desc2_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
Status = DftiCommitDescriptor(My_Desc2_Handle)
Status = DftiComputeForward(My_Desc2_Handle, Y_in, Y_out)
Status = DftiFreeDescriptor(My_Desc2_Handle)
! result is given by Y_out in CCS format.

Y_out is complex, yes?

Re: FFT examples wrong?

Dmitry_B_Intel — Fri, 08 May 2009 03:24:42 GMT

Hi,

In this real-to-complex transform the resulting data, which is32/2+1=17 complex-valued items, is laid out to 34 real-valued items of array Y_out, according to the default layout for the conjugate-even domain (which is DFTI_COMPLEX_REALwith CCS packing). The complex values are stored as ,,,,... and adding a stride will separate real and imaginary part of each complex number.

The result may really be stored as complex numbers, [RI],[RI],... so that adding a stride will keep the complex numbers. For this one should use CCE storage for conjugate-even domain, by calling
DftiSetValue(hand,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX) and using complex Y_out(17).

In the case stride is 1 the two abovementioned layouts, i.e. CCS andCCE, are equivalent.

Internally the actual type of the data arguments in DftiComputeForward is disregarded (starting with MKL 10.1 Update 1) and the data is laid out according to the configuration of the descriptor, so it should make no difference whether Y_out is declared ascomplex orreal array.

Thanks,
Dima

Re: FFT examples wrong?

scrognoid — Fri, 08 May 2009 16:02:21 GMT

So, for me, it sounds like the example will work and I can simply change Y-out to complex(17), or if I want to spell it out for clarity, I would add the line you indicate:

! Fortran example.
! 1D real to conjugate-even
Use MKL_DFTI
Real:: Y_in(32)
Complex:: Y_out(17)
type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle, My_Desc2_Handle
Integer :: Status
! Perform a real to complex conjugate-even transform
Status = DftiCreateDescriptor(My_Desc2_Handle, DFTI_SINGLE,
DFTI_REAL, 1, 32)
Status = DftiSetValue( My_Desc2_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
Status = DftiSetValue(My_Desc2_Handle,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX)
Status = DftiCommitDescriptor(My_Desc2_Handle)
Status = DftiComputeForward(My_Desc2_Handle, Y_in, Y_out)
Status = DftiFreeDescriptor(My_Desc2_Handle)
! result is given by Y_out in CCE format.

Do I have that right?

Re: FFT examples wrong?

Dmitry_B_Intel — Tue, 12 May 2009 01:36:55 GMT

Yes, that's right.