详解eigen中的表达式模板_eigen 模板元编程-CSDN博客

本文链接：https://blog.csdn.net/weixin_44604698/article/details/136479441

1.起因

试想一种这样的情况：

    Matrix<double, 4> a;
    Matrix<double, 4> b;
    Matrix<double, 4> c;
    Matrix<double, 4> d;
    for(int i = 0; i < 4; i++){
        a[i] = i;
        b[i] = i;
        c[i] = i + 2;
    }
    d = a + b - c;

常规解法是：先计算a+b，然后存到一个临时Matrix中，然后在计算减去c，最后将该值赋给d。这种解法的代价是：

多构造一个临时空间，存储a+b的值；
要做两遍for循环，加法和减法各一遍；

可不可以实现：

for(int i = 0; i < 4; i++){
    d[i] = a[i] + b[i] - c[i];
}

这就用到了我们今天要提到了表达式模板，我为了方便大家理解抽象eigen中的核心组件；具体代码位置为：https://github.com/Beichen-Wang/HPCTest/tree/master/SystemOptimizer/ExpressionTemplate

2.核心思想

核心思想是利用模板元编程，在编译时展开，重载operator=，在赋值时进行计算。

3.基本组件

主要组件就是这5个class:

Matrix：基础的矩阵类，主要用于存储矩阵元素，其中getR和get分别返回对应位置的元素；
BinaryOps：主要用于存储加法和减法的返回类，不直接计算，在evaluator中计算；
MatrixBase：核心基类，重载operator+、operator-和operator=，加法和减法不进行计算，最后的返回值为BinaryOps，在operator=的时候具体计算；
evaluator<Matrix>：通过getR得到Matrix的引用用于计算；

evaluator<BinaryOps>::get(intdex)：计算加法和减法的位置，最后得到具体的计算值；

    typename Trait<TL>::Scalar evaluator<BinaryOps<BinaryOP, TL, TR>>::get(int Index) {
        return m_data.func()(m_data.left().get(Index), m_data.right().get(Index));
    }

3.1operator=

重点介绍一下operator=，核心是调用了call_assignment，我这里面以引用的方式将evaluator作为src和dst，在general_assignment进行了赋值；

template <typename Dst, typename Src>
inline void call_assignment(Dst & dst, Src & src){
    using EDstType = evaluator<Dst>;
    using ESrcType = evaluator<Src>;
    EDstType EDst(dst);
    ESrcType ESrc(src);
    using Kernel = general_assignment<EDstType, ESrcType>;
    Kernel kernel(EDst, ESrc);
    Assignment<Kernel, 0, Trait<Dst>::Size>::run(kernel);
}

这里面值得注意的是：我运用到了模板类，这样就可以将赋值过程循环展开，也就是变成了：

d[0] = a[0] + b[0] - c[0];
d[1] = a[1] + b[1] - c[1];
d[2] = a[2] + b[2] - c[2];
d[3] = a[3] + b[3] - c[3];

算是实现了一点小的优化，当然eigen中关于这部分的优化更为完备，通过自己的成本模型判断出不同的unrolling type；

enum UnrollingType {
  /** \internal Do not unroll loops. */
  NoUnrolling,
  /** \internal Unroll only the inner loop, but not the outer loop. */
  InnerUnrolling,
  /** \internal Unroll both the inner and the outer loop. If there is only one loop, 
    * because linear traversal is used, then unroll that loop. */
  CompleteUnrolling
};

3.2evaluator

我刚开始很好奇，为什么eigen要搞这么麻烦，单独抽象一个evaluator类，像一般的表达式模版，直接特化operator[]不是来的简单易懂，后来我发现单独抽象evaluator类，可以实现优化，尤其涉及到矩阵乘法，试想这样一种情况：scalar * ( A * B )，理论上这种情况就需要临时变量存储A*B，但是如果遇到这种情况，改成(A*scalar) * B就不需要了。通过特化evaluator就可以实现这种方式，具体实现方式可以查看：https://github.com/libigl/eigen/blob/1f05f51517ec4fd91eed711e0f89e97a7c028c0e/Eigen/src/Core/ProductEvaluators.h#L40

4.总结

通过参考eigen，实现了自己的一个表达式模板，目前可以完成加法和减法运算，并预留了其他运算的接口。具体链接为：https://github.com/Beichen-Wang/HPCTest/tree/master/SystemOptimizer/ExpressionTemplate。其中用到的优化技巧有：