Working with FFTW 2.1.5 MPI

This is how you do it in FFTW 2.1.5 MPI for an array with real-space dimensions $[n_{x},n_{y},n_{z}]$, where nx is the ``fast'' index for a ``row-major'' array as described in Sec. 1.4.3 (we'll also assume n_fields=1).

1. Following the usual call to MPI_Init(), you call

    fplan = fftw3d_mpi_create_plan(mpi_comm,nx,ny,nz,
      FFTW_FORWARD,flags);
    if (output_order == FFTW_NORMAL_ORDER) {
      iplan = fftw3d_mpi_create_plan(mpi_comm,nx,ny,nz,
        FFTW_INVERSE,flags);
    } else {
      iplan = fftw3d_mpi_create_plan(mpi_comm,ny,nz,nx,
        FFTW_INVERSE,flags);
    }
   

2. You next call the routine

    fftwnd_mpi_local_sizes(fplan,&local_nz,&local_z_start,
      &local_ny_after_transpose,&local_y_start_after_transpose,
      &total_local_size);
   

   to get information on the part of the data that you will be working with. You must then use malloc() to create total_local_size pixels of complex data (we'll call this sub-array local_data), and an equal amount of workspace array data (we'll call this sub-array workspace).
3. You can then load your part of the data. Let's say you have a function f(i) which provides the initial value for each pixel in a 3D real-space array with i=ix+iy*nx+iz*nx*ny. You would load data into this array in real space with

    for (iz=0; iz<local_nz; iz++)
      for (iy=0; iy<ny; iy++)
        for (ix=0; ix<nx; ix++)
          i = ix+iy*nx+(iz+local_z_start)*nx*ny;
          c_re(local_data,i) = f(i);
   

4. You then do a transform on the data with

    fftwnd_mpi(plan,n_fields,local_data,workspace,output_order);
   

   You can then refer to Fourier space data in the output as

    if (output_order == FFTW_NORMAL_ORDER) {
      for (iz=0; iz<local_nz; iz++)
        for (iy=0; iy<ny; iy++)
          for (ix=0; ix<nx; ix++)
            i = ix+iy*nx+(iz+local_z_start)*nx*ny;
            this_real = c_re(local_data,i);
    } else {
      for (iy=0; iy<local_ny_after_transpose; iy++)
        for (iz=0; iz<nz; iz++)
          for (ix=0; ix<nx; ix++)
            i = ix+iz*nx+(iy+local_y_start_after_transpose)*nx*nz;
            this_real = c_re(local_data,i);
    }
   

5. At the end of the program, you should call fftwnd_mpi_destroy_plan() on fplan and iplan, and MPI_Finalize().