FFT solution 1

C

# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <fftw3.h>

int main ( void )

{
  ptrdiff_t i;
  const ptrdiff_t n = 1024;
  fftw_complex *in;
  fftw_complex *out;
  fftw_complex *newout;
  fftw_plan plan_backward;
  fftw_plan plan_forward;
  double cdiff, prec, diffr, diffi;
/* Create arrays. */
   in = fftw_malloc ( sizeof ( fftw_complex ) * n );
   out = fftw_malloc ( sizeof ( fftw_complex ) * n );
   newout = fftw_malloc ( sizeof ( fftw_complex ) * n );
/* Initialize data */
   for ( i = 0; i < n; i++ )
   { 
      if (i >= (n/4-1) && (3*n/4-1)) 
        {
          in[i][0] = 1.;
          in[i][1] = 0.;
        } 
      else
        {
          in[i][0] = 0.;
          in[i][1] = 0.;
        }       
   }
/* Create plans. */
   plan_forward = fftw_plan_dft_1d ( n, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
   plan_backward = fftw_plan_dft_1d ( n, out, newout, FFTW_BACKWARD, FFTW_ESTIMATE );
/* Compute transform (as many times as desired) */
   fftw_execute ( plan_forward );
/* Normalization */
   for ( i = 0; i < n; i++ ) { 
   out[i][0] = out[i][0]/n;
   out[i][1] = out[i][1]/n;
   }
/* Compute anti-transform */
   fftw_execute ( plan_backward );
/* Check results */
   prec = 1.0e-7;
   cdiff = 0.;
   for( i = 0; i < n; i++ ) {
      diffr = fabs(in[i][0] - newout[i][0]);
      diffi = fabs(in[i][1] - newout[i][1]);
      if ( cdiff < diffr) {
         cdiff = diffr;
      }
      if( cdiff < diffi) {
         cdiff = diffi;
      }
   }
   printf("max diff =%g\n",cdiff);
   if(cdiff > prec)
      printf("Results are incorrect\n");
   else 
      printf("Results are correct\n");

/* Print results */
   for ( i = 0; i < n; i++ )
   { 
    printf("i = %d, in = (%f), out = (%f,%f), newout = (%f)\n",i,in[i][0],out[i][0],out[i][1],newout[i][0]);
   }
/* deallocate and destroy plans */
   fftw_destroy_plan ( plan_forward );
   fftw_destroy_plan ( plan_backward );
   fftw_free ( in );
   fftw_free ( newout );
   fftw_free ( out );

  return 0;
}

FORTRAN

      program FFTW1D
        use, intrinsic :: iso_c_binding
        implicit none
        include 'fftw3.f03'
        integer(C_INTPTR_T):: L = 1024 
        integer(C_INT) :: LL
        type(C_PTR) :: plan1, plan2
        type(C_PTR) :: p_idata, p_odata, p_newdata
        real, parameter :: prec = 1.0e-7
        real :: cdiff, diff
        complex(C_DOUBLE_COMPLEX), dimension(:), pointer :: odata
        complex(C_DOUBLE), dimension(:), pointer :: idata, newdata
        integer :: i
        character(len=41), parameter :: filename='serial_data_optim.txt'
!! Allocate
           LL = int(L,C_INT)
           p_idata = fftw_alloc_complex(L)
           p_odata = fftw_alloc_complex(L)
           p_newdata = fftw_alloc_complex(L)
           call c_f_pointer(p_idata,idata,(/L/))
           call c_f_pointer(p_odata,odata,(/L/))
           call c_f_pointer(p_newdata,newdata,(/L/))
!!   create MPI plan for in-place forward DF
           plan1 = fftw_plan_dft_1d(LL, idata, odata, FFTW_FORWARD, FFTW_ESTIMATE)
           plan2 = fftw_plan_dft_1d(LL, odata, newdata, FFTW_BACKWARD,FFTW_ESTIMATE)
!! initialize data 
          do i = 1, L 
              if ( (i .ge. (L/4)) .and. (i .le. (3*L/4)) ) then
                  idata(i) = (1.,0.)
              else
                  idata(i) = (0.,0.)
              endif
          end do
!! compute transform (as many times as desired)
          call fftw_execute_dft(plan1, idata, odata)
!! Normalizzation
          odata = odata/L
!! Compute anti-transform
          call fftw_execute_dft(plan2, odata, newdata)
!! Check data
          cdiff =0.
          do i = 1, L
             diff = abs(idata(i) - newdata(i))
             if(cdiff .lt. diff) then
                cdiff = diff
             endif
          end do
          write(*,*) 'Cdiff_max =', cdiff
          if(cdiff .gt. (prec)) then
            write(*,*) 'Results are incorrect:'
          else
            write(*,*) 'Results are correct!'
          endif

!! print data
          OPEN(UNIT=45, FILE=filename, ACCESS='SEQUENTIAL')
              do i = 1, L
                   write(45,*) i, idata(i),odata(i),newdata(i)
              end do
          CLOSE(45)
!! deallocate and destroy plans
          call fftw_destroy_plan(plan1)
          call fftw_destroy_plan(plan2)
          call fftw_free(p_idata)
          call fftw_free(p_odata)
          call fftw_free(p_newdata)
      end

For compilation (on FERMI):

#!/bin/bash
# rm job.* MPI_Test

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

# fortran
bgxlf2003_r -O3 -qstrict -qarch=qp -qtune=qp -I$FFTW3_INCLUDE Es_1.f90 -L$FFTW3_LIB -lfftw3f -lfftw3 -lm -o Es-1

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

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
ifort -O3 -fpp -I. -I$FFTW_INC es1.f90 -L$FFTW_LIB -lfftw3 -o Es_1

# C
icc  -O3 -I. -I$FFTW_INC es1.c -L$FFTW_LIB -lfftw3 -o Es_1