Unroll Loops Optimization with Template Meta Programming
作者：王鹏下载源代码
简介 逆风编程技术
在《C++ Templates: The Complete Guide》一书中（以下简称书），提出了模板元编程最早的实际应用之一：在数值运算中进行解循环优化。
而本文的标题是噱头！本文的真正目的是指出这种优化措施在增加复杂性的同时，并不一定能明显改善效率。应当谨慎使用该技术——默认不使用该技术，在点积计算确实是效率瓶颈时考虑采用该技术，并认真测试该技术是否真能提高效率。
背景
数值运算库中常常需要提供向量点积（dot_product）运算。其定义用C++代码描述也许更清楚～
template<typename T>
T dot_product(int dim,const T v1[],const T v2[]) {
T result=0;
for (int i=0;i<dim;++i)
result += v1[i]*v2[i];
return result;
}
我们可以使用这个函数，求2个向量的点积
int v1[] = {1,2,3};
int v2[] = {4,5,6};
int r1 = dot_product(3,v1,v2);
得到r1=32
书中指出：“这个结果是正确的，但是在要求高效率的应用中，它耗费了太多的时间”，对于这类特殊的问题，“简单的把循环展开”，如：r1=v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]
“反而会好得多”。
如何便捷的展开循环？将每个dot_product(dim,...) 手工重写成展开式？还是设计一整套函数: dot_product_dim_2，dot_product_dim_3，... ？
无疑这是一个冗长乏味易出错的方案。
书中提出了一种使用模板元编程解决该问题的方案，让我们来看看他是怎么做的。
模板元编程——解循环
// 首先是递归模板
template<int DIM,typename T>
struct DotProduct {
static T execute(const T v1[],const T v2[]);
}
template<int DIM,typename T>
T DotProduct<DIM,T>::execute(const T v1[],const T v2[]) {
return v1[0]*v2[0] + DotProduct<DIM-1,T>::execute(v1+1,v2+1);
}
// 递归边界模板
template<typename T>
struct DotProduct<1,T> {
static T execute(const T v1[],const T v2[]);
}
template<typename T>
T DotProduct<1,T>::execute(const T v1[],const T v2[]) {
return v1[0]*v2[0];
}
我们可以使用DotProduct 来进行点积计算了：
int v1[] = {1,2,3}; int v2[] = {4,5,6};
int r2 = DotProduct<3,int>::execute(v1,v2);
计算r2的函数，在DotProduct<3,int>::execute 被实例化时就确定了。
编译器将会实例化
int DotProduct<3,int>::execute(const int v1[],const int v2[]) {
return v1[0]*v2[0] + DotProduct<2,int>::execute(v1+1,v2+1);
}
然后编译器继续实例化
int DotProduct<2,int>::execute(const int v1[],const int v2[]) {
return v1[0]*v2[0] + DotProduct<1,int>::execute(v1+1,v2+1);
}
这里，我们有一个 DotProduct<1,T> 的偏特化版本，计算函数将被实例化为
int DotProduct<1,int>::execute(const int v1[],const int v2[]) {
return v1[0]*v2[0];
}
而这3个函数都足够短小，编译器通常会为它们做inline工作，使得r2的计算被展开为
r2 = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
DotProduct<3,int>::execute(v1,v2) 的调用语法不是很友好。
书中还设计了一个helper函数。
template<int DIM,typename T>
T dot_product(const T v1[],const T v2[]) {
return DotProduct<DIM,T>::execute(v1,v2);
}
如同STL的make_pair、bind2nd，这个模板函数将模板参数T的推导工作交给编译器，使DotProduct的界面更友好，更不容易出错。当然，维度DIM必须指定。
现在可以按如下方式使用 ：
int r2 = dot_product<3>(v1,v2); // 是不是和 int r1 = dot_product(3,v1,v2); 十分相似？
本文以下部分将从编译器生成的机器码层次来验证这一技术，同时与普通循环点积计算做比较。
测试环境 ：Windows Xp sp2，Microsoft Visual Studio 2005 Team Suite，release默认设置
验证1：
该技术真的能和预想的一样，实例化多个函数，并将它们inline，达到解循环的效果吗？
让我们写个短小程序验证一下：
int main() {
int v1[] = {1,2,3};
int v2[] = {4,5,6};
int r1 = dot_product(3,v1,v2); //循环版本
int r2 = dot_product<3>(v1,v2); //解循环版本
cout<<r1<<endl;
cout<<r2<<endl;
return 0;
}
程序输出当然是2个32，让我们进入反汇编看看实际工作情况：
令人吃惊的是，前4行代码根本没有生成任何机器码
生成机器码的语句只有2条输出语句和return 0;
cout<<r1<<endl; // 生成7条指令
mov eax,dword ptr [__imp_std::endl (402040h)
mov ecx,dword ptr [__imp_std::cout (402044h)]
push eax
push 20h
call dword ptr [__imp_std::basic_ostream<char,std::char_traits<char> >::operator<< (40203Ch)]
mov ecx,eax
call dword ptr [__imp_std::basic_ostream<char,std::char_traits<char> >::operator<< (402038h)]
这些指令的的含义是什么呢？
我们知道 cout<<r1<<endl; 的完整形式是
(cout<<r1)<<endl; 即是 cout.operator<<(r1).operator<<(endl);
cout.opeator<<(...) 返回的是自身引用，所以可以进行链式输出
前4条指令执行后：
堆栈应该是：
std::endl
20h
ECX 保存的是 std::cout
std::cout 的类型是std::basic_ostream<char,std::char_traits<char> >
所以第5条指令正是对cout的成员函数“ opeator<<”进行调用
(this指针 std::cout,已经存入ECX，而且可以猜想这个成员函数正是 operator(int) 重载)
而参数，对的，是call指令前，堆栈栈顶“20h”—— 32的16进制。
第5条指令将执行 cout.operator<<(20h); ！
让我们把 r1被编译器偷偷替换成20h的事放在一边，继续解析后2条指令：
cout.operator<< 使用__thiscall调用约定，函数自己弹出堆栈中参数20h， 返回自身引用，返回值存入EAX。
所以第5条指令执行后：
堆栈应该是：
std::endl;
EAX 保存的是返回值，即是std::cout
第6条指令 mov ecx,eax 将返回值重新放入ECX，设置好this指针。
第7条指令 继续调用 operator<< 的另一个重载版本。（这时 std::endl 在栈顶）
完成语句 cout.operator<<(endl);
cout<<r2<<endl; 生成了相似的7条指令。
这个原本就短小的程序，被编译器“篡改”为更短小的程序：
int main() {
cout<<32<<endl;cout<<32<<endl;
return 0;
}
如果改用 printf 输出：
printf(“%d\n”,r1); printf(“%d\n”,r2);
从反汇编可以更清楚的看出，编译器实际做出的代码是：
printf(“%d\n”,0x20); printf(“%d\n”,0x20);
（详细代码见 sample1）
比较合理的解释是，编译器知道了太多的上下文，做出了足够的优化，将本来应该运行时得出的结果—— 一个不变的结果 ——放入的可执行文件中，直接输出。
为了验证模板元解循环技术，我们要另外想办法。
验证2：
我们需要“阻挠”编译器，让它不能猜测出运行结果，从而生成真正的计算代码。
这里使用了一个比较简单的方案：
template<class OutIt>
void __stdcall MyGenerate(OutIt first,OutIt last,const char *name) {
cout<<"generating "<<name<<" with clock\n";
generate(first,last,clock);
ostream_iterator<int> oit(cout," ");
copy(first,last,oit);
cout<<endl;
}
用clock 函数返回值填充一个序列，即使编译器“胆大包天”，也不敢猜测clock返回结果了～
多余的 name 参数和一些输出语句是因为在release模式下，watch窗口中看不到值。
而__stdcall 是让 MyGenerate后不生成多余的平衡堆栈代码。
按照如下方式使用 ：
MyGenerate(v1,v1+3,”v1”);
MyGenerate(v2,v2+3,”v2”);
int r1 = dot_product(3,v1,v2);
MyGenerate(v1,v1+3,”v1”);
MyGenerate(v2,v2+3,”v2”);
int r2 = dot_product<3>(v1,v2);
cout<<r1<<endl;
cout<<r2<<endl;
仍然不够， 编译器会将r2的计算，推迟到 cout<<r2<<endl;中，使得r2的计算不够清晰。
所以我们再加入一个函数
void ForceCalcResult(int result) { }
强制编译器立刻计算点积。
同时，使用函数指针调用，避免编译器对该函数inline。
程序主干如下 ：
int main() {
const int dim = 3;
int v1[dim];
int v2[dim];
void ( *const ForceCalcResult_non_inline)(int) = ForceCalcResult;
MyGenerate(v1,v1+dim,"v1");
MyGenerate(v2,v2+dim,"v2");
int r1 = dot_product(dim,v1,v2);
ForceCalcResult_non_inline(r1);
MyGenerate(v1,v1+dim,"v1");
MyGenerate(v2,v2+dim,"v2");
int r2 = dot_product<dim>(v1,v2);
ForceCalcResult_non_inline(r2);
cout<<r1<<endl;
cout<<r2<<endl;
return 0;
}
这样，编译器就屈服了，乖乖的把r1和r2的计算代码展现出来。
//MyGenerate(v1,v1+dim,"v1");
mov eax,offset string "v1" (402134h)
lea edi,[esp+24h]
lea ebx,[esp+18h]
call MyGenerate<int *> (401220h)
//MyGenerate(v2,v2+dim,"v2");
mov eax,offset string "v2" (402138h)
lea edi,[esp+18h]
lea ebx,[esp+0Ch]
call MyGenerate<int *> (401220h)
从这4条指令，我们可以看出v1和v2的地址：
v1[0]:esp+18h, v1[1]:esp+1Ch, v1[2]:esp+20h, v1[dim]:esp+24h
v2[0]:esp+0Ch, v2[1]:esp+10h, v2[2]:esp+14h, v2[dim]:esp+18h
让我们先看模板解循环版本：
//int r2 = dot_product<dim>(v1,v2);
mov edi,dword ptr [esp+10h]
mov edx,dword ptr [esp+14h]
imul edi,dword ptr [esp+1Ch]
imul edx,dword ptr [esp+20h]
// edi = v2[1]*v1[1]
// edx = v2[2]*v1[2]
mov eax,dword ptr [esp+0Ch]
imul eax,dword ptr [esp+18h]
// eax = v2[0]*v1[0]
add edi,edx
add edi,eax
// edi = edi+edx+eax = v2[1]*v1[1]+v2[2]*v1[2]+v2[0]*v1[0]
循环被解开了！利用3个寄存器edi,edx,eax，最终结果应该是保存在edi中。
再看接下来的代码
//ForceCalcResult_non_inline(r2);
push edi
call ForceCalcResult (401000h)
结果确实是存放在edi中的。
我们接着看普通“循环”版本：
//int r1 = dot_product(dim,v1,v2);
mov esi,dword ptr [esp+10h]
mov eax,dword ptr [esp+14h]
imul esi,dword ptr [esp+1Ch]
imul eax,dword ptr [esp+20h]
// esi = v2[1]*v1[1]
// eax = v2[2]*v1[2]
mov ecx,dword ptr [esp+0Ch]
imul ecx,dword ptr [esp+18h]
// ecx = v2[0]*v1[0]
add esi,eax
add esi,ecx
// esi = esi+eax+ecx = v2[1]*v1[1] + v2[2]*v1[2] + v2[0]*v1[0]
//ForceCalcResult_non_inline(r1);
push esi
call ForceCalcResult (401000h)
几乎相同的代码，使用的是esi,ecx,eax，结果保存在esi中。
循环同样被解开了！
编译器在我们的重重算计下，被迫在运行时计算结果，同时将运算方法展现出来；但同时，它依然不屈不饶的向我们展示它的强大优化能力！
（详细代码见sample2）
验证2+：
在验证2中，编译器知道了 const int dim = 3; 这一上下文。从而将普通循环版本的点积函数进行展开。
让我们来看看，dim取更大值时，编译器会如何。
dim=10：
普通“循环”点积计算被展开，使用esi,eax,ecx,edx，10条mov与imul指令，9条add指令，最终将计算结果存入esi
模板元点击计算也被展开，使用edi,eax,ecx,edx，10条mov与imul指令，9条add指令，最终将计算结果存入edi
dim=11：
循环点积计算发生了变化，代码如下：
//MyGenerate(v1,v1+dim,"v1");
00401017 mov eax,offset string "v1" (402134h)
0040101C lea edi,[esp+68h]
00401020 lea ebx,[esp+3Ch]
00401024 call MyGenerate<int *> (401280h)
//MyGenerate(v2,v2+dim,"v2");
00401029 mov eax,offset string "v2" (402138h)
0040102E lea edi,[esp+3Ch]
00401032 lea ebx,[esp+10h]
00401036 call MyGenerate<int *> (401280h)
// v1[0]:esp+3Ch, v2[0]:esp+10h
//int r1 = dot_product(dim,v1,v2);
0040103B xor ebp,ebp
// int result = ebp=0;
0040103D xor eax,eax
// eax=0;
0040103F nop
// align
// i=eax;
// loops:
00401040 mov ecx,dword ptr [esp+eax+10h]
// ecx = *(v2+i);
00401044 imul ecx,dword ptr [esp+eax+3Ch]
// ecx *= *(v1+i);
00401049 add eax,4
// ++i;
0040104C add ebp,ecx
// result+=ecx
0040104E cmp eax,2Ch
// if (i<11)
00401051 jl main+30h (401040h)
// goto loops;
编译器已经不能“忍受”这样长度的循环展开，将其老老实实的翻译成真正的循环。
但是，编译器仍然对函数做了inline工作。
对于模板元解循环，因为代码要求11层函数调用，编译器能做的只有将11层调用inline，得到的是就地展开的计算。
dim=12：
循环点积计算： 编译器甚至连inline也不做了，进行正式的函数调用。
模板元点击计算：依然就地展开。
dim=32：
循环点积计算：不展开，调用
模板元点积计算：
比较有趣的是，编译器也没有将函数就地展开，而是分成3步。就地计算前9个乘法然后与DotProduct<23,int>的返回值相加。DP<23>也没有就地计算全部，而是计算前11个，然后与DP<12>的返回值相加。
3步计算当中都没有循环。
值得注意的是，循环点积计算函数也没有很老实的按照源代码的方式进行计算，而是进行了不完全的循环展开，类似于如下代码：
template< typename T>
T dot_product(int DIM,const T v1[],const T v2[]) {
const int step = complierKnow;
T results[step] = {0};
int i=0;
for (;i+step-1<DIM;i+=step) {
results[0] += v1[i]*v2[i];
results[1] += v1[i+1]*v2[i+1];
...
results[step-1] += v1[i+step-1]*v2[i+step-1];
}
for (;i<DIM;++i)
results[0] += v1[i]*v2[i];
return results[0]+result[1]+...result[step-1];
}
DIM和step似乎没有什么固定的规律。
（详细代码见sample2，取不同的dim值即可。）
验证3：
sample2中，编译器对被计算的向量的上下文依然知情：v1,v2在main函数的栈中。
sample3做了小小的改动：dot_product_loop和dot_product_unloop将无从得知传递给它们的向量在何处。（当然，在main函数中，仍然使用函数指针来调用这2个函数，使得编译器不做inline）
sample3得到的结果和sample2基本相同，只是对向量寻址由esp+v1_offset,esp+v2_offset变为[esp+4],[esp+8]
（详细代码见sampl3）
验证4：
在sample3的基础上，通过控制台输入读取dim，使得编译对其不知情。
当然，在这种情况下，是无法使用模板元解循环的，因为它要求dim是编译时常量。
在这种情况下，循环点积计算终于被编译器翻译成真正的“循环”了。
总结：
循环版本有更大的灵活性：在dim为编译时常量时，编译器会根据其大小，进行完全解循环或部分解循环。同时也支持dim为运行时的值。
模板元版本必须使用编译时常量作为dim，而且总是完全解开循环。
循环版本与模板元版本最终都会使用相同次数的乘法与运算，区别在于乘法指令和跳转指令的数目。
模板元版本在dim较大的时候，可能会使用跳转指令，将部分工作交给更小维度的点积计算。但是乘法指令数目与维数一样。
循环版本可能会使用更多的跳转指令，使得乘法指令数目大大减少。
多出的跳转指令会占用运行时间，但也能减少目标代码体积。权衡应该使用那一种版本与权衡某个函数是否应该inline十分相似。
书中17章最后部分也提醒读者，不计后果的解循环并不总是能优化运行时间。
借用大师的观点来结束本文 —— 优化的第一原则是：不要优化。优化的第二原则（仅适用于专家）是：还是不要优化。再三测试，而后优化。
《C++编程规范》 第8条
| 
您的位置：逆风者
→ VC++
→ 正文
本文章所属分类：首页
→ VC++