Strange behaviour when doing a fftw3 MPI Fourier transform from complex to real

Consider a problem where you want to transform a spectral representation of f(x,y) ~ cos(x) to coordinate space. So exp(i*x) + exp(-i*x) ----> f(x,y), where f(x,y) = some_factor * cos(x) .

In column-major layout (my example is written below in Fortran 2003/8) the spectral array would be initialised as,

src=(0.,0.)
src(2,1)=(1.,0.)  .

After the transform one should get a tgt array of reals (doubles), where all columns are identical and each column behaves as a cos(x) function.

However , using the code given below, I get a tgt array where only the odd columns behave correctly, while even are all zero.

What, if anything am I doing wrong?

My fftw3 version is 3.4.4 .

    program c2r
  ! COMPLEX --> REAL 2D MPI TRANSFORM

  ! compiled using: mpif90 -O0 -g -I/usr/include c2rnot.F90 -o c2rnot -L /usr/lib/x86_64-linux-gnu -lfftw3_mpi -lfftw3
  use, intrinsic :: iso_c_binding
  use MPI

  implicit none
  include 'fftw3-mpi.f03'
  integer :: ierr,id

  ! logical dimensions of the transform
  integer(C_INTPTR_T), parameter :: N1 = 32
  integer(C_INTPTR_T), parameter :: N2 = 32

  ! helper variables: storage size, extents along the slow dim, offsets
  integer(C_INTPTR_T) :: alloc_local, slicec,slicer,offc,offr

  ! pointers to planner, real (image) and complex (original) arrays
  type(C_PTR) :: plan, ctgt, csrc

  ! fortran array representation of csrc ...
  complex(8), pointer :: src(:,:)
  ! ..., ctgt and a total array to hold the mpi-gathered result
  real(8), pointer :: tgt(:,:),total(:,:)


  integer(C_INTPTR_T) :: i,y0

  ! init mpi and fftw ...
  call mpi_init(ierr)

  call mpi_comm_rank(MPI_COMM_WORLD,id,ierr)

  call fftw_mpi_init()
  ! done with init


  ! calculate local complex array storage requirements and layout
  alloc_local = fftw_mpi_local_size_2d(N2,N1/2+1, &
        & MPI_COMM_WORLD,slicec,offc)
  ! allocate local complex array
  csrc = fftw_alloc_complex(alloc_local)
  ! fortran representation
  call c_f_pointer(csrc, src, [N1/2+1,slicec])

  ! allocate local real array
  ctgt = fftw_alloc_real(2*alloc_local)
  ! for non-transposed arrays slices are the same
  slicer=slicec
  offr=offc
  ! fortran representation
  call c_f_pointer(ctgt, tgt, [2*(N1/2+1),slicer])

  ! allocate global container (holds gathered tgt arrays)
  if (id==0) allocate(total(N1,N2))

  ! N1 - fast, N2 - slow dimension, in: src, out: tgt
  plan =  fftw_mpi_plan_dft_c2r_2d(N2,N1,src,tgt, MPI_COMM_WORLD, FFTW_MEASURE)


  src=(0.D0,0.D0)

  ! ~ exp(i*x) + exp(-i*x) which should turn into ~ cos(x)
  src(2,1)=(1.D0,0.D0)

  tgt=0.D0
  ! do transform
  call fftw_mpi_execute_dft_c2r(plan,src,tgt)

  ! collect tgt-s into total in process 0
  call mpi_gather(tgt(1:n1,1:slicer),slicer*n1,mpi_double,total,slicer*n1,mpi_double,0,mpi_comm_world,ierr)

  ! The scalar field total has x: 1..N1, y:1..N2 layout. Print
  ! behaviour along the x-axis for a given, constant y. The result
  ! should be ~ cos(x) function for any y0 in [1..N2]. However, I get
  ! such result only for odd y0. For even, total(i,y0) is 0.
  y0=3
  if (id == 0) then
     do i=1,n1
        print *, i, total(i,y0)
     enddo
  endif

  call mpi_finalize(ierr)
end program c2r

Answers


Remember that all processus execute the code in mpi.

So as you write src(2,1)=(1.D0,0.D0) you are setting size frequencies, one for each process.

Therefore, your code seems to work if a single process is used mpiexec -np 1 c2rnot. But if you run many processus mpiexec -np 2 c2rnot, you get something else...

Could you try if(id==0) src(2,1)=(1.D0,0.D0) ? id is the rank, as in your code.

You may also use MPI_REAL8 in mpi_gather() : my compiler complained about mpi_double having no IMPLICIT type.

Good luck ! Bye,


Need Your Help

How to average all the frames of a video file in which objects are not moving in OpenCV?

c++ c opencv

I have different frames of a video file. Now, on observing each frames separately, I noticed there are many frames in which objects has not moved. I need to do averaging of all those frames and mak...

how to train the stanford LexicalizedParser to recognize new words as nouns?

parsing nlp stanford-nlp

I am trying to figure out how to train the stanford LexicalizedParser

About UNIX Resources Network

Original, collect and organize Developers related documents, information and materials, contains jQuery, Html, CSS, MySQL, .NET, ASP.NET, SQL, objective-c, iPhone, Ruby on Rails, C, SQL Server, Ruby, Arrays, Regex, ASP.NET MVC, WPF, XML, Ajax, DataBase, and so on.