通过fftw_mpi_r2c_2d和fftw_mpi_c2r_2d输出不正确 [英] Incorrect output via fftw_mpi_r2c_2d and fftw_mpi_c2r_2d

问题描述

我编写了一个简单的测试程序，以便在2d域(使用Fortran)中使用MPI实现FFTW.该域的宽度为"Ny x Nx"，并在第二个('x')索引中进行了分区.

I wrote a simple test program in order to implement the FFTW with MPI in a 2d domain (with Fortran). The domain is 'Ny x Nx' wide and partitioned in the second ('x') index.

在正确(我相信吗?)声明和分配变量和计划之后，我调用了fftw_mpi r2c_2d函数，接下来我用fftw_mpi c2r_2d转换回其输出，以检查是否获得了原始输入. r2c_2d部分似乎工作正常.但是，在使用c2r_2d函数将输出转换回(部分归一化)后，我没有得到原始输入:结果向量在索引(:，j)上显示零"，其中j对应于"Ny/2"的倍数.我究竟做错了什么?谢谢！

After proper (I believe?) declaration and allocation of variables and plans, I call the fftw_mpi r2c_2d function and next I transform back its output with the fftw_mpi c2r_2d, in order to check if I get the original input. The r2c_2d part seems to work fine. However, I don't get the original input after transforming back the output (apart normalization) with the c2r_2d function: the resulting vector displays 'zeros' at the indices (:,j) with j corresponding to multiples of 'Ny/2'. What am I doing wrong? Thanks!

以下是代码的摘录:

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

推荐答案

该问题归因于

尽管实际数据在概念上为n0×n1×n2×…×nd-1，但实际存储为n0×n1×n2×…×[2(nd-1/2 + 1)]数组，其中最后一个尺寸已被填充，使其与复杂输出的尺寸相同.这与就地串行r2c/c2r接口(请参阅实数据的多维DFT)非常相似，除了在MPI中，即使对于就地数据也需要填充.

Although the real data is conceptually n0 × n1 × n2 × … × nd-1 , it is physically stored as an n0 × n1 × n2 × … × [2 (nd-1/2 + 1)] array, where the last dimension has been padded to make it the same size as the complex output. This is much like the in-place serial r2c/c2r interface (see Multi-Dimensional DFTs of Real Data), except that in MPI the padding is required even for out-of-place data.

因此，用于16x16变换的输入数组因此是16x18数组.每行末尾的额外两个数字的值在实际空间中是没有意义的.但是，在将c指针强制转换为fortran 2D数组时，不要忘记这些额外的数字:

Hence, the input array for a 16x16 transform is therefore a 16x18 array. The value of the extra two numbers at the end of each rows are meaningless in the real space. Yet, these extra numbers must not be forgotten as the c pointer is cast to a fortran 2D array:

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

多余的数字仍会打印在每一行的末尾.可以对数组进行切片以避免打印出这些毫无价值的值:

The extra numbers are still printed at the end of each row. The array can be sliced to avoid printing these worthless values:

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

这是完整的代码，基于您的代码和是否可以使用2D转换?可以由mpif90 main.f90 -o main -I/usr/include -L/usr/lib -lfftw3_mpi -lfftw3 -lm编译并由mpirun -np 2 main运行.

Here is the complete code, based on yours and the one of How to do a fftw3 MPI "transposed" 2D transform if possible at all? It can be compiled by mpif90 main.f90 -o main -I/usr/include -L/usr/lib -lfftw3_mpi -lfftw3 -lm and ran by 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

这篇关于通过fftw_mpi_r2c_2d和fftw_mpi_c2r_2d输出不正确的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！