Chinaunix首页 | 论坛 | 博客
  • 博客访问: 537166
  • 博文数量: 576
  • 博客积分: 40000
  • 博客等级: 大将
  • 技术积分: 5020
  • 用 户 组: 普通用户
  • 注册时间: 2008-10-13 14:47
文章分类

全部博文(576)

文章存档

2011年(1)

2008年(575)

我的朋友

分类:

2008-10-14 14:55:08


利用模板元编程实现解循环优化
Unroll Loops Optimization with Template Meta Programming

作者:王鹏

简介

在《C++ Templates: The Complete Guide》一书中(以下简称书),提出了模板元编程最早的实际应用之一:在数值运算中进行解循环优化。
而本文的标题是噱头!本文的真正目的是指出这种优化措施在增加复杂性的同时,并不一定能明显改善效率。应当谨慎使用该技术——默认不使用该技术,在点积计算确实是效率瓶颈时考虑采用该技术,并认真测试该技术是否真能提高效率。

背景

数值运算库中常常需要提供向量点积(dot_product)运算。其定义用C++代码描述也许更清楚~

template
T dot_product(int dim,const T v1[],const T v2[]) {
	T result=0;
	for (int i=0;i我们可以使用这个函数,求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
struct DotProduct {
	static T execute(const T v1[],const T v2[]);
}
template
T DotProduct::execute(const T v1[],const T v2[]) {
	return v1[0]*v2[0] + DotProduct::execute(v1+1,v2+1);
}

// 递归边界模板
template
struct DotProduct<1,T> {
	static T execute(const T v1[],const T v2[]);
}

template
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
T dot_product(const T v1[],const T v2[]) {
	return DotProduct::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<

程序输出当然是2个32,让我们进入反汇编看看实际工作情况:
令人吃惊的是,前4行代码根本没有生成任何机器码
生成机器码的语句只有2条输出语句和return 0;
cout< >::operator<< (40203Ch)]
mov	ecx,eax
call	dword ptr [__imp_std::basic_ostream >::operator<< (402038h)]
这些指令的的含义是什么呢?

我们知道 cout< (cout< cout.opeator<<(...) 返回的是自身引用,所以可以进行链式输出

前4条指令执行后:
堆栈应该是:
std::endl
20h
ECX 保存的是 std::cout

std::cout 的类型是std::basic_ostream >
所以第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<

这个原本就短小的程序,被编译器“篡改”为更短小的程序:

int main() {
	cout<<32<如果改用 printf 输出:
printf(“%d\n”,r1); printf(“%d\n”,r2);
从反汇编可以更清楚的看出,编译器实际做出的代码是:
printf(“%d\n”,0x20); printf(“%d\n”,0x20);

(详细代码见 sample1)

比较合理的解释是,编译器知道了太多的上下文,做出了足够的优化,将本来应该运行时得出的结果—— 一个不变的结果 ——放入的可执行文件中,直接输出。

为了验证模板元解循环技术,我们要另外想办法。

验证2:

我们需要“阻挠”编译器,让它不能猜测出运行结果,从而生成真正的计算代码。
这里使用了一个比较简单的方案:

template
void __stdcall MyGenerate(OutIt first,OutIt last,const char *name) {
	cout<<"generating "< oit(cout," ");
	copy(first,last,oit);
	cout<
  

用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<
  

仍然不够, 编译器会将r2的计算,推迟到 cout<

所以我们再加入一个函数

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(v1,v2);
	ForceCalcResult_non_inline(r2);

	cout<这样,编译器就屈服了,乖乖的把r1和r2的计算代码展现出来。
//MyGenerate(v1,v1+dim,"v1");
mov		eax,offset string "v1" (402134h) 
lea		edi,[esp+24h] 
lea		ebx,[esp+18h] 
call		MyGenerate (401220h) 
//MyGenerate(v2,v2+dim,"v2");
mov		eax,offset string "v2" (402138h) 
lea		edi,[esp+18h] 
lea		ebx,[esp+0Ch] 
call		MyGenerate (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(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 (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 (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和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条


--------------------next---------------------

阅读(226) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~