c++ - 遍历 Eigen 矩阵的最有效方法

我正在创建一些函数来执行诸如负数和正数的“分隔和”，kahan，成对和其他东西，这些东西与从矩阵中获取元素的顺序无关紧要，例如:

template <typename T, int R, int C>
inline T sum(const Eigen::Matrix<T,R,C>& xs)
{
  T sumP(0);
  T sumN(0);
  for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nRows; ++i)
   for (size_t j = 0; j < nCols; ++j)
   {
        if (xs(i,j)>0)
          sumP += xs(i,j);
        else if (xs(i,j)<0) //ignore 0 elements: improvement for sparse matrices I think
          sumN += xs(i,j);
   }
 return sumP+sumN;
}

现在，我想使它尽可能高效，所以我的问题是，像上述那样遍历每行的每一列会更好，还是像下面这样做相反的事情:

for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nCols; ++i)
  for (size_t j = 0; j < nRows; ++j)

(我想这取决于矩阵元素在内存中分配的顺序，但是我在Eigen的手册中找不到)。

此外，还有其他替代方法，例如使用迭代器(它们存在于Eigen中吗？)可能会更快一些？

最佳答案

默认情况下，Eigen按列优先(Fortran)顺序分配矩阵(documentation)。

迭代矩阵的最快方法是按照存储顺序进行，以错误的方式进行迭代将增加高速缓存未命中的次数(如果矩阵不适合L1，则将占用您的计算时间，因此，读取会增加计算时间)减少了缓存行/elemsize的因子(可能是64/8 = 8)。

如果您的矩阵适合L1缓存，则不会有什么区别，但是好的编译器应该能够对循环进行矢量化处理，如果启用AVX(在 Shiny 的新内核i7上)，则可以使速度提高4倍。 (256位/64位)。

最后，不要指望Eigen的任何内置函数都能为您提速(我不认为有迭代器，但我可能会误会)，它们只是会为您提供相同的功能(非常简单) ) 代码。

TLDR:交换迭代顺序，您需要最快速地更改行索引。