C++科学计算技巧

第一章

1.1 C++ 编译器

1.1.1 安抚子标准的编译器

所有的C++编译器仅支持ISO/ANSI C++标准的不同子集。试图提供可以运行在不同平台上的C++库或程序是一个让人头疼的事情。

编译器中最大的兼容性问题是

  1. Template
  2. STL

    如果你需要写可以在不同平台下工作的代码,大部分时间你会在煎熬中度过。在Blitz++库中,我采用了扩展的预处理Kludges来解决不同的编译器。例如下面这段代码:


00001: #include <cmath>
00002: 
00003: namespace blitz{
00004:     double pow(double x,double y)
00005:     {return std::pow(x,y);}
00006: };

这段简单的代码在Blitz namespace(名字空间)中定义了一个函数pow,它调用了标准的pow函数。一些编译器不支持<cmath>头文件,所以需要替换为<math.h>。一些编译器支持,但是C math函数不在std的namespace中。甚至一些编译器根本就不支持namespace。我使用的解决方案如下:


00001: #ifdef BZ_MATH_FN_IN_NAMESPACE_STD
00002:   #include <cmath>
00003: #else
00004:   #include <math.h>
00005: #endif
00006: 
00007: #ifdef BZ_NAMESPACES
00008: double pow(double x,double y)
00009: {
00010:   return BZ_MATHFN_SCOPE(pow)(x,y);
00011: }
00012: #endif

这种解决方法很难看。不幸的是,这种情形仍会持续很多年,直到所有的编译器都满足标准。

Blitz++自带有一个脚本(bzconfig),它会运行许多小的测试程序来侦测出当前的编译器需要什么样的kludges。Luc Maisonobe已经将这些程序加入了Autoconf工具中。你也可以修改Blitz++中的bzconfig脚本来使用它。

另一种方法来避免多个平台支持的方法时仅支持gcc编译器。因为它在支持标准方面要好过很多编译器,并且在标准库实现上有较高的提升。可是它在优化上没有KAI c++编译器和一些特定厂商的编译器强。但对于大型的out-of-cache程序,gcc做的很好。(Out-out-cache 程序对于优化不敏感,因为大多数的cpu时间用在了内存上(Stalling for memory).

1.1.2 编译器一览

多数的c++编译器使用由EDG(Edison Design Group)写的前端。这个前端是一个很接近标准,健壮可靠,可以提供很有意义的错误信息。你可以在任何平台上找到EDG的身影。其中的一些是:SGI C++,KAI C++,DEC cxx,Gray C++,Intel C++,Portland Group C++,Tera C++。

大部分的c++编译器(Borland/Inprise,Microsoft,etc)不符合标准的程度令人惊讶。它们都关注在COM,DCOM,WFC等上。一点都不关心标准的支持。Microsoft's 6.0 C++编译器在编译简单的成员模版时都会经常的崩溃。在和Microsoft开发人员对话的过程中,他们说这是因为他们的客户不会关心是否符合标准,而是在关心一些WFC/COM等东西。也许这是对的。但是如果你在自己最喜欢的编译器上由于标准问题而受挫,那就因该向开发商投诉,否则他们是不会考虑这些的。

1.1.3 特定的C++优化

大部分的编译器在优化C++程序上令人失望。对于在c++中经常出现的嵌套的聚集和零时变量,常规的优化技巧需要修改。但是KAI编译器是一个例外,它的优化做的很好。另外一个是SGI编译器(虽然对于KAI,目前在零时变量消除上还做的不够好)

KAI在优化技术上有美国专利,导致其他厂商不能实现相同的优化技术。尽管如此,它还是给EPC++编译器提供了优化引擎的使用许可证。

1.2 编译时间

对Blitz++抱怨最多的是编译时间过长(以及其他的template C++ 库)。下面小节将介绍为何C++程序需要这么长的时间编译。

1.2.1 头文件

头文件是一个大问题。C++没有一个很有效的模块系统。导致大量的头文件在每一个编译模块都需要被解析和分析。例如, 包括<iostream>在有些编译器上会拉入成千上万行的代码。

Template 库的情况更加严重。因为所有的模版聚现都在编译期进行。所以除了函数声明,头文件还包括了函数定义。

例如,在我的桌面工作站用gcc编译Blitz++的 example/array.cpp例子需要37秒。如果你除去其中所有的代码,只留下#include<blitz/array.h>这一行,编译竟然还需要22秒。也就是是说,60%的时间都被用来编译头文件了!

预先编译头文件技术可以在一定程度上缓解这个问题,但是有时候他们也会让问题更加严重。在Blitz++中开启Precompiled Header可以产生大约64Mb的dumps,读取它会比它取代的原头文件华更多的时间。

本节告诉你要尽量保持你的Header越小越好。当开发一个库时,提供很多个小的头文件来在相应的情况下分别包含,而不应该仅仅提供一个包含所有程序库的超大头文件。

1.2.2 Prelinking

模版给编译器导致了一个严重的问题。什么时候他们才需要被聚现化?这通常有两种方式:

1.2.3 程序数据库方法-visual age c++

IBM的Visual Age C++ 编译器对于分别编译有一个有趣的技术。它把所有的代码保存在一个数据库中。由于完整的源代码都是可以在编译期间取得的(而不是仅仅一个编译单元的代码),VAC++ 不需要做prelinking或者担心分别编译问题。编译器做增量编译,但是仅针对代码改变的情况。这种技术可以大幅度的降低编译template c++ 库的时间。

1.2.4 二次/三次 模版算法

像Blitz++的库使用了大量的template。一些编译器在模版聚现上有自己的算法,但适应性不是很好。例如下面的代码:


00001: template<int N>
00002: void foo()
00003: {
00004:   for<N-1>();
00005: }
00006: 
00007: template<>
00008: void foo<0>()
00009: {}
00010: 

如果聚现foo<N>,那么foo<N-1>被聚现,同样foo<N-1>,...,直到foo<0>。编译需要多少时间?聚现的函数个数是N个,所以编译时间应该随N线性增长,对吗?

不幸的是对一些编译器并不是这样。简单的说,许多编译器自身保存一个链表,确定一个template是否被聚现过需要线性搜索。所以编译时间为O(N^2).这对于很少的template的程序不成问题,但对于像Blitz++那样有上千个template需要聚现的程序库就很致命。

同样的问题可以发生在优化算法上。比如别名检查。我在一些主流C++编译器上发现了O(N^2),O(N^3)的优化算法。

1.3 静态多态

1.3.1 虚函数是邪恶的吗?

在C++的世界里,多态意味着虚函数。这种多态是很有作用的设计方法。但是使用虚函数会有大的性能损失。


00001: class Matrix{
00002:   public:virtual double operator()(int i,int j)=0;
00003: };
00004: 
00005: class SymmetricMatrix:public Matrix{
00006:   public:virtual double operator()(int i ,int j);
00007: };
00008: 
00009: class UpperTriangularMatrix:public Matrix{
00010:   public:virtual double operator()(int i,int j);
00011: };

虚函数分派operator()来从矩阵中读取元素,这会破坏任何一种矩阵算法的性能。

1.3.2 解决方法A: 简单引擎(simple engines)

"engines"来自在LANL的POOMA团队。设计模式社区也有一个同样的名称。

想法是把动态多态(虚函数)替换为静态多态(template parameters)。这是一个采用了这种技术的Matrix的骨架。


00001: class Symmetric{
00002: //封装针对对称矩阵的存储信息
00003: };
00004: 
00005: class UpperTriangular{
00006: //封装针对三角矩阵的存储信息
00007: };
00008: 
00009: tempalte<class T_engine>
00010: class Matrix{
00011: private:
00012:  T_engine engine;
00013: };
00014: 
00015: //采用任何矩阵结构的函数
00016: template<class T_engine>
00017: double sum(Matrix<T_engine>& A);
00018: 
00019: //使用例子...
00020: Matrix<Symmetric> A;
00021: sum(A);
00022: 

在上例中,矩阵结构的不同变化隐藏在engines Symmetric和UpperTriangular中。Matrix类把engine作为template parameter,并且包含这个engine的实例作为成员变量。

用户使用的符号变了:Matrix<Symmetric> 而不是原来的SymmetricMatrix,不过这个不是一个大问题。我们可以使用typedefine来隐藏这个细节。同样,在Matrix中,大部分操作要被代理给engines来处理。

这种解决方法最严重的问题是所有的Matrix子类型都必须有同样的成员函数。例如,如果我们想要Matrix<Symmetric>有一个isSymmetricPositiveDefinite()的方法,典型的实现是:


00001: template<class T_engine>
00002: class Matrix{
00003: public:
00004:   //Delegate to the engine
00005:   bool isSymmetricPositiveDefinite()
00006:   {return engine.isSymmetricPositiveDefinite();}
00007: 
00008: private:
00009:   T_engine engine;
00010: };
00011: 
00012: class Symmetric{
00013: public:
00014:   bool isSymmetricPositiveDefinite()
00015:   {
00016:    ...
00017:   }
00018: };

但是这并没有使Matrix<UpperTriangular>有理由包含isSymmetricPositiveDefinite()函数,因为upper triangular矩阵不能为对称。从上例可以看出这种方法使基类包含了所有子类型所拥有的函数。所以导致产生了很多在engine中会抛出异常的成员函数。


00001: class UpperTriangular{
00002: public:
00003:   bool isSymmetricPositiveDefinite()
00004:   {
00005:     throw makeError("Method isSymmetricPositiveDefinite() is not defined for UpperTriangular矩阵");
00006:     return false;
00007:   }
00008: 
00009: };
00010: 

另一种方法:忽略这些函数定义会把你的用户推进火坑。因为当用户试图使用这些未定义的函数的时候,会得到一大堆奇怪的模版聚现错误。

可喜的是,有一种更好的解决方法。

1.3.3 解决方案B: Barton 和Nackman 技术

这种技术称为"Barton and Nackman Trick",因为 他们在他们的杰出著作<<Scientific and Engineering C++>>中使用了这种技术。Geoff取了一个比较有意义的名字"Curiously Recursive Template Pattern"。

如果你没有见过这种技巧,那么它一定会让你开始抓头的。是的,它是合法的。


00001: //基类有一个模版参数。这个参数定义了该类应该从哪里继承。
00002: template<class T_leaftype>
00003: class Matrix{
00004: public:
00005:  T_leaftype& asLeaf()
00006: {return static_cast<T_leaftype&>(*this);}
00007: 
00008: double operator()(int i,int j)
00009: {return asLeaf()(i,j);} //代理给leaf
00010: 
00011: }
00012: 
00013: class SymmetricMatirx : public Matrix<SymmetricMatrix>{
00014: ...
00015: };
00016: 
00017: class UpperTriMatrix: public Matrix<UpperTriMatrix>{
00018: ...
00019: };
00020: 
00021: // 使用例子,可以接受任何矩阵类
00022: template<class T_leaftype>
00023: double sum(Matrix<T_leaftype>& A);
00024: 
00025: //例子
00026: SymmetricMatrix A;
00027: sum(A);

这种技术可以实现一般化的继承体系结构。关键点在于基类有一个模版参数,该参数为后代(叶子)类的类型。保证了在编译期间该类型信息是可以得到的-不再需要虚函数的分派。

在继承类中可以任意选择性的进行定制(例如,基类可以实现默认功能,继承类重载相应的功能)

与上节解决方法不同,继承类可以有其特殊的,仅属于自己的(基类没有定义)成员函数。

这种方法,基类仍然要将一切调用代理给继承类。

使你的继承树多于一层是可行的,但需要多加考虑。

1.4 Callback inlining 技术

Callback 经常会出现在程序中。一些典型的例子如下


00001: double integrate(double a,double b,int numSamplePoints,double(*f)(double){
00002: double delta=(b-a)/numSamplePoints-1);
00003: double sum = 0.0;
00004: 
00005: for(int i=0;i<numSamplePoints;++i)
00006:   sum+=f(a+i*delta);
00007: 
00008: return sum*(b-a)/numSamplePoints;
00009: }

这种方法工作的很好,但是对性能有影响:在内部循环中的函数调用会导致较大的性能损失。一些优化编译器也许会执行inter-procedural和指针分析来确定f指向的确切函数,但是大多数编译器不会做。

1.4.1 Callbacks: C++方式

在C++库中,callbacks一般采用虚函数。新一点的方式是定义一个基类,如Ietegrator,定义一个虚函数,在继承子类中定义实际需要被积分的函数。


00001: class Integrator{
00002: public:
00003: 
00004: double integrate(double a,double b,int numSamplePoints)
00005: {
00006: double delta=(b-a)/..;
00007: ...
00008: 
00009: for(int i=0;i<numSamplePoints;++i)
00010:   sum+=functionToIntegrate(a+i*delta);
00011: reutrn sum*(b-a)/numSamplePoints;
00012: }
00013: 
00014: virtual double functionToIntegrate(double x)=0;
00015: };
00016: 
00017: class IntegrateMyFunction:Integrator{
00018: public:
00019:   virtual double functionToIntegrate(double x)
00020:   {
00021:      return cos(x);
00022:   }
00023: };

同C类似,虚函数分派在内部循环会导致性能下降(如果该函数执行时间比较长,则虚函数的损失可以忽略不计)。对于性能,需要在积分函数中采用inline the callback函数。

下面就介绍三种方法实现inlined callbacks。

  1. Expression templates

    这是一种精细,复杂的解决方法。它扩展和维护起来比较难。一般不必浪费你的精力去选择这种方法。

  2. stl风格的函数对象

    Callback函数可以被封装在一个类中,从而创建一个函数对象或者叫仿函数。


00001: class MyFunction{
00002: public:
00003: 
00004:    //注意:为inline 成员函数
00005:   double operator()(double x)
00006:   {return 1.0/(1.0+x);}
00007: };
00008: 
00009: template<class T_function>
00010: double integrate(double a,double b,int numSamplePoints,T_function f)
00011: {
00012:   double delta=(b-a)/(numSamplePoints-a);
00013:   double sum=0;
00014: 
00015:   for(int i=0;i<numSamplePoints;++i)
00016:     sum+=f(a+i*delta);
00017: 
00018:   return sum*(b-a)/numSamplePoints;
00019: 
00020: }
00021: 
00022: //用法
00023: integrate(0.3,0.5,30,MyFunction())
00024: 

一个MyFunction的对象作为参数传入integrate函数。C++编译器会聚现化integrate, 用MyFunction来替换T_function;这样会导致MyFunction::operator()被inlined进integrate()的内部循环中。

这种方法的缺点是用户必须将所有他想调用的函数封装在类中。同样,如果Callback函数需要参数,则必须提供给MyFunction的构造函数,并保存为成员变量。

注意,你可以向上例的integrate传入函数指针。但这和C风格的callback函数一样,不会在循环中将你传入的函数指针进行inline化。

  1. 函数指针作为模版参数 第三种方法很简单:

00001: double function1(double x)
00002: {
00003:   return 1.0/(1.0+x);
00004: }
00005: 
00006: template<double T_function(double)>
00007: double integrate(double a,double b,int numSamplePoints)
00008: {
00009: 
00010:   double delta=(b-a)/(numSamplePoints-1);
00011:   double sum=0.0;
00012: 
00013:   for(int i=0;i<numSamplePoints;++i)
00014:      sum+=T_function(a+i*delta);
00015: 
00016:   return sum*(b-a)/numSamplePoints;
00017: 
00018: }
00019: 
00020: //..
00021: integrate<function1>(1.0,2.0,100);

integrate()函数将指向callback函数的函数指针作为模版参数。由于这个指针在编译期间就可以确定,所以function1被inline进入内部循环。参考oon-list档案的Sumclass algorithms条目。

1.5 控制代码膨胀

代码膨胀对于基于模版的库来说是一个问题。例如,你使用queue的一些实例:


00001:   queue<float>
00002:   queue<int>
00003:   queue<Event>

你的程序会有三个不同版本的queue<>。这对于使用基于虚函数多态的老风格c++库并不是问题。

但是,实际中编译器会生成一些荒谬的多余代码使得问题更加严重。

1.5.1 避免kitchen-sink模版参数

将类的每个属性都作为模版参数是很有诱惑力的。下面就是一个夸张的例子。


00001: template<typename number_type,int N,intM,bool sparse,typename allocator,bool symmetric,bool banded,int bandwidthIfBanded>
00002: class Matrix;

每个多余的模版参数都为代码膨胀提供了机会。如果你使用了很多种"kitchen-sink"的模版类,你应该考虑通过虚函数,成员模版等减少一些模版参数。

例如,我们有一个类list<>。我们想要用户有能力提供自己的内存分配器,所以我们可以将Alloctor作为一个模板参数:


00001: list<double ,allocator<double>> A;
00002: list<double , my_alloctor<double>>B;

这样会为A和B生成不同的list版本,尽管两者大部分都是一模一样的(比如你要遍历list,你就不需要考虑它在内存中是怎么分配的)。采用alloactor作为模版参数有什么好处吗?几乎没有-内存分配是很慢的,所以附加的虚函数损失对他来说是可以接受的。下面的方法可以获得同样灵活的设计:


00001: list<double> A(allocator<double>());
00002: list<double> B(my_allocator<double>());

其中allocator和my_allocator是多态的。

1.5.2 将函数定义在模版外

你可以通过将函数定义移动至类定义的外面来避免创建不需要的重复成员函数。例如,下面是一个Array3D类的框架。


00001: template<class T_numtype>
00002: class Array3D{
00003: public:
00004:   //确保下面的array有同样的大小
00005:   template<class T_numtype2>
00006:   bool shapeCheck(Array3D<T_numtype2>&);
00007: 
00008: private:
00009:   int base[3];              //起始x,y,z的index值
00010:   int extent[3];            //每个x,y,z的长度
00011: };
00012: 
00013: template<class T_numtype> template<class T_numtype2>
00014: bool Array3D<T_numtype>::shapeCheck(Array3D<T_numtype2>&B)
00015: {
00016:   for(int i=0;i<3;++i)
00017:     if((base[i]!=B.base(i))||(extent[i]!=B.extent(i)))
00018:       return false;
00019: 
00020:    return true;
00021: }
00022: 

shapeCheck()函数确定两个array是否有同样的起始index和长度。如果你有很多种不同的array(Array3D<float>,Array3D<int>,...),你会最终拥有很多个shapeCheck()函数的实例。但是,其实这些实例是完全一样的,因为shapeCheck()并不使用参数T_numtype和T_numtype2。

为了避免产生多重的实例,你可以将shapeCheck()移动至没有模版参数的基类,或者移至全局函数,如下:


00001: class Array3D_base{
00002: 
00003: public:
00004:   bool shapeCheck(Array3D_base&);
00005: 
00006: private:
00007: 
00008:   int base[3];
00009:   int extent[3];
00010: };
00011: 
00012: template<class T_numtype>
00013: class Array3D:public Array3D_base{
00014:   ...
00015: };

由于shapeCheck()现在不在模版类中,所以只有一个它的实例。如果你正好是C++编译器的作者,读到这里,你应该注意,这个问题可以通过自动的进行模版参数有效性检测来避免。

1.5.3 Inlining 级别

一般情况,我趋向于在Blitz++中将很多东西inline化,因为这样子会协助编译器优化(对于那些不做interprocedural分析和强化inlining的编译器)。如果有一天,我感觉Blitz++中的inline太多了,这有一种策略可以处理这个问题。


00001: #if BZ_INLINING_LEVEL==0
00002:   //NO INLINING
00003:   #define _bz_inline1
00004:   #define _bz_inline2
00005: #elif BZ_INLINING_LEVEL==1
00006:   #define _bz_inline1 inline
00007:   #define _bz_inline2
00008: #elif BZ_INLINING_LEVEL==2
00009:   #define _bz_inline1 inline
00010:   #define _bz_iiline2 inline
00011: #endif
00012: 
00013: _bz_inline void foo(Array<int>&A)
00014: {
00015:   //...
00016: }

用户如果不想要任何的inline,可以-DBZ_INLINING_LEVEL=0来编译,还可以控制其为1,2来控制INLINE的程度。

1.6 容器

容器类可以大致分为两个部分:STL风格(基于iterator的)和Data/View风格的。

1.6.1 风格容器

1.6.2 Data/View 容器

1.7 Aliasig和限制

Aliasing对C++科学计算库来说是影响性能的一大问题。例如,考虑计算1阶矩阵A <- A+xxT:


00001: void rank1Update(Matrix&A,const Vector&x)
00002: {
00003:   for(int i=0;i<A.rows();++i)
00004:      for(int j=0;j<A.cols();++j)
00005:         A(i,j)+=x(i)*x(j);
00006: }

对于好的性能,编译器需要保持x(i)在寄存器中而不是在内部循环里重复的载入它。尽管如此,编译器不能确定的是在A中的数据不会覆盖x(i)中的数据。这对于矩阵和向量类的设计来说是不可能采用的,但是编译器不知道这些。Aliasing的可能性导致编译器产生低效率的代码。

另一个Aliasing的结果是中断循环的软件流水线。 例如,一个DAXPY操作y <- y + a*x


00001:   void daxpy(int n,double a,double*x,double *y)
00002: {
00003:   for(int i=0;i<n;++i)
00004:     y[i]=y[i]+a*x[i];
00005: }
00006: 

为了取得好的效率,一个优化编译器会部分的展开循环,然后进行软件流水线来找到有效的内部循环方法。设计这个方法经常需要改变循环中的载入和存储顺序,例如:


00001:   //虚拟的流水化DAXPY循环
00002:   double* xt=x;
00003:   double* yt=y;
00004:   for(;i<n;i+=4)
00005: 
00006:      double t0=yt[0]+a*xt[0];
00007:      double t1=yt[1]+a*xt[1];
00008:      yt[0]=t0;
00009:      double t2=yt[2]+a*xt[2];
00010:      double t3=yt[3]+a*xt[3];
00011:      yt[1]=t1
00012:      double t3=yt[3]+a*xt[3];
00013:      i+=4;
00014:      xt+=4;
00015:      yt[2]=t2;
00016:      yt+=4;
00017:      yt[3]=t3;
00018:   }

不幸的是,编译器无法改变载入和存储的顺序因为x,y指针可能在内存中重迭。在Fortran77/90中,编译器可以假设数组地址没有重叠。对于在cache中的数组这会导致20-50%的性能差异。这依赖于特定的循环/架构。Aliasing对于out-of-cache数组来说不是大问题。

优秀的c/c++编译器会做Aliasing分析,目的是确定是否可以认为指针可被认为是指向不重叠内存的。Alias分析可以减少许多的Aliasing可能性,但是他不是完美的。例如,分析像daxpy函数的程序必须做全局分析。编译器必须在同一时间考虑你执行的所有函数。除非你的编译器做了全局alias分析(SGI会),否则你就会有aliasing问题。

NCEG(Numerical C Extensions Group)设计了一个关键字restrict,可以被一些编译器识别(KAI,Cray,SGi,Intel,gcc)。restrict是一个可以告诉编译器该变量没有Aliasing的修饰符。例如:


00001:   double * restrict a;  //a[0],a[1]..没有aliasing
00002:   int & restrict b ;    //b没有aliases

在C++中restrict的明确定义并不是很nailed down。但是当restrict使用后,编译器可以自由的优化而不去考虑装载/存储顺序。

不幸的是将restrict加入iso/ansi c++标准的提议被否决了。我们最终得到的是一个能力有争议的valarray。valarray数组类设计为没有aliasing,因此可以允许编译器对其上的操作进行优化。

尽管如此,restrict已经被ansi/iso c列为标准,因此有可能在下一次中加入c++标准。当可以使用restrict时,开发程序可以使用预先编译kludge。例如,blitz++定义了如下的宏:


00001:  #ifdef BZ_HAS_RESTRICT
00002:   #define _bz_restrict restrict
00003:  #else
00004:   #define _bz_restrict
00005:  #endif
00006: 
00007: 
00008: class Vector{
00009: double * _bz_restrict data_;
00010: 
00011: };

在一些平台上,编译器有选项可以允许编译器假设没有aliasing。这可以导致好的效果,但是必须给予注意。这些选项会改变你所有程序中的指针语意义semantics of pointers。如果你对指针做了写特殊处理,这可能会导致问题。

1.8 Traits萃取器

Traits技术是由Nathan Myers在C++ Report 1995年June发表的《Traits:a new and useful template technique》中最新提出的。 Trait类保存一些映射关系,将一些概念映射到另一些上去。你可以从以下东西发出映射:

1.8.1 例子: average()

考虑下面的代码,它计算一个数组的平均值:


00001: template<class T>
00002: T average(const*T data,int numElements)
00003: {
00004:   T sum=0;
00005:   for (int i=0;i<numElements;++i)
00006:       sum+=data[i];
00007:   return sum/numElements;
00008: }
00009: 

如果T是一个浮点数,它会工作的很好,但是如果是一下的情况,就不好说了:

Valid XHTML 1.0 :: Valid CSS :: Made with Emacs-Muse