Complex to real IFFT using MKL

karagiannis__dionysi · ‎06-18-2015

Hi ,

I need to perform an IDFT of some complex data to real in FORTRAN. I have read online a lot of how to do that with MKL but it confused me more. So I need to ask some questions on the configurations I need to make to so some a transformation.

My complex data come from a simulation of N^3 grid, where N=512. They are all inside a complex array Z and are conjugate symmetric in the sense that each line corresponds to a vector k of coordinates (i,j,l) in such a way that N*N*i+N*j+l is the line number. The very first element of the array is real, I think is the mean value of the field (DC), and Z(N/2) element is again real. After that Z(N/2+l) = conjg(Z(N/2-l)). So in general I have the same packing for all dimensions: Z(N*N*(N/2+i)+N*(N/2+j)+N/2+l) = conjg(Z(N*N*(N/2-i)+N*(N/2-j)+N/2-l)).

My first question is do I need to perform a 3D or a 1D FFT? Can I specify the structure of my complex input array with a stride so MKL will know where are the other dimensions? If I have understood well how to make the strides I need to do something like this:

stat = DftiCreateDescriptor( desc_handle, DFTI_DOUBLE, DFTI_REAL, 1,)
stat = DftiSetValue(desc_handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)

stat = DftiSetValue(desc_handle, DFTI_INPUT_STRIDES,(/0,1,N,N**2/))

Then I need to perform the backward transformation but I do not know in what format should my input complex array should be. In the examples of intel it shows that the default storage scheme is the CCE, where I need to save half of them for each dimension. But then I saw that there is also a CCS storage scheme. In which of the two schemes must my input array be?

Due to the conjugate symmetry mathematically my output array is expected to be a complex fortran array with zero imaginary part. I have tried a test code out but I dont get that. My code runs perfect but the results are not those expected. What have I've done wrong? Thanks in advance for your help and time.

Regards,

Dionisis

Ying_H_Intel · ‎06-19-2015

Hi Dionisis,

Just quick comment,

Regarding 3D or 1D FFT, it may depend on what you need to do. 3D FFT doesn't equail of several 1D FFTs.

Regarding the IFFT, MKL provide some fortran code. You may build them and see if it can help.

MKLexample\dftf\source\basic_dp_real_dft_3d.f90. It include real to complex, complex to real IFFT.

Best Regards,

Ying

print *,"Example basic_dp_real_dft_3d"
print *,"Forward-Backward double-precision real out-of-place 3D FFT"
print *,"Configuration parameters:"
print *,"DFTI_PRECISION = DFTI_DOUBLE"
print *,"DFTI_FORWARD_DOMAIN = DFTI_REAL"
print *,"DFTI_DIMENSION = 3"
print '(" DFTI_LENGTHS = /"I0","I0","I0"/" )', N1, N2, N3
print *,"DFTI_PLCEMENT = DFTI_NOT_INPLACE"
print *,"DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX"

print *,"Create DFTI descriptor for real transform"
status = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_REAL, 3, [N1,N2,N3])
if (0 /= status) goto 999

print *,"Set out-of-place"
status = DftiSetValue(hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
if (0 /= status) goto 999

print *,"Set CCE storage"
status = DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE, &
DFTI_COMPLEX_COMPLEX)
if (0 /= status) goto 999

cstrides = [0, 1, N1/2+1, N2*(N1/2+1)]
rstrides = [0, 1, N1, N2*N1]

print '(" Set input strides = "4(I0:", "))', rstrides
status = DftiSetValue(hand, DFTI_INPUT_STRIDES, rstrides)
if (0 /= status) goto 999

print '(" Set output strides = "4(I0:", "))', cstrides
status = DftiSetValue(hand, DFTI_OUTPUT_STRIDES, cstrides)
if (0 /= status) goto 999

print *,"Commit DFTI descriptor"
status = DftiCommitDescriptor(hand)
if (0 /= status) goto 999

print *,"Allocate data arrays"
allocate(x_real(N1, N2, N3))
allocate(x_cmplx(N1/2+1, N2, N3))

print *,"Initialize data for real-to-complex FFT"
call init_r(x_real, N1, N2, N3, H1, H2, H3)

print *,"Compute forward transform"
status = DftiComputeForward(hand, x_real(:, 1, 1), x_cmplx(:, 1, 1))
if (0 /= status) goto 999

print *,"Verify the result"
status = verify_c(x_cmplx, N1, N2, N3, H1, H2, H3)
if (0 /= status) goto 999

print *,"Reconfigure DFTI descriptor for backward transform"

print '(" Set input strides = "4(I0:", "))', cstrides
status = DftiSetValue(hand, DFTI_INPUT_STRIDES, cstrides)
if (0 /= status) goto 999

print '(" Set output strides = "4(I0:", "))', rstrides
status = DftiSetValue(hand, DFTI_OUTPUT_STRIDES, rstrides)
if (0 /= status) goto 999

print *,"Recommit DFTI descriptor"
status = DftiCommitDescriptor(hand)
if (0 /= status) goto 999

print *,"Initialize data for complex-to-real FFT"
call init_c(x_cmplx, N1, N2, N3, H1, H2, H3)

print *,"Compute backward transform"
status = DftiComputeBackward(hand, x_cmplx(:, 1, 1), x_real(:, 1, 1))
if (0 /= status) goto 999

print *,"Verify the result"
status = verify_r(x_real, N1, N2, N3, H1, H2, H3)
if (0 /= status) goto 999

karagiannis__dionysi · ‎06-19-2015

Hi Ying,

Thank you for your answer and help. I need to do a 3D transformation since the complex coefficients correspond to a 3D Fourier space vector.

I have seen these examples but they haven't helped much. What I want to know is if I keep the CCE storage format which is the defaults for the DFTI_CONJUGATE_EVEN_STORAGE= DFTI_COMPLEX_COMPLEX how should my input complex data look like?

My original data are a 1D fortran array which follow this:

Suppose you have 3D array A of size NxNxN, where:

1st dimension has size N and has stride 1
2nd dimension has size N and has stride N
3rd dimension has size N and has stride N*N

So element A should be accessed as A[k*N*N + j*N + i].

So the stride should be: cstride=(/0,1,N,N*N/) right or do I need to put only N/2+1 elements for each dimension since my data have a conjugate symmetry for every dimension. Then how should the stride and the input data should look like so the MKL FFT routines can read and use this symmetry properly and give the real results I need?

Kind Regards,

Dionisis

karagiannis__dionysi · ‎06-19-2015

Hi again,

I order to simplify my problem I am just testing a 1D IFFT. The input is a complex array with z(1) and z(N/2+1) being real. The total size of the array is N=256. If I prepare the IFFT like the following:

stat = DftiCreateDescriptor( desc_handle, DFTI_DOUBLE, DFTI_REAL, 1,256)
stat = DftiSetValue(desc_handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)

stat = DftiSetValue(desc_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)

stat = DftiComputeBackward(desc_handle, X_in, M_out)

,where X_in has the conjugate symmetry even described above. That means that M_out is expected to be a real array in the sense that if:

real(dp) :: M_out(N+2)

then the expected real fortran array with have every other element equal to zero.

complex(dp) :: M_out(N/2)

then the expected complex fortran array will have N/2 size with the imaginary part being zero.

However the results I get are not real after doing the above. It is like the routines do not understand that the input complex fortran array does not have the conjugate even symmetry. Why is that? Do I have to add any other preference parameters to ensure the structure of the input is read correctly?

Regards,

Dionisis

Ying_H_Intel · ‎06-23-2015

Hi Dionisis,

You may refer to the sample, basic_dp_real_dft_1d.f90. It is exactly sample for

DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX and input for IDFT is x_cmplx(N/2+1)

print *,"Example basic_dp_real_dft_1d"
print *,"Forward-Backward double-precision real-to-complex", &
" out-of-place 1D transform"
print *,"Configuration parameters:"
print *,"DFTI_PRECISION = DFTI_DOUBLE"
print *,"DFTI_FORWARD_DOMAIN = DFTI_REAL"
print *,"DFTI_DIMENSION = 1"
print '(" DFTI_LENGTHS = /"I0"/" )', N
print *,"DFTI_PLACEMENT = DFTI_NOT_INPLACE"
print *,"DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX"

allocate(x_real(N))
allocate(x_cmplx(N/2+1))

If you don't need the Forward part, you can delete them. i modify the sample and attach here.

Not sure , how your input look likes (is a complex array with z(1) and z(N/2+1) being real ? , where X_in has the conjugate symmetry even described above)

If you setDFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX. Then your X_in should be (x_cmplx(N/2+1)) here. and I mark it as bold. They are half of the FFT 's result.

Initialize data for real-to-complex FFT
102.000000000000 55.0000000000000 89.0000000000000
12.0000000000000 3.00000000000000 45.0000000000000
9.00000000000000 8.00000000000000
Compute forward transform
Verify the result
(323.000000000000,0.000000000000000E+000) (103.242640687119,-89.8994949366
(7.00000000000000,-80.0000000000000) (94.7573593128807,70.1005050633883)
(83.0000000000000,0.000000000000000E+000)
323.000000000000 103.242640687119 7.00000000000000
94.7573593128807 83.0000000000000
Create DFTI descriptor
Set out-of-place
Set CCE storage
Commit DFTI descriptor
Compute backward transform
Verify the result
102.000000000000 55.0000000000000 89.0000000000000
12.0000000000000 3.00000000000000 45.0000000000000
9.00000000000000 8.00000000000000
Release the DFTI descriptor
Release the DFTI_I descriptor
Deallocate data arrays
TEST PASSED

If you set DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_real, and DFTI_PACKED_FORMAT = DFTI_CCS_FORMAT.

Then the x_in should be N+2. or (equivalence a complex array with z(1) and z(N/2+1) being real, same above blod ) and output is x_real(N).

Best Regards,

Ying

Ying_H_Intel · ‎06-23-2015

Attach the fortran code file.

The packed format like CCS, park, perm are described in MKl reference manual: https://software.intel.com/zh-cn/node/521971 =>DFTI_PACKED_FORMAT .

Regards,

Ying

karagiannis__dionysi · ‎06-24-2015

Hi Ying,

Thanks for your help it did the trick. If I want to perform a 3D IFFT then I have to put half of them again right? So my X_in will be X_in((N/2+1)*N*N) correct?

Thanks again for your time and help.

Regards,

Dionisis

Ying_H_Intel · ‎06-24-2015

Hi Dionisis,

Right, if 3D IFFT, you have to put half of them and your X_in for IDFT will be X_in((N/2+1)*N*N).

the sample basic_dp_real_dft_3d.f90 code under <MKL install directory>/sample/mkl_core.zip=> unzip, \dftf\source.