FFTW advanced interface with MPI in C - Stack Overflow

admin管理员组

文章数量:1279011

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 ***

本文标签： FFTW advanced interface with MPI in CStack Overflow