FFT solution 2

C
# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <mpi.h>
# include <fftw3-mpi.h>
     
     int main(int argc, char **argv)
     {
         const ptrdiff_t L = 1024, M = 1024;
         fftw_plan plan;
         fftw_complex *data ;
         ptrdiff_t alloc_local, local_L, local_L_start, i, j, ii;
         int proc_id, nproc;
         double xx, yy, rr, r2, t0, t1, t2, t3, tplan, texec;
         const double amp = 0.25;
         /* Initialize */
         MPI_Init(&argc, &argv);
         fftw_mpi_init();
         MPI_Comm_size(MPI_COMM_WORLD,&nproc);
         MPI_Comm_rank(MPI_COMM_WORLD,&proc_id);
     
         /* get local data size and allocate */
         alloc_local = fftw_mpi_local_size_2d(L, M, MPI_COMM_WORLD, &local_L, &local_L_start);
         data = fftw_alloc_complex(alloc_local);
         /* create plan for in-place forward DFT */
         t0 = MPI_Wtime();
         plan = fftw_mpi_plan_dft_2d(L, M, data, data, MPI_COMM_WORLD, FFTW_FORWARD, FFTW_ESTIMATE);
         t1 = MPI_Wtime();
         /* initialize data to some function my_function(x,y) */
         for (i = 0; i < local_L; ++i) for (j = 0; j < M; ++j)
         {
	    ii = i + local_L_start;
            xx = ( (double) ii - (double)L/2 )/(double)L;
            yy = ( (double)j - (double)M/2 )/(double)M;
            r2 = pow(xx, 2) + pow(yy, 2);
            rr = sqrt(r2);
            if (rr <= amp)
            {
              data[i*M + j][0] = 1.;   
              data[i*M + j][1] = 0.; 
            }
            else
            {
              data[i*local_L + j][0] = 0.;   
              data[i*local_L + j][1] = 1.; 
            }
         }
         /* compute transforms, in-place, as many times as desired */
         t2 = MPI_Wtime();
         fftw_execute(plan);
         t3 = MPI_Wtime();
         /* Print results */
         tplan = t1 - t0;
         texec = t2 - t1;
         if (proc_id == 0) printf(" T_plan = %f, T_exec = %f \n",tplan,texec);
         /* deallocate and destroy plans */ 
         fftw_destroy_plan(plan);
         fftw_mpi_cleanup();
         fftw_free ( data );
         MPI_Finalize();
     }

FORTRAN

        program FFT_MPI_3D
        use, intrinsic :: iso_c_binding
        implicit none
	include 'mpif.h'
        include 'fftw3-mpi.f03'
        integer(C_INTPTR_T), parameter :: L = 1024
        integer(C_INTPTR_T), parameter :: M = 1024
        type(C_PTR) :: plan, cdata
        complex(C_DOUBLE_COMPLEX), pointer :: fdata(:,:)
        real(C_DOUBLE) :: t1, t2, t3, t4, tplan, texec
        integer(C_INTPTR_T) :: alloc_local, local_M, local_j_offset
        integer(C_INTPTR_T) :: i, j
        complex(C_DOUBLE_COMPLEX) :: fout
        integer :: ierr, myid, nproc
!
! Initialize
	  call mpi_init(ierr)
	  call MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr)
	  call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
	  call fftw_mpi_init()
!
!   get local data size and allocate (note dimension reversal)
	  alloc_local = fftw_mpi_local_size_2d(M, L, &
     &                  MPI_COMM_WORLD, local_M, local_j_offset)
	  cdata = fftw_alloc_complex(alloc_local)
	  call c_f_pointer(cdata, fdata, [L,local_M])
!
!   create MPI plan for in-place forward DFT (note dimension reversal) 
          t1 = MPI_wtime()
	  plan = fftw_mpi_plan_dft_2d(M, L, fdata, fdata, &
     &         MPI_COMM_WORLD, FFTW_FORWARD, FFTW_MEASURE)
          t2 = MPI_wtime()
!
! initialize data to some function my_function(i,j) 
	  do j = 1, local_M
	    do i = 1, L
	      call initial(i, (j + local_j_offset), L, M, fout)
	      fdata(i, j) = fout
	    end do
	  end do
!
! compute transform (as many times as desired) 
          t3 = MPI_wtime()
	  call fftw_mpi_execute_dft(plan, fdata, fdata)
          t4 = MPI_wtime()
!
! print solutinos
          tplan = t2-t1
          texec =t4-t3
          if (myid.eq.0) print*,'Tplan=',tplan,'   Texec=',texec
!
! deallocate and destroy plans
	  call fftw_destroy_plan(plan)
	  call fftw_mpi_cleanup()
	  call fftw_free(cdata)
	  call mpi_finalize(ierr)
!
      end
!
! ***** Subroutines ****************************************************
!
      subroutine initial(i, j, L, M, fout)
        use, intrinsic :: iso_c_binding
	integer(C_INTPTR_T), intent(in) :: i, j, L, M
        complex(C_DOUBLE_COMPLEX), intent(out) :: fout
	real(C_DOUBLE), parameter :: amp = 0.25
	real(C_DOUBLE) :: xx, yy, LL, MM, r1

          xx = real(i, C_DOUBLE) - real((L+1)/2, C_DOUBLE)
	  yy = real(j, C_DOUBLE) - real((M+1)/2, C_DOUBLE)
	  LL =  real(L, C_DOUBLE)
	  MM = real(M, C_DOUBLE)
        
          r1 = sqrt(((xx/LL)**2.) + ((yy/MM)**2.))
	  if (r1 .le. amp) then
	    fout = CMPLX(1., 0. , C_DOUBLE_COMPLEX)
          else
            fout = CMPLX(0., 0. , C_DOUBLE_COMPLEX)
          endif

      return
      end
! **********************************************************************

For compilation (on FERMI):

module purge
module load autoload fftw/3.3.2--bgq-xl--1.0 bgq-xl/1.0

# fortran
mpixlf2003_r -O3 -qstrict -qarch=qp -qtune=qp -I$FFTW3_INCLUDE Es_2.f90 -L$FFTW3_LIB -lfftw3_mpi -lfftw3f -lfftw3 -lm -o Es_2

# C
mpixlc_r -O3 -qstrict -qarch=qp -qtune=qp -I$FFTW3_INCLUDE Es_2.c -L$FFTW3_LIB -lfftw3_mpi -lfftw3f -lfftw3 -lm -o Es_2

For compilation (on EURORA):

module purge
module load autoload profile/advanced
module load fftw/3.3.3--openmpi--1.6.4--intel--cs-xe-2013--binary

# fortran
mpif90 -O3 -fpp -I. -I$FFTW_INC es2.f90 -L$FFTW_LIB -lfftw3_mpi -lfftw3 -o Es_2

# C
mpicc -O3 -I. -I$FFTW_INC es2.c -L$FFTW_LIB -lfftw3_mpi -lfftw3 -o Es_2