使用FFTW3与OpenMPI创建场向量时遇内存释放错误求助

我正在结合FFTW3与OpenMPI开发，想要创建一个可执行FFTW变换的场向量，但运行代码时在执行fftw_free(kVect[dim]);时，mpiexec抛出“double free or corruption”错误。注释掉三层循环的初始化部分后程序不再崩溃，怀疑内存分配方式有问题，求助正确的内存分配方法。

对应的代码如下：

#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;
}

崩溃的核心原因是内存越界访问，并非内存分配方式错误：

• 当N2=1时，N2/2 +1 = 1，alloc_local对应的是local_n0 × N1 × 1的复数元素总量。
• 但你在循环中计算索引时用了(i * N1 + j) * (2 * (N2 / 2 + 1)) + k，相当于把数组当成local_n0 × N1 × 2的大小访问，直接超出了分配的内存范围。
• 越界写操作会破坏FFTW内部的内存管理元数据，导致后续fftw_free触发双重释放或内存损坏错误。

修正步骤：

1. 索引计算要与fftw_mpi_local_size_3d的参数匹配，直接用(i * N1 + j) * (N2 / 2 + 1) + k作为复数数组的索引（fftw_complex本身已包含实部和虚部两个double，无需额外乘2）。
2. 内存分配流程本身是正确的：用malloc分配指针数组，fftw_alloc_complex分配各维度的复数内存，最后分别用fftw_free和free释放即可。

修正后的循环代码：

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

额外注意：后续创建FFT变换计划时，需确保维度参数与内存分配时的N0,N1,N2/2+1一致，避免再次出现不匹配问题。

内容的提问来源于stack exchange，提问作者Sinai