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 閱讀(900) 評論(0)  編輯 收藏 引用 所屬分類: 數值計算