Dynamic Allocated Matrix and FFT routine in MKL

[pietro_vinettigmail_]

Hi everybody,

I need to perform a 2D FFT on a dynamically allocated matrix of MKL_Complex8.
Let us assume that Dim is the variable containing the dimension of the square matrix I want to process, to handle the matrix DimxDim I use a double pointer structure such as the following code:

DFTI_DESCRIPOR* a;
int dimension[2];
dimension[0]=Dim;
dimension[1]=Dim;

MKL_Complex8** A,**B;
A=(MKL_Complex8*)malloc(Dim*sizeof(MKL_Complex8*));
B=(MKL_Complex8*)malloc(Dim*sizeof(MKL_Complex8*));

for (it k=0;k A=(MKL_Complex8*)malloc(Dim*sizeof(MKL_Complex8);
B=(MKL_Complex8*)malloc(Dim*sizeof(MKL_Complex8);
for (j=0;j A=...//Initialization of Akj element
B[]j=0;
}
}


Then I perform a 2D FFT with the following code:

DftiCreateDescriptor( &a, DFTI_SINGLE, DFTI_COMPLEX, 2,dimension );
status = DftiSetValue( a, DFTI_PLACEMENT, DFTI_NOT_INPLACE );
DftiCommitDescriptor( a );
status = DftiComputeForward ( a, A,B);

THe code is compiled correctly and no errors appear but when I check the B matrix after the processing I do not get the correct result and the strange thing is that B is not anymore a matrix (if I try to access to B

after FFT I get an error).

How can I manage dinamically allocated matrix with MKL FFT? Where is the issue in my approach?
Thanks in advance for any help.

Bye,

Pietro

[Dmitry_B_Intel]

Hi Pietro,

I would not recommend you double-subscript syntax for dynamically allocated arrays. This syntax is good for static arrays, where compiler knows the size of the array and can insert offset computation automatically for you. Instead, I would suggest you to compute the offsets explicitly, and optionally use preprocessor define to hide the offset computation and make the code more readable. An example:

MKL_Complex8 *A, *B; // not **A, **B
MKL_LONG dimension[2]; // MKL_LONG, not int!
A = malloc( sizeof(A[0]) * dimension[0] * dimension[1] );
B = /* similarly */;
#define A(i,j) A[(i)* dimension[1] + (j)]
#define B(i,j) B[(i)* dimension[1] + (j)]
for(...) { A(i,j) = ... ; }
DftiComitDescriptor(..., 2, dimension);
DftiComputeForward(..., A,B); // this will work, because both A and B point to actual data, not to the tables of pointers. When B is a table of pointers, FFT overwrites them and thus produces errors you have seen.

Thanks
Dima

[pietro_vinettigmail_]

With your suggestions, now everything works properly!

Many thanks Dima.

Bye,

Pietro