Blitz++与MTL两大数值计算程序库(C++)的简介 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, 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

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, dense<>, row_major>::type Matrix;

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

MatrixA(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
{

for (Matrix::size_type j=0; 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, dense, 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, dense<>, row_major>::type LUMatrix;

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

mtl::copy(A, LU);

typedef dense1D 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

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, dense, 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, dense<>, row_major>::type CMatrix;

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

mtl::copy(A, LU);

typedef dense1D 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.
【Blitz++与MTL两大数值计算程序库(C++)的简介】