I am using the advance interface of FFTW combined with MPI in order to perform several 1D forward and backward transforms at a time. When the transform is complex to complex (both input and output have size N) my code works fine, but when I change the transform to be real to complex, and vice-versa, my code does not work anymore. Of course I adapted the plans and execute functions.
I suspect that the problem comes from the number of elements allocated for the complex array that should be roughly half of the elements allocated for the real array. Despite I tried N/2 and N/2+1 this still does not work. Below you can find my code and error message. By the way, isn't that weird that fftw_execute_r2r appears in the error message while I am using fftw_execute_r2c ?
Does someone knows what I am doing wrong ?
Thank you in advance
Code:
#include <stdio.h> #include <stdlib.h> #include <mpi.h> #include <unistd.h> #include <complex.h> #include <fftw3-mpi.h> int main(int argc, char *argv[]) { /* --- MPI Init --------------------------------------------------- */ int rank, size; int Nproc_x; int Nx, Ny; int nx = 3; MPI_Init(&argc, &argv); MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &Nproc_x); Nx = Nproc_x*nx; Ny = 40; /* --- FFTW INITIALISATION ----------------------------------------------- */ int i, N; int fftw_rank; ptrdiff_t fftw_dims[1]; ptrdiff_t howmany, block0; /* local data size */ ptrdiff_t fftw_alloc_local; ptrdiff_t fftw_local_nx; ptrdiff_t fftw_local_x_start; /* Complex arrays */ double * rho_real; fftw_complex * rho_cmplx; /* Plans for FFTW transforms */ static fftw_plan r2c_plan; static fftw_plan c2r_plan; /* Other variables */ static double FFT_norm; fftw_mpi_init(); fftw_rank = 1; fftw_dims[0] = Nx; howmany = Ny; block0 = nx; /* Array size in radial dir on current processor */ FFT_norm = (double) 1/Nx; /* get local data size and allocate data array */ fftw_alloc_local = fftw_mpi_local_size_many(fftw_rank, fftw_dims, howmany, block0, MPI_COMM_WORLD, /* Transform in azimuth */ &fftw_local_nx, &fftw_local_x_start); N = fftw_alloc_local; rho_real = fftw_alloc_real(N); rho_cmplx = fftw_alloc_complex(N/2+1); /* plan for forward transformation: real density to complex density */ r2c_plan = fftw_mpi_plan_many_dft_r2c(fftw_rank, fftw_dims, howmany, block0, FFTW_MPI_DEFAULT_BLOCK, rho_real, rho_cmplx, MPI_COMM_WORLD, FFTW_ESTIMATE); c2r_plan = fftw_mpi_plan_many_dft_c2r(fftw_rank, fftw_dims, howmany, FFTW_MPI_DEFAULT_BLOCK, block0, rho_cmplx, rho_real, MPI_COMM_WORLD, FFTW_ESTIMATE); /* Populate real array */ for (int i = 0; i < fftw_local_nx; ++i) { for (int j = 0; j < Ny; ++j) { rho_real[i * Ny + j] = i * Ny + j ; } } fftw_mpi_execute_dft_r2c(r2c_plan, rho_real, rho_cmplx); fftw_mpi_execute_dft_c2r(c2r_plan, rho_cmplx, rho_real); printf("INFO(gravBessel): fftw execute success ! \n"); MPI_Finalize(); return 0; } Error message:
[euler:12536] *** Process received signal *** [euler:12536] Signal: Segmentation fault (11) [euler:12536] Signal code: Address not mapped (1) [euler:12536] Failing at address: (nil) [euler:12536] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x43090)[0x7f2f7ff04090] [euler:12536] [ 1] /lib/x86_64-linux-gnu/libfftw3.so.3(fftw_execute_r2r+0x4)[0x7f2f802e59d4] [euler:12536] [ 2] my_program(+0x1509)[0x55da6283a509] [euler:12536] [ 3] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf3)[0x7f2f7fee5083] [euler:12536] [ 4] my_program(+0x122e)[0x55da6283a22e] [euler:12536] *** End of error message *** 
(nil). I suggest to build a debug version of your program to see where exactly the error occurs. Running it withvalgrindor build a version with the compiler's address sanitizer might find more problems. See airbus-seclab.github.io/c-compiler-security/#gcc-tldr Simply trying different array sizes without understanding is not a good approach.NULLwhen the rank is equal to 1, which is precisely my case.