Trouble computing a 2D iFFT with the MKL using DftiComputeBackward

[CSA_O]

Hello,

I'm trying to use a 2D backward FFT with the MKL.

I can compute a 1D or a 2D by One-dimensional transform but I have trouble to expand it to a 2D.

My guess would be that I don't fill my input data correctly or that the strides are not done correctly. I've found somewhere that the input and output strides have to be switched when using DftiComputeBackward, but I'm not 100% sure of that.

I import my input data from two files I filled using fft() in Matlab, one with the real part and the other one with the imaginary part. The complex data is stored in tfNewComp.

The final code should be able to work in 3D but for the time beeing I'm just trying to make it work in 2D, that's why there is a dimension (nX) equal to 1.

I only get 0 and -0 from the DftiComputeBackward function.

My sandbox is pretty straight forward so I can only assume I'm missing something in the configuration before calling the function.

Thank you for your help.

void testMatlab()
{
    // Size of the three dimensions of the input matrix
    int nX = 1;
    int nY = 300;
    int nZ = 70;

    // Transformed parameters
    double *tfNewRe = (double*)malloc(nX*nY*nZ*sizeof(double));
    double *tfNewIm = (double*)malloc(nX*nY*nZ*sizeof(double));
    MKL_Complex16 *tfNewComp = (MKL_Complex16*)malloc(nX*nY*nZ*sizeof(MKL_Complex16));

    // Parameter to be retrieved by the backward transform
    //double *data = (double*)malloc(nX*nY*nZ*sizeof(double));
    MKL_Complex16 *data = (MKL_Complex16*)malloc(nX*nY*nZ*sizeof(MKL_Complex16));

    // Misc parameters
    char line[128];
    int indice = 0;
    double value = 0;
    MKL_LONG status;

    // Reading the real input data
    FILE *file;
    file = fopen("TFNewRe.txt", "r");
    if (file == NULL)
    {
        printf("Error with TFNewRe.txt\n");
        exit(1);
    }

    while (fscanf(file, "%lf", &value) == 1 && indice < nX*nY*nZ) 
    {
        tfNewRe[indice] = value; 
        indice++;
    }

    indice = 0;

    fclose(file);

    // Reading the imaginary input data
    file = fopen("tfNewIm.txt", "r");
    if (file == NULL)
    {
        printf("Error with tfNewIm.txt\n");
        exit(1);
    }

    while (fscanf(file, "%lf", &value) == 1 && indice < nX*nY*nZ) 
    {
        tfNewIm[indice] = value; 
        indice++;
    }

    // Seting the input data into a complex array
    for(int i = 0; i < nX*nY*nZ; i++)
    {
        tfNewComp[i].real = tfNewRe[i];
        tfNewComp[i].imag = tfNewIm[i];
    }

    // Create a descriptor handle
    DFTI_DESCRIPTOR_HANDLE desc;

    // Create the descriptor
    MKL_LONG dim[2] = {nY, nZ};
    status = DftiCreateDescriptor(&desc, DFTI_DOUBLE, DFTI_COMPLEX, (MKL_LONG)2, dim);

    // Set the distance between two lines
    DftiSetValue(desc, DFTI_INPUT_DISTANCE, nZ);
    DftiSetValue(desc, DFTI_OUTPUT_DISTANCE, nZ);

    // Set the strides
    /* MKL_LONG inputStrides[3] = {0, 2*(nZ/2+1), 1};
    MKL_LONG outputStrides[3] = {0, nZ/2+1, 1}; */
    MKL_LONG outputStrides[3] = {0, nZ, 1};
    MKL_LONG inputStrides[3] = {0, nZ/2+1, 1};
    DftiSetValue(desc, DFTI_INPUT_STRIDES, inputStrides);
    DftiSetValue(desc, DFTI_OUTPUT_STRIDES, outputStrides);

    // Initialize the placement to not in place
    status += DftiSetValue(desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE);

    status += DftiSetValue(desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);

    // Set the number of FFT to compute
    status += DftiSetValue(desc, DFTI_NUMBER_OF_TRANSFORMS, 1);

    // Set the normalization on the backward transform as Matlab does
    status += DftiSetValue(desc, DFTI_BACKWARD_SCALE, nY*nZ);

    // Commit the descriptor
    status += DftiCommitDescriptor(desc);

    if(status != 0)
    {
        printf("Erreur dans la configuration\n");
    }

    // Perform the backward Fourier transform
    status = DftiComputeBackward(desc, tfNewComp, data);

    // Free the descriptor
    DftiFreeDescriptor(&desc);

    // Write the output
    file = fopen("sortieC.txt", "w");
    if (file == NULL)
    {
        printf("Erreur lors de l'ouverture du fichier sortieC.txt\n");
        exit(1);
    }

    for (indice = 0; indice < nX*nY*nZ; indice++)
    {
        fprintf(file, "%f\n", data[indice].real);
    }

    fclose(file);

    // Free the memory allocation
    free(tfNewRe);
    free(tfNewIm);
    free(tfNewComp);
    free(data);
}

[AryanK_Intel]

Hi Côme Saint-André,

Thanks for posting in Intel Communities.

Could you please provide us the dependency files like "TFNewRe.txt", "tfNewIm.txt", "sortieC.txt", so that we can reproduce the issue and investigate further?

Best regards,

Sri Raj Aryan.K

[CSA_O]

Here attached TFNewRe.txt file with the real parts of the transform and TFNewIm.txt file with the imaginary parts.
sortieC.txt is where I try to write the result of the backward transform.

The 2 files are way bigger than what I need for this simple case as my final goal is to perform a 300X300X70 backward FFT.

Thanks,

[AryanK_Intel]

Hi Côme Saint-André,

Thanks for the details. We are looking into your issue internally and we will update you soon and also, Could you please refer to the examples on dft in the below path:

For Windows: C:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\examples\examples_core_c.zip\c\dft

For Linux: /opt/intel/oneapi/mkl/latest/examples/examples_core_c.tgz/c/dft

Thanks and Regards,

Sri Raj Aryan.K

[AryanK_Intel]

Hi Côme Saint-André,

We would like to know what exactly you are trying to do. We feel the following part of the code needs to be relooked.

// Create the descriptor

  MKL_LONG dim[2] = { nY, nZ};

  status = DftiCreateDescriptor(&desc, DFTI_DOUBLE, DFTI_COMPLEX, (MKL_LONG)2, dim);

  // Set the distance between two lines

  DftiSetValue(desc, DFTI_INPUT_DISTANCE, nZ);

  DftiSetValue(desc, DFTI_OUTPUT_DISTANCE, nZ);

From the above code given by you, DFTI_INPUT_DISTANCE and DFTI_OUTPUT_DISTANCE are the same but it should be the distance between successive 2D transforms. So the number needs to be at least nY*nZ. 

 // Create the descriptor

  MKL_LONG dim[2] = { nY, nZ};

  status = DftiCreateDescriptor(&desc, DFTI_DOUBLE, DFTI_COMPLEX, (MKL_LONG)2, dim);

  ...

  // Set the strides

  MKL_LONG outputStrides[3] = {0, nZ, 1};

  MKL_LONG inputStrides[3] = {0, nZ/2+1, 1};

  DftiSetValue(desc, DFTI_INPUT_STRIDES, inputStrides);

  DftiSetValue(desc, DFTI_OUTPUT_STRIDES, outputStrides);


  ...


  status += DftiSetValue(desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);

We need some clarification on these n/2+1 style inputStrides and setting of DFTI_CONJUGATE_EVEN_STORAGE, 

both of which would indicate the intention to compute a real backward transform, but the descriptor is set up to do complex transforms. 

Could you please let us know whether this helps and if our understanding is correct?

Thanks and Regards,

Sri Raj Aryan.K

[CSA_O]

Hi,

Thank you for getting back to me.

I'm trying to translate a Matlab code into C. There is a use of the Matlab ifftn function in the code and I want to use dftiComputeBackward to get a similar output. Only the real part of the result is used but I was trying to get there step by step to make sure there is no mistake in what I do, that's why I was going for a complex to complex transform.

On that simple step I am looking to compute the backward transform of the first 2D matrix even though in the end I should compute the nX transforms of my 3D input.

You are right about the distance, in my case it should have been nY*nZ.

From what I understood from the basic_dp_complex_dft_2d.c example I don't need to set the stride in my case.

There was also a mistake with the scale I had written as I had written nY*nZ instead of 1/(nY*nZ).

Here is how I modified the code :

 // Create a descriptor handle
    DFTI_DESCRIPTOR_HANDLE desc;

    // Create the descriptor
    MKL_LONG dim[2] = {nY, nZ};
    status = DftiCreateDescriptor(&desc, DFTI_DOUBLE, DFTI_COMPLEX, (MKL_LONG)2, dim);

    // Set the distance between two 2D FFT
    DftiSetValue(desc, DFTI_INPUT_DISTANCE, nY*nZ);
    DftiSetValue(desc, DFTI_OUTPUT_DISTANCE, nY*nZ);

    // Initialize the placement to not in place
    status += DftiSetValue(desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE);

    status += DftiSetValue(desc, DFTI_COMPLEX_STORAGE, DFTI_COMPLEX_COMPLEX);

    // Set the number of FFT to compute
    status += DftiSetValue(desc, DFTI_NUMBER_OF_TRANSFORMS, 1);

    // Set the normalization on the backward transform as Matlab does
    status += DftiSetValue(desc, DFTI_BACKWARD_SCALE, 1/(nY*nZ));

    // Commit the descriptor
    status += DftiCommitDescriptor(desc);

    if(status != 0)
    {
        printf("Erreur dans la configuration\n");
    }

    // Perform the backward Fourier transform
    status = DftiComputeBackward(desc, tfNewComp, data);

    // Free the descriptor
    DftiFreeDescriptor(&desc);

    // Write the output
    file = fopen("sortieC.txt", "w");
    if (file == NULL)
    {
        printf("Error with file sortieC.txt\n");
        exit(1);
    }

    for (indice = 0; indice < nX*nY*nZ; indice++)
    {
        fprintf(file, "%f\n", data[indice].real);
    }

    fclose(file);

But it still looks like my understanding is not right as my output file is filled with 0.

[AryanK_Intel]

Hi Côme Saint-André,

We are looking into your issue internally, we will get back to you soon with an update.

Thanks and Regards,

Sri Raj Aryan.K

[AryanK_Intel · Accepted solution]

Hi Côme Saint-André,

Could you please try changing your code?

From this line of code

status += DftiSetValue(desc, DFTI_BACKWARD_SCALE, 1/(nY*nZ));

To

status += DftiSetValue(desc, DFTI_BACKWARD_SCALE, 1.0/(nY*nZ)); 

As the scale value is set to zero by integer division, all the outputs are zero. So, we changed the line of code to fix this.

And we observed that there was a change in the output and also attaching the output log file. 

We also request you to refer to the below documentation, since it advises you to use a value with the same precision as the descriptor, but an integer input will be implicitly converted to double and still function.

https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/dfti-forward-scale-dfti-backward-scale.html#GUID-2F357ED2-F4DE-4E1D-AA21-51B8E3E2F7C4

And, could you please confirm whether it is the expected output? If not, could you please provide the details and issues you are facing?

Thanks and Regards,

Aryan 

[CSA_O]

Hello,

I've tested with the correction and I get the result I was looking for.

Thank you for your help.

[AryanK_Intel]

Hi Côme Saint-André,

It’s great to know that the issue has been resolved, in case you run into any other issues, please feel free to create a new thread.

Have a Good Day!

Thanks & Regards,

Aryan