Main Page | Namespace List | Class Hierarchy | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages
Base

TMatrix.h

Go to the documentation of this file.
00001 #ifndef PTMATRIX_TMATRIX_H
00002 #define PTMATRIX_TMATRIX_H
00003 /** 
00004     @page PTM ポータブル テンプレート 行列クラスライブラリ
00005     
00006     @author 長谷川 晶一
00007 
00008     @date 2001年6月10日,2003年10月20日更新
00009 
00010     @section introM はじめに
00011         このドキュメントはポータブル テンプレート 行列クラスライブラリ
00012         のドキュメントです．
00013         ポータブル テンプレート 行列クラスライブラリは，
00014         テンプレートによるM×N行列のクラスライブラリです．
00015         行列のサイズをテンプレートで持つバージョンと変数で持つバージョンがあります．
00016         
00017     @section specM このライブラリの特徴
00018     @subsection tmpM テンプレートで解決
00019         テンプレート版は，全てをテンプレートで静的に解決しています．
00020         行列のサイズやポインタなどを保持するための管理領域を持ちません．
00021         そのため，
00022         <ul>
00023             <li> 行列の一部(部分行列，ベクトル，行ベクトル，列ベクトル)などを
00024                 直接参照することができる(代入も可能)．
00025             <li> 配列をキャストして行列として使用することも可能．
00026             <li> 行列の掛け算などで，行列のサイズが合わない場合，コンパイル時に
00027                 コンパイラがエラーを出力する．
00028             <li> 行列のサイズを動的に変更することができない．
00029         </ul>
00030         といった特徴をもちます．
00031     @subsection portM 移植性
00032         処理系依存部分の多いテンプレート機能を活用していながら，全機能が，
00033         3つのコンパイラで使用できます．サポートしているコンパイラは，
00034         <ul>
00035             <li> CL (MS Visual C++ 6.0)
00036             <li> bcc32(Borland C++ 5.5.1)   ごめんなさいまだ未確認です．多分動作しません．
00037             <li> gcc(GNU c compiler 2.95.3-5)
00038         </ul>
00039         です．
00040     @subsection why 新たにクラスライブラリを作成した理由
00041         すでに多くの行列ライブラリがありながら，新たに作成した理由は，
00042         - TNT，MTL の行列は管理領域をメモリに持つため，配列をキャストして
00043             行列として使用することができない．
00044         - Blitz は TinyMatrix, TinyVector を持つが，Visual C++ で
00045             使用できない．
00046         - 私の知る限り，部分行列，部分ベクトルへの参照を返す行列ライブラリ
00047             は存在しない．
00048 
00049     からです．
00050     
00051     @section ptm_usage< 使い方
00052     ポータブル テンプレート 行列クラスライブラリは，ヘッダファイルだけからなる
00053     クラスライブラリなので, TMatrix.h, TMatrixUtility.h, TVector.h
00054     を同じフォルダに入れておき，.cppファイルからヘッダをインクルードするだけで
00055     使用できます．
00056     @subsection sampleM サンプル
00057         簡単なサンプルです．適当なファイル名(たとえば sample.cpp) で保存して
00058         コンパイルしてください．コンパイルするためには，
00059         <DL>
00060         <DT> visual C++ の場合
00061         <DD> cl -GX sample.cpp
00062         <DT> gccの場合
00063         <DD> g++ sample.cpp
00064         <DT> bcc の場合
00065         <DD> bcc32 sample.cpp
00066         </DL>
00067     としてください．
00068     @verbatim
00069 #include "TMatrix.h"    //  行列ライブラリのインクルードする．
00070 #include <iostream>
00071 using namespace PTM;    //  行列クラスはPTM名前空間の中で宣言されている．
00072 void main(){
00073     TMatrixRow<2,2,float> mat;     //  2行2列の行列を宣言
00074     mat[0][0] = 1;  mat[0][1] = 2;
00075     mat[1][0] = 3;  mat[1][1] = 4;
00076     TVector<2,float> vec;       //  2次元のベクトルを宣言
00077     vec[0] = 1; vec[1] = 0;
00078     std::cout << mat;
00079     std::cout << vec << std::endl;
00080     std::cout << mat * vec << std::endl;    //  掛け算
00081     std::cout << mat + mat << std::endl;    //  足し算
00082     std::cout << mat - mat << std::endl;    //  引き算
00083     std::cout << mat.trans() << std::endl;  //  転置
00084     std::cout << mat.inv() << std::endl;    //  逆行列
00085 }
00086 @endverbatim
00087     @subsection vecfunc ベクトルの機能
00088     次の演算ができます．
00089     <ul>
00090     <li> +:和, -:差, *:内積/定数倍, /:定数分の1
00091     <li> ==:比較, =:代入
00092     <li> <<:出力, >>:入力
00093     <li> %:外積(2・3次元のみ)
00094     </ul>
00095     次のメンバ関数を持ちます．
00096     <ul>
00097     <li> unit(): 向きが等しい単位ベクトルを返す．
00098     <li> norm(): ベクトルの大きさ(ノルム)を返す．
00099     <li> sub_vector(): 部分ベクトルへの参照を返す．
00100     </ul>
00101     変数版はサイズの変更ができます．
00102     <ul>
00103     <li> resize(int h, int w):  サイズの変更
00104     <li> height():  高さ＝行列の行数
00105     <li> width():   幅＝行列の列数
00106     </ul>
00107     @subsection matfunc 行列の機能
00108     次の演算ができます．
00109     <ul>
00110         <li> +:和, -:差, *:積/定数倍, /:定数分の1
00111         <li> ==:比較, =:代入
00112         <li> <<:出力, >>:入力
00113     </ul>
00114     次のメンバ関数を持ちます．
00115     <ul>
00116         <li> det(): 行列式を返す．
00117         <li> inv(): 逆行列を返す．
00118         <li> gauss(): ガウスの消去法で方程式を解く．
00119         <li> sub_matrix(): 部分行列への参照を返す．
00120         <li> row(): 行ベクトルへの参照を返す．
00121         <li> col(): 列ベクトルへの参照を返す．
00122     </ul>
00123         sub_matrix()やrow()には値を代入することもできます．
00124         @verbatim
00125     TMatrixRow<3,3,float> mat; TVector<3, float> vec;
00126     mat.row() = vec;@endverbatim
00127     @section pub 再配布の条件
00128     再配布する場合は，原著作者の署名・連絡先を改変・削除しないでください．
00129     閉鎖的な会社でソフトウエアを書く人にも自由に使用できるようにしたい
00130     ので，GPL・LGPLにしませんでした．
00131     もちろんGPL・LGPLに改変して再配布してくださっても結構です．
00132     @section support サポート
00133     バグ，不具合，意味不明なコンパイルエラーなどを見つけた場合は，
00134     長谷川 晶一 (hase@hi.pi.titech.ac.jp) までご連絡ください．
00135     できる限りサポートし，より良いライブラリにしていくつもりです．<br>
00136     特にこのライブラリはテンプレートクラスライブラリなので，使用時に，
00137     内部の動作が分らないと意味のわからないコンパイルエラーに出会うことも
00138     あると思います．そのような問題には対応するつもりですので，まずはご連絡ください．
00139     @section thanksM 謝辞
00140     LU分解，逆行列，ガウス消去法などの行列計算アルゴリズムは，<br>
00141     「『Ｃ言語による最新アルゴリズム事典』全ソースコード」<br>
00142     ftp://ftp.matsusaka-u.ac.jp/pub/algorithms<br>
00143     奥村 晴彦 Haruhiko Okumura<br>
00144     を改変して流用させていただきました．
00145     自由にコードを使わせてくださって，ありがとうございます．    */
00146 //-----------------------------------------------------------------------------
00147 /** @file TMatrix.h
00148     テンプレートによるN×M行列型の定義.
00149     要素の型とサイズをテンプレートの引数にすることで，
00150     管理情報をメモリに持たずに，一般の行列を扱う．
00151     配列をキャストして行列にすることもできる．
00152     部分行列や行や列を行列やベクトルとして取り出すこともできる．
00153     sparse matrix には対応していない．                                      */
00154 //------------------------------------------------------------------------------
00155 
00156 #include "TVector.h"
00157 
00158 /// Portable Template Matrixライブラリの名前空間
00159 namespace PTM{;
00160 
00161 #ifdef _WIN32
00162  #pragma pack(push, 4)
00163  #ifdef _DEBUG
00164   #pragma optimize ("awgity", on)
00165   #pragma auto_inline(on)
00166   #pragma inline_recursion(on)
00167  #endif
00168 #endif
00169 
00170 /** 部分行列型作成のためのユーティリティークラス.
00171     TSubMatrixDim<top, left, height, width> と次元を指定できる。*/
00172 template <size_t T, size_t L, size_t H, size_t W>
00173 class TSubMatrixDim{
00174 public:
00175     DIMDEF(DIMENC(T),       TOP);
00176     DIMDEF(DIMENC(L),       LEFT);
00177     DIMDEF(DIMENC(H),       HEIGHT);
00178     DIMDEF(DIMENC(W),       WIDTH);
00179 };
00180 /** 部分行列型作成のためのユーティリティークラス.
00181     TMatDim<height, width> と次元を指定できる。 */
00182 template <size_t H, size_t W>
00183 class TMatDim{
00184 public:
00185     DIMDEF(DIMENC(H),       HEIGHT);
00186     DIMDEF(DIMENC(W),       WIDTH);
00187 };
00188 
00189 template <size_t H, size_t W, class OD> class TSubMatrixRow;
00190 template <class T, class Z=T, class U=Z> class ESubMatrixRow;
00191 template <size_t H, size_t W, class OD> class TSubMatrixCol;
00192 template <class T, class Z=T, class U=Z> class ESubMatrixCol;
00193 /// 行ベクトルの基本型  サイズ：テンプレート
00194 template <class desc>
00195 class TMakeSubMatrixRow{
00196 public:
00197     typedef TYPENAME desc::element_type element_type;
00198     typedef TYPENAME desc::exp_type exp_type;
00199     exp_type& exp(){ return *(exp_type*)this;}
00200     const exp_type& exp() const { return *(exp_type*)this;}
00201     /// @name 部分行列
00202     //@{
00203     /// 部分行列（テンプレート版）
00204     template <class SUB>
00205         TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(SUB){
00206         return (TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(DIMDEC(SUB::TOP), DIMDEC(SUB::LEFT));
00207     }
00208     template <class SUB>
00209         const TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(SUB) const {
00210         return (TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(DIMDEC(SUB::TOP), DIMDEC(SUB::LEFT));
00211     }
00212     /// 部分行列（サイズだけテンプレート版）
00213     template <class SUB>
00214         TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(size_t t, size_t l, SUB){
00215         return (TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(t,l);
00216     }
00217     template <class SUB>
00218         const TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(size_t t, size_t l, SUB) const {
00219         return (TSubMatrixRow<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(t,l);
00220     }
00221     /// 部分行列（変数版）
00222     ESubMatrixRow<element_type> vsub_matrix(int t, int l, int h, int w){
00223         return ESubMatrixRow<element_type>(h, w, exp().stride(), &exp().item(t,l));
00224     }
00225     //@}
00226 };
00227 /// 列ベクトルの基本型  サイズ：テンプレート
00228 template <class desc>
00229 class TMakeSubMatrixCol{
00230 public:
00231     typedef TYPENAME desc::element_type element_type;
00232     typedef TYPENAME desc::exp_type exp_type;
00233     exp_type& exp(){ return *(exp_type*)this;}
00234     const exp_type& exp() const { return *(exp_type*)this;}
00235     /// @name 部分行列
00236     //@{
00237     /// 部分行列（テンプレート版）
00238     template <class SUB>
00239         TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(SUB){
00240         return (TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(DIMDEC(SUB::TOP), DIMDEC(SUB::LEFT));
00241     }
00242     template <class SUB>
00243         const TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(SUB) const {
00244         return (TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(DIMDEC(SUB::TOP), DIMDEC(SUB::LEFT));
00245     }
00246     /// 部分行列（サイズだけテンプレート版）
00247     template <class SUB>
00248         TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(size_t t, size_t l, SUB){
00249         return (TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(t,l);
00250     }
00251     template <class SUB>
00252         const TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>& sub_matrix(size_t t, size_t l, SUB) const {
00253         return (TSubMatrixCol<DIMDEC(SUB::HEIGHT), DIMDEC(SUB::WIDTH), desc>&)exp().item(t,l);
00254     }
00255     /// 部分行列（変数版）
00256     ESubMatrixCol<element_type> vsub_matrix(int t, int l, int h, int w){
00257         return ESubMatrixCol<element_type>(h, w, exp().stride(), &exp().item(t,l));
00258     }
00259     //@}
00260 };
00261 /// 行ベクトルの基本型  サイズ：変数
00262 template <class desc>
00263 class EMakeSubMatrixRow{
00264 public:
00265     typedef TYPENAME desc::element_type element_type;
00266     typedef TYPENAME desc::exp_type exp_type;
00267     exp_type& exp(){ return *(exp_type*)this;}
00268     const exp_type& exp() const { return *(exp_type*)this;}
00269     /// @name 部分行列
00270     //@{
00271     /// 部分行列（変数版）
00272     ESubMatrixRow<element_type> vsub_matrix(int t, int l, int h, int w){
00273         return ESubMatrixRow<element_type>(h, w, exp().stride(), &exp().item(t,l));
00274     }
00275     //@}
00276 };
00277 /// 列ベクトルの基本型  サイズ：変数
00278 template <class desc>
00279 class EMakeSubMatrixCol{
00280 public:
00281     typedef TYPENAME desc::element_type element_type;
00282     typedef TYPENAME desc::exp_type exp_type;
00283     exp_type& exp(){ return *(exp_type*)this;}
00284     const exp_type& exp() const { return *(exp_type*)this;}
00285     /// @name 部分行列
00286     //@{
00287     /// 部分行列（変数版）
00288     ESubMatrixCol<element_type> vsub_matrix(int t, int l, int h, int w){
00289         return ESubMatrixCol<element_type>(h, w, exp().stride(), &exp().item(t,l));
00290     }
00291     //@}
00292 };
00293 
00294 
00295 //----------------------------------------------------------------------------
00296 /// @name 行列演算の実装
00297 //@{
00298 template <class DESC> class MatrixImp;
00299 template <DIMTYPE H, DIMTYPE W, class D> class TMatrixBaseBase;
00300 template <DIMTYPE H, DIMTYPE W, class D> class TMatrixBase;
00301 template <class DESC> class EMatrixBase;
00302 
00303 /** 代入(*this = b).
00304     @param b 同じサイズの行列.  */
00305 template <class AD, class BD>
00306 void assign(MatrixImp<AD>& a, const MatrixImp<BD>& b) {
00307     a.resize(b.height(), b.width());
00308     a.size_assert(b);
00309     for(size_t i=0; i<a.height(); ++i)
00310         for(size_t j=0; j<a.width(); ++j)
00311             a.item(i,j) = (TYPENAME AD::element_type)b.item(i,j);
00312 }
00313 /** 代入(*this = b).
00314     @param b 同じサイズのベクトル.  */
00315 template <class AD>
00316 void assign(MatrixImp<AD>& a, const TYPENAME AD::element_type* b) {
00317     for(size_t i=0; i<a.height(); ++i)
00318         for(size_t j=0; j<a.width(); ++j)
00319             a.item(i,j) = b[i*a.width()+j];
00320 }
00321 /// 比較
00322 template <class AD, class BD>
00323 bool equal(const MatrixImp<AD>& a, const MatrixImp<BD>& b){
00324     if (!a.size_check(b)) return false;
00325     for(size_t i=0; i<a.height(); ++i)
00326         for(size_t j=0; j<a.width(); ++j)
00327             if (a.item(i,j) != b.item(i,j)) return false;
00328     return true;
00329 }
00330 /// 加算
00331 template <class AD, class BD>
00332 void add(MatrixImp<AD>& a, const MatrixImp<BD>& b){
00333     a.size_assert(b);
00334     for(size_t i=0; i<a.height(); ++i)
00335         for(size_t j=0; j<a.width(); ++j)
00336             a.item(i,j) += b.item(i,j);
00337 }
00338 /// 減算
00339 template <class AD, class BD>
00340 void sub(MatrixImp<AD>& a, const MatrixImp<BD>& b){
00341     a.size_assert(b);
00342     for(size_t i=0; i<a.height(); ++i)
00343         for(size_t j=0; j<a.width(); ++j)
00344             a.item(i,j) -= b.item(i,j);
00345 }
00346 /// 行列とスカラーの掛け算
00347 template <class AD>
00348 void multi(MatrixImp<AD>& a, TYPENAME AD::element_type b){
00349     for(size_t i=0; i<a.height(); ++i)
00350         for(size_t j=0; j<a.width(); ++j)
00351             a.item(i,j) *= b;
00352 }
00353 /// 行列とベクトルの掛け算
00354 template <class RD, class AD, class BD>
00355 void multi(VectorImp<RD>& r, const MatrixImp<AD>& a, const VectorImp<BD>& b){
00356     r.resize(a.height());
00357     for(size_t n=0; n<r.size(); ++n){
00358         r.item(n) = a.row(n) * b;
00359     }
00360 }
00361 /// 行列とベクトルの掛け算  :   3x3
00362 template <class RD, class AD, class BD>
00363 void multi(TVectorBase<DIMENC(3), RD>& r, const TMatrixBase<DIMENC(3),DIMENC(3),AD>& a, const TVectorBase<DIMENC(3), BD>& b){
00364     r(0) = a.row(0) * b;
00365     r(1) = a.row(1) * b;
00366     r(2) = a.row(2) * b;
00367 }
00368 /// 行列の掛け算
00369 template <class RD, class AD, class BD>
00370 void multi(MatrixImp<RD>& r, const MatrixImp<AD>& a, const MatrixImp<BD>& b){
00371     typedef TYPENAME RD::zero zero;
00372     assert(a.width()==b.height());
00373     r.resize(a.height(), b.width());
00374     for(size_t i=0; i<a.height(); ++i){
00375         for(size_t j=0; j<b.width(); ++j){
00376             r.item(i,j) = zero(0);
00377             for(size_t k=0; k<a.width(); ++k){
00378                 r.item(i,j) += a.item(i,k) * b.item(k,j);
00379             }
00380         }
00381     }
00382 }
00383 /// 行列の掛け算    サイズ：3x3
00384 #ifndef __BORLANDC__
00385 template <class RD, class AD, class BD>
00386 void multi(TMatrixBase<DIMENC(3), DIMENC(3), RD>& r, const TMatrixBase<DIMENC(3), DIMENC(3), AD>& a, const TMatrixBase<DIMENC(3), DIMENC(3), BD>& b){
00387         typedef TYPENAME RD::element_type ET;
00388 #define CALC(i,j) r.item(i,j) = ET( a.item(i,0)*b.item(0,j) + a.item(i,1)*b.item(1,j) + a.item(i,2)*b.item(2,j) )
00389     CALC(0,0);  CALC(0,1);  CALC(0,2);
00390     CALC(1,0);  CALC(1,1);  CALC(1,2);
00391     CALC(2,0);  CALC(2,1);  CALC(2,2);
00392 #undef CALC
00393 }
00394 #endif
00395 /// 行列の掛け算    サイズ：4x4
00396 template <class RD, class AD, class BD>
00397 void multi(TMatrixBase<DIMENC(4), DIMENC(4), RD>& r, const TMatrixBase<DIMENC(4), DIMENC(4), AD>& a, const TMatrixBase<DIMENC(4), DIMENC(4), BD>& b){
00398         typedef TYPENAME RD::element_type ET;
00399 #define CALC(i,j) r.item(i,j) = ET( a.item(i,0)*b.item(0,j) + a.item(i,1)*b.item(1,j) + a.item(i,2)*b.item(2,j) + a.item(i,3)*b.item(3,j) )
00400     CALC(0,0);  CALC(0,1);  CALC(0,2);  CALC(0,3);
00401     CALC(1,0);  CALC(1,1);  CALC(1,2);  CALC(1,3);
00402     CALC(2,0);  CALC(2,1);  CALC(2,2);  CALC(2,3);
00403     CALC(3,0);  CALC(3,1);  CALC(3,2);  CALC(3,3);
00404 #undef CALC
00405 }
00406 
00407 /// 行列式  サイズ：変数
00408 template <class AD>
00409 TYPENAME AD::element_type det(const MatrixImp<AD>& a){
00410     TYPENAME AD::ret_type tmp(a);
00411     VVector<int> ip;
00412     VVector<TYPENAME AD::element_type> w;
00413     ip.resize(a.height());
00414     w.resize(a.height());
00415     return lu(tmp, ip, w);
00416 }
00417 /// 行列式  サイズ：テンプレート
00418 template <DIMTYPE H, DIMTYPE W, class AD>
00419 TYPENAME AD::element_type det(const TMatrixBaseBase<H,W,AD>& a){
00420     TYPENAME AD::ret_type tmp(a);
00421     TVector<DIMDEC(H), int> ip;
00422     TVector<DIMDEC(H), TYPENAME AD::element_type> w;
00423     return lu(tmp, ip, w);
00424 }
00425 #ifndef __BORLANDC__
00426 /// 行列式  サイズ：2x2
00427 template <class AD>
00428 TYPENAME AD::element_type det(const TMatrixBase<DIMENC(2),DIMENC(2),AD>& a){
00429     return a.item(0,0) * a.item(1,1) - a.item(0,1) * a.item(1,0);
00430 }
00431 #endif
00432 /// 行列式  サイズ：3x3
00433 template <class AD>
00434 TYPENAME AD::element_type det(const TMatrixBase<DIMENC(3),DIMENC(3),AD>& a){
00435     return 
00436         ( a.item(0,0) * a.item(1,1) * a.item(2,2)  +  a.item(1,0) * a.item(2,1) * a.item(0,2)  +  a.item(2,0) * a.item(0,1) * a.item(1,2) ) -
00437         ( a.item(2,0) * a.item(1,1) * a.item(0,2)  +  a.item(0,0) * a.item(2,1) * a.item(1,2)  +  a.item(1,0) * a.item(0,1) * a.item(2,2) );
00438 }
00439 /// LU分解を行う。a を書き換える。行列式を返す。
00440 template <class AD>
00441 TYPENAME AD::element_type lu(MatrixImp<AD>& a, int* ip, TYPENAME AD::element_type* weight){
00442 #define ABS_LU_MATRIX(a)    ((a)>0 ? (a) : -(a))
00443     assert(a.width() == a.height());
00444     int i, j, k, ii, ik;
00445     int n = a.height();
00446     TYPENAME AD::element_type t, u, det_;
00447     
00448     det_ = 0;                   // 行列式
00449     for (k = 0; k < n; k++) {  // 各行について
00450         ip[k] = k;             // 行交換情報の初期値
00451         u = 0;                 // その行の絶対値最大の要素を求める
00452         for (j = 0; j < n; j++) {
00453             t = ABS_LU_MATRIX(a.item(k,j));
00454             if (t > u) u = t;
00455         }
00456         if (u == 0) goto PTM_EXIT; // 0 なら行列はLU分解できない
00457         weight[k] = 1 / u;     // 最大絶対値の逆数
00458     }
00459     det_ = 1;                   // 行列式の初期値
00460     for (k = 0; k < n; k++) {  // 各行について
00461         u = -1;
00462         for (i = k; i < n; i++) {  // より下の各行について
00463             ii = ip[i];            // 重み×絶対値 が最大の行を見つける
00464             t = ABS_LU_MATRIX(a.item(ii, k)) * weight[ii];
00465             if (t > u) {  u = t;  j = i;  }
00466         }
00467         ik = ip[j];
00468         if (j != k) {
00469             ip[j] = ip[k];  ip[k] = ik;  // 行番号を交換
00470             det_ = -det_;  // 行を交換すれば行列式の符号が変わる
00471         }
00472         u = a.item(ik, k);  det_ *= u;  // 対角成分
00473         if (u == 0) goto PTM_EXIT;    // 0 なら行列はLU分解できない
00474         for (i = k + 1; i < n; i++) {  // Gauss消去法
00475             ii = ip[i];
00476             t = (a.item(ii, k) /= u);
00477             for (j = k + 1; j < n; j++)
00478                 a.item(ii, j) -= t * a.item(ik, j);
00479         }
00480     }
00481     PTM_EXIT:
00482     return det_;           // 戻り値は行列式
00483 }
00484 //  a x + b = 0 の1次方程式を解く．LU分解済みの必要あり．
00485 template <class AD, class XD, class BD>
00486 void solve(MatrixImp<AD>& a, VectorImp<XD>& x, const VectorImp<BD>& b, int* ip){
00487     int i, j, ii;
00488     TYPENAME XD::element_type t;
00489     const int n = a.height();
00490     for (i = 0; i < n; i++) {       // Gauss消去法の残り
00491         ii = ip[i];  t = b[ii];
00492         for (j = 0; j < i; j++) t -= a.item(ii, j) * x[j];
00493         x[i] = t;
00494     }
00495     for (i = n - 1; i >= 0; i--) {  // 後退代入
00496         t = x[i];  ii = ip[i];
00497         for (j = i + 1; j < n; j++) t -= a.item(ii, j) * x[j];
00498         x[i] = t / a.item(ii, i);
00499     }
00500 }
00501 /// コレスキー法．a,sを書き換える．
00502 template <class AD, class BD>
00503 void cholesky(MatrixImp<AD>& a, VectorImp<BD>& s){
00504     int i,j,k;
00505     int num = a.height();
00506     
00507     //reduction  foreward
00508     a.item(0,0) = sqrt(a.item(0,0));
00509     s.item(0) /= a.item(0, 0);
00510     for(i=1;i<num;i++) {
00511         a.item(0,i) /= a.item(0,0);
00512     }
00513 
00514     for(i=1;i<num;i++){
00515         for(k=0;k<i;k++){
00516             a.item(i,i) -= a.item(k,i)*a.item(k,i);
00517             s.item(i) -= a.item(k,i)*s.item(k);
00518         }
00519         a.item(i,i) = sqrt(a.item(i,i));
00520         for(j=i+1; j<num; j++){
00521             for (k=0; k<i; k++){
00522                 a.item(i,j) -= a.item(k,i) * a.item(k,j);
00523             }
00524             a.item(i,j) /= a.item(i,i);
00525         }
00526         s.item(i) /= a.item(i,i);
00527     }
00528 
00529     // backwark substitution 
00530     for(i=num-1; i>=0; i--){
00531         for(j=i+1; j<num; j++){         
00532             s.item(i) -= a.item(i,j) * s.item(j);
00533         }
00534         s.item(i) /= a.item(i,i);
00535     }
00536 }
00537 /// ガウスの消去法，作業領域(行交換の記録)として， int ip[height()];  が必要．
00538 template <class AD, class XD, class BD>
00539 TYPENAME AD::element_type gauss(MatrixImp<AD>& a, VectorImp<XD>& x, const VectorImp<BD>& b, int* ip){
00540     TYPENAME AD::element_type det_;     // 行列式
00541     AD::col_vector_type::ret_type w;
00542     det_ = lu(a, ip, w);                // LU分解
00543     if (det_ != 0) solve(a, x, b, ip);  // LU分解の結果を使って連立方程式を解く
00544     return det_;                        // 戻り値は行列式
00545 }
00546 
00547 /** 逆行列を求める．
00548     @param a        元の行列(LU分解される)
00549     @param b        逆行列
00550     @param ip       作業領域(行交換の記録)
00551     @param weight   作業領域(行の重み付け)
00552 */
00553 template <class RD, class AD>
00554 TYPENAME AD::element_type inv(MatrixImp<RD>& r, MatrixImp<AD>& a, int* ip, TYPENAME AD::element_type* weight) {
00555     assert(a.height() == a.width());
00556     r.resize(a.height(), a.width());
00557     int i, j, k, ii;
00558     int n = a.height();
00559     TYPENAME AD::element_type t, det;
00560 
00561     det = a.lu(ip, weight);
00562     if (det != 0){
00563         for (k = 0; k < n; k++) {
00564             for (i = 0; i < n; i++) {
00565                 ii = ip[i];  t = (ii == k);
00566                 for (j = 0; j < i; j++)
00567                     t -= a.item(ii, j) * r.item(j, k);
00568                 r.item(i, k) = t;
00569             }
00570             for (i = n - 1; i >= 0; i--) {
00571                 t = r.item(i, k);  ii = ip[i];
00572                 for (j = i + 1; j < n; j++)
00573                     t -= a.item(ii, j) * r.item(j, k);
00574                 r.item(i, k) = t / a.item(ii, i);
00575             }
00576         }
00577     }
00578     return det;
00579 }
00580 
00581 ///  逆行列を返す．
00582 template <class AD>
00583 TYPENAME AD::ret_type inv(const MatrixImp<AD>& a){
00584     typedef TYPENAME AD::ret_type ret_type;
00585     ret_type r, tmp(a);
00586     VVector<int> ip;
00587     ip.resize(a.height());
00588     VVector<TYPENAME AD::element_type> w;
00589     w.resize(a.height());
00590     inv(r, tmp, (int*)ip, (TYPENAME AD::element_type*)w);
00591     return r;
00592 }
00593 ///  逆行列を返す． サイズ：テンプレート
00594 template <class AD, DIMTYPE H, DIMTYPE W>
00595 TYPENAME AD::ret_type inv(const TMatrixBaseBase<H,W,AD>& a){
00596     TYPENAME AD::ret_type r, tmp(a);
00597     TVector<DIMDEC(H), int> ip;
00598     TVector<DIMDEC(H), TYPENAME AD::element_type> w;
00599     inv(r, tmp, (int*)ip, (TYPENAME AD::element_type*)w);
00600     return r;
00601 }
00602 #ifndef __BORLANDC__
00603 /// 逆行列を返す
00604 template <class AD>
00605 TYPENAME AD::ret_type inv(const TMatrixBase<DIMENC(2), DIMENC(2), AD>& a){
00606     TYPENAME AD::element_type d = a.det();
00607     TYPENAME AD::ret_type rv;
00608     rv.item(0,0) =  a.item(1,1) / d;
00609     rv.item(0,1) = -a.item(0,1) / d;
00610     rv.item(1,0) = -a.item(1,0) / d;
00611     rv.item(1,1) =  a.item(0,0) / d;
00612     return rv;
00613 }
00614 #endif
00615 /// 逆行列
00616 template <class AD>
00617 TYPENAME AD::ret_type inv(const TMatrixBase<DIMENC(3), DIMENC(3), AD>& a){
00618 #define DET2_INV_TMATRIXBASE(a,b,c,d)   (a*d - b*c)
00619     TYPENAME AD::ret_type rtv;
00620     TYPENAME AD::element_type det_ = 1 / a.det();
00621     rtv.item(0,0) = DET2_INV_TMATRIXBASE(a.item(1,1), a.item(1,2), a.item(2,1), a.item(2,2)) * det_;
00622     rtv.item(1,0) = DET2_INV_TMATRIXBASE(a.item(1,2), a.item(1,0), a.item(2,2), a.item(2,0)) * det_;
00623     rtv.item(2,0) = DET2_INV_TMATRIXBASE(a.item(1,0), a.item(1,1), a.item(2,0), a.item(2,1)) * det_;
00624         
00625     rtv.item(0,1) = DET2_INV_TMATRIXBASE(a.item(2,1), a.item(2,2), a.item(0,1), a.item(0,2)) * det_;
00626     rtv.item(1,1) = DET2_INV_TMATRIXBASE(a.item(2,2), a.item(2,0), a.item(0,2), a.item(0,0)) * det_;
00627     rtv.item(2,1) = DET2_INV_TMATRIXBASE(a.item(2,0), a.item(2,1), a.item(0,0), a.item(0,1)) * det_;
00628     
00629     rtv.item(0,2) = DET2_INV_TMATRIXBASE(a.item(0,1), a.item(0,2), a.item(1,1), a.item(1,2)) * det_;
00630     rtv.item(1,2) = DET2_INV_TMATRIXBASE(a.item(0,2), a.item(0,0), a.item(1,2), a.item(1,0)) * det_;
00631     rtv.item(2,2) = DET2_INV_TMATRIXBASE(a.item(0,0), a.item(0,1), a.item(1,0), a.item(1,1)) * det_;
00632     return rtv;
00633 #undef DET2_INV_TMATRIXBASE
00634 }
00635 //----------------------------------------------------------------------------
00636 /** 行列のインタフェース ．行列クラスの基本クラスとなり，行列へのアクセスを提供する．
00637     行列の実体は，TMatrix / VMatrix / EMatrix があり，
00638     テンプレート版，変数版，変数版外部バッファとなっている．
00639     行列の添え字とサイズの意味は以下の通り．
00640     @verbatim
00641     (0,0) (0,1)....................(0,m)...  ^
00642     (1,0) (1,1)                        :...  | 
00643       :                                :... height()
00644       :                                :...  |
00645     (n,0) (n,1)....................(n,m)...  V
00646     <---------  width()  -------------->
00647     <---------  stride()  ---------------->@endverbatim
00648     メモリのイメージとしては，複数の行で構成される行列(???Row)と
00649     複数の列で構成される行列(???Col)がある．
00650 */
00651 /// 行列計算の実装
00652 template <class DESC>
00653 class MatrixImp: public DESC::make_sub_matrix{
00654 public:
00655     typedef DESC desc;
00656     typedef TYPENAME desc::exp_type             exp_type;
00657     typedef TYPENAME desc::ret_type             ret_type;
00658     typedef TYPENAME desc::element_type         element_type;
00659     typedef TYPENAME desc::row_vector_ref       row_vector_ref;
00660     typedef TYPENAME desc::const_row_vector_ref const_row_vector_ref;
00661     typedef TYPENAME desc::col_vector_ref       col_vector_ref;
00662     typedef TYPENAME desc::const_col_vector_ref const_col_vector_ref;
00663     typedef TYPENAME desc::trans_ref                trans_ref;
00664     typedef TYPENAME desc::const_trans_ref      const_trans_ref;
00665     typedef TYPENAME desc::zero                 zero;
00666     typedef TYPENAME desc::unit                 unit;
00667     /// コンストラクタ
00668     MatrixImp(){}
00669 
00670     ///@name 行列の参照
00671     //@{
00672     /// 実体の取得
00673     exp_type& exp(){ return *(exp_type*)this; }
00674     const exp_type& exp() const { return *(const exp_type*)this; }
00675     /// n番目の要素を返す(基数は0).
00676     element_type& item(size_t n, size_t m){ return exp().item_impl(n,m); }
00677     const element_type& item(size_t n, size_t m) const { return exp().item_impl(n,m); }
00678     /// 行ベクトル
00679     row_vector_ref row(size_t n){ return exp().row_impl(n) ;}
00680     const_row_vector_ref row(size_t n) const { return exp().row_impl(n) ;}
00681     template <class I>
00682     row_vector_ref operator [] (I n){ return exp().row_impl(n) ;}
00683     template <class I>
00684     const_row_vector_ref operator [] (I n) const { return exp().row_impl(n) ;}
00685     /// 列ベクトル
00686     col_vector_ref col(size_t m){ return exp().col_impl(m) ;}
00687     const_col_vector_ref col(size_t m) const { return exp().col_impl(m) ;}
00688     /// 転置行列
00689     trans_ref trans() { return exp().trans_impl(); }
00690     const_trans_ref trans() const { return exp().trans_impl(); }
00691     /// element_type* へ変換
00692     operator element_type*(){ return &item(0,0); }
00693     operator const element_type*() const { return &item(0,0); }
00694     
00695     /// 行数の取得
00696     size_t height() const { return exp().height_impl(); }
00697     /// 列数の取得
00698     size_t width() const { return exp().width_impl(); }
00699     /// サイズの設定
00700     void resize(size_t h, size_t w){ exp().resize_impl(h,w); }
00701     /// stride の取得
00702     size_t stride() const { return exp().stride_impl(); }
00703     /** ベクトルの要素を全てvにする.
00704         @param v 要素型の値.    */
00705     void clear(const element_type v=zero(0)){
00706         for(size_t i=0; i<height(); ++i)
00707             for(size_t j=0; j<width(); ++j)
00708                 item(i,j) = v;
00709     }
00710     //@}
00711 
00712     /// @name チェック
00713     //@{
00714     template <class B> void size_assert(const MatrixImp<B>& b) const {
00715         assert(height() == b.height() && width() == b.width());
00716     }
00717     template <class B> bool size_check(const MatrixImp<B>& b) const {
00718         return height() == b.height() && width() == b.width();
00719     }
00720     //@}
00721 
00722     /// @name 演算
00723     //@{
00724     /** 代入(*this = b) @param b 同じサイズの行列.  */
00725     template <class BD> void assign(const MatrixImp<BD>& b) { PTM::assign(exp(), b); }
00726     /** 代入(*this = b).@param b 同じサイズのベクトル.  */
00727     void assign(const element_type* b) { PTM::assign(exp(), b); }
00728     /// 比較
00729     template <class BD> bool equal(const MatrixImp<BD>& b) const { return PTM::equal(exp(), b); }
00730     /// 加算
00731     template <class BD> void add(const MatrixImp<BD>& b){ PTM::add(exp(), b); }
00732     /// 減算
00733     template <class BD> void sub(const MatrixImp<BD>& b){ PTM::sub(exp(), b); }
00734     /// スカラーの掛け算
00735     void multi(element_type b){ PTM::multi(exp(), b); }
00736     /// LU分解を行う。thisを書き換える。行列式を返す。
00737     element_type lu(int* ip, element_type* weight){ return PTM::lu(exp(), ip, weight); }
00738     //  (*this) x + b = 0 の1次方程式を解く．LU分解済みの行列でないとだめ．
00739     template <class XD, class BD> void solve(VectorImp<XD>& x, const VectorImp<BD>& b, int* ip){ PTM::solve(exp(), x, b, ip); }
00740     /// コレスキー法
00741     template <class VBASE> void cholesky(VectorImp<VBASE>& s){ PTM::cholesky(exp(), s); }
00742     /// ガウスの消去法，作業領域(行交換の記録)として， int ip[height()];  が必要．
00743     template <class XD, class BD> element_type gauss(VectorImp<XD>& x, const VectorImp<BD>& b, int* ip){ return PTM::gauss(exp(), x, b, ip); }
00744     /// 逆行列を求める。
00745     template <class B> element_type inv(MatrixImp<B>& a_inv, int* ip, element_type* weight) { return PTM::inv(a_inv, exp(), ip, weight); }
00746     /// 逆行列を求める。
00747     ret_type inv() const { return PTM::inv(exp()); }
00748     /// 行列式
00749     element_type det() const { return PTM::det(exp()); }
00750     //@}
00751     ///@name    ストリーム入出力
00752     //@{
00753     /// 出力
00754     void print(std::ostream& os, char* sep="( )") const {
00755 //      os << "sz:" << height() << "," << width() << std::endl;
00756         int w = os.width();
00757         os.width(0);
00758         for(size_t i=0; i<height(); ++i){
00759             for(int j=0; j<w; ++j){
00760                 if (sep[1]) os << sep[1];
00761             }
00762             if (i == 0){
00763                 if(sep[0]) os << sep[0];
00764             }else{
00765                 if(sep[1]) os << sep[1];
00766             }
00767             row(i).print(os);
00768             if (i==height()-1){
00769                 if(sep[2]) os << sep[2];
00770             }
00771             os << std::endl;
00772         }
00773         os.width(w);
00774     }
00775     /// 入力
00776     void input(std::istream& is){
00777         char ch;
00778         is >> ch;
00779         for(int i=0; i<height(); ++i) is >> row(i);
00780         is >> ch;
00781     }
00782     //@}
00783 protected:
00784     /// バッファの初期化．
00785     void init_buffer(){}
00786     /// デフォルトコンストラクタが呼ぶ．
00787     void set_default(){}
00788 };
00789 
00790 /** MatrixImp型派生クラスに必要なメンバの定義.
00791     派生クラスを作るたびに，このマクロを使ってメンバを作る.
00792     @param  THIS    新たに宣言する派生クラスの型名.
00793     @see    MatrixImp   */
00794 
00795 #define DEF_MATRIXD_BASIC_MEMBER(THIS)                                      \
00796     typedef THIS                        this_type;                          \
00797     typedef TYPENAME desc::exp_type             exp_type;                   \
00798     typedef TYPENAME desc::ret_type             ret_type;                   \
00799     typedef TYPENAME desc::row_vector_ref       row_vector_ref;             \
00800     typedef TYPENAME desc::const_row_vector_ref const_row_vector_ref;       \
00801     typedef TYPENAME desc::col_vector_ref       col_vector_ref;             \
00802     typedef TYPENAME desc::const_col_vector_ref const_col_vector_ref;       \
00803     typedef TYPENAME desc::element_type         element_type;               \
00804     typedef TYPENAME desc::zero                 zero;                       \
00805     typedef TYPENAME desc::unit                 unit;                       \
00806     typedef TYPENAME desc::trans_ref            trans_ref;                  \
00807     typedef TYPENAME desc::const_trans_ref      const_trans_ref;            \
00808     /*  行列 b を代入   */                                                  \
00809     template <class B>                                                      \
00810     THIS& operator =(const PTM::MatrixImp<B>& b){                           \
00811         assign(b); return *this;                                            \
00812     }                                                                       \
00813     THIS& operator =(const THIS& b){                                        \
00814         assign(b); return *this;                                            \
00815     }                                                                       \
00816     /** +=演算子(*this = *this + b).                                        \
00817         @param  b   次元が等しい行列    */                                  \
00818     template <class B>                                                      \
00819     this_type& operator +=(const PTM::MatrixImp<B>& b){                     \
00820         add(b); return *this;                                               \
00821     }                                                                       \
00822     /** -=演算子(*this = *this - b). @param b   次元が等しい行列    */      \
00823     template <class B>                                                      \
00824     this_type& operator -=(const PTM::MatrixImp<B>& b){                     \
00825         sub(b); return *this;                                               \
00826     }                                                                       \
00827     /** - 演算子 (return -*this).   */                                      \
00828     ret_type operator- () { ret_type r(*this); r*=-1; return r; }           \
00829     /** *=演算子(*this = b * *this). @param b   要素型  */                  \
00830     this_type operator*= (element_type b){                                  \
00831         multi(b);                                                           \
00832         return *this;                                                       \
00833     }                                                                       \
00834     /** /=演算子(*this = *this / b). @param b   要素型  */                  \
00835     this_type operator/= (element_type b){                                  \
00836         div(b);                                                             \
00837         return *this;                                                       \
00838     }                                                                       \
00839 
00840 #define DEF_MATRIX_BASIC_MEMBER(THIS)                                       \
00841     DEF_MATRIXD_BASIC_MEMBER(THIS)                                          \
00842     /*  デフォルトコンストラクタ    */                                      \
00843     THIS(){ init_buffer(); set_default();}                                  \
00844     /*  行列 b による初期化     */                                          \
00845     template <class B>                                                      \
00846     THIS(const PTM::MatrixImp<B>& b){init_buffer(); assign(b);}             \
00847 
00848 //----------------------------------------------------------------------------
00849 //  次元をテンプレートで持つベクトル    T???Matrix
00850 //
00851 /// 次元をテンプレートで持つベクトルの基本型
00852 template<DIMTYPE H, DIMTYPE W, class D>
00853 class TMatrixBaseBase: public MatrixImp<D> {
00854 protected:
00855     /// バッファの初期化は不用
00856     void init_buffer(){};
00857     /// 継承専用
00858     TMatrixBaseBase(){} 
00859 public:
00860     DIMDEF(H, HEIGHT);          ///<    行数
00861     DIMDEF(W, WIDTH);           ///<    列数
00862     DIMDEF(D::STRIDE, STRIDE);  ///<    列(行)の間隔
00863     typedef D desc;
00864     typedef MatrixImp<desc> base_type;
00865     DEF_MATRIXD_BASIC_MEMBER(TMatrixBaseBase);
00866 
00867     /// 行数
00868     size_t height_impl() const { return DIMDEC(H); }
00869     /// 列数
00870     size_t width_impl() const { return DIMDEC(W); }
00871     size_t stride_impl() const { return DIMDEC(D::STRIDE); }
00872     void resize_impl(size_t h, size_t w) { assert(h==height() && w==width()); }
00873     /// 行ベクトル
00874     row_vector_ref row_impl(size_t n){ return (row_vector_ref)item(n,0); }
00875     const_row_vector_ref row_impl(size_t n) const { return (row_vector_ref)item(n,0); }
00876     /// 列ベクトル
00877     col_vector_ref col_impl(size_t m){ return (col_vector_ref)item(0,m); }
00878     const_col_vector_ref col_impl(size_t m) const { return (col_vector_ref)item(0,m); }
00879     /// 転置
00880     trans_ref trans_impl() { return (trans_ref)item(0,0); }
00881     const_trans_ref trans_impl() const { return (const_trans_ref)item(0,0); }
00882 };
00883 template<DIMTYPE H, DIMTYPE W, class D>
00884 class TMatrixBase: public TMatrixBaseBase<H,W,D> {
00885 public:
00886     DIMDEF(H, HEIGHT);          ///<    行数
00887     DIMDEF(W, WIDTH);           ///<    列数
00888     DIMDEF(D::STRIDE, STRIDE);  ///<    列(行)の間隔
00889     typedef D desc;
00890     typedef TMatrixBaseBase<H,W,D> base_type;
00891     DEF_MATRIXD_BASIC_MEMBER(TMatrixBase);
00892 };
00893 
00894 template <class EXP, class TRANS, size_t H, size_t W, size_t STR, class T, class Z=T, class U=Z>
00895 class TMatrixDescBase{
00896 public:
00897     DIMDEF(DIMENC(STR),             STRIDE);                ///<    行列バッファの幅
00898     typedef EXP                     exp_type;               ///<    実体
00899     typedef exp_type                ret_type;               ///<    返り値型
00900     typedef T                       element_type;           ///<    要素の型
00901     typedef Z                       zero;                   ///<    zero(0)が 0 を返す型
00902     typedef U                       unit;                   ///<    unit(1)が 1 を返す型
00903     typedef TRANS                   trans_type;             ///<    転地行列の型
00904     typedef trans_type&             trans_ref;              ///<    転地行列の参照型
00905     typedef const trans_type&       const_trans_ref;        ///<    転地行列の参照型
00906 };
00907 template <class EXP, class TRANS, size_t H, size_t W, size_t STR, class T, class Z=T, class U=Z>
00908 class TMatrixDescRow: public TMatrixDescBase<EXP,TRANS,H,W,STR,T,Z,U>{
00909 public:
00910     typedef TMakeSubMatrixRow< TMatrixDescRow<EXP,TRANS,H,W,STR,T,Z,U> > make_sub_matrix;
00911     typedef TVector<W,T>            row_vector_type;        ///<    行ベクトル型
00912     typedef row_vector_type&        row_vector_ref;         ///<    行ベクトルの参照
00913     typedef const row_vector_type&  const_row_vector_ref;   ///<    行ベクトルの参照
00914     typedef TVectorSlice<H,STR,TVector<H*STR,T> >
00915                                     col_vector_type;        ///<    列ベクトル型
00916     typedef col_vector_type&        col_vector_ref;         ///<    列ベクトルの参照
00917     typedef const col_vector_type&  const_col_vector_ref;   ///<    列ベクトルの参照
00918 };
00919 template <class EXP, class TRANS, size_t H, size_t W, size_t STR, class T, class Z=T, class U=Z>
00920 class TMatrixDescCol: public TMatrixDescBase<EXP,TRANS,H,W,STR,T,Z,U>{
00921 public:
00922     typedef TMakeSubMatrixCol< TMatrixDescCol<EXP,TRANS,H,W,STR,T,Z,U> > make_sub_matrix;
00923     typedef TVectorSlice<W,STR,TVector<W*STR,T> >
00924                                     row_vector_type;        ///<    行ベクトル型
00925     typedef row_vector_type&        row_vector_ref;         ///<    行ベクトルの参照
00926     typedef const row_vector_type&  const_row_vector_ref;   ///<    行ベクトルの参照
00927     typedef TVector<H,T>            col_vector_type;        ///<    列ベクトル型
00928     typedef col_vector_type&        col_vector_ref;         ///<    列ベクトルの参照
00929     typedef const col_vector_type&  const_col_vector_ref;   ///<    列ベクトルの参照
00930 };
00931 template <class EXP, class TRANS, size_t H, size_t W, class OD>
00932 class TSubMatrixDescRow: public TMatrixDescRow<EXP, TRANS, H, W, DIMDEC(OD::STRIDE), TYPENAME OD::element_type, TYPENAME OD::zero, TYPENAME OD::unit>{
00933 public:
00934     typedef TMakeSubMatrixRow< TSubMatrixDescRow<EXP,TRANS,H,W,OD> > make_sub_matrix;
00935 };
00936 template <class EXP, class TRANS, size_t H, size_t W, class OD>
00937 class TSubMatrixDescCol: public TMatrixDescCol<EXP, TRANS, H, W, DIMDEC(OD::STRIDE), TYPENAME OD::element_type, TYPENAME OD::zero, TYPENAME OD::unit>{
00938 public:
00939     typedef TMakeSubMatrixCol< TSubMatrixDescCol<EXP,TRANS,H,W,OD> > make_sub_matrix;
00940 };
00941 
00942 template <size_t H, size_t W, class T, class Z=T, class U=Z> class TMatrixCol;
00943 /** 行列型. TMatrixRow<3,3, float> m; のように使う
00944     @param  H   行数．
00945     @param  W   列数．
00946     @param  T   要素の型.   */
00947 template <size_t H, size_t W, class T, class Z=T, class U=Z>
00948 class TMatrixRow:public TMatrixBase<DIMENC(H), DIMENC(W), TMatrixDescRow<TMatrixRow<H,W,T,Z,U>, TMatrixCol<W,H,T,Z,U>, H, W, W, T, Z, U> >{
00949 public:
00950     /// 型情報
00951     typedef TMatrixDescRow<TMatrixRow<H,W,T,Z,U>, TMatrixCol<W,H,T,Z,U>, H, W, W, T, Z, U> desc;
00952     typedef TMatrixBase<DIMENC(H),DIMENC(W),desc> base_type;        ///<    基本クラス型
00953     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
00954     DEF_MATRIX_BASIC_MEMBER(TMatrixRow);
00955 
00956 public:
00957     /// 要素のアクセス
00958     element_type& item_impl(size_t i, size_t j){ return data[i][j]; }
00959     const element_type& item_impl(size_t i, size_t j) const { return data[i][j]; }
00960 private:
00961     element_type data[H][W];                                        ///<    データ
00962 };
00963 
00964 /** 列行列型. TMatrixCol<3,3, float> m; のように使う
00965     @param  H   行数．
00966     @param  W   列数．
00967     @param  T   要素の型.   */
00968 template <size_t H, size_t W, class T, class Z, class U>
00969 class TMatrixCol:public TMatrixBase<DIMENC(H), DIMENC(W), TMatrixDescCol<TMatrixCol<H,W,T,Z,U>, TMatrixRow<W,H,T,Z,U>, H,W,H,T,Z,U> >{
00970 public:
00971     /// 型情報
00972     typedef TMatrixDescCol<TMatrixCol<H,W,T,Z,U>, TMatrixRow<H,W,T,Z,U>, H, W, H, T, Z, U> desc;
00973     typedef TMatrixBase<DIMENC(H),DIMENC(W),desc> base_type;        ///<    基本クラス型
00974     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
00975     DEF_MATRIX_BASIC_MEMBER(TMatrixCol);
00976 
00977 public:
00978     /// 要素のアクセス
00979     element_type& item_impl(size_t i, size_t j){ return data[j][i]; }
00980     const element_type& item_impl(size_t i, size_t j) const { return data[j][i]; }
00981 private:
00982     element_type data[W][H];                ///<    データ
00983 };
00984 
00985 template <size_t H, size_t W, class OD> class TSubMatrixCol;
00986 
00987 /// 部分行列(テンプレート版)
00988 template <size_t H, size_t W, class OD>
00989 class TSubMatrixRow:public TMatrixBase<DIMENC(H),DIMENC(W),TSubMatrixDescRow<TSubMatrixRow<H,W,OD>,TSubMatrixCol<W,H,OD>, H,W,OD> >{
00990 public:
00991     typedef TSubMatrixDescRow<TSubMatrixRow<H,W,OD>,TSubMatrixCol<W,H,OD>, H,W,OD> desc;
00992     typedef TMatrixBase<DIMENC(H),DIMENC(W),desc> base_type;
00993     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
00994     DEF_MATRIX_BASIC_MEMBER(TSubMatrixRow);
00995     DIMDEF(base_type::HEIGHT, HEIGHT);
00996     DIMDEF(base_type::WIDTH, WIDTH);
00997 
00998     /// 要素のアクセス
00999     element_type& item_impl(size_t i, size_t j){ return data[i][j]; }
01000     const element_type& item_impl(size_t i, size_t j) const { return data[i][j]; }
01001 protected:
01002     element_type data[DIMDEC(HEIGHT)][DIMDEC(STRIDE)];
01003 };
01004 
01005 /// 部分行列(テンプレート版)
01006 template <size_t H, size_t W, class OD>
01007 class TSubMatrixCol:public TMatrixBase<DIMENC(H),DIMENC(W),TSubMatrixDescCol<TSubMatrixCol<H,W,OD>,TSubMatrixRow<W,H,OD>, H,W,OD> >{
01008 public:
01009     typedef TSubMatrixDescCol<TSubMatrixCol<H,W,OD>,TSubMatrixRow<W,H,OD>, H,W,OD> desc;
01010     typedef TMatrixBase<DIMENC(H),DIMENC(W),desc> base_type;
01011     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
01012     DEF_MATRIX_BASIC_MEMBER(TSubMatrixCol);
01013     DIMDEF(base_type::HEIGHT, HEIGHT);
01014     DIMDEF(base_type::WIDTH, WIDTH);
01015 
01016     /// 要素のアクセス
01017     element_type& item_impl(size_t i, size_t j){ return data[j][i]; }
01018     const element_type& item_impl(size_t i, size_t j) const { return data[j][i]; }
01019 protected:
01020     element_type data[DIMDEC(WIDTH)][DIMDEC(STRIDE)];
01021 };
01022 
01023 //-----------------------------------------------------------------------------
01024 //  サイズを変数で持つベクトル
01025 //
01026 
01027 template <class T> class EMatrixRow;
01028 template <class T> class VMatrixRow;
01029 template <class T> class EMatrixCol;
01030 template <class T> class VMatrixCol;
01031 template <class EXP, class TRANS, class T, class Z=T, class U=Z> class EMatrixDescBase{
01032 public:
01033     typedef EXP                         exp_type;               ///<    実体
01034     typedef T                           element_type;           ///<    要素の型
01035     typedef Z                           zero;                   ///<    zero(0)が 0 を返す型
01036     typedef U                           unit;                   ///<    unit(1)が 1 を返す型
01037     typedef TRANS                       trans_type;             ///<    転地行列
01038     typedef trans_type                  trans_ref;              ///<    転地行列の参照型
01039     typedef trans_type                  const_trans_ref;        ///<    転地行列の参照型
01040 };
01041 template <class EXP, class TRANS, class T, class Z=T, class U=Z>
01042 class EMatrixDescRow: public EMatrixDescBase<EXP, TRANS, T, Z, U>{
01043 public:
01044     typedef VMatrixRow<T>                           ret_type;               ///<    返り値型
01045     typedef EVector<T>                              row_vector_type;        ///<    行ベクトル型
01046     typedef row_vector_type                         row_vector_ref;         ///<    列ベクトルの参照
01047     typedef TYPENAME row_vector_type::const_type    const_row_vector_ref;   ///<    行ベクトルの参照
01048     typedef EVectorSlice<T>                         col_vector_type;        ///<    列ベクトル型
01049     typedef col_vector_type                         col_vector_ref;         ///<    行ベクトルの参照
01050     typedef TYPENAME col_vector_type::const_type    const_col_vector_ref;   ///<    列ベクトルの参照
01051     typedef EMakeSubMatrixRow< EMatrixDescRow<EXP, TRANS, T, Z, U> > make_sub_matrix;
01052 };
01053 template <class EXP, class TRANS, class T, class Z=T, class U=Z>
01054 class EMatrixDescCol: public EMatrixDescBase<EXP, TRANS, T, Z, U>{
01055 public:
01056     typedef VMatrixCol<T>                           ret_type;               ///<    返り値型
01057     typedef EVectorSlice<T>                         row_vector_type;        ///<    行ベクトル型
01058     typedef row_vector_type                         row_vector_ref;         ///<    列ベクトルの参照
01059     typedef TYPENAME row_vector_type::const_type    const_row_vector_ref;   ///<    行ベクトルの参照
01060     typedef EVector<T>                              col_vector_type;        ///<    列ベクトル型
01061     typedef col_vector_type                         col_vector_ref;         ///<    行ベクトルの参照
01062     typedef TYPENAME col_vector_type::const_type    const_col_vector_ref;   ///<    列ベクトルの参照
01063     typedef EMakeSubMatrixCol< EMatrixDescCol<EXP, TRANS, T, Z, U> > make_sub_matrix;
01064 };
01065 
01066 /// サイズを変数でもつ行列の基本型
01067 template <class D>
01068 class EMatrixBase:public MatrixImp<D>{
01069 public:
01070     /// 型情報
01071     typedef D desc;
01072     typedef MatrixImp<desc> base_type;      ///<    基本クラス型
01073     DEF_MATRIXD_BASIC_MEMBER(EMatrixBase);
01074 
01075     size_t height_impl() const { return height_; }
01076     size_t width_impl() const { return width_; }
01077         
01078     /// 行ベクトル
01079     row_vector_ref row_impl(size_t n){
01080         return row_vector_ref(width(), 1, &item(n,0)); 
01081     }
01082     const_row_vector_ref row_impl(size_t n) const {
01083         return const_row_vector_ref(width(), 1, &item(n,0));
01084     }
01085     /// 列ベクトル
01086     col_vector_ref col_impl(size_t m){
01087         return col_vector_ref(height(), stride(), &item(0,m));
01088     }
01089     const_col_vector_ref col_impl(size_t m) const {
01090         return const_col_vector_ref(height(), stride(), &item(0,m));
01091     }
01092     /// 転置
01093     trans_ref trans_impl(){
01094         return trans_ref(width(), height(), stride(), &item(0,0));
01095     }
01096     const_trans_ref trans_impl() const {
01097         return const_trans_ref(width(), height(), stride(), &item(0,0));
01098     }
01099     ///
01100     void resize_impl(size_t h, size_t w) { assert(h==height_ && w==width_);}
01101 
01102 protected:
01103     size_t height_;
01104     size_t width_;
01105     element_type* data;                     ///<    データ
01106     void init_buffer(){height_=0; width_=0; data=0; }
01107     EMatrixBase(size_t h, size_t w, const element_type* d):height_(h), width_(w), data((element_type*)d){}
01108     EMatrixBase():height_(0), width_(0), data(0){}
01109 };
01110 
01111 /// サイズを変数でもつ行列の基本型
01112 template <class D>
01113 class EMatrixBaseRow:public EMatrixBase<D>{
01114 protected:
01115     EMatrixBaseRow(){}
01116     EMatrixBaseRow(const EMatrixBaseRow& m): EMatrixBase<D>(m){}
01117 public:
01118     typedef D desc;
01119     typedef EMatrixBase<D> base_type;
01120     DEF_MATRIXD_BASIC_MEMBER(EMatrixBaseRow)
01121     EMatrixBaseRow(size_t h, size_t w, const element_type* d):EMatrixBase<D>(h,w,d){}
01122     /// 要素のアクセス
01123     element_type& item_impl(size_t i, size_t j){ return data[i*stride()+j]; }
01124     const element_type& item_impl(size_t i, size_t j) const { return data[i*stride()+j]; }
01125     /// ストライド
01126     size_t stride_impl() const { return width_; }
01127     /// 行ベクトル
01128     row_vector_ref row_impl(size_t n){
01129         return row_vector_ref(width(), 1, &item(n,0));
01130     }
01131     const_row_vector_ref row_impl(size_t n) const {
01132         return const_row_vector_ref(width(), 1, &item(n,0));
01133     }
01134     /// 列ベクトル
01135     col_vector_ref col_impl(size_t m){
01136         return col_vector_ref(height(), stride(), &item(0,m));
01137     }
01138     const_col_vector_ref col_impl(size_t m) const {
01139         return const_col_vector_ref(height(), stride(), &item(0,m));
01140     }
01141 };
01142 
01143 /// サイズを変数でもつ行列の基本型
01144 template <class D>
01145 class EMatrixBaseCol:public EMatrixBase<D>{
01146 protected:
01147     EMatrixBaseCol(){}
01148     EMatrixBaseCol(const EMatrixBaseCol& m): EMatrixBase<D>(m){}
01149 public:
01150     typedef D desc;
01151     typedef EMatrixBase<D> base_type;
01152     DEF_MATRIXD_BASIC_MEMBER(EMatrixBaseCol)
01153     EMatrixBaseCol(size_t h, size_t w, const element_type* d):EMatrixBase<D>(h,w,d){}
01154     /// 要素のアクセス
01155     element_type& item_impl(size_t i, size_t j){ return data[j*stride()+i]; }
01156     const element_type& item_impl(size_t i, size_t j) const { return data[j*stride()+i]; }
01157     /// ストライド
01158     size_t stride_impl() const { return height_; }
01159     /// 行ベクトル
01160     row_vector_ref row_impl(size_t n){
01161         return row_vector_ref(width(), stride(), &item(n,0));
01162     }
01163     const_row_vector_ref row_impl(size_t n) const {
01164         return const_row_vector_ref(width(), stride(), &item(n,0));
01165     }
01166     /// 列ベクトル
01167     col_vector_ref col_impl(size_t m){
01168         return col_vector_ref(height(), 1, &item(0,m));
01169     }
01170     const_col_vector_ref col_impl(size_t m) const {
01171         return const_col_vector_ref(height(), 1, &item(0,m));
01172     }
01173 };
01174 
01175 /** 行列型. EMatrixRow<float> m(3,3,buf); のように使う
01176     @param  T   要素の型.   */
01177 template <class T>
01178 class EMatrixRow:public EMatrixBaseRow< EMatrixDescRow<EMatrixRow<T>, EMatrixCol<T>, T> >{
01179 public:
01180     /// 型情報
01181     typedef EMatrixDescRow<EMatrixRow<T>, EMatrixCol<T>, T> desc;
01182     typedef EMatrixBaseRow<desc> base_type;     ///<    基本クラス型
01183     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
01184     DEF_MATRIX_BASIC_MEMBER(EMatrixRow);
01185     EMatrixRow(const EMatrixRow& m):base_type(m.height_, m.width_,m.data){}
01186     EMatrixRow(size_t h, size_t w, size_t str, const element_type* d):base_type(h,w,d){ assert(str == w); }
01187 };
01188 
01189 /** 行列型. VMatrixRow<float> m(3,3); のように使う
01190     @param  T   要素の型.   */
01191 template <class T>
01192 class VMatrixRow:public EMatrixBaseRow< EMatrixDescRow<VMatrixRow<T>, EMatrixCol<T>, T> >{
01193 public:
01194     /// 型情報
01195     typedef EMatrixDescRow<VMatrixRow<T>, EMatrixCol<T>, T> desc;
01196     typedef EMatrixBaseRow<desc> base_type;     ///<    基本クラス型
01197     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
01198     DEF_MATRIX_BASIC_MEMBER(VMatrixRow);
01199     VMatrixRow(const VMatrixRow& m){ init_buffer(); assign(m); }
01200     ~VMatrixRow(){ delete [] data; }
01201     ///
01202     void resize_impl(size_t h, size_t w) {
01203         if (height()*width()<h*w){
01204             delete [] data;
01205             data = new T[h*w];
01206             height_ = h;
01207             width_ = w;
01208         }
01209     }
01210 };
01211 
01212 /** 行列型. EMatrixCol<float> m(3,3,buf); のように使う
01213     @param  T   要素の型.   */
01214 template <class T>
01215 class EMatrixCol:public EMatrixBaseCol< EMatrixDescCol<EMatrixCol<T>, EMatrixRow<T>, T> >{
01216 public:
01217     /// 型情報
01218     typedef EMatrixDescCol<EMatrixCol<T>, EMatrixRow<T>, T> desc;
01219     typedef EMatrixBaseCol<desc> base_type;     ///<    基本クラス型
01220     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
01221     DEF_MATRIX_BASIC_MEMBER(EMatrixCol);
01222     EMatrixCol(const EMatrixCol& m):base_type(m.height_, m.width_, m.data){}
01223     EMatrixCol(size_t h, size_t w, size_t str, const element_type* d):base_type(h,w,d){ assert(str == h); }
01224 };
01225 
01226 /** 行列型. VMatrixCol<float> m(3,3); のように使う
01227     @param  T   要素の型.   */
01228 template <class T>
01229 class VMatrixCol:public EMatrixBaseCol< EMatrixDescCol<VMatrixCol<T>, EMatrixRow<T>, T> >{
01230 public:
01231     /// 型情報
01232     typedef EMatrixDescCol<VMatrixCol<T>, EMatrixRow<T>, T> desc;
01233     typedef EMatrixBaseCol<desc> base_type;     ///<    基本クラス型
01234     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
01235     DEF_MATRIX_BASIC_MEMBER(VMatrixCol);
01236     VMatrixCol(const VMatrixCol& m){ init_buffer(); assign(m); }
01237     ~VMatrixCol(){ delete [] data; }
01238     ///
01239     void resize_impl(size_t h, size_t w) {
01240         if (height()*width()<h*w){
01241             delete [] data;
01242             data = new T[h*w];
01243         }
01244         height_ = h;
01245         width_ = w;
01246     }
01247 };
01248 
01249 template <class T, class Z, class U> class ESubMatrixCol;
01250 /// 部分行列(変数版)
01251 template <class T, class Z, class U>
01252 class ESubMatrixRow:public EMatrixBaseRow< EMatrixDescRow<ESubMatrixRow<T>, ESubMatrixCol<T>, T,Z,U> >{
01253 public:
01254     typedef EMatrixDescRow<ESubMatrixRow<T>, ESubMatrixCol<T>, T,Z,U> desc;
01255     typedef EMatrixBaseRow<desc> base_type;
01256     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
01257     DEF_MATRIX_BASIC_MEMBER(ESubMatrixRow);
01258     ///
01259     ESubMatrixRow(size_t h, size_t w, size_t str, const element_type* d):base_type(h, w, d), stride_(str){}
01260     ///
01261     ESubMatrixRow(const ESubMatrixRow& m):base_type(m.height_, m.width_, m.data), stride_(m.stride_){}
01262     ///
01263     size_t stride_impl() const { return stride_; }
01264 protected:
01265     size_t stride_;
01266     void init_buffer() { height_=0; width_=0; stride_=0; data=0; }
01267 };
01268 
01269 /// 部分行列(変数版)
01270 template <class T, class Z, class U>
01271 class ESubMatrixCol:public EMatrixBaseCol< EMatrixDescCol<ESubMatrixCol<T>, ESubMatrixRow<T>, T,Z,U> >{
01272 public:
01273     typedef EMatrixDescCol<ESubMatrixCol<T>, ESubMatrixRow<T>, T,Z,U> desc;
01274     typedef EMatrixBaseCol<desc> base_type;
01275     /// 継承されない基本的なメンバの定義. @see ::DEF_MATRIX_BASIC_MEMBER
01276     DEF_MATRIX_BASIC_MEMBER(ESubMatrixCol);
01277     ///
01278     ESubMatrixCol(size_t h, size_t w, size_t str, const element_type* d):base_type(h,w,d), stride_(str){}
01279     ///
01280     ESubMatrixCol(const ESubMatrixCol& m):base_type(m.height_, m.width_, m.data), stride_(m.stride_){}
01281     ///
01282     size_t stride_impl() const { return stride_; }
01283 protected:
01284     size_t stride_;
01285     void init_buffer() { height_=0; width_=0; stride_=0; data=0; }
01286 };
01287 
01288 
01289 //----------------------------------------------------------------------------
01290 /// @name 行列の演算子
01291 //@{
01292 /// 行列の比較
01293 template <class AD, class BD>
01294 bool operator == (const MatrixImp<AD>& a, const MatrixImp<BD>& b){
01295     return a.equal(b);
01296 }
01297 /// 行列の比較
01298 template <class AD, class BD>
01299 bool operator != (const MatrixImp<AD>& a, const MatrixImp<BD>& b){
01300     return !a.equal(b);
01301 }
01302 /// 行列の和
01303 template <class AD, class BD>
01304 TYPENAME AD::ret_type operator + (const MatrixImp<AD>& a, const MatrixImp<BD>& b){
01305     TYPENAME AD::ret_type r(a);
01306     r.add(b);
01307     return r;
01308 }
01309 /// 行列の差
01310 template <class AD, class BD>
01311 TYPENAME AD::ret_type operator - (const MatrixImp<AD>& a, const MatrixImp<BD>& b){
01312     TYPENAME AD::ret_type r(a);
01313     r.sub(b);
01314     return r;
01315 }
01316 /// 行列の掛け算    サイズ：テンプレート
01317 template <DIMTYPE AH, DIMTYPE AW, class AD, DIMTYPE BW, class BD>
01318 TMatrixCol<DIMDEC(AH), DIMDEC(BW), TYPENAME AD::element_type> operator * (
01319     const TMatrixBase<AH, AW, AD>& a,
01320     const TMatrixBase<AW, BW, BD>& b){
01321     TMatrixCol<DIMDEC(AH), DIMDEC(BW), TYPENAME AD::element_type> r;
01322     multi(r, a, b);
01323     return r;
01324 }
01325 /// 行列の掛け算    サイズ：変数
01326 template <class AD, class BD>
01327 VMatrixCol<TYPENAME AD::element_type> operator * (const MatrixImp<AD>& a, const MatrixImp<BD>& b){
01328     VMatrixCol<TYPENAME AD::element_type> r;
01329     multi(r, a, b);
01330     return r;
01331 }
01332 
01333 //  ベクトルと行列の掛け算
01334 /// ベクトルと行列の掛け算  サイズ：テンプレート
01335 template <DIMTYPE AH, DIMTYPE AW, class AD, class BD>
01336 TVector<DIMDEC(AH), TYPENAME AD::element_type> operator * (const TMatrixBase<AH, AW, AD>& a, const VectorImp<BD>& b){
01337     TVector<DIMDEC(AH), TYPENAME AD::element_type> r;
01338     multi(r, a, b);
01339     return r;
01340 }
01341 /// ベクトルと行列の掛け算  サイズ：変数
01342 template <class AD, class BD>
01343 VVector<TYPENAME BD::element_type> operator * (const MatrixImp<AD>& a, const VectorImp<BD>& b){
01344     VVector<TYPENAME BD::element_type> r;
01345     multi(r, a, b);
01346     return r;
01347 }
01348 /// ベクトルと行列の掛け算  サイズ：テンプレート
01349 template <class AD, DIMTYPE BH, DIMTYPE BW, class BD>
01350 TVector<DIMDEC(BW), TYPENAME AD::element_type> operator * (const VectorImp<AD>& a, const TMatrixBase<BH, BW, BD>& b){
01351     TVector<DIMDEC(BW), TYPENAME AD::element_type> r;
01352     multi(r, b.trans(), a);
01353     return r;
01354 }
01355 /// ベクトルと行列の掛け算  サイズ：変数
01356 template <class AD, class BD>
01357 VVector<TYPENAME AD::element_type> operator * (const VectorImp<AD>& a, const MatrixImp<BD>& b){
01358     VVector<TYPENAME AD::element_type> r;
01359     multi(r, b.trans(), a);
01360     return r;
01361 }
01362 /// 行列の定数倍
01363 template <class D>
01364 TYPENAME D::ret_type operator * (const MatrixImp<D>& a, TYPENAME D::element_type b){
01365     TYPENAME D::ret_type r(a);
01366     r.multi(b);
01367     return r;
01368 }
01369 /// 行列の定数倍
01370 template <class D>
01371 TYPENAME D::ret_type operator * (TYPENAME D::element_type a, const MatrixImp<D>& b){
01372     TYPENAME D::ret_type r(b);
01373     r.multi(a);
01374     return r;
01375 }
01376 
01377 /// 表示
01378 template <class D>
01379 std::ostream& operator << (std::ostream& os, const MatrixImp<D>& m){
01380     m.print(os);
01381     return os;
01382 }
01383 /// 入力
01384 template <class D>
01385 std::istream& operator >> (std::istream& is, MatrixImp<D>& m){
01386     m.input(is);
01387     return is;
01388 }
01389 
01390 #ifdef _WIN32
01391  #ifdef _DEBUG
01392   #pragma optimize ("", on)
01393   #pragma auto_inline(off)
01394   #pragma inline_recursion(off)
01395  #endif
01396  #pragma pack(pop)
01397 #endif
01398 
01399 }   //  namespace PTM
01400 #endif
Generated on Sun Apr 16 02:07:13 2006 for Springhead by  doxygen 1.4.1