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

Jozsef — Wed, 22 Sep 2010 00:16:17 GMT

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

3D FFT examples: forward->backward works, backward->forward doe

Dmitry_B_Intel — Wed, 22 Sep 2010 07:36:05 GMT

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

3D FFT examples: forward->backward works, backward->forward doe

Jozsef — Fri, 29 Oct 2010 17:21:58 GMT

But I thought the conjugate-symmetric part is not stored. Is it?
J

3D FFT examples: forward->backward works, backward->forward doe

Dmitry_B_Intel — Sat, 30 Oct 2010 16:17:36 GMT

Most of it is not. However, some symmetry should be present (e.g. Im(Y(0,0,0)==0).
Dima

topic 3D FFT examples: forward->backward works, backward->forward doe in Intel® oneAPI Math Kernel Library

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

3D FFT examples: forward->backward works, backward->forward doe

3D FFT examples: forward->backward works, backward->forward doe

3D FFT examples: forward->backward works, backward->forward doe