[SOLVED] How can I create a field vector with FFTW3 in C++?

How can I create a field vector with FFTW3 in C++?

I am using FFTW3 with OpenMPI and I want a field vector on which I will be able to perform FFTW transforms. With the code it will be clearer:

#include <complex>
#include <fftw3-mpi.h>

int main(int argc, char **argv)
{

MPI_Init(&argc, &argv);
fftw_mpi_init();

ptrdiff_t N0=50, N1=10, N2=1;
ptrdiff_t alloc_local, local_n0, local_0_start;

alloc_local = fftw_mpi_local_size_3d(N0, N1, N2 / 2 + 1, MPI_COMM_WORLD, &local_n0,&local_0_start);

fftw_complex **kVect;
kVect = (fftw_complex **)malloc(3 * sizeof(fftw_complex *));
for (int dim = 0; dim < 3; dim++)
    kVect[dim] = fftw_alloc_complex(alloc_local);

ptrdiff_t i, j, k;

for (i = 0; i < local_n0; i++)
    for (j = 0; j < N1; j++)
        for (k = 0; k < N2; k++)
        {
            kVect[0][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][0] = 1.0 / (2.0 * M_PI);
            kVect[0][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][1] = 0.0;
            kVect[1][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][0] = 1.0 / (2.0 * M_PI);
            kVect[1][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][1] = 0.0;
            kVect[2][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][0] = 1.0 / (2.0 * M_PI);
            kVect[2][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][1] = 0.0;
        }

for (int dim = 0; dim < 3; dim++)
    fftw_free(kVect[dim]);

free(kVect);

MPI_Finalize();

return 0;
}

I thought with this code it should work. However, at the free step fftw_free(kVect[dim]); mpiexec send an error message double free or corruption.

I tried to comment the initialisation part with the 3 for loop and the code does not crash. I suppose the way I use malloc is not suitable considering FFTW3 but I have no idea how to allocate memory for this vector.

Solution

I found the solution of this problem. It comes from the definition of arrays by the FFTW. To make FFT, the complex and the real arrays must have the same dimension, however the complex one is double she size of the real. Therefore, to allocate memory to a real field, one uses fftw_alloc_real(2*alloc_local). And to access the field at position (i,j,k), one uses the equation (i * N1 + j) * (2 * (N2 / 2 + 1)) + k as mention in the example: https://www.fftw.org/fftw3_doc/Multi_002ddimensional-MPI-DFTs-of-Real-Data.html. For the complex field, one should use k + j*N2 + i*N1*N2 instead, otherwise it will go beyond the boundaries of the field and trigger this error.