健康，快乐，勇敢的宁帅！！

Blitz++ 与 MTL 两大数值计算程序库 (C++) 的简介

Blitz++ 提供了一个 N 维（ 1—10 ）的 Array 类 , 这个 Array 类以 reference counting 技术实现，支持任意的存储序 (row-major 的 C-style 数组， column-major 的 Fortran-style 数组 ) ，数组的切割 (slicing), 子数组的提取 (subarray), 灵活的 Array 相关表达式处理。另外提供了可以产生不同分布的随机数 (F,Beta,Chi-Square ，正态，均匀分布等 ) 的类也是很有特色的。

MTL 专注于线性代数相关的计算任务，如各种形式矩阵的生成 ( 对角，共轭，稀疏，对称等 ) ，相关的计算，变换，以及与一维向量的运算。

两个程序库对于从 Matlab 导入导出数据都有不错的支持。

本文主要介绍如何在 Visual C++7.1 编译器下运用这两个程序库。

以前的 VC6 编译器由于对 ISO C++98 标准的支持不够，特别是在 template 方面，以至于很难编译这种完全用 template 技术构造起来的程序库。 Blitz++ 是完全不支持 VC6 的。

到了 VC7.1 ，由于对于 ISO 标准的支持达到了 98% ，使得我们可以很轻松的编译使用这两个程序库。

不过这两个程序库的文档不是那么友好，特别是 MTL ，仅仅提供了类似于 reference 的文档，对于具体的使用方法则不作介绍。 Blitz++ 相对来说好一些，还提供一份介绍性的入门文档 。所以使用这两个程序库阅读其源代码往往是必要的。当然了，两个程序库都是 template 代码，源代码必定是全开放的。

先来介绍一下配置吧 。

1，  Blitz++, 目前最高版本是 0.7 ， Blitz++ 已经成为 SourceForge 的项目了，所以可以在 SourceForge.net 下载到。下载后解压缩，你会看到 \Blitz++-0.7\blitz 和 \Blitz++-0.7\random 两个文件夹，这是 blitz 的源代码所在处。 \Blitz++-0.7\manual 是文档所在文件夹。 \Blitz++-0.7\benchmarks,\Blitz++-0.7\examples 和 \Blitz++-0.7\testsuite 中都有很多好的使用实例可供参考。

现在将 VC++ 的 IDE 的 Include 设置为 \Blitz++-0.7 ，因为 blitz 源码中都有这样形式的 #include ,#include 。或者就干脆把两个源码文件夹整个得 copy 到 include 文件夹内。然后将 blitz 文件夹下的 config.h 改为其它名字，而将 config-VS.NET2003.h 的名字改为 config.h 。 OK, 现在你就可以编译所有的 testsuite 和 benchmarks 了。

1，  MTL 的配置相对来说麻烦一点，现在 http://www.osl.iu.edu/research/mtl/ 这里下载一个 VC++7 的，不过还不能马上用。由于 VC++7.1 对标准的支持更近了一步，同时对于某些语法细节的检查更为严格 ( 主要是对于 typename 和 template partial specialization ），我们要对代码做一些小小地修改，特别是 mtl/mtl_config.h 这个文件。有一些地方要加入 typename 。另外有两个模板偏特化的情况需要修改，加上 template <> 。在这里 http://newsuppy.go.nease.net/mtl.zip 我提供了一个修改完成的版本，不过我不保证我的修改可能引入的新的 bugs ，所以请谨慎使用。 MTL 的内部使用一定数量的 STL 组件和算法。 MTL 的源代码都在 mtl 文件夹内，由于 mtl 内部的 include 都是 #include “…” 的形式，使用时把 mtl 文件夹复制到当前 project 下就可以。如果要设 VC++ 的 Include 目录，则应该先把所有的 #include “…” 改为 #include  <…> 这样的形式。

不过刚开始使用 MTL 还是有一些不太容易让人接受的地方。比如 mtl::matrix 这个模板类并不能够产生实际的矩阵对象，而要通过它的 type 成员产生一个对应模板参数的类型，再通过这个类型来实例化对象。

比如 typedef mtl::matrix< SPAN>, rectangle<>, dense<>, row_major >::type Matrix; Matrix A;

这里的 A 才是真正的矩阵对象，而 Matrix 则是一个元素为 float ，矩形，密集，行主 (C-style) 的矩阵类。

       下面我提供三个简单的入门例子解释 MTL 的使用。分别有矩阵的加法，乘法，转置，求逆以及一个线性方程组求解的例子。

       另外 mtl 的 test 和 contrib 文件夹下也有很多不错的示例代码可以查阅。

MTL 使用示例 1 ，矩阵的加法，乘法和转置。  

#include 

#include 

#include  " mtl/mtl.h "                                                                                         

#include 

using   namespace  std; 

using   namespace  mtl; 

                                                                                                                           

template  <   class  Matrix >  

void  print_matrix(Matrix &  mat, const   string &  description) 

{ 

       std::cout  <<  description; 

  

       std::cout  << ' [ ' ; 

        for  (Matrix::iterator i  =  mat.begin(); i != mat.end(); ++ i) 

        { 

          for  (Matrix::OneD::iterator j  = ( * i).begin(); j != ( * i).end(); ++ j) 

          { 

                std::cout  << ' \t ' <<* j; 

         }  

  

         std::cout  << ((i + 1 ==  mat.end()) ? " \t]\n " :   " \n " ); 

       }  

}  

  

int  main( int  argc, char *  argv[]) 

{ 

  typedef matrix < float , rectangle <> , dense <> , row_major > ::type Matrix; 

  

   const  Matrix::size_type MAX_ROW  = 3 , MAX_COL  = 3 ; 

  

  Matrix  A(MAX_ROW,MAX_COL),B(MAX_ROW,MAX_COL),C(MAX_ROW,MAX_COL); 

  

   //  fill Matrix A with the index syntax 

   for  (Matrix::size_type i = 0 ; i < MAX_ROW; ++ i) 

   { 

          for  (Matrix::size_type j = 0 ; j < MAX_COL; ++ j) 

          { 

                A(i, j) =  Matrix::value_type(rand() % 50 ); 

         }  

  }  

  

   //  fill Matrix B with the iterator syntax 

   for  (Matrix::iterator i = B.begin(); i != B.end(); ++ i) 

   { 

          for  (Matrix::OneD::iterator j = ( * i).begin(); j != ( * i).end(); ++ j) 

          { 

                 * j  =  Matrix::value_type(rand() % 50 ); 

         }  

  }  

  

  print_matrix(A, " A=\n " ); 

  print_matrix(B, " B=\n " ); 

  

   //  Matrix C = A + B 

  add(A, C); 

  add(B,C); 

  print_matrix(C, " C = A + B  \n " ); 

  

   //  Matrix C = A * B^T,  B^T: transpose of B 

  transpose(B); 

  print_matrix(B, " B^T=\n " ); 

  zero_matrix(C);         //  because mult(A, B, C): C += A*B
  mult(A,B,C); 


  print_matrix(C, " C = A * B^T\n " ); 

   return   0  ; 

}  

#include 

#include 

#include 

#include  " mtl/mtl.h "

#include  " mtl/lu.h "

#include 

using   namespace  std;

using   namespace  mtl;

 

int  main( int  argc,  char *  argv[])

{

  typedef matrix < float , rectangle <> , dense < external > , row_major > ::type Matrix;

   //  dense : data copy from a float array,not generate them with yourself

 

   const  Matrix::size_type MAX_ROW  =   3 , MAX_COL  =   3 ;

 

   //  solve the equation Ax=b

   //  { 4x - y + z = 7

   //     4x - 8y + z= -21

   //    -2x + y + 5z = 15 }

   //  A = [ 4 -1 1

   //           4 -8 1

   //           -2 1 5 ]

   //  b = [7 - 21 15]^T

   float  a[]  =   { 4.0f ,  - 1.0f ,  1.0f ,  4.0f ,  - 8.0f ,  1.0f ,  - 2.0f ,  1.0f ,  5.0f } ;

  Matrix A(a, MAX_ROW, MAX_COL);

    

  typedef matrix < float , rectangle <> , dense <> , row_major > ::type LUMatrix;

  LUMatrix LU(A.nrows(), A.ncols());

  mtl::copy(A, LU);

 

  typedef dense1D < float >  Vector;

  Vector pvector(A.nrows());

  lu_factor(LU, pvector);

 

  Vector b(A.nrows()), x(A.nrows());

  b[ 0 ]  =   7.0f , b[ 1 ]  =   - 21.0f , b[ 2 ]  =   15.0f ;

  lu_solve(LU, pvector, b, x);

 

   for  (Vector::iterator i = x.begin(); i != x.end();  ++ i)

         cout  <<   * i  <<   ' \t ' ;

 

  system( " pause " );

   return   0 ;

}

3 ，矩阵求逆

#include 

#include 

#include 

#include  " mtl/mtl.h "

#include  " mtl/lu.h "

#include 

using   namespace  std;

using   namespace  mtl;

 

template  < class  Matrix >

void  print_matrix(Matrix &  mat,  const   string &  description)

{                                                                                             

       std::cout  <<  description;

 

       std::cout  <<   ' [ ' ;

        for  (Matrix::iterator i  =  mat.begin(); i != mat.end();  ++ i)

        {

          for  (Matrix::OneD::iterator j  =  ( * i).begin(); j != ( * i).end();  ++ j)

          {

                std::cout  <<   ' \t '   <<   * j;

         }

 

         std::cout  <<  ((i + 1   ==  mat.end())  ?   " \t]\n "  :   " \n " );

       }

}

 

int  main( int  argc,  char *  argv[])

{

  typedef matrix < float , rectangle <> , dense < external > , row_major > ::type Matrix;

   //  dense : data copy from a float array,not generate them with yourself

 

   const  Matrix::size_type MAX_ROW  =   3 , MAX_COL  =   3 ;

 

   //  inverse matrix A

   //  A = [ 4 -1 1

   //           4 -8 1

   //           -2 1 5 ]

   float  a[]  =   { 4.0f ,  - 1.0f ,  1.0f ,  4.0f ,  - 8.0f ,  1.0f ,  - 2.0f ,  1.0f ,  5.0f } ;

  Matrix A(a, MAX_ROW, MAX_COL);

    

  typedef matrix < float , rectangle <> , dense <> , row_major > ::type CMatrix;

  CMatrix LU(A.nrows(), A.ncols());

  mtl::copy(A, LU);

 

  typedef dense1D < float >  Vector;

  Vector pvector(A.nrows());

  lu_factor(LU, pvector);

 

  CMatrix InvA(A.nrows(), A.ncols()); 

  lu_inverse(LU, pvector, InvA);

  

  print_matrix(A,  " A = \n " );

  print_matrix(InvA,  " A^(-1) = \n " );

  system( " pause " );

   return   0 ;

}

posted on 2006-11-22 23:33 ningfangli 阅读(811) 评论(0)  编辑 收藏 引用 所属分类: 数值计算