Неверный вывод через fftw_mpi_r2c_2d и fftw_mpi_c2r_2d

Question

Я написал простую тестовую программу для реализации FFTW с MPI в домене 2d (с помощью Fortran). Область «Ny x Nx» широта и разделена во втором («x») индексе.

После правильного (я верю?) Объявления и распределения переменных и планов я вызываю функцию fftw_mpi r2c_2d, а затем я преобразовываю его вывод с помощью fftw_mpi c2r_2d, чтобы проверить, не получаю ли я исходный ввод. Кажется, что часть r2c_2d работает нормально. Тем не менее, я не получаю исходный ввод после преобразования назад (отдельно нормализации) с помощью функции c2r_2d: результирующий вектор отображает «нули» по индексам (:, j) с j, соответствующими кратным «Ny/2», , Что я делаю не так? Благодаря!

Вот выдержка из кода:

Program TEST 

use, intrinsic :: iso_c_binding 

Implicit none 

include 'mpif.h' 
include 'fftw3-mpi.f03' 

Integer*8,parameter :: nx=16, ny=16 

!MPI 
integer*8 :: ipe,npe 
integer*8 ::mpi_realtype,icomm=mpi_comm_world,istat(mpi_status_size),ierr 

! FFTW VARIABLES DECLARATION 
type(C_PTR)   :: p1, p2, cdatar, cdatac 
integer(C_INTPTR_T) :: alloc_local, local_L, local_L_offset, local_M, local_M_offset 
real(C_DOUBLE), pointer :: faux(:,:) ! real input 2d function 
complex(C_DOUBLE), pointer :: gaux(:,:) ! complex output of 2d FFTW (transposed) 


! MPI initialization 
call mpi_init(ierr) 

call mpi_comm_rank(icomm,ipe,ierr) 
call mpi_comm_size(icomm,npe,ierr) 


! FFTW ALLOCATIONS AND PLANS 

call fftw_mpi_init() 


alloc_local = fftw_mpi_local_size_2d(ny/2+1,nx & 
    ,MPI_COMM_WORLD, local_L, local_L_offset) 


cdatac = fftw_alloc_complex(alloc_local) 

call c_f_pointer(cdatac, gaux, [nx,local_L]) !transposed 

alloc_local = fftw_mpi_local_size_2d(nx,ny/2+1, MPI_COMM_WORLD, & 
    local_M, local_M_offset) 


cdatar = fftw_alloc_real(2*alloc_local) 

call c_f_pointer(cdatar, faux, [ny,local_M]) 

! Create plans 

p1 = fftw_mpi_plan_dft_r2c_2d(nx,ny,faux,gaux, MPI_COMM_WORLD, & 
     ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT)) 


p2 = fftw_mpi_plan_dft_c2r_2d(nx,ny,gaux,faux, MPI_COMM_WORLD, & 
     ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN)) 

! EXECUTE FFTW 

call random_number(faux) 


print *, "real input:", real(faux(1,:)) 

call fftw_mpi_execute_dft_r2c(p1,faux,gaux) 

call fftw_mpi_execute_dft_c2r(p2, gaux, faux) 

print *, "real output:", real(faux(1,:))/(nx*ny) 

call fftw_destroy_plan(p1) 
call fftw_destroy_plan(p2) 


call mpi_finalize(ierr) 


End Program TEST 

Answer 1

Проблема обусловлена ​​padding необходимой FFTW:

Хотя реальные данные концептуально n0 × n1 × n2 × ... × nd-1, он физически хранится как массив n0 × n1 × n2 × ... × [2 (nd-1/2 + 1)], где последнее измерение было дополнено, чтобы сделать его того же размера, что и сложный вывод. Это очень похоже на интерфейс последовательного интерфейса r2c/c2r на месте (см. Многомерные DFT-данные реальных данных), за исключением того, что в MPI отступы требуются даже для данных вне места.

Следовательно, входной массив для преобразования 16x16 является, следовательно, массивом 16x18. Значение дополнительных двух чисел в конце каждого ряда не имеет смысла в реальном пространстве. Тем не менее, эти дополнительные цифры не следует забывать, как указатель с отливают фортрановский 2D массива:

call c_f_pointer(cdatar, faux, [2*(ny/2+1),local_M]) 

Лишние цифры по-прежнему печатаются в конце каждой строки. Массив может быть нарезано, чтобы избежать печати этих бесполезных значения:

print *, "real input:", real(faux(1:ny,:)) 
... 
print *, "real output:", real(faux(1:ny,:))/(nx*ny) 

Вот полный код, основанный на вашем и один из How to do a fftw3 MPI "transposed" 2D transform if possible at all? Он может быть скомпилирован с помощью mpif90 main.f90 -o main -I/usr/include -L/usr/lib -lfftw3_mpi -lfftw3 -lm и выбежала на mpirun -np 2 main.

Program TEST 

use, intrinsic :: iso_c_binding 

Implicit none 

include 'mpif.h' 
include 'fftw3-mpi.f03' 

Integer*8,parameter :: nx=4, ny=8 

!MPI 
integer*8 :: ipe,npe 
integer*8 ::mpi_realtype,icomm=mpi_comm_world,istat(mpi_status_size),ierr 

! FFTW VARIABLES DECLARATION 
type(C_PTR)   :: p1, p2, cdatar, cdatac 
integer(C_INTPTR_T) :: alloc_local, local_L, local_L_offset, local_M, local_M_offset 
real(C_DOUBLE), pointer :: faux(:,:) ! real input 2d function 
complex(C_DOUBLE), pointer :: gaux(:,:) ! complex output of 2d FFTW (transposed) 


! MPI initialization 
call mpi_init(ierr) 

call mpi_comm_rank(icomm,ipe,ierr) 
call mpi_comm_size(icomm,npe,ierr) 


! FFTW ALLOCATIONS AND PLANS 

call fftw_mpi_init() 


alloc_local = fftw_mpi_local_size_2d(ny/2+1,nx & 
    ,MPI_COMM_WORLD, local_L, local_L_offset) 


cdatac = fftw_alloc_complex(alloc_local) 

call c_f_pointer(cdatac, gaux, [nx,local_L]) !transposed 

alloc_local = fftw_mpi_local_size_2d(nx,ny/2+1, MPI_COMM_WORLD, & 
    local_M, local_M_offset) 


cdatar = fftw_alloc_real(2*alloc_local) 

call c_f_pointer(cdatar, faux, [2*(ny/2+1),local_M]) 

! Create plans 

p1 = fftw_mpi_plan_dft_r2c_2d(nx,ny,faux,gaux, MPI_COMM_WORLD, & 
     ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT)) 


p2 = fftw_mpi_plan_dft_c2r_2d(nx,ny,gaux,faux, MPI_COMM_WORLD, & 
     ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN)) 

! EXECUTE FFTW 

call random_number(faux) 


print *, "real input:", real(faux(1:ny,:)) 

call fftw_mpi_execute_dft_r2c(p1,faux,gaux) 

call fftw_mpi_execute_dft_c2r(p2, gaux, faux) 

print *, "real output:", real(faux(1:ny,:))/(nx*ny) 

call fftw_destroy_plan(p1) 
call fftw_destroy_plan(p2) 


call mpi_finalize(ierr) 


End Program TEST