# Introduction
## 基本クラス
本ライブラリを利用する際に重要なクラスは以下の通りです。
* SparseMat: 実数疎行列のクラス
* SparseMatC: 複素実数疎行列のクラス
* MatSolvers: 線形ソルバクラス
## 疎行列の構築方法
疎行列は、まず行数を指定して実体を作成します。この段階ではc++のmap[]を用いて、一時保存状態になります。
疎行列の値は、addメソッドを用いて代入していきます。FEMなどでは各位置の値が一度に決まることはないと思いますので、addを使って順次足していってください。
仕様上、最初に指定した以上の行に代入するとエラーになるので注意してください。
疎行列の構築が終わったら、fixメソッドを使うことで疎行列の配置を確定させます。具体的には、Eigenの疎行列を作成します。
行列の確定後（fixを呼び出した後）は、四則演算をEigenの機能を用いることで高速に行えます。
（なお、確定後もaddを使って行列の更新は可能です。）
### 疎行列の構築サンプル
疎行列の構築サンプルです。
```
double matA[5][5] ={
{ 61, 0, 3.1, -0.1, 1},
{ 0, 42.5, 2.0, 0.5, 4},
{ 3.1, 2.0, 10.5, 0, 0},
{-0.1, 0.5, 0, 80.4, 1.5},
{ 1, 4, 0, 1.5, 1000.0}
};
size = 5
//疎行列を定義
SRLfem::SparseMat matAs(size);
/* 行列A,B,Cの値をセット */
for(int i = 0 ; i < size ; i++){
for(int j = 0 ; j < size ; j++){
if( fabs(matA[i][j]) > 1.0e-12 ){
matAs.add(i, j, matA[i][j]);
}
}
}
//疎行列位置を確定
matAs1.fix();
```
なお、fixを呼ぶと、現在の要素位置に従ってEigen疎行列を作成します。
つまり、現在、長方形になっている場合はそのまま長方形で確定させます。
数値計算上、正方形にしたい場合は、fixに引数trueを追加してください。
```
matAs1.fix(true);
```
こうすることで、正方形になるようにゼロを代入してEigen疎行列を構築します。
### 疎行列の四則演算
fixで確定後は、疎行列の演算を実行可能です。
（なお、fix前の行列に演算しようと思うと未定義動作となりバグります。）
```
SRLfem::SparseMat matAs(5);
SRLfem::SparseMat matBs(5);
......
SRLfem::SparseMat matCs1 = matAs1 + matBs;
SRLfem::SparseMat matCs2 = matAs1 * matBs;
///...etc
//行列・ベクトル積も対応可能
double vec[5] ={1,2,3,4,5};
double* vec2 = matAs*vec;
vector vec2{1,2,3,4,5};
vector ans;
ans = matAs*vec2;
Eigen::VectorXd vec3
vec3 << 1,2,3,4,5;
Eigen::VectorXd ans2 = matAs*vec3;
```
ほかにも、転置、部分抽出、上三角、下三角取得、など、数値計算に必要そうな処理は大体実装済みです。
SparseMat.hppの、メソッド定義をご覧ください。
## 線形ソルバの利用
線形方程式Ax=bを解くソルバを複数実装しています。
特に電磁気のFEMで有効とされる加速係数付きICCG法を実装している、数少ないOSSと思います。
本OSS独自実装のソルバは下記のとおりです。
* 加速係数付きICCG法（実数・複素）
* 加速係数付きMRTR法（実数・複素）
* Eisenstatの方法を導入した前処理付きMRTR法（実数・複素）(SGS-MRTR)
なお、現状は、**大規模並列化には対応していません**。Eigenの並列計算機能を用いた小規模並列のみです。数コア程度なら並列化で多少早くなりますが、4,6コア以上を使ってもほとんど効果がないと思います。
デフォルトではEigenの並列化がonになっていると思いますので、Eigen機能を使って並列数を調整ください。
不要な場合はEigenの並列化機能をOFFすればソルバ内の並列化もOFFになります。
```
Eigen::setNbThreads(4);
int n = Eigen::nbThreads( );
cout << "Thred Eigen " << n << endl;
```
### ソルバの実行サンプル
独自実装のソルバはMatSolversクラスを使って実行できます。
収束判定を設定するために内部状態を持っているので、インスタンス化して利用してください。
以下は、ICCG法の実行例です。
```
SRLfem::MatSolvers solvers;
//@param1: 行数
//@param2: 収束判定
//@param3: 最大反復数
//@param4: 加速係数（1以上だとその値で固定。負の値にすると、その絶対値を超えない範囲で自動計算）
//@param5: 疎行列A
//@param6: 右辺ベクトルb
//@param7: 解ベクトルx
//@param8: ソルバ開始前の解ベクトルゼロ初期化有無。trueでxをゼロで初期化。falseだと何もしない(省略時はfalse)。
bool bl = solvers.solveICCG(size, 1.0e-8, 10000, 1.05, matAs, vec, results, true);
```
MatSolverクラスは、収束判定などにいくつかの内部変数を持っています。
代表的なものは以下。
* **setSaveLog(true)**: 残差のログを保存します。デフォルトはOFF。
* **setSaveBest(true)**: 最良解を保存するか。デフォルトはOFF。ONのとき、発散時には保存している最良解を解ベクトルxに代入します。
* **setConvNormalizeType(int x)**: 収束判定の正規化タイプ(0:右辺、1:初期残差、2:外部指定)
* **setDiagScale(true)**: 対角スケーリングの有無をセット。デフォルトはOFFです。
残差を記録している場合は、以下のように取得できます。
```
SRLfem::MatSolvers solvers;
vector log1;
solvers.getResidualLog(log1);
for(auto x : log1){
cout << x << endl;
}
```
ICCG以外のソルバの実行例です。
```
SRLfem::MatSolvers solvers;
// [IC-MRTR]
//@param1: 行数
//@param2: 収束判定
//@param3: 最大反復数
//@param4: 加速係数（1以上だとその値で固定。負の値にすると、その絶対値を超えない範囲で自動計算）
//@param5: 疎行列A
//@param6: 右辺ベクトルb
//@param7: 解ベクトルx
//@param8: ソルバ開始前の解ベクトルゼロ初期化有無。trueでxをゼロで初期化。falseだと何もしない(省略時はfalse)。
bool bl = solvers.solveICMRTR(size, 1.0e-8, 10000, 1.05, matAs, vec, results, true);
// [SGS-MRTR]: 加速係数不要。内部で対角スケーリング必須なので、フラグがOFFでも対角スケーリングは実行されます
//@param1: 行数
//@param2: 収束判定
//@param3: 最大反復数
//@param5: 疎行列A
//@param6: 右辺ベクトルb
//@param7: 解ベクトルx
//@param8: ソルバ開始前の解ベクトルゼロ初期化有無。trueでxをゼロで初期化。falseだと何もしない(省略時はfalse)。
bool bl = solvers.solveSGSMRTR(size, 1.0e-8, 10000, matAs, vec, results, true);
```