魔法师の高塔

表达式模板

表达式模板是诸如Eigen、GSL、Mshadow等C++矩阵库所使用的主要技术。

很普遍的一个问题,如何计算两个向量相加?$\vec c = \vec a + \vec b$
在C++里,可以用数组来表示这两个向量,那么我们可以写出如下代码。

1
2
3
4
5
void add(const float * a, const float * b,float* c){
for (size_t i = 0; i < N; i++) {
c[i] = a[i] + b[i];
}
}

自然地,如果有一个Vec类来描述这个向量,使用重载符号等技术,我们可以将其改为

1
c = a+b;

示例代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
struct Vec {
int len;
float* dptr;
Vec(int len) : len(len) {
std::cout<<__func__<<std::endl;
dptr = new float[len];
}
Vec(float * dptr_,int len) : len(len) {
std::cout<<__func__<<std::endl;
dptr = new float[len];
memcpy(dptr, dptr_, sizeof(float) * len);
}
Vec(const Vec& src) : len(src.len) {
std::cout<<__func__<<std::endl;
dptr = new float[len];
memcpy(dptr, src.dptr, sizeof(float) * len );
}
inline Vec &operator=(const Vec &src) {
dptr = new float[src.len];
for (int i = 0; i < len; ++i) {
dptr[i] = src.dptr[i];
}
return * this;
}
~Vec(void) {
delete [] dptr;
std::cout<<__func__<<std::endl;
}
};
inline Vec operator+(const Vec &lhs, const Vec &rhs) {
Vec res(lhs.len);
for (int i = 0; i < lhs.len; ++i) {
res.dptr[i] = lhs.dptr[i] + rhs.dptr[i];
}
return res;
}

main函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
const int n = 3;
int main(){
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);
A = B + C + C;
for (int i=0; i<n; ++i) {
std::cout<<i<<" "<<A.dptr[i]<<"="<<B.dptr[i]<<"+"<<C.dptr[i]<<"+"<<C.dptr[i]<<std::endl;
}
return 0;
}

可以看到,在式子A=B+C+C这个计算过程中,产生了两个中间变量,生成这两个变量,需要分配空间,消耗资源,在要求效率的前提下,这样显然不能令人满意。
一个可行的解决方法,是让计算延迟到需要求值时。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
struct Exp{
virtual float Eval(int i) const = 0;
};
template <typename T1, typename T2>
struct BinaryAddExp : public Exp {
const T1 & lhs;
const T2 & rhs;
BinaryAddExp(const T1 & lhs, const T2 & rhs):lhs(lhs),rhs(rhs){}
inline float Eval(int i) const{
return lhs.Eval(i)+rhs.Eval(i);
}
};
struct Vec : public Exp {
int len;
float * dptr;
Vec() = default;
Vec(float * dptr, int len):len(len),dptr(dptr){
std::cout<<__func__<<std::endl;
}
inline float Eval(int i) const{
return dptr[i];
}
inline Vec & operator = (const Exp& src_) {
for (int i = 0; i < len; ++i) {
dptr[i] = src_.Eval(i);
}
return * this;
}
~Vec(){
std::cout<<__func__<<std::endl;
}
};
template <typename T1, typename T2>
inline BinaryAddExp<T1,T2> operator+(const T1 & lhs, const T2 & rhs){
return BinaryAddExp<T1,T2>(lhs,rhs);
}

这就是表达式模板。如果考虑伸缩复用,基类Exp不应该直接定义Eval函数,可以考虑支持一个返回子类类型的函数。
可以修改为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
template <typename T>
struct Exp{
inline const T& self() const {
return * static_cast<const T*>(this);
}
};
template <typename T1, typename T2>
struct BinaryAddExp : public Exp<BinaryAddExp<T1, T2>> {
const T1 & lhs;
const T2 & rhs;
BinaryAddExp(const T1 & lhs, const T2 & rhs):lhs(lhs),rhs(rhs){}
inline float Eval(int i) const{
return lhs.Eval(i)+rhs.Eval(i);
}
};
struct Vec : public Exp<Vec> {
int len;
float * dptr;
Vec() = default;
Vec(float * dptr, int len):len(len),dptr(dptr){
std::cout<<__func__<<std::endl;
}
inline float Eval(int i) const{
return dptr[i];
}
template <typename T>
inline Vec& operator = (const Exp<T>& src_) {
const T & src = src_.self();
for (int i = 0; i < len; ++i) {
dptr[i] = src.Eval(i);
}
return * this;
}
~Vec(){
std::cout<<__func__<<std::endl;
}
};
template <typename T1, typename T2>
inline BinaryAddExp<T1,T2> operator+(const T1 & lhs, const T2 & rhs){
return BinaryAddExp<T1,T2>(lhs,rhs);
}

考虑更多计算结构
BinaryAddExp进一步泛化为BinaryMapExp,

1
2
3
4
5
6
7
8
9
template <typename Op, typename T1, typename T2>
struct BinaryMapExp : public Exp<BinaryMapExp<Op,T1, T2>> {
const T1 & lhs;
const T2 & rhs;
BinaryMapExp(const T1 & lhs, const T2 & rhs):lhs(lhs),rhs(rhs){}
inline float Eval(int i) const{
return Op::Map(lhs.Eval(i),rhs.Eval(i));
}
};

设计Op

1
2
3
4
5
6
7
8
9
10
struct add{
inline static float Map(float a, float b){
return a+b;
}
};
struct mul{
inline static float Map(float a, float b){
return a*b;
}
};

重载符号

1
2
3
4
5
6
7
8
template <typename T1, typename T2>
inline BinaryMapExp<add,T1,T2> operator+(const T1 & lhs, const T2 & rhs){
return BinaryMapExp<add,T1,T2>(lhs,rhs);
}
template <typename T1, typename T2>
inline BinaryMapExp<mul,T1,T2> operator*(const T1 & lhs, const T2 & rhs){
return BinaryMapExp<mul,T1,T2>(lhs,rhs);
}

这就是表达式模板的基础知识。
进一步优化,我们可以加入标量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#include <iostream>
template <typename T>
struct Trait{
typedef T const & ExprRef;
};
template <typename T>
struct Exp{
inline const T& self() const {
return * static_cast<const T*>(this);
}
};
template <typename T>
struct Scalar :Exp<Scalar<T>>{
T s;
Scalar(T v):s(v){}
inline T operator[](size_t idx) const{
return s;
}
};
template <>
struct Trait<float> {
typedef Scalar<float> ExprRef;
};
template <>
struct Trait<int> {
typedef Scalar<int> ExprRef;
};
template <>
struct Trait<double> {
typedef Scalar<double> ExprRef;
};
template <typename Op,typename T1, typename T2>
struct BinaryMapExp : public Exp<BinaryMapExp<Op,T1, T2>> {
typename Trait<T1>::ExprRef lhs;
typename Trait<T2>::ExprRef rhs;
BinaryMapExp(const T1 & lhs, const T2 & rhs):lhs(lhs),rhs(rhs){
std::cout<<__func__<<std::endl;
}
inline float operator[](size_t i) const{
return Op::Map(lhs[i],rhs[i]);
}
};
struct Vec : public Exp<Vec> {
size_t len;
float * dptr;
Vec() = default;
Vec(float * dptr, size_t len):len(len),dptr(dptr){
std::cout<<__func__<<std::endl;
}
Vec(size_t len):len(len),dptr(new float(len)){}
inline float operator[](size_t i) const{
return dptr[i];
}
template <typename T>
inline Vec& operator = (const T& src_) {
for (int i = 0; i < len; ++i) {
dptr[i] = src_[i];
}
return * this;
}
~Vec(){
std::cout<<__func__<<std::endl;
}
};
struct add {
inline static float Map(float a,float b){
return a+b;
}
};
struct mul {
inline static float Map(float a,float b){
return a*b;
}
};
template <typename T1, typename T2>
inline BinaryMapExp<add,T1,T2> operator+(const T1 & lhs, const T2 & rhs){
return BinaryMapExp<add,T1,T2>(lhs,rhs);
}
template <typename T1, typename T2>
inline BinaryMapExp<mul,T1,T2> operator*(const T1 & lhs, const T2 & rhs){
return BinaryMapExp<mul,T1,T2>(lhs,rhs);
}
const int n = 3;
int main(){
float sb[n] = {2,3,4};
float sc[n] = {3,4,5};
Vec A(n),B(sb,n),C(sc,n);
A = B + C * 3.0 + B*C;
for (int i=0; i<3; i++) {
std::cout<<A[i]<<std::endl;
}
return 0;
}

在DL中,模板表达式的用途就是展开整个计算图,然后可以利用一些手段对其优化。这方面接触不多,还得继续学习。

学习使我忘忧:-D