Program rupak
use, intrinsic :: iso_c_binding
implicit none
include 'mpif.h'
integer (C_INTPTR_T), parameter :: nx = 7
integer (C_INTPTR_T), parameter :: ny = 5
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) :: in,in2
complex(C_DOUBLE_COMPLEX), dimension(nh,ny) :: out
integer :: ierr, myid, nproc, seed
type(C_PTR) :: plan, plan1, cdatar, cdatac
integer(C_INTPTR_T) :: alloc_local, local_ny, local_j_offset
real(C_DOUBLE), pointer :: rdata(:,:)
complex(C_DOUBLE_complex), pointer :: cdata(:,:)
seed = 123456789
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, nx
do j = 1, ny
in(i,j) = ran(seed)
seed=in(i,j)*seed
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_ny_offset)
cdatar = fftw_alloc_real(2 * alloc_local)
call c_f_pointer(cdatar, rdata, [2*nh,local_ny])
alloc_local = fftw_mpi_local_size_2d(nh,ny, MPI_COMM_WORLD, &
local_nx, local_nx_offset)
cdatac = fftw_alloc_complex(alloc_local)
call c_f_pointer(cdatac, cdata, [ny,local_nx])
! create MPI plan for out-of-place DFT (note dimension reversal)
plan1 = fftw_mpi_plan_dft_r2c_2d(ny,nx,rdata,cdata, MPI_COMM_WORLD, &
ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT))
plan2 = fftw_mpi_plan_dft_c2r_2d(ny,nx,cdata,rdata2, MPI_COMM_WORLD, &
ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN))
! initialize data to some function my_function(i,j)
do i = 1, nx
do j = 1, local_ny
rdata(i, j) = in(i, j + local_j_offset)
write ( *, '(2x,i4,2x,i4,2x,2g14.6)' ) i,j+local_j_offset,rdata(i,j),local_j_offset!,myid
end do
end do
! compute transform (as many times as desired)
call fftw_mpi_execute_dft_r2c(plan, rdata, cdata)
! COPY TO OUTPUT ARRAY
do i = 1, nh
do j = 1, local_ny
out(i,j+local_j_offset) = cdata(i,j)
write ( *, '(2x,i4,2x,i4,2x,2g14.6)' ) i, j+local_j_offset, cdata(i,j)!,out(i,j+local_j_offset)!, local_j_offset
end do
end do
call fftw_mpi_execute_dft_c2r(plan1, cdata, rdata)
! COPY TO OUTPUT ARRAY
do i = 1, nx
do j = 1, local_ny
in2(i,j+local_j_offset) = rdata(i,j)
write ( *, '(2x,i4,2x,i4,2x,2g14.6)' ) i, j+local_j_offset, rdata(i,j)/(ny*nx)!,in2(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, nh
do j = 1, ny
!write ( *, '(2x,i4,2x,i4,2x,2g14.6)' ) i, j, out(i,j)
end do
end do
do i = 1, nx
do j = 1, ny
!write ( *, '(2x,i4,2x,i4,2x,2g14.6)' ) i, j, in2(i,j)/(ny*nx)
end do
end do
call mpi_finalize(ierr)
end