FFTW (2-D MPI) (Complex to Real to Complex)
Program rupak
use, intrinsic :: iso_c_binding
implicit none
include 'mpif.h'
integer (C_INTPTR_T), parameter :: nx = 5
integer (C_INTPTR_T), parameter :: ny = 3
integer (C_INTPTR_T), parameter :: nh = ( nx / 2 ) + 1
include "fftw3-mpi.f03"
integer(C_INTPTR_T) :: i, j
real(C_DOUBLE), dimension(nx,ny) :: ux
complex(C_DOUBLE_COMPLEX), dimension(nh,ny) :: ukx,ukx_dum
integer :: ierr, myid, nproc
type(C_PTR) :: plan, plan1, cdatar, cdatac
integer(C_INTPTR_T) :: alloc_local, local_ny, local_j_offset
real(C_DOUBLE), pointer :: idata(:,:)
complex(C_DOUBLE_complex), pointer :: odata(:,:)
integer,parameter :: seed = 99999999
call srand(seed)
call mpi_init(ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
do i = 1, nh
do j = 1, ny
ukx(i,j) = (1.0d0,0.0d0)*rand() + (0.0d0,1.0d0)*rand()
end do
end do
call fftw_mpi_init()
! get local data size and allocate (note dimension reversal)
alloc_local = fftw_mpi_local_size_2d(ny, nh, MPI_COMM_WORLD, &
local_ny, local_j_offset)
cdatar = fftw_alloc_real(2 * alloc_local)
cdatac = fftw_alloc_complex(alloc_local)
call c_f_pointer(cdatar, idata, [2*nh,local_ny])
call c_f_pointer(cdatac, odata, [nh,local_ny])
! create MPI plan for out-of-place DFT (note dimension reversal)
plan = fftw_mpi_plan_dft_c2r_2d(ny, nx, odata, idata, MPI_COMM_WORLD, &
FFTW_MEASURE)
plan1 = fftw_mpi_plan_dft_r2c_2d(ny, nx, idata, odata, MPI_COMM_WORLD, &
FFTW_MEASURE)
! initialize data to some function my_function(i,j)
do i = 1, nh
do j = 1, local_ny
odata(i, j) = ukx(i, j + local_j_offset)
!write ( *, * ) i,j+local_j_offset,odata(i,j),local_j_offset!,myid
end do
end do
! compute transform (as many times as desired)
call fftw_mpi_execute_dft_c2r(plan, odata, idata)
! COPY TO OUTPUT ARRAY
do i = 1, nx
do j = 1, local_ny
idata(i,j) = idata(i,j)/(float(Nx)*float(Ny))
ux(i,j+local_j_offset) = idata(i,j)
write ( *, * ) i, j+local_j_offset, idata(i,j)!,ux(i,j+local_j_offset)!, local_j_offset
end do
end do
call fftw_mpi_execute_dft_r2c(plan1, idata, odata)
! COPY TO OUTPUT ARRAY
do i = 1, nh
do j = 1, local_ny
ukx_dum(i,j+local_j_offset) = odata(i,j)
!write ( *, * ) i, j+local_j_offset, odata(i,j)!,ukx_dum(i,j+local_j_offset)!, local_j_offset
end do
end do
! deallocate and destroy plans
call fftw_destroy_plan(plan)
call fftw_mpi_cleanup()
call fftw_free(cdatar)
call fftw_free(cdatac)
do i = 1, nx
do j = 1, ny
!write ( *, '(2x,i4,2x,i4,2x,2g14.6)' ) i, j, ux(i,j)
end do
end do
do i = 1, nh
do j = 1, ny
!write ( *, '(2x,i4,2x,i4,2x,2g14.6)' ) i, j, ukx_dum(i,j)
end do
end do
call mpi_finalize(ierr)
end