如何在 C++ 中使用 FFTW3 创建场向量？

Question

我将 FFTW3 与 OpenMPI 结合使用，并且我想要一个能够执行 FFTW 变换的场向量。有了代码就更清楚了：

#include <complex>
#include <fftw3-mpi.h>

int main(int argc, char **argv)
{

MPI_Init(&argc, &argv);
fftw_mpi_init();

ptrdiff_t N0=50, N1=10, N2=1;
ptrdiff_t alloc_local, local_n0, local_0_start;

alloc_local = fftw_mpi_local_size_3d(N0, N1, N2 / 2 + 1, MPI_COMM_WORLD, &local_n0,&local_0_start);

fftw_complex **kVect;
kVect = (fftw_complex **)malloc(3 * sizeof(fftw_complex *));
for (int dim = 0; dim < 3; dim++)
    kVect[dim] = fftw_alloc_complex(alloc_local);

ptrdiff_t i, j, k;

for (i = 0; i < local_n0; i++)
    for (j = 0; j < N1; j++)
        for (k = 0; k < N2; k++)
        {
            kVect[0][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][0] = 1.0 / (2.0 * M_PI);
            kVect[0][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][1] = 0.0;
            kVect[1][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][0] = 1.0 / (2.0 * M_PI);
            kVect[1][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][1] = 0.0;
            kVect[2][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][0] = 1.0 / (2.0 * M_PI);
            kVect[2][(i * N1 + j) * (2 * (N2 / 2 + 1)) + k][1] = 0.0;
        }

for (int dim = 0; dim < 3; dim++)
    fftw_free(kVect[dim]);

free(kVect);

MPI_Finalize();

return 0;
}

我认为这段代码应该可以工作。但是，在空闲步骤

fftw_free(kVect[dim]);

mpiexec 发送错误消息

double free or corruption

。

我尝试用 3

for

循环注释初始化部分，并且代码不会崩溃。考虑到 FFTW3，我想我使用 malloc 的方式不适合，但我不知道如何为这个向量分配内存。

Answer 1

我找到了这个问题的解决方案。它来自 FFTW 对数组的定义。为了进行 FFT，复数数组和实数数组必须具有相同的维度，但复数数组的大小是实数数组的两倍。因此，要将内存分配给实际字段，可以使用

fftw_alloc_real(2*alloc_local)

。 要访问位置

(i,j,k)

处的字段，请使用示例中提到的方程

(i * N1 + j) * (2 * (N2 / 2 + 1)) + k

：https://www.fftw.org/fftw3_doc/Multi_002dDimension-MPI-DFTs-of-Real-Data.html 。对于复杂的字段，应该使用

k + j*N2 + i*N1*N2

代替，否则会超出字段的边界并触发此错误。