TMatrixで行列の行を配列として取得する

TMatrixDの使い方で理解不足のところがあったのでメモ書き程度に。 TMatrixDTMatrixT<double>と、double型の行列要素を扱えるクラスとして定義されているものです。

まず、TMatrixDの定義やら要素の入力などを。

root [] TMatrixD mat(2,2); // 2 x 2行列の生成
// TMatrixTには(rown,coln)というoperatorが定義されているので、この例のように行列要素にアクセス出来ます
root [] mat(0,0) = 1.;
root [] mat(0,1) = 2.;
root [] mat(1,0) = 0.;
root [] mat(1,1) = 2.;
// mat.Print()で、行列を並べて表示される
root [] mat.Print()

2x2 matrix is as follows

     |      0    |      1    |
 -------------------------------
   0 |          1           2
   1 |          0           2
		   

ここで、以下のようにして、行列の列を配列として取り出せます。これで詰まった。

root [] double *ptr;
root [] ptr = mat[0].GetPtr(); // これで行列の0行目をdoubleの配列として取得出来る
root [] ptr[0]
(double) 1.000000
root [] ptr[1]
(double) 2.000000

原理的には、TMatrixDには TMatrixTRow<double> という型の戻り値を持つOperator [] が定義されていて、double* を戻り値に持つTMatrixTRow::GetPtr()を呼び出している事になっています。 また、戻り値はdouble*型であるので、std::inner_product() を使用してベクトルの内積を取ることも出来ます。

root [] std::vector<double> v = {1.,2.};
root [] std::inner_product(v.begin(),v.end(),mat[0].GetPtr(),0.)
(double) 5.000000

あとは基本的な行列計算の方法を書いておきます。行列演算はもとの行列が入れ替わる(関数の戻り値はvoidです)ので注意が必要なのかな

root [] mat.Determinant() // 行列式を求める
(double) -2.0000000
root [] mat.T() // 転置行列に置き換える
root [] mat.Print();

2x2 matrix is as follows

     |      0    |      1    |
 -------------------------------
   0 |          1           0
   1 |          2           2
		   
root [] mat.Invert(); //逆行列に置き換える(4次以上の場合はInvertFast()を推奨されているようです)
root [] mat.Print();

2x2 matrix is as follows

     |      0    |      1    |
 -------------------------------
   0 |          1           0
   1 |         -1         0.5
		   

積に関しては、その積に合うサイズの行列のオブジェクトを定義し、TMatrixT::Mult(TMatrixD &a, TMatrixD &b)関数を用いて、解を代入するような形になります(日本語が下手過ぎるので、実際に例を示します)

root [] TMatrixD mat2(2.2); //空の2 x 2行列を用意
root [] mat2.Mult(mat,mat); //ひとまずmatを二乗して上で用意したmat2に積を詰めます
root [] mat2.Print();

2x2 matrix is as follows

     |      0    |      1    |
 -------------------------------
   0 |          1           0
   1 |       -1.5        0.25
		   

同様に、TMatrixD::TMult(A,B)で $ A^TB $,TMatrixD::MultT(A,B)で $ AB^T $が計算されます。

Shoichiro Masuoka

CNS, the Univ. of Tokyo. Dcotoral student