Странное поведение при выполнении fftw3 MPI Преобразование Фурье от сложного к реальному

Question

Рассмотрим проблему, в которой вы хотите преобразовать спектральное представление f (x, y) ~ cos (x) в координатное пространство. Таким образом, exp (i * x) + exp (-i * x) ----> f (x, y), где f (x, y) = some_factor * cos (x).

В макете столбцов (мой пример написан ниже в Fortran 2003/8) спектральный массив будет инициализирован, как,

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

После преобразования необходимо получить TGT массива чисел (double), где все столбцы идентичны, и каждый столбец ведет себя как функция cos (x).

Однако, используя приведенный ниже код, я получаю массив tgt, где только нечетные столбцы ведут себя корректно, хотя и равны нулю.

Что, если что-то я делаю неправильно?

Моя версия fftw3 - 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 

Answer 1

Помните, что весь процесс выполняет код в mpi.

Так как вы пишете src(2,1)=(1.D0,0.D0), вы устанавливаете частоты size, по одному для каждого процесса.

Следовательно, ваш код работает, если используется один процесс mpiexec -np 1 c2rnot. Но если вы запустите много processus mpiexec -np 2 c2rnot, вы получите что-то еще ...

Не могли бы вы попробовать if(id==0) src(2,1)=(1.D0,0.D0)? id - это ранг, как в вашем коде.

Вы также можете использовать MPI_REAL8 в mpi_gather(): мой компилятор пожаловался на mpi_double, не имея типа IMPLICIT.

Удачи вам! Bye,