Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6956 Discussions

3D FFT examples: forward->backward works, backward->forward doesn't

Jozsef
Beginner
335 Views
Hi all,

I am doing a very similar test to the dfti-examples.

Take example real_3d_cce_double_ex2.c:

(1) It does a forward and a backward fft and it checks the error on the real domain to see if we get back the original data.

Now I'm trying to do the opposite with the same fft settings:

(2) I first do a backward and then a forward fft and I check the error on the complex domain to see if I can get back the original
(complex) data.

However, (1) works (i.e. I get error ~1.0e-15), but in case (2) the error is large.

Details:

The code below uses the following constants:
[cpp]const MKL_LONG N = 32;
const MKL_LONG LENGTHS[3] = {N,N,N};
const MKL_LONG PS[4] = {0,N*N,N,1};
const MKL_LONG SS[4] = {0,N*(N/2+1),N/2+1,1};
const MKL_LONG PO = N*N*N;
const MKL_LONG SO = N*N*(N/2+1);
const MKL_LONG n = 1;[/cpp]

My settings for (both forward and backward) fft are as follows:

[cpp]if ( status = DftiCreateDescriptor(fft, DFTI_DOUBLE, DFTI_REAL, 3, LENGTHS) )
    ERR(DftiErrorMessage(status));

  if ( status = DftiSetValue(*fft, DFTI_PLACEMENT, DFTI_NOT_INPLACE) )
    ERR(DftiErrorMessage(status));

  if ( status = DftiSetValue(*fft, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX) )
    ERR(DftiErrorMessage(status));

  if ( status = DftiSetValue(*fft, DFTI_PACKED_FORMAT, DFTI_CCE_FORMAT) )
    ERR(DftiErrorMessage(status));

  if ( status = DftiSetValue(*fft, DFTI_NUMBER_OF_TRANSFORMS, n) )
    ERR(DftiErrorMessage(status));

  if ( status = DftiSetValue(*fft, DFTI_FORWARD_SCALE, 1.0) )
    ERR(DftiErrorMessage(status));

  if ( status = DftiSetValue(*fft, DFTI_BACKWARD_SCALE, 1.0/(double)PO) )
    ERR(DftiErrorMessage(status));[/cpp]

Settings for only forward fft:

[cpp] // specifically initialize for forward (physical to spectral) transform
    if ( status = DftiSetValue(fft, DFTI_INPUT_DISTANCE, PO) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiSetValue(fft, DFTI_OUTPUT_DISTANCE, SO) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiSetValue(fft, DFTI_INPUT_STRIDES, PS) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiSetValue(fft, DFTI_OUTPUT_STRIDES, SS) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiCommitDescriptor(fft) )
      ERR(DftiErrorMessage(status));

    // perform forward FFT
    if ( status = DftiComputeForward(fft,u,U) )
      ERR(DftiErrorMessage(status));[/cpp]

Settings for only backward fft:

[cpp] // specifically initialize for backward (spectral to physical) transform
    if ( status = DftiSetValue(fft, DFTI_INPUT_DISTANCE, SO) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiSetValue(fft, DFTI_OUTPUT_DISTANCE, PO) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiSetValue(fft, DFTI_INPUT_STRIDES, SS) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiSetValue(fft, DFTI_OUTPUT_STRIDES, PS) )
      ERR(DftiErrorMessage(status));

    if ( status = DftiCommitDescriptor(fft) )
      ERR(DftiErrorMessage(status));

    // perform backward FFT
    if ( status = DftiComputeBackward(fft,U,u) )
      ERR(DftiErrorMessage(status));[/cpp]

The code I use to initialize, transform and then test the complex data is as follows:

[cpp]  MKL_LONG ki;
  double maxerr;
  double complex *P, *expP;
  double *p;


  // allocate memory
  if ( !(P = (double complex*)malloc(SO*sizeof(double complex))) )
    ERR("Can't allocate memory!");
  if ( !(exP = (double complex*)malloc(SO*sizeof(double complex))) )
    ERR("Can't allocate memory!");
  if ( !(p = (double*)malloc(PO*sizeof(double))) )
    ERR("Can't allocate memory!");

  // initialize complex data
  for ( ki = 0; ki < SO; ki++ )
    P[ki] = some complex random number

  // inverse & forward fft  
  do backward FFT P -> p
  do forward FFT p -> exP

  // test error
  maxerr = 0.0;
  for ( ki = 0; ki < SO; ki++ )
  {
    if (fabs(creal(P[ki]-exP[ki])) > maxerr)
      maxerr = fabs(creal(P[ki]-exP[ki]));
    if (fabs(cimag(P[ki]-exP[ki])) > maxerr)
      maxerr = fabs(cimag(P[ki]-exP[ki]));
  }
  printf("maxerr = %g\n",maxerr);[/cpp]

QUESTION: Are there any particular settings that need to be done as a setup to the FFT routines if one wants to apply backward fft first? What am I missing?

Thanks,
Jozsef
0 Kudos
3 Replies
Dmitry_B_Intel
Employee
335 Views
Hi Jozsef,

Probably, it is input data issue. Complex-to-real transform expects the input is conjugate-even [that isP(n1,n2,n3)==conj(P(N-n1,N-n2,N-n3)) ]. Ifthe input is not conjugate-even,the transform produces junk, because it doesn't attempt to 'symmetrize' the input.

Thanks
Dima
0 Kudos
Jozsef
Beginner
335 Views
But I thought the conjugate-symmetric part is not stored. Is it?
J
0 Kudos
Dmitry_B_Intel
Employee
335 Views
Most of it is not. However, some symmetry should be present (e.g. Im(Y(0,0,0)==0).
Dima
0 Kudos
Reply