Monday, 15 April 2013

c - How to do a fftw3 MPI "transposed" 2D transform if possible at all? -



c - How to do a fftw3 MPI "transposed" 2D transform if possible at all? -

consider 2d transform of form l x m (column major setup), complex array src real array tgt. or , in fortranese,

complex(c_double_complex), pointer :: src(:,:) real(8), pointer :: tgt(:,:) .

corresponding pointers are

type(c_ptr) :: csrc,ctgt .

i allocate them in next manner:

! complex array first alloc_local = fftw_mpi_local_size_2d(m,l/2+1,mpi_comm_world,local_m,local_offset1) csrc = fftw_alloc_complex(alloc_local) phone call c_f_pointer(csrc, src, [l/2,local_m]) ! real array alloc_local = fftw_mpi_local_size_2d(2*(l/2+1),m, & mpi_comm_world,local_l,local_offset2) ctgt = fftw_alloc_real(alloc_local) phone call c_f_pointer(ctgt, tgt, [m,local_l])

now, plan created as:

! create c-->r transform 1 transposition left out plan = fftw_mpi_plan_dft_c2r_2d(m,l,src,tgt, mpi_comm_world, & ior(fftw_measure,fftw_mpi_transposed_out))

finally, transform performed as:

call fftw_mpi_execute_dft_c2r(plan, src, tgt)

however, prescription not work. lastly phone call causes segmentation fault. @ first, thought might have how allocate src , tgt arrays, playing different amount of memory allocated tgt did not give result. so, either doing silly, or not possible @ all.

edit : minimalistic compileable example

program trashingfftw use, intrinsic :: iso_c_binding utilize mpi implicit none include 'fftw3-mpi.f03' integer(c_intptr_t), parameter :: l = 256 integer(c_intptr_t), parameter :: m = 256 type(c_ptr) :: plan, ctgt, csrc complex(c_double_complex), pointer :: src(:,:) real(8), pointer :: tgt(:,:) integer(c_intptr_t) :: alloc_local, local_m, & & local_l,local_offset1,local_offset2 integer :: ierr,id phone call mpi_init(ierr) phone call mpi_comm_rank(mpi_comm_world,id,ierr) phone call fftw_mpi_init() alloc_local = fftw_mpi_local_size_2d(m,l/2+1, mpi_comm_world, & local_m, local_offset1) csrc = fftw_alloc_complex(alloc_local) phone call c_f_pointer(csrc, src, [l/2,local_m]) alloc_local = fftw_mpi_local_size_2d(2*(l/2+1),m, mpi_comm_world, & & local_l, local_offset2) ctgt = fftw_alloc_real(alloc_local) phone call c_f_pointer(ctgt, tgt, [m,local_l]) plan = fftw_mpi_plan_dft_c2r_2d(m,l,src,tgt, mpi_comm_world, & ior(fftw_measure, fftw_mpi_transposed_out)) phone call fftw_mpi_execute_dft_c2r(plan, src, tgt) phone call mpi_finalize(ierr) end programme trashingfftw

and reply is:

for mpi real transforms, there 2 allowed combinations of transpositions , directions:

real complex transform , fftw_mpi_transposed_out complex real transform , fftw_mpi_transposed_in

i have found while digging within fftw3 ver. 3.3.4 code, file "rdft2-problem.c", comment on line 120.

edit:

minimal compilable and working example:

program trashingfftw use, intrinsic :: iso_c_binding utilize mpi implicit none include 'fftw3-mpi.f03' integer(c_intptr_t), parameter :: l = 256 integer(c_intptr_t), parameter :: m = 256 type(c_ptr) :: plan, ctgt, csrc complex(c_double_complex), pointer :: src(:,:) real(8), pointer :: tgt(:,:) integer(c_intptr_t) :: alloc_local, local_m, & & local_l,local_offset1,local_offset2 integer :: ierr,id phone call mpi_init(ierr) phone call mpi_comm_rank(mpi_comm_world,id,ierr) phone call fftw_mpi_init() alloc_local = fftw_mpi_local_size_2d(l/2+1,m, mpi_comm_world, & local_l, local_offset1) print *, id, "alloc complex=",alloc_local, local_l csrc = fftw_alloc_complex(alloc_local) phone call c_f_pointer(csrc, src, [m,local_l]) !caveat: must partition real storage according complex layout, why ! using m , l/2+1 instead of m, 2*(l/2+1) done in original post alloc_local = fftw_mpi_local_size_2d(m,l/2+1, mpi_comm_world, & & local_m, local_offset2) print *, id, "alloc real=",alloc_local, local_m ! 2 reals per complex ctgt = fftw_alloc_real(2*alloc_local) ! first l relevant, rest dangling space (see fftw3 docs) phone call c_f_pointer(ctgt, tgt, [l,local_m]) plan = fftw_mpi_plan_dft_c2r_2d(m,l,src,tgt, mpi_comm_world, & ior(fftw_measure, fftw_mpi_transposed_in)) ! should non-null print *, 'plan:', plan src(3,2)=(1.,0) phone call fftw_mpi_execute_dft_c2r(plan, src, tgt) phone call mpi_finalize(ierr) end programme thrashingfftw

c fortran mpi fftw

No comments:

Post a Comment