通过 fftw_mpi_r2c_2d 和 fftw_mpi_c2r_2d 输出不正确

我编写了一个简单的测试程序，以便在 2d 域（使用 Fortran）中使用 MPI 实现 FFTW。该域的宽度为“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

该问题是由于padding http://fftw.org/doc/Multi_002ddimensional-MPI-DFTs-of-Real-Data.htmlfftw 需要：

虽然真实数据在概念上是 n0 × n1 × n2 × … × nd-1 ，但它在物理上存储为 n0 × n1 × n2 × … × [2 (nd-1/2 + 1)] 数组，其中最后一个维度已被填充以使其与复数输出的大小相同。这很像就地串行 r2c/c2r 接口（请参阅实际数据的多维 DFT），只不过在 MPI 中，即使对于异地数据也需要填充。

因此，16x16 变换的输入数组是 16x18 数组。每行末尾额外的两个数字的值在真实空间中是没有意义的。然而，当 c 指针被转换为 Fortran 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)

这是完整的代码，基于您的代码和以下之一如果可能的话，如何进行 fftw3 MPI“转置”2D 变换？ https://stackoverflow.com/questions/24290014/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