I am trying to make the 2D MPI FFTW example from http://www.fftw.org/doc/2d-MPI-example.html#g_t2d-MPI-example work.
The example code is not complete so I have to write some extra lines to make it compile and test (below is the code).
For some reason, the output (in place FFT) is just zeros as far as I can tell.
The output after transform is simply 0 0
(all zeros).
Am I using the MPI FFTW in the wrong way? The example code is simple enough.
// compile with: mpicc simple_mpi_example.c -Wl,-rpath=/usr/local/lib -lfftw3_mpi -lfftw3 -o simple_mpi_example */
#include <fftw3-mpi.h>
int main(int argc, char **argv){
const ptrdiff_t N0 = 1000, N1 = 1000;
fftw_plan plan;
fftw_complex *data; //local data of course
ptrdiff_t alloc_local, local_n0, local_0_start, i, j;
MPI_Init(&argc, &argv);
fftw_mpi_init();
/* get local data size and allocate */
alloc_local = fftw_mpi_local_size_2d(N0, N1, MPI_COMM_WORLD,
&local_n0, &local_0_start);
data = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * alloc_local);
MPI_Barrier(MPI_COMM_WORLD);
printf("%i %i\n", local_0_start, local_n0);
MPI_Barrier(MPI_COMM_WORLD);
/* create plan for forward DFT */
plan = fftw_mpi_plan_dft_2d(N0, N1, data, data, MPI_COMM_WORLD,
FFTW_FORWARD, FFTW_ESTIMATE);
/* initialize data to some function my_function(x,y) */
for (i = 0; i < local_n0; ++i) for (j = 0; j < N1; ++j){
data[i*N1 + j][0]=local_0_start + i;
data[i*N1 + j][1]=i;
}
MPI_Barrier(MPI_COMM_WORLD);
printf("%f %f\n", data[10*N1 + 10][0], data[10*N1 + 10][1]);
MPI_Barrier(MPI_COMM_WORLD);
/* compute transforms, in-place, as many times as desired */
fftw_execute(plan);
printf("%f %f\n", data[10*N1 + 10][0], data[10*N1 + 10][1]);
fftw_destroy_plan(plan);
fftw_free(data);
MPI_Finalize();
printf("finalize\n");
return 0;
}