C++矩阵优化算法

写在前面的话：此矩阵优化算法并非原创，笔者只是转述了一种智慧。 


在C++中，定义一个矩阵通常是这样的:   
class   MyMatrix   
{   
        ........   
        public:   
        ........   
        float   data[50000];   
} 
              
这里为了方便说明问题，使用固定大小的数组，实际使用中更多的是动态分配。   在对MyMatrix重载operator   +和-后，我们就可以进行如下计算了:   
MyMatrix   m1,m2,m3,m4;   
.....   
m4=m1+m2-m3;//表达式   
......   


C++编译器将表达式解释为：先把m1+m2计算好赋给一个临时MyMatrix类变量(tmp1)，tmp1-m3后生成新的临时变量tmp2,然后才把tmp2赋值给m4。程序执行过程中，会产生临时变量tmp1,tmp2（有些C++编译器可以优化掉其中一个或全部）,由于data一般较大，分配内存会占用时间和空间。这就是为什么在数值计算方面强大的C++比Fortran慢的重要原因之一。但如果放弃MyMatrix的operator   +和-重载,添加operator   []取data[]的重载后，写如下代码:   


struct   plus;   struct   minus;   
template   <class   L,   class   OpTag,   class   R>   
struct   Expression   
{   
        Expression(L   const&   l,   R   const&   r)   
                    :   l(l),   r(r)   {} 


        float   operator[](unsigned   index)   const;   


        L   const&   l;   
        R   const&   r;   
};   


template   <class   L,   class   R>   
Expression <L,plus,R>   operator+(L   const&   l,   R   const&   r)   
{   
        return   Expression <L,plus,R> (l,   r);   
}   
template   <class   L,   class   R>   
Expression <L,minus,R>   operator-(L   const&   l,   R   const&   r)   
{   
        return   Expression <L,minus,R> (l,   r);   
}   
struct   plus   
{   
        static   float   apply(float   a,   float   b)   
        {   return   a   +   b;   }   
};   


struct   minus   
{   
        static   float   apply(float   a,   float   b)   
        {   return   a   -   b;   }   
};   


对MyMatrix添加=重载:   
template   <class   Expr>   
MyMatrix   &MyMatrix::operator=(Expr   const&   x)   
{   
        for   (unsigned   i   =   0;   i   <   50000;   i++)   
                (*this)   =   x;   
        return   *this;   
}   


然后计算:   
......   
m4=m1+m2-m3;   




这时生成的临时变量类型是Expression <L,plus,R> ,和Expression <Expression <L,plus,R>   ,minus,R> ,   它们占用的内存远比MyMatrix小,Expression直到operator=MyMatrix   的时候才展开计算，实际上是把矩阵运算变成了加法，所以大大加快了速度。这正是数值计算所需要的。   
这种方法叫Expression   templates优化。更多细节参考 < <C++   template> > 、Blitz++库、boost::ublas库和MTL库等。 


后话：编程语言之争一直是热门话题。因为笔者所接触的有限元计算程序涉及大规模的数值计算，所以笔者也曾花了2年时间搞腾Fortran。最后，笔者还是放弃古老的Fortran，转投C++的阵营了——毕竟连微软都早放弃Fortran了:) 


Fortran在数值计算方面的确有些优势，但Fortran在软件工程方面的表现得实在让人汗颜。须知，面向对象并不是一个炒作出来的概念，面向对象对改进软件工程方面起了里程碑式的作用，大大提高了软件开发效率。 


从本质上看，可执行程序都编译成二进制代码了，怎么会存在Fortran比C++快的道理？其实最主要原因还是开发人员对C++的数值计算不熟悉所致。从上面的例子也可以看出，Fortran只是在某些方面做了一些适合数值计算的技术处理，C++并非不能做到，只是一般人平时没留心罢了。 


Fortran计算速度快还有一个原因是大量使用了全局变量。Fortran有一个公共变量块的定义，可以在那里设置全局变量。如果你愿意，你也可以在C++程序中大量使用全局变量提高程序运行速度。问题是这样做的后果是在软件工程方面带来无尽的后遗症，软件维护变得异常困难。 


C++比Fortran慢有一个不好克服的原因就是使用了面向对象。虚函数的重载要查找虚函数表，必然会降低效率。所以在听说了Fortran要支持面向对象以后，彻底让笔者转投到C++的阵营去了——如果Fortran支持面向对象，他在数值计算方面的优势将荡然无存。 


Fortran，过时的工具了，还是放弃吧。 


我的博客:   http://chenkong0500.blog.163.com/ 
欢迎交流