3D DFT

[jed_p_]

I am trying to do a 3D FFT in MKL 11.3, but its not working for me... most likely I am doing something wrong, but I am at a loss.

I get an error when I attempt to call DftiCommitDescriptor.  I tried two different cases and got two different errors.  My code is shown below... Any help would be appreciated!

#include <mkl.h>
#include <stdio.h>

int main()

{

  DFTI_DESCRIPTOR_HANDLE dh;

  float *data;
  MKL_LONG sz[3];
  MKL_LONG dim = 3;
  int val;
  int error_code;

  sz[0] = 128;
  sz[1] = 128;
  sz[2] = 128;

  dim = 2;
  error_code = DftiCreateDescriptor(&dh, DFTI_SINGLE, DFTI_REAL, dim, sz);
  error_code = DftiCommitDescriptor(dh);
  printf("error_code: %d\n", error_code);

  dim = 3;
  error_code = DftiCreateDescriptor(&dh, DFTI_SINGLE, DFTI_REAL, dim, sz);
  error_code = DftiCommitDescriptor(dh);
  printf("error_code: %d\n", error_code);

  printf("%s\n", DftiErrorMessage(error_code));

  dim = 3;
  error_code = DftiCreateDescriptor(&dh, DFTI_SINGLE, DFTI_REAL, dim, sz);
  error_code = DftiSetValue(dh, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
  error_code = DftiCommitDescriptor(dh);
  printf("error_code: %d\n", error_code);

  printf("%s\n", DftiErrorMessage(error_code));

}

Note that the first call works (dim is set to 2).

The second call fails (only difference is that dim is 3), indicating that the functionality is not implemented.

I read in the documentation that I should use DFTI_CONJUGATE_EVEN_STORAGE <- DFTI_COMPLEX_COMPLEX, so I tried this as well (third call), but it failed with the error "Inconsistent configuration parameters". 

The output I get is as follows:

error_code: 0
error_code: 6
Intel MKL DFTI ERROR: Functionality is not implemented
error_code: 3
Intel MKL DFTI ERROR: Inconsistent configuration parameters

[jed_p_]

OK... I figured it out... I am posting what I had done wrong here in case anyone else does the same thing.  Also, I still have one thing that bugs me (see bottom of this comment):

I found in the documentation that the default settings of the strides is not valid for 3D FFTs.  They must be set explicitly.  The third section of my code now looks like this:

  dim = 3;
  error_code = DftiCreateDescriptor(&dh, DFTI_DOUBLE, DFTI_REAL, dim, sz);
  error_code = DftiSetValue(dh, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
  rstrides[0] = 0;
  rstrides[1] = sz[1]*(sz[2]/2+1)*2;
  rstrides[2] = (sz[2]/2+1)*2;
  rstrides[3] = 1;
  error_code = DftiSetValue(dh, DFTI_INPUT_STRIDES, rstrides);
  cstrides[0] = 0;
  cstrides[1] = sz[1]*(sz[2]/2+1);
  cstrides[2] = (sz[2]/2+1);
  cstrides[3] = 1;
  error_code = DftiSetValue(dh, DFTI_OUTPUT_STRIDES, cstrides);
  error_code = DftiCommitDescriptor(dh);
  printf("error_code: %d\n", error_code);

There is no error now.  :-)

One thing still troubles me though:  It seems like if I want to do a 3D FFT, I need to ensure that the stride in the first (fastest changing) dimension is bigger (by one or two) than the size of that dimension.  That means that the input data needs to have some padding (gaps).  Is there any way around this?  It seems to me that refactoring my data to conform to this may take a significant amount of time, though I have not timed it yet... maybe it is negligible compared to the time of the FFT.

[Dmitry_B_Intel · Accepted solution]

Hi jed,

Padding is necessary for 3D real-to-complex FFT. MKL does not support 'packed storage' for 3D. It also does not support real-to-real FFT.

If you rearrange your data allocation so as to allow padding, you will solve the problem. And you may also get some small performance increase if you set the leading dimension to a 64-byte multiple size, e.g. instead of "rstrides[2] = (sz[2]/2+1)*2"   set something like "rstrides[2] = UP_TO_MUL8((sz[2]/2+1)*2  )" with UP_TO_MUL8(n) defined as  ((((n)+7)/8)*8). Of course rstrides[1] should then be defined as "rstrides[1] = sz[1]*rstrides[2]" and cstrides set accordingly.

Thanks,
Dima

[jed_p_]

Dima,

Thanks so much for the great response.