牵着老婆满街逛

严以律己,宽以待人. 三思而后行.
GMail/GTalk: yanglinbo#google.com;
MSN/Email: tx7do#yahoo.com.cn;
QQ: 3 0 3 3 9 6 9 2 0 .

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

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

Blitz++ MTL 都是基于 C++ template 高效数值计算程序库,不过他们专注于不同的方向。

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 ;

}

 

2 ,下面是一个线性方程组的解法

#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 ;

}

 

 

参考: 1 ,数值方法 (Matlab ) 3rd

           John H.Mathews, Kurtis D.Fink 著, 陈渝,周璐,钱方 等译

       2 Matlab 6.5 的文档     The MathWorks, Inc.


作者Blog: http://blog.csdn.net/newsuppy/

posted on 2006-07-01 13:26 杨粼波 阅读(3063) 评论(1)  编辑 收藏 引用

评论

# re: Blitz++与MTL两大数值计算程序库(C++)的简介 2008-08-20 16:20 littlewater

不知道MTL71版本能不能更新一下,下载不能t.T  回复  更多评论   


只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理