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

Question

是否可以实现一个接收 C 风格指针作为模板参数并以某种方式解析为静态特征矩阵但使用提供的内存的类？假设一个声明看起来像这样:

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

我确实了解映射，但我在下面提供的示例代码已证明，与静态特征矩阵相比，它们的速度要慢 20%。

前提是:

• 我需要提供我自己的 C 指针。这样我就可以高效地重用 C 代码而不会产生拷贝。
• 生成的矩阵对于 Eigen 而言应该是静态的，以便 Eigen 可以像在编译时使用静态数组一样进行优化。看看我上面的例子，在编译时我会提供矩阵大小(静态)和 C 指针。
• CMatrix 应该退回到 Eigen::Matrix。如果未提供 C 数组的附加模板参数，我将获得正常的特征矩阵。
• 我不打算进行完整的 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;
}

Answer 1

有点离题，但既然你强调了性能:

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];                                                                                                                                                                              
  }
}

您可能会通过交错加载和 mul/adds 获得一些额外的改进 - 我尽量保持简单。

结果如下:

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