Expression Template
这类技术整体上是Lazy Evaluation的思想, 可以看作在C++语法的层次上通过模板类型推导，静态地构建出待求值的矩阵运算表达式的表达式树，然后等到整个表达式构建完成，需要赋值给一个Tensor的时候，再对表达式树惰性求值。
基本原理
MShadow的官方guide提供了一个非常赞的教程，这里以自己的思路直接分析一下。
:::C++
template<typename SubType>
struct Exp {
// returns const reference of the actual type of this expression
inline const SubType& self(void) const {
return *static_cast<const SubType*>(this);
}
};
首先我们构建一个表达式的基类Exp，它是一个模板基类，唯一的功能就是提供一个方法，返回实际子类类型的引用，就是一个基类指针自动转子类指针的功能。
:::C++
template<typename OP, typename TLhs, typename TRhs>
struct BinaryMapExp: public Exp<BinaryMapExp<OP, TLhs, TRhs>> {
const TLhs& lhs;
const TRhs& rhs;
BinaryMapExp(const TLhs& lhs, const TRhs& rhs)
:lhs(lhs), rhs(rhs) {}
// evaluation function, evaluate this expression at position i
inline float Eval(int i) const {
return OP::Map(lhs.Eval(i), rhs.Eval(i));
}
};
然后我们定义二元运算符的类型，还是按照模板泛型的套路，这个模板接受三个类型参数，分别为当前二元运算的操作OP（加减乘除maxmin）和左右operand。然后它的实现实际上会调用OP::Map来进行二元运算。
那么比如乘法类型的运算就可以这样来定义：
:::C++
struct mul {
inline static float Map(float a, float b) {
return a * b;
}
};
此时我们就有足够的工具来定义抽象的二元操作了。但此时用上述的模板类，我们只能够定义抽象的表达式树，却还无法对它进行求值，因为Exp可并没有实现Eval方法，他只是抽象的表达式概念，无法拿来做求值的。所以此时我们引入定义了这个方法的可求值的实体Vec对象：
:::C++
struct Vec: public Exp<Vec> {
int len;
float* dptr;
Vec(void) {}
Vec(float *dptr, int len)
: len(len), dptr(dptr) {}
// here is where evaluation happens
template<typename EType>
inline Vec& operator=(const Exp<EType>& src_) {
const EType &src = src_.self();
for (int i = 0; i < len; ++i) {
dptr[i] = src.Eval(i);
}
return *this;
}
// evaluation function, evaluate this expression at position i
inline float Eval(int i) const {
return dptr[i];
}
};
增加辅助模板函数和运算符重载来作为语法糖，让表达式写起来更加自然一些：
:::C++
template<typename OP, typename TLhs, typename TRhs>
inline BinaryMapExp<OP, TLhs, TRhs>
F(const Exp<TLhs>& lhs, const Exp<TRhs>& rhs) {
return BinaryMapExp<OP, TLhs, TRhs>(lhs.self(), rhs.self());
}
template<typename TLhs, typename TRhs>
inline BinaryMapExp<mul, TLhs, TRhs>
operator*(const Exp<TLhs>& lhs, const Exp<TRhs>& rhs) {
return F<mul>(lhs, rhs);
}
最后我们就可以用如下的简洁语法来书写向量运算表达式，并且完全不用分配额外的空间，因为对表达式的解析和运算合并都在静态推导的过程中完成了。
:::C++
const int n = 3;
int main(void) {
float sa[n] = {1, 2, 3};
float sb[n] = {2, 3, 4};
float sc[n] = {3, 4, 5};
Vec A(sa, n), B(sb, n), C(sc, n);
// run expression, this expression is longer:)
A = B * F<maximum>(C, B);
for (int i = 0; i < n; ++i) {
printf("%d:%f == %f * max(%f, %f)\n",
i, A.dptr[i], B.dptr[i], C.dptr[i], B.dptr[i]);
}
return 0;
}
开了优化的话，编译出来的二进制代码不包含上述的对象构造之类的东西，就约等于是一个for循环把结果算出来而已，甚至循环都会直接被展开。祭出反汇编器，我们会发现上述代码实际上被编译成了如下的伪代码：
:::C++
int _main() {
xmm0 = intrinsic_movsd(xmm0, *0x100000f60);
xmm2 = intrinsic_movsd(xmm2, *0x100000f68);
xmm1 = intrinsic_movsd(xmm1, *0x100000f70);
xmm3 = intrinsic_movapd(xmm3, xmm1);
printf("%d:%f == %f * max(%f, %f)\n", 0x0, rdx, rcx, r8, r9);
xmm0 = intrinsic_movsd(xmm0, *0x100000f78);
xmm2 = intrinsic_movsd(xmm2, *0x100000f80);
xmm1 = intrinsic_movsd(xmm1, *0x100000f68);
xmm3 = intrinsic_movapd(xmm3, xmm1);
printf("%d:%f == %f * max(%f, %f)\n", 0x1, rdx, rcx, r8, r9);
intrinsic_movsd(xmm0, *0x100000f88);
intrinsic_movsd(xmm2, *0x100000f90);
intrinsic_movapd(xmm3, intrinsic_movsd(xmm1, *0x100000f80));
printf("%d:%f == %f * max(%f, %f)\n", 0x2, rdx, rcx, r8, r9);
return 0x0;
}