Community
cancel
Showing results for 
Search instead for 
Did you mean: 
158 Views

Complex to real IFFT using MKL

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

 

0 Kudos
7 Replies
Ying_H_Intel
Employee
158 Views

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

 

158 Views

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

158 Views

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
Employee
158 Views

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
Employee
158 Views

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

 

158 Views

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
Employee
158 Views

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. 

Best Regards,

Ying 

 

 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))