C++中Eigen::Map映射FFTW数组到Eigen矩阵结果错误问题

使用Eigen的Map特性将已有内存映射为Eigen矩阵时，映射FFTW的C++复数数组仅返回部分正确切片，未覆盖完整内存块。

问题复现代码 #

static const int nx = 10;
static const int ny = 10; 
static const int nyk = ny/2 + 1;
static const int nxk = nx/2 + 1;
static const int ncomp = 2;

fftw_complex *uhk; // this is pert Te
uhk= (fftw_complex*) fftw_malloc((((nx)*(ny+1))*nyk)* sizeof(fftw_complex)); 
memset(uhk, 42, (((nx))*nyk)* sizeof(fftw_complex));

for (int i = 0; i < nx; i++){
    for (int j = 0; j < nyk; j++){
        for (int k = 0; k < ncomp; k++){
            uhk[i + nyk*j][k] = //taking fft of some expression
        }
    }
}

Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned> uhOut(reinterpret_cast<std::complex<double>*>(uhkOut),nyk,nx);
std::cout << uhOut<< '\n';

异常输出 #

(0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00) (0.58,-0.28) (0.95,-0.46) (0.95,-0.46) (0.58,-0.28) (0.00,-0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)
 (0.00,0.00) (0.59,-0.09) (0.95,-0.14) (0.95,-0.14) (0.59,-0.09) (0.00,-0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)  (0.00,0.00)

预期输出 #

(0.00,0.00)     (5.14,0.00)     (8.32,0.00)     (8.32,0.00)     (5.14,0.00)     (0.00,0.00)    (-5.14,0.00)    (-8.32,0.00)    (-8.32,0.00)    (-5.14,0.00)
     (0.00,0.00)    (2.43,-2.35)    (3.93,-3.81)    (3.93,-3.81)    (2.43,-2.35)    (0.00,-0.00)    (-2.43,2.35)    (-3.93,3.81)    (-3.93,3.81)    (-2.43,2.35)
     (0.00,0.00)    (0.52,-1.04)    (0.85,-1.68)    (0.85,-1.68)    (0.52,-1.04)    (0.00,-0.00)    (-0.52,1.04)    (-0.85,1.68)    (-0.85,1.68)    (-0.52,1.04)
     (0.00,0.00)    (0.57,-0.55)    (0.93,-0.89)    (0.93,-0.89)    (0.57,-0.55)    (0.00,-0.00)    (-0.57,0.55)    (-0.93,0.89)    (-0.93,0.89)    (-0.57,0.55)
     (0.00,0.00)    (0.58,-0.28)    (0.95,-0.46)    (0.95,-0.46)    (0.58,-0.28)    (0.00,-0.00)    (-0.58,0.28)    (-0.95,0.46)    (-0.95,0.46)    (-0.58,0.28)
     (0.00,0.00)    (0.59,-0.09)    (0.95,-0.14)    (0.95,-0.14)    (0.59,-0.09)    (0.00,-0.00)    (-0.59,0.09)    (-0.95,0.14)    (-0.95,0.14)    (-0.59,0.09)

和Map的功能本身无关，代码存在三处明确错误：

1. 数组索引计算完全错误：写入uhk时使用的偏移i + nyk*j不符合线性内存排布逻辑。要映射的是nyk行、nx列的列优先Eigen矩阵，矩阵第i列、第j行元素的内存偏移应为j + nyk * i。错误索引导致大部分FFT计算结果被写到了错误的内存位置，仅当i、j取循环上限附近的数值时，偏移刚好和正确位置重合，因此仅最后两行的部分值和预期匹配，其余位置全是memset填充的初始值，看起来像被截断。
2. Map传入野指针：Map构造时传入的指针是uhkOut，但实际申请内存、写入数据的指针是uhk，访问未初始化指针属于未定义行为。
3. 内存操作长度不匹配：fftw_malloc申请了nx*(ny+1)*nyk个复数长度的内存，但memset仅初始化了nx*nyk个复数长度，多余的申请内存没有实际用途，还容易干扰索引判断。

补充说明：fftw_complex和std::complex<double>内存布局完全兼容，reinterpret_cast的类型转换本身没有问题，不需要怀疑类型匹配性。

• 修正循环内的数组索引，匹配Eigen列优先的内存排布规则
• 把Map构造函数的传入指针改为实际持有数据的uhk
• 内存申请和初始化长度保持一致，不需要申请多余内存

修正后的核心代码如下：

static const int nx = 10;
static const int ny = 10; 
static const int nyk = ny/2 + 1;
static const int ncomp = 2;

fftw_complex *uhk;
// 申请实际需要的内存长度
uhk= (fftw_complex*) fftw_malloc(nx * nyk * sizeof(fftw_complex)); 
memset(uhk, 0, nx * nyk * sizeof(fftw_complex));

for (int i = 0; i < nx; i++){
    for (int j = 0; j < nyk; j++){
        for (int k = 0; k < ncomp; k++){
            // 修正索引，匹配列优先内存偏移
            uhk[j + nyk * i][k] = // 写入FFT计算结果
        }
    }
}
// 传入正确的指针uhk
Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned> uhOut(reinterpret_cast<std::complex<double>*>(uhk), nyk, nx);
std::cout << uhOut<< '\n';

// 用完记得释放fftw申请的内存
fftw_free(uhk);

如果写入数据时习惯按行优先顺序排布，可以在Map时指定行优先存储标志，不需要强行调整索引，写法如下：

Eigen::Map<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>, Eigen::Unaligned> 
uhOut(reinterpret_cast<std::complex<double>*>(uhk), nyk, nx);

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