2D DFT через FFTW

Question

Я пытаюсь переместить 2D-DFT на FFTW, но у меня проблемы.

В моем старом коде я выполнил 1D C2C вдоль z, а затем 1D C2R вдоль направления x. Сейчас я пытаюсь использовать 2D C2R от FFTW, но мне чего-то не хватает, и результат другой ..

Вот что я пробовал:

#include <mpi.h>
#include <fftw3-mpi.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define N 4

int main (int argv, char **argc)
{
int nz=N*2+1, nx=N+1;
int nxd=8, nzd=12;
int rank, size;
double *real;   
fftw_complex *complx, *data;
double t1, t2;
fftw_plan plan;
ptrdiff_t alloc_local, local_n0, local_0_start, alloc_local_data, local_n0_data, local_0_start_data;

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

MPI_Comm_rank (MPI_COMM_WORLD, &rank);
MPI_Comm_size (MPI_COMM_WORLD, &size);

alloc_local_data = fftw_mpi_local_size_2d (nz, nx, 
                                      MPI_COMM_WORLD, 
                                      &local_n0_data,  
                                      &local_0_start_data);
alloc_local = fftw_mpi_local_size_2d (nzd, nxd, 
                                      MPI_COMM_WORLD, 
                                      &local_n0,  
                                      &local_0_start);
data = fftw_alloc_complex(alloc_local);
real = (double*) fftw_malloc(sizeof(double) * 2*alloc_local);
complx = fftw_alloc_complex(alloc_local);

plan = fftw_mpi_plan_dft_c2r_2d (nzd,(nxd-1)*2, complx, real, 
                             MPI_COMM_WORLD, 
                                 FFTW_MEASURE);
/*plan_f = fftw_mpi_plan_dft_r2c_2d (nzd, nxd, complx, real,
                             MPI_COMM_WORLD, 
                                 FFTW_MEASURE);*/
// realitialize realput with some numbers   
for (int i = 0; i < local_n0_data; i++) {
    for (int k = 0; k < nz; k++) {
        data[k+i*nz][0] = 2*(k+(i+local_0_start_data)*(nz));
        data[k+i*nz][1] = 2*(k+(i+local_0_start_data)*(nz))+1;
    //  printf("[%d]complex[%d]= %f+I*%f\n",rank,k+i*(nz), data[k+i*(nz)][0], data[k+i*(nz)][1]  );
    }
    //printf("-----------------\n");
}

/*Move data inside nxd=8xnzd=12 array*/
for(int i=0; i < 5; i++){
    for(int k=0; k < 5; k++) {
        complx[i*nzd+k][0]=data[(k+nz/2)+i*nz][0];
        complx[i*nzd+k][1]=data[(k+nz/2)+i*nz][1];
        //printf("[%d]complex[%d]= %f+I*%f\n",rank,k+i*(nzd), complx[k+i*(nzd)][0], complx[k+i*(nzd)][1]  );
    }   
}
for(int i=0; i < 5; i++){
    for(int k=5; k < 8; k++) {
        complx[i*nzd+k][0]=0;
        complx[i*nzd+k][1]=0;
        //printf("[%d]complex[%d]= %f+I*%f\n",rank,k+i*(nzd), complx[k+i*(nzd)][0], complx[k+i*(nzd)][1]  );
    }   
}
for(int i=0; i < 5; i++){
    for(int k=8; k < nzd; k++) {
        complx[i*nzd+k][0]=data[(k-nz+1)+i*nz][0];
        complx[i*nzd+k][1]=data[(k-nz+1)+i*nz][1];
        //printf("[%d]complex[%d]= %f+I*%f\n",rank,k+i*(nzd), complx[k+i*(nzd)][0], complx[k+i*(nzd)][1]  );
    }   
}
for(int i=5; i < nxd;i++){
    for(int k=0; k<nzd;k++){
        complx[i*nzd+k][0]=0;
        complx[i*nzd+k][1]=0;
    }
}
/*for(int i=0;i<nxd;i++){
    for(int k=0; k< nzd;k++){
        printf("complx[%d]=%f+I*%f\n",i*nzd+k,complx[i*nzd+k][0],complx[i*nzd+k][1]);
    }
    printf("----------------\n");
}*/

/*transform*/
fftw_execute(plan);

for(int i=0;i<2*(nxd-1);i++){
    for(int k=0; k< nzd;k++){
        printf("real[%d]=%f\n",i*nzd+k,real[i*nzd+k]);
    }
    printf("----------------\n");
}

// Clean up and get complx
//free(real);   free(complx);   fftw_free(data);
fftw_destroy_plan (plan);
MPI_Finalize ();
return 0;
}

Как видите, я использую FFTW с поддержкой MPI, но на данный момент я работаю в одноядерном режиме, поэтому не беспокойтесь о возможных ошибках привязки массива.

Рабочий массив 8x12 (nx, nz), который заполняется значениями, которые поступают из данных массива, то есть 5x9.

Я следовал парам примеров из руководства, здесь и здесь

Любой совет будет приветствоваться!