c++ - 如何使用 Eigen 避免稀疏表达式中的内存分配

Question

我有一个稀疏模式恒定的应用程序。可以说我的计算采用以下形式

sm3 = sm1 + sm2

但是，即使我将所有这些操作数中的稀疏模式设置为相同，我的分析器也显示大部分时间都花在分配和取消分配结果矩阵上。

这是我的 MWE:

#include <eigen3/Eigen/Sparse>
#include <iostream>

int main(int argc, char *argv[])
{
  using namespace Eigen;

  SparseMatrix<double> sm1(2, 2), sm2(2, 2), sm3(2, 2);

  // Populate sm1 and sm2
  sm1.insert(0,0) = 2.0;

  sm2.insert(1,1) = 3.0;

  // Compute the result pattern
  sm3 = sm1 + sm2;

  // Copy the augmented pattern into the operands
  sm1 = sm2 = sm3;

  // This loop triggers a lot of new[] and delete[] calls
  for(int i = 0; i < 1000; i++)
    sm3 = sm2 + sm1;
}

这些分配操作可以避免吗？

Answer 1

目前这是不可能的，因为默认情况下稀疏矩阵被假定为别名。例如，如果您这样做:

m3 = m3 + m1;

如果 m1 模式未完全包含在 m3 模式中，则不可能直接在 m3 中计算表达式。在 future ，我们可以使用如下语法强制重用目标内存:

m3.noalias() = m1 + m2;

同时，由于您的矩阵很小，您可以通过强制 m1 和 m2 的模式与模式相同来解决问题，甚至获得更高的性能通过添加一些显式零来计算m3。然后，使用 Eigen 3.3，您可以将稀疏加法转换为密集 vector 加法:

m3.coeffs() = m1.coeffs() + m2.coeffs();

即使 m1 和 m2 之间的交集很小，您也会获得非常高的加速(可能是一个数量级)，因为您消除了内存间接寻址，并从矢量化中受益(不要忘记启用 AVX，例如 -mavx)。