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.