c++ - 如何在 Eigen 中对子矩阵求和

Question

我有一些矩阵定义为:

Eigen::MatrixXd DPCint = Eigen::MatrixXd::Zero(p.szZ*(p.na-1),p.szX);

\\ perform some computations and fill every sub-matrix of size [p.szZ,p.szX] with some values
#pragma omp parallel for
for (int i=0; i < p.na-1; i++)
{
...
DPCint(Eigen::seq(i*p.szZ,(i+1)*p.szZ-1),Eigen::all) = ....;
}

\\ Now sum every p.szZ rows to get a matrix that is [p.szZ,p.szX]

在 Matlab 中，此操作快速而简单。如果我想将循环与 OpenMP 并行化，我不能简单地在这里执行 += 操作。类似地，我可以遍历每组 p.szZ 行并对其求和，但该循环无法并行化，因为每个线程都会输出相同的数据。有没有一些有效的方法可以使用 Eigen 的索引操作来求和子矩阵？这似乎是一个简单的操作，我觉得我遗漏了一些东西，但我已经有一段时间没能找到解决方案了。

澄清

基本上，在上面的循环之后，我想在一行中完成:

for (int i = 0; i < p.na-1; i++)
{
DPC += DPCint(Eigen::seq(i*p.szZ,(i+1)*p.szZ-1),Eigen::all);
}

在 matlab 中，我可以简单地将矩阵 reshape 为 3D 矩阵并沿三维求和。我对Eigen的tensor库不熟悉，希望这个操作不用借助tensor库也可以实现。但是，我的首要任务是速度和效率，所以我愿意接受任何建议。

Answer 1

在基于 na 的轴上执行并行归约效率不高。事实上，这个维度对于多线程来说已经非常小了，但它也(几乎)强制线程在效率低下的临时矩阵上运行(这是内存限制的，所以它不能很好地扩展)。

另一种解决方案是并行化szZ 维度。每个线程都可以在一个切片上工作，并在没有临时矩阵的情况下执行局部缩减。此外，这种方法还应该改进 CPU 缓存的使用(因为每个线程计算的 DPC 部分更有可能适合缓存，因此它们不会从 RAM 中重新加载)。这是一个(未经测试的)示例:

// All thread will execute the following loops (all iterations but on different data blocks)
#pragma omp parallel
for (int i = 0; i < p.na-1; i++)
{
    // "nowait" avoid a synchronization but this require a 
    // static schedule which is a good idea to use here anyway.
    #pragma omp for schedule(static) nowait
    for (int j = 0; j < p.szZ; j++)
        DPC(j, Eigen::all) += DPCint(i*p.szZ+j, Eigen::all);
}

正如@chtz 所指出的，最好避免使用临时DPCint 矩阵，因为内存吞吐量是一种非常有限的资源(尤其是在并行代码中)。

编辑: 我假设矩阵存储在行优先存储顺序中，默认情况下并非如此。这可以修改(参见 the doc )，实际上它会使第一个和第二个循环缓存高效。然而，混合存储顺序通常容易出错，并且使用行优先排序会迫使您重新定义基本类型。 @Homer512 的解决方案是一种替代实现，当然更适合列主矩阵。