Eigen:高效实现C++数组矩阵

5

是否可以实现一个类,该类将C风格的指针作为模板参数并以提供的内存解析为静态Eigen矩阵?声明可能如下所示:

EIGEN_ALIGN16 double array[9];
CMatrix<double,3,3,array> :: m;

我对地图有所了解,但是下面提供的示例代码已经证明相比静态Eigen矩阵,使用它们会慢20%。

下面是前提条件:

  • 我需要提供自己的C指针,这样我就可以高效地重用C代码而不需要进行拷贝。
  • 结果矩阵应该对Eigen来说看起来是静态的,这样Eigen就可以像在编译时处理静态数组一样进行优化。请看我的上面的示例,在编译时我会同时提供矩阵大小(静态)和C指针。
  • CMatrix 应该回退到 Eigen::Matrix。当未提供C数组的附加模板参数时,我将获得普通的Eigen矩阵。
  • 我不打算制作完整的Eigen扩展。我的意思是,我不关心提供所有类型的检查以为其他用户提供一个整洁的扩展。我只想找到这个问题的最有效的解决方案。

通过添加一个新的构造函数实现一个解决方案是否可能?比如:

EIGEN_ALIGN16 double data[9];
Eigen::Matrix<double,3,3> m(data); //where data is NOT copied but used to replace the static allocation used by default.

在下面找到了我的用于基准测试Map和Matrix效率的代码。它是自包含的,您可以使用以下方式进行编译:
g++ -Ofast -DNDEBUG -DEIGEN_NO_MALLOC -I/path_to_my_Eigen benchmark1.cpp -o benchmark1 -lrt

这里是代码:
#include <Eigen/Eigen>
#include <bench/BenchTimer.h>

#include <iostream>

using namespace Eigen;
using namespace std;

//#define CLASSIC_METHOD
#define USE_MAPS

EIGEN_DONT_INLINE void classic(double VO[4], double AT[4][4], double VI[4])
{
  for (int ii=0; ii<4; ii++)
    {
      VO[ii] = AT[ii][0] * VI[0] +
               AT[ii][1] * VI[1] +
               AT[ii][2] * VI[2] +
               AT[ii][3] * VI[3];
    }
};

template <typename OutputType, typename MatrixType, typename VectorType>
EIGEN_DONT_INLINE void modern(MatrixBase<OutputType>& VOE, const MatrixBase<MatrixType>& A44, const MatrixBase<VectorType>& VIE)
{
  VOE.noalias() = A44.transpose()*VIE;
};

int main()
{
  EIGEN_ALIGN16 double AT[4][4] = {0.1, 0.2, 0.3, 2.0, 0.2, 0.3, 0.4, 3.0, 0.3, 0.4, 0.5, 4.0, 0.0, 0.0, 0.0, 1.0};
  EIGEN_ALIGN16 double VI[4] = {1, 2, 3, 4};
  EIGEN_ALIGN16 double VO[4];

//Eigen matrices

#ifndef USE_MAPS
  Matrix4d A44 = Matrix4d::MapAligned(AT[0]);
      Vector4d VIE = Vector4d::MapAligned(VI);
  Vector4d VOE(0,0,0,0);
#else
  Map<Matrix4d,Aligned> A44(AT[0]);
  Map<Vector4d,Aligned> VIE(VI);
  Map<Vector4d,Aligned> VOE(VO);

  // Map<Matrix4d> A44(AT[0]);                                                                                                                                                                                                                                      
  // Map<Vector4d> VIE(VI);                                                                                                                                                                                                                                           
  // Map<Vector4d> VOE(VO);

#endif

#ifdef EIGEN_VECTORIZE
  cout << "EIGEN_VECTORIZE defined" << endl;
#else
    cout << "EIGEN_VECTORIZE NOT defined" << endl;
#endif

  cout << "VIE:" << endl;
  cout << VIE << endl;

  VI[0] = 3.14;
  cout << "VIE:" << endl;
  cout << VIE << endl;

  BenchTimer timer;

  const int num_tries = 5;
  const int num_repetitions = 200000000;

#ifdef CLASSIC_METHOD
  BENCH(timer, num_tries, num_repetitions, classic(VO, AT, VI));
  std::cout << Vector4d::MapAligned(VO) << std::endl;
#else
  BENCH(timer, num_tries, num_repetitions, modern(VOE, A44, VIE));
  std::cout << VOE << std::endl;
#endif

  double elapsed = timer.best();
  std::cout << "elapsed time: " << elapsed*1000.0 << " ms" << std::endl;

  return 0;
}

使用 std::size_t; 模板 <typename T, size_t Rows, size_t Cols, T (&Var)[Rows][Cols]> 结构体 MyEfficientVector { 使用 matrix_type = decltype(Var); matrix_type m = Var; } ; - Massa
@gha.st。你认为应该如何实现这个呢? - Alejandro
@Alejandro 这是一个非常好的问题。我不知道是否可能,因为 Eigen::MatrixBase 的数据在 Eigen::DenseStorage 中。也许你应该先看一下 Eigen/DenseStorage.h,了解你应该如何替换其默认存储模式为你设计的存储模式。无论如何,如果你想扩展 Eigen 的功能,你的“快速而肮脏”的计划都失败了,因为你肯定需要制作一个“正确的”Eigen扩展... - Massa
@gha.st。谢谢。虽然调用data()对我来说不够。想象一下这样的情况,我想要矩阵a和矩阵b,但我希望它们两个在内部都存储相同的数组。 - Alejandro
@gha.st。嗯,我在网上发现这种做法通常是不被鼓励的。对我来说,这看起来像是一个恶心的把戏。虽然我曾经考虑过这个选项,但我会尽量避免它。 - Alejandro
显示剩余6条评论
1个回答

5

虽然有些跑题,但由于您强调性能:

Eigen汇编并不总是最优的 - 由于寄存器复用差和写回内存的开销(这绝对不是Eigen的错 - 在通用模板中做到这一点是不可能的)。

如果您的核心代码比较简单(QCD?),我建议手动编写汇编代码(使用指令集)。

以下是使用指令集重写的经典核心代码,比Eigen版本更快,对于Map/Matrix类型也是如此(因此您无需发明自己的类型)。

EIGEN_DONT_INLINE void classic(double * __restrict__ VO, const double * __restrict__ AT, const double * __restrict__ VI) {
  __m128d vi01 = _mm_load_pd(VI+0);
  __m128d vi23 = _mm_load_pd(VI+2);
  for (int i = 0; i < 4; i += 2) {
    __m128d v00, v11;
    // v[i+0,i+0]                                                                                                                                                                                                   
    {
      int ii = i*4;
      __m128d at01 = _mm_load_pd(&AT[ii + 0]);
      __m128d at23 = _mm_load_pd(&AT[ii + 2]);
      v00 = _mm_mul_pd(at01, vi01);
      v00 = _mm_add_pd(v00, _mm_mul_pd(at23, vi23));
    }
    // v[i+1,i+1]                                                                                                                                                                                                   
    {
      int ii = i*4 + 4;
      __m128d at01 = _mm_load_pd(&AT[ii + 0]);
      __m128d at23 = _mm_load_pd(&AT[ii + 2]);
      v11 = _mm_mul_pd(at01, vi01);
      v11 = _mm_add_pd(v11, _mm_mul_pd(at23, vi23));
    }

    __m128d v = _mm_hadd_pd(v00, v11);
    // v is now [v00[0] + v00[1], v11[0] + v11[1]]                                                                                                                                                                               
    _mm_store_pd(VO+i, v);
    // VO[i] += AT[j+0 + i*4]*VI[j+0];                                                                                                                                                                              
    // VO[i] += AT[j+1 + i*4]*VI[j+1];                                                                                                                                                                              
  }
}

通过交错加载和乘加运算,您可以获得一些额外的改进-我尽力保持简单。

以下是结果:

g++ -Ofast -DNDEBUG -DEIGEN_NO_MALLOC -DCLASSIC_METHOD -I /usr/local/eigen benchmark1.cpp -o benchmark1 -lrt -msse4; ./benchmark1 
elapsed time: 611.397 ms

g++ -Ofast -DNDEBUG -DEIGEN_NO_MALLOC -DCLASSIC_METHOD -DUSE_MAPS -I /usr/local/eigen benchmark1.cpp -o benchmark1 -lrt -msse4; ./benchmark1
elapsed time: 615.541 ms

g++ -Ofast -DNDEBUG -DEIGEN_NO_MALLOC -DUSE_MAPS -I /usr/local/eigen benchmark1.cpp -o benchmark1 -lrt -msse4; ./benchmark1
elapsed time: 981.941 ms

g++ -Ofast -DNDEBUG -DEIGEN_NO_MALLOC -I /usr/local/eigen benchmark1.cpp -o benchmark1 -lrt -msse4; ./benchmark1 
elapsed time: 838.852 ms

进一步介绍,如果您的矩阵被转置,则可能编写更好的simd内核-水平加法(_mm_hadd_pd)非常昂贵。
在评论中增加讨论:将Eigen映射移动到函数内会消除映射和矩阵参数之间的时间差异。
EIGEN_DONT_INLINE void mapped(double (&VO)[4], const double (&AT)[4][4], const double (&VI)[4]) {
  Map<const Matrix4d,Aligned> A44(&AT[0][0]);
  Map<const Vector4d,Aligned> VIE(VI);
  Map<Vector4d,Aligned> VOE(VO);
  VOE.noalias() = A44.transpose()*VIE;
}

当将Map传递给函数(未内联函数)时,这是汇编代码的顶部。

    movq    (%rsi), %rcx
    movq    (%rdx), %rax
    movq    (%rdi), %rdx
    movapd  (%rcx), %xmm0
    movapd  16(%rcx), %xmm1
    mulpd   (%rax), %xmm0
    mulpd   16(%rax), %xmm1

与传递数组引用(及其中的映射)或矩阵相比。
    movapd  (%rsi), %xmm0
    movapd  16(%rsi), %xmm1
    mulpd   (%rdx), %xmm0
    mulpd   16(%rdx), %xmm1

哇!这是一项惊人的工作!谢谢!然而,为每个操作编写汇编语言超出了我有限的时间范围。您知道Eigen是否在后台执行此操作?如果不是,您是否知道其他任何库可以实现此功能?您对我提出的将C数组作为附加模板参数传递的想法有何看法?您是否知道这是否可行?尝试在编译时欺骗Eigen,即创建一个静态数组,但实际上是由用户传递的内存。 - Alejandro
我的猜测是Map和Matrix之间的差异是由于no-inline语句引起的 - Map需要进行一些额外的操作来从函数参数中提取数组的地址。 除此之外,它们完全做同样的事情。将map移动到函数内部并将数组传递给它 - 结果将是相同的。我正在添加代码示例以演示我的意思。 - Anycorn
哇,是的,你说得对!!我尝试了你的代码,确实看到了区别。所以似乎将Map作为函数参数传递是一个不好的主意,除非该函数被内联? 我仍然不明白你的声明“Map需要一些额外的操作来从函数参数中提取数组的地址”的含义。你能详细解释一下吗?谢谢! - Alejandro
我宁愿补丁编译器,这会更有用。遗憾的是,编译器仍然无法正确利用可用的众多寄存器。 - ggael
嗨,关于汇编函数,我只想简短地评论一下:它依赖于使用_mm_load_pd对齐的数据,并且当数据未对齐时可能会导致崩溃,因此在创建用作函数参数的变量时要小心,并使用例如align16(...)来创建它们。 - gilgamash
显示剩余2条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接