如何用标量/整数对Eigen张量做除法？FFT归一化遇零值问题

为什么会得到全零结果？ #

• 整数除法陷阱：1/(nx*ny*nz)中，nx/ny/nz都是整数，C++里整数间除法会执行整数除法，结果直接取整截断。当nx*ny*nz > 1时，这个表达式的结果就是0，后续和张量相乘自然得到全零。
• 错误调用data()：rArr.data()返回的是张量底层数据的指针，不是张量对象，直接和数值相乘是完全错误的操作（本质是把指针地址和数值做乘法，和张量元素毫无关系）。

高效的归一化实现方式 #

不需要额外创建dummy张量，直接在rArr上执行原地运算即可，既节省内存开销又提升执行效率：

// 用1.0触发浮点数除法，规避整数除法问题
rArr *= 1.0 / (nx*ny*nz);

如果担心nx*ny*nz整数溢出（比如三个维度数值都很大），可以拆分转换为浮点数再计算，更安全：

const double scale = 1.0 / static_cast<double>(nx) / ny / nz;
rArr *= scale;

修正后的完整代码片段 #

void c2rfft3d(Eigen::Tensor<std::complex<double>, 3>& cArr, Eigen::Tensor<double, 3>& rArr){    
    fftw_complex *input_array = (fftw_complex*) fftw_malloc(nx*ny*nz * sizeof(fftw_complex));
    memcpy(input_array, cArr.data(), nx*ny*nz * sizeof(fftw_complex));
    
    fftw_complex *output_array = (fftw_complex*) fftw_malloc(nx*ny*nz * sizeof(fftw_complex));
    
    fftw_plan backward = fftw_plan_dft_3d(nx, ny, nz, input_array, output_array, FFTW_BACKWARD, FFTW_ESTIMATE);
    
    fftw_execute(backward);
    fftw_destroy_plan(backward);
    fftw_cleanup();

    memcpy(rArr.data(), output_array, nx*ny*nz * sizeof(double));
    
    // 原地归一化，高效且避免错误
    const double scale = 1.0 / static_cast<double>(nx*ny*nz);
    rArr *= scale;

    // 释放fftw分配的内存，避免泄漏
    fftw_free(input_array);
    fftw_free(output_array);
}

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