计算科学 - 为什么 Eigen 分配一个临时值来并行评估 A.noalias() = B.transpose()*C？ - 吾爱随笔录

我编写了一个程序，它使用矩阵乘法迭代地转换数据。为了尽量减少大内存分配的数量，我使用两个大致相等的std::vector<double>'sdata和temp_data分别存储当前数据和新数据。对于矩阵运算，我将每个缓冲区包装成一个Eigen::Map<Eigen::MatrixXd>, 分别称为data_matrix和temp_data_matrix。然后，我需要计算以下操作：

temp_data_matrix.noalias() = data_matrix.transpose()*factor;

wherefactor只是一个Eigen::MatrixXd远小于data_matrixortemp_data_matrix的。多亏了.noalias()，这条线通常不会分配任何大型临时文件。然而，当多线程程序突然分配一个临时（大致）等于temp_data！在实践中，这可以将内存使用量增加约 50%，所以我对以下几点感到好奇：

为什么只有在多线程时才会发生这种情况？什么算法方面需要在输出缓冲区顶部的临时缓冲区？
为什么这条线不会发生这种情况temp_data_matrix.noalias() = factor*data_matrix？毕竟，我的操作似乎更高效缓存：因为我们转置data_matrix，我们只需要从data_matrix和计算每对列的标量积，factor考虑到两个矩阵都是列优先的，这应该很快。但是，我没有考虑阻塞的影响，我知道这对性能非常重要。
data在不分配or顺序的任何临时对象的情况下，实现上述行的最快方法是什么temp_data？如果并行算法需要一些“临时空间”，它不能temp_data_matrix用于那个吗？我不关心复制factor，因为它在实践中相对较小。目前最快的解决方案似乎是这条线切换回单线程，但我希望这可以改进。我尝试使用lazyProduct，但结果比单线程解决方案慢得多。

附录：上面的观察是在我的原始程序和下面的示例代码中使用 Massif (valgrind) 进行的。我使用命令编译了示例程序g++ -O3 -DNDEBUG -I /usr/include/eigen3 -fopenmp eigen_memory_test.cpp -o eigen_memory_test。

#include <iostream>
#include <Eigen/Dense>
#include <algorithm>
#include <random>
#include <cstring>

int main(int argc, char** argv) {
    
    bool parallel = (argc > 1 && strcmp(argv[1], "parallel") == 0);
    std::cout << "parallel: " << (parallel  ? "true" : "false") << std::endl;
    bool lazyproduct = (argc > 2 && strcmp(argv[2], "lazyproduct") == 0);
    std::cout << "lazyproduct: " << (lazyproduct  ? "true" : "false") << std::endl;
    bool test = true;
    
    // Initialize parallellism
    int threads = 4;
    Eigen::setNbThreads(threads);
    
    // Initialization and memory allocation
    size_t columns = 100000;
    std::vector<size_t> rows_per_step = {100, 90, 80, 70, 60, 50};
    size_t steps = rows_per_step.size() - 1;
    std::vector<double> data(rows_per_step[0]*columns);
    std::vector<double> temp_data(rows_per_step[0]*columns);
    
    // Generate random data
    std::random_device rnd_device;
    std::mt19937 mersenne_engine {rnd_device()};
    std::uniform_real_distribution<double> distribution(-1.0, 1.0);
    auto generator = std::bind(distribution, std::ref(mersenne_engine));
    std::generate(data.begin(), data.end(), generator);
    
    // If testing, print first elements to verify later
    if (test)
        std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;
    
    // Process steps
    for (size_t step = 0; step < steps; step++) {
        
        size_t initial_rows = rows_per_step[step];
        size_t new_rows = rows_per_step[step + 1];
        
        // Choose factor
        Eigen::MatrixXd factor = (test ? 
            Eigen::MatrixXd(Eigen::MatrixXd::Identity(initial_rows, new_rows)) : 
            Eigen::MatrixXd(Eigen::MatrixXd::Random(initial_rows, new_rows))
        );
        
        // Transform data
        Eigen::Map<Eigen::MatrixXd> data_matrix(data.data(), initial_rows, columns);
        Eigen::Map<Eigen::MatrixXd> temp_data_matrix(temp_data.data(), columns, new_rows);
        if (!parallel)
            Eigen::setNbThreads(1);
        // In my real problem setting, the transpose in the following matrix multiplication is actually necessary
        if (lazyproduct)
            temp_data_matrix.noalias() = data_matrix.transpose().lazyProduct(factor);
        else
            temp_data_matrix.noalias() = data_matrix.transpose()*factor;
        Eigen::setNbThreads(threads);
        
        // Transpose back to original data
        // This doesn't happen in the real problem, but I put this here to keep this example simpler
        Eigen::Map<Eigen::MatrixXd> transposed_data_matrix(data.data(), new_rows, columns);
        transposed_data_matrix.noalias() = temp_data_matrix.transpose();
        
    }
    
    // Print first elements to make sure matrix operations aren't optimized away
    std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;
    
    return 0;
    
}
```