Matrixクラス 行列/ベクトル計算クラス
FRONT PAGE

コレは何?/What is this?

C++上で行列/ベクトルの計算を超絶簡単にするクラス。 ヘッダをインクルードするだけで使用可能!
これで,C++でMATLABと同じように気軽に行列計算ができます。
(「Eigen」の車輪の再開発のようなもの。機能の充実性と超大規模な行列の計算速度を求めるならば,このクラスの代わりに「Eigen」の導入をオススメしておく。 なんでこんなん作ってるのかというと,権利関係を完全にクリアにしたかったので。)

Matrixクラスの主な機能

ダウンロード

ソースコード一式は以下からダウンロードできます。

gcc 4.8 以上を推奨。とりあえずMSYS2+MinGW64,Scientific Linux 6.5で動作確認済み。ARCS上での動作も確認。
注意!C++11じゃないと動きません!!(一部,templateとconstexprを使用しているので)

ファイル一覧

コンパイルの前に…

コンパイルする前に,自分の処理系に合わせてMakefileの中のCFLAGSをいじって下さい。デフォルトではMSYS2+MINGW64用になってます。

コンパイル&実行

make でコンパイル。その後に ./MatrixTest を起動すると,LU分解の計算例が表示されます。

使用可能な演算子

行列/ベクトル/スカラー間の,+(加算),-(減算),*(乗算),+=(加算代入),-=(減算代入),^(冪乗) の演算子が使えます。

関数一覧

各々の関数における計算結果はMATLABと比較して合っていることを確認している。
ただし,LU分解などの答えが一意に求まらない関数では,当然,MATLABとは異なる結果を出力する。

void SetElem(unsigned int m, unsigned int n, double val);// 指定した要素番号に値を設定する関数 m行目,n列目にvalを代入
double GetElem(unsigned int m, unsigned int n);// m行n列目の値を返す関数
void FillAll(double u);// すべての要素を指定した値で埋める関数
unsigned int GetRowLength(void);// 行の長さを返す関数 紛らわしいので名称変更
unsigned int GetColumnLength(void);// 列の長さを返す関数 紛らわしいので名称変更
unsigned int GetWidthLength(void);// 行列の幅(列数)を返す関数
unsigned int GetHeightLength(void);// 行列の高さ(行数)を返す関数
void _PrintMatSize(const Matrix& u, const std::string varname);// 行列のサイズの表示 u : 表示する行列 この関数はマクロを介して呼ばれることを想定している。
void _PrintMatrix(const Matrix& u, const std::string format, const std::string varname);// 行列の要素を表示 format : 表示形式 (%1.3e とか %5.3f とか printfと同じ) u : 表示する行列 この関数はマクロを介して呼ばれることを想定している。
Matrix tp(const Matrix& u);// 転置行列を返す関数
double tr(const Matrix& u);// 行列のトレースを返す関数
double prod(const Matrix& u);// 行列の対角要素の総積を返す関数
Matrix diag(const Matrix& u);// 行列の対角要素を返す関数
double det(const Matrix& u);// 行列式の値を返す関数
Matrix inv(const Matrix& A);// 逆行列を返す関数 (LU分解とAx=bの解から求めるので高速なはず)
Matrix lpinv(const Matrix& A);// 左擬似逆行列を返す関数 (Aが縦長行列の場合)
Matrix rpinv(const Matrix& A);// 右擬似逆行列を返す関数 (Aが横長行列の場合)
Matrix abs(const Matrix& U);// 行列要素の絶対値を返す関数
Matrix sumrow(const Matrix& U);// 行方向へ加算して列ベクトルを出力する関数
Matrix sumcolumn(const Matrix& U);// 列方向へ加算して行ベクトルを出力する関数
double max(const Matrix& u);// ベクトル要素の最大値を返す関数
double absmax(const Matrix& u);// ベクトル要素の絶対値の最大値を返す関数
unsigned int maxidx(const Matrix& u);// ベクトル要素の最大値の要素番号を返す関数
unsigned int absmaxidx(const Matrix& u);// ベクトル要素の絶対値の最大値の要素番号を返す関数
double infnorm(const Matrix& U);// 行列の絶対値ノルムを返す関数
Matrix getrow(const Matrix& u, unsigned int m);// 指定した行を抽出する関数
void setrow(Matrix& u, const Matrix& v, unsigned int m);// 指定した行を行ベクトルで上書きする関数
void swaprow(Matrix& U, unsigned int m1, unsigned int m2);// 指定した行と行を入れ替える関数
void fillrow(Matrix& U, double a, unsigned int m, unsigned int n1, unsigned int n2);// m行目のn1列目からm2列目を数値aで埋める関数  (n1 <= n2 であること)
Matrix getcolumn(const Matrix& u, unsigned int n);// 指定した列を抽出する関数
void setcolumn(Matrix& u, const Matrix& v, unsigned int n);// 指定した列を列ベクトルで上書きする関数
void swapcolumn(Matrix& U, unsigned int n1, unsigned int n2);// 指定した列と列を入れ替える関数
void fillcolumn(Matrix& U, double a, unsigned int n, unsigned int m1, unsigned int m2);// n列目のm1行目からm2行目を数値aで埋める関数  (m1 <= m2 であること)
Matrix orderrow(const Matrix& U, const Matrix& v);// 並び替え記憶列ベクトルvに従って,入力行列Uの行を並び替える関数
Matrix reorderrow(const Matrix& U, const Matrix& v);// 並び替え記憶列ベクトルvが昇順になるように,入力行列Uの行を並び替えて元に戻す関数
Matrix solve(const Matrix& A, const Matrix& b);// Ax = b の形の線形連立1次方程式をxについて解く関数
Matrix expm(const Matrix& A, unsigned int M);// 行列指数関数 e^(A) A : 入力行列,M : パデ近似の次数
Matrix integral_expm(const Matrix& U, const double T, const unsigned long N, const unsigned int P);// 指数行列の数値定積分[0,T)をする関数 U:入力行列,T:積分範囲の終わり,N:分割数,P:パデ近似の次数
Matrix::LUperm LU(const Matrix& A, Matrix& L, Matrix& U, Matrix& v);// LU分解  A:入力行列,L:下三角行列,U:上三角行列,v:並べ替え記憶列ベクトル 並べ替えが偶数回/奇数回発生したことを返す。返す値は Matrix::ODD か Matrix::EVEN のどちらか
void Set(const T1& u1, const T2&... u2);// 行列要素に値を設定する関数 (非スレッドセーフなので注意)
Matrix zeros(unsigned int n, unsigned int m);// m行n列の零行列を返す関数
Matrix ones(unsigned int n, unsigned int m);// m行n列の要素がすべて1の行列を返す関数
Matrix ident(unsigned int n);// n行n列の単位行列を返す関数
Matrix eye(unsigned int n);// n行n列の単位行列を返す関数

定義一覧

#define PrintMatSize(a) (_PrintMatSize(a,#a))    // 行列サイズ表示マクロ
#define PrintMatrix(a,b) (_PrintMatrix(a,b,#a))    // 行列要素表示マクロ

// LU分解の際の並べ替えが偶数回/奇数回発生したことを返すための定義
enum LUperm {
    ODD,    // 奇数
    EVEN    // 偶数
};

使用例

<行列要素への代入>

コード例 実行結果
Matrix A(4,4);
A.Set(
    2.6,  1.8,  8.7,  8.5,
    8,    2.6,  5.8,  6.2,
    4.3,  1.5,  5.5,  3.5,
    9.1,  1.4,  1.4,  5.1
);
PrintMatrix(A,"%6.3f");
A =
[  2.600  1.800  8.700  8.500 ]
[  8.000  2.600  5.800  6.200 ]
[  4.300  1.500  5.500  3.500 ]
[  9.100  1.400  1.400  5.100 ]

<加減乗と冪>

コード例 実行結果 (一連の計算は続けて実行)
A = A*2;
PrintMatrix(A,"%6.3f");
A =
[  5.200  3.600 17.400 17.000 ]
[ 16.000  5.200 11.600 12.400 ]
[  8.600  3.000 11.000  7.000 ]
[ 18.200  2.800  2.800 10.200 ]
A = 0.5*A;
PrintMatrix(A,"%6.3f");
A =
[  2.600  1.800  8.700  8.500 ]
[  8.000  2.600  5.800  6.200 ]
[  4.300  1.500  5.500  3.500 ]
[  9.100  1.400  1.400  5.100 ]
A = A + 1;
PrintMatrix(A,"%6.3f");
A =
[  3.600  2.800  9.700  9.500 ]
[  9.000  3.600  6.800  7.200 ]
[  5.300  2.500  6.500  4.500 ]
[ 10.100  2.400  2.400  6.100 ]
A = A - 1;
PrintMatrix(A,"%6.3f");
A =
[  2.600  1.800  8.700  8.500 ]
[  8.000  2.600  5.800  6.200 ]
[  4.300  1.500  5.500  3.500 ]
[  9.100  1.400  1.400  5.100 ]
A = 1.2 + A;
PrintMatrix(A,"%6.3f");
A =
[  3.800  3.000  9.900  9.700 ]
[  9.200  3.800  7.000  7.400 ]
[  5.500  2.700  6.700  4.700 ]
[ 10.300  2.600  2.600  6.300 ]
A = 1.2 - A;
PrintMatrix(A,"%6.3f");
A =
[ -2.600 -1.800 -8.700 -8.500 ]
[ -8.000 -2.600 -5.800 -6.200 ]
[ -4.300 -1.500 -5.500 -3.500 ]
[ -9.100 -1.400 -1.400 -5.100 ]
A = +A;
PrintMatrix(A,"%6.3f");
A =
[ -2.600 -1.800 -8.700 -8.500 ]
[ -8.000 -2.600 -5.800 -6.200 ]
[ -4.300 -1.500 -5.500 -3.500 ]
[ -9.100 -1.400 -1.400 -5.100 ]
A = -A;
PrintMatrix(A,"%6.3f");
A =
[  2.600  1.800  8.700  8.500 ]
[  8.000  2.600  5.800  6.200 ]
[  4.300  1.500  5.500  3.500 ]
[  9.100  1.400  1.400  5.100 ]
A += A;
PrintMatrix(A,"%6.3f");
A =
[  5.200  3.600 17.400 17.000 ]
[ 16.000  5.200 11.600 12.400 ]
[  8.600  3.000 11.000  7.000 ]
[ 18.200  2.800  2.800 10.200 ]
A -= A*0.5;
PrintMatrix(A,"%6.3f");
A =
[  2.600  1.800  8.700  8.500 ]
[  8.000  2.600  5.800  6.200 ]
[  4.300  1.500  5.500  3.500 ]
[  9.100  1.400  1.400  5.100 ]
A += 3.14;
PrintMatrix(A,"%6.3f");
A =
[  5.740  4.940 11.840 11.640 ]
[ 11.140  5.740  8.940  9.340 ]
[  7.440  4.640  8.640  6.640 ]
[ 12.240  4.540  4.540  8.240 ]
A -= 3.14;
PrintMatrix(A,"%6.3f");
A =
[  2.600  1.800  8.700  8.500 ]
[  8.000  2.600  5.800  6.200 ]
[  4.300  1.500  5.500  3.500 ]
[  9.100  1.400  1.400  5.100 ]
Matrix B = A*A;
PrintMatrix(B,"%6.3f");
B =
[ 135.920 34.310 92.810 107.060 ]
[ 122.960 38.540 125.260 136.040 ]
[ 78.680 24.790 81.260 82.950 ]
[ 87.290 29.260 102.130 116.940 ]
B = A^2;
PrintMatrix(B,"%6.3f");
B =
[ 135.920 34.310 92.810 107.060 ]
[ 122.960 38.540 125.260 136.040 ]
[ 78.680 24.790 81.260 82.950 ]
[ 87.290 29.260 102.130 116.940 ]
ソースコードの例 → MatrixTest1.cc

<LU分解>

コード例 実行結果
Matrix A(4,4);
A.Set(
    26,    18,    87,    85,
    80,    26,    58,    62,
    43,    15,    55,    35,
    91,    14,    14,    51
);
PrintMatrix(A,"%6.3f");

Matrix L(4,4), U(4,4), v(1,4);

LU(A,L,U,v);

PrintMatrix(L,"%6.3f");
PrintMatrix(U,"%6.3f");
A =
[ 26.000 18.000 87.000 85.000 ]
[ 80.000 26.000 58.000 62.000 ]
[ 43.000 15.000 55.000 35.000 ]
[ 91.000 14.000 14.000 51.000 ]

L =
[  1.000  0.000  0.000  0.000 ]
[  0.286  1.000  0.000  0.000 ]
[  0.879  0.978  1.000  0.000 ]
[  0.473  0.599  0.037  1.000 ]

U =
[ 91.000 14.000 14.000 51.000 ]
[  0.000 14.000 83.000 70.429 ]
[  0.000  0.000 -35.484 -51.716 ]
[  0.000  0.000  0.000 -29.349 ]

<線形3元連立1次方程式の解 (Ax=bのxを求める)>

コード例 実行結果
Matrix A(3,3), b(1,3), x(1,3);
A.Set(
	1,  1,  1,
	2,  3, -2,
	3, -1,  1
);
b.Set(
	9,
	5,
	7
);
PrintMatrix(A,"%6.4f");
PrintMatrix(b,"%6.4f");

x = solve(A,b);
PrintMatrix(x,"%6.15f");
A =
[ 1.0000 1.0000 1.0000 ]
[ 2.0000 3.0000 -2.0000 ]
[ 3.0000 -1.0000 1.0000 ]

b =
[ 9.0000 ]
[ 5.0000 ]
[ 7.0000 ]

x =
[ 2.000000000000000 ]
[ 3.000000000000000 ]
[ 4.000000000000000 ]

<行列指数関数 (指数行列)>

コード例 実行結果
Matrix A(3,3);
A.Set(
	1,  1,  1,
	2,  3, -2,
	3, -1,  1
);
PrintMatrix(A,"%6.3f");

Matrix Y = expm(A,6);
PrintMatrix(Y,"%1.15e");

printf("det(Y)     = %1.15e\n",det(Y));
printf("exp(tr(A)) = %1.15e\n",exp(tr(A)));
A =
[  1.000  1.000  1.000 ]
[  2.000  3.000 -2.000 ]
[  3.000 -1.000  1.000 ]

Y =
[ 1.020275609460608e+01 1.172847841880298e+01 -1.845697590221385e+00 ]
[ 5.358055492240175e+00 3.818443826855350e+01 -2.345695683760592e+01 ]
[ 1.256180857470163e+01 -2.679027746120072e+00 1.020275609460605e+01 ]

det(Y)     = 1.484131591025755e+02
exp(tr(A)) = 1.484131591025766e+02
det(expm(A)) = exp(tr(A)) の公式により検算でき,14桁の精度があることが分かる。 上記の例ではパデ近似の次数を6としているが,これ以上次数を上げても計算精度はあんまり変わらないみたいだ。 (det(A)が極端に小さい時は未検証)
ちなみにMATLABの計算結果は以下の通り。
expm(A) = 1.020275609460608e+01     1.172847841880297e+01    -1.845697590221378e+00
          5.358055492240140e+00     3.818443826855344e+01    -2.345695683760592e+01
          1.256180857470166e+01    -2.679027746120060e+00     1.020275609460608e+01

det(expm(A))  = 1.484131591025765e+02
exp(trace(A)) = 1.484131591025766e+02
MATLABと比較しても14桁ぐらいまでは合ってるようだ。 これで状態遷移行列が計算できるので,状態方程式の離散化がC++上で コンパイル時 or オンライン で実行できてかなり美味しい。
上記のexpm関数を利用して,MATLABのc2d関数に当たるものを作りました→Discretize離散化関数
その他の例は後々増やしていきます。。。

今後の予定

ライセンス/Licenses

Copyright (C) 2011-2015 Yuki YOKOKURA
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details . Besides, you can negotiate about other options of licenses instead of GPL. If you would like to get other licenses, please contact us.





- 261 -

研究室の横の倉庫 - Side Warehouse of Laboratory
Copyright(C), Side Warehouse, All rights reserved.