反復ソルバー

Gmm++ で提供されているソルバーのほとんどは,ITLから少しだけ変更を加えたものです(gresが最適化され,複雑な行列に適用されています).ファイル gmm/gmm_iter_solvers.h をインクルードします.

イテレーション

Gmm++ の反復オブジェクトはITLのイテレーションオブジェクトを変更したものです.これはITLのテンプレート型ではありません.

最も簡単な初期化は次のとおりです.

gmm::iteration iter(2.0E-10);

ここで, 2.0E-10 は収束を得るために得られる(相対的)残差です.以下のようにいくつかの呼び出しが可能です.

iter.set_noisy(n) // n = 0 : no output
                  // n = 1 : output of iterations on the standard output
                  // n = 2 : output of iterations and sub-iterations
                  //         on the standard output
                  // ...
iter.get_iteration() // after a computation, gives the number of
                     // iterations made.
iter.converged()     // true if the method converged.
iter.set_maxiter(n)  // Set the maximum of iterations.
                     // A solver stops if the maximum of iteration is
                     // reached, iter.converged() is then false.

線形ソルバー

使用可能な線形ソルバのリストを次に示します.

gmm::row_matrix< std::vector<double> > A(10, 10);  // The matrix
std::vector<double> B(10); // Right hand side
std::vector<double> X(10); // Unknown
gmm::identity_matrix PS;   // Optional scalar product for cg
gmm::identity_matrix PR;   // Optional preconditioner
...
gmm::iteration iter(10E-9);// Iteration object with the max residu
size_t restart = 50;       // restart parameter for GMRES

gmm::cg(A, X, B, PS, PR, iter); // Conjugate gradient

gmm::bicgstab(A, X, B, PR, iter); // BICGSTAB BiConjugate Gradient Stabilized

gmm::gmres(A, X, B, PR, restart, iter) // GMRES generalized minimum residual

gmm::qmr(A, X, B, PR, iter) // Quasi-Minimal Residual method.

gmm::least_squares_cg(A, X, B, iter) // unpreconditionned least square CG.

ソルバ gmm::constrained_cg(A, C, X, B, PS, PR, iter); は線形拘束でシステムを解きます , C は拘束を表す行列です.しかし,まだ実験段階です.

(バージョン1.7)ソルバー gmm::bfgs(F, GRAD, X, restart, iter) はBFGSの準Newtonアルゴリズムであり,Wolfe式を用いて大規模な問題を探索します.制約なしに関数 F を最小化し,勾配 GRAD を与えます. restart は保存されている更新ベクトルの最大数です.

前処理行列

次の前処理行列は,線形ソルバーで使用できます.

gmm::identity_matrix P;   // No preconditioner

gmm::diagonal_precond<matrix_type> P(SM); // diagonal preconditioner

gmm::mr_approx_inverse_precond<matrix_type> P(SM, 10, 10E-17);
                                             // preconditioner based on MR
                                             // iterations

gmm::ildlt_precond<matrix_type> P(SM); // incomplete (level 0) ldlt
                                      // preconditioner. Fast to be
                                      // computed but less efficient than
                                      // gmm::ildltt_precond.

// incomplete ldlt with k fill-in and threshold preconditioner.
// Efficient but could be costly.
gmm::ildltt_precond<matrix_type> P(SM, k, threshold);

gmm::ilu_precond<matrix_type> P(SM);  // incomplete (level 0) ilu
                                      // preconditioner. Very fast to be
                                      // computed but less efficient than
                                      // gmm::ilut_precond.


// incomplete LU with k fill-in and threshold preconditioner.
// Efficient but could be costly.
gmm::ilut_precond<matrix_type> P(SM, k, threshold);

// incomplete LU with k fill-in, threshold and column pivoting preconditioner.
// Try it when ilut encounter too small pivots.
gmm::ilutp_precond<matrix_type> P(SM, k, threshold);

これらの前処理行列は, ildltt\_precond を除き,すべてITLから提供されています. ilut_precond は最適化されて単純化されています.また,安定性の理由から,不完全なLDLT 前処理行列で修正および変換されています(同様に, cholesky_precond を追加します.これは,実際には閾値前処理を持つ不完全なLDLT前処理行列です).もちろん, ildlt\_precondildltt_precond は対称的な実数あるいはHermitian複素行列を主にcgで使うために設計されています.

加算Schwarz法

加算Schwarz法は巨大線形システム(手法の原理については, [SCHADD] を参照してください)を解くことを可能にする分解領域法です.

現時点では,手法は並列化されていません(これは完了させる必要があります...).呼び出しは次のとおりです.

gmm::sequential_additive_schwarz(A, u, f, P, vB, iter, local_solver, global_solver)

A は線形システムの行列です. u は未知のベクトルである. f は右辺です. P はローカルソルバーの最終的な前処理行列です. vB は矩形疎行列(of type const std::vector<vBMatrix> ( vBMatrix は疎行列型))のベクトルであり,これらの行列の各々はサイズ \(N \times N_i\) であり,\(N\)A のサイズであり, \(N_i\)\(i^{th}\) サブドメイン内の変数の数です.行列の各列は, \(i^{th}\) サブドメインを表すサブ空間の基本ベクトルである. iter はイテレーションオブジェクトです.SuperLuがインストールされている場合は, local_solver はリスト gmm::using_gmres(), gmm::using_bicgstab(), gmm::using_cg(), gmm::using_qmr() から選択する必要があります. global_solver は,リスト gmm::using_gmres(), gmm::using_bicgstab(), gmm::using_cg(), gmm::using_qmr() から選択する必要があります.

テストプログラム schwarz_additive.C はGetFEMのディレクトリ tests にあり,より良い前処理(つまり,サブドメインの1つは,実際には粗いメッシュを表します)を作るために粗いメッシュの使用による弾性静的問題の付加的Schwarz法による求解の例です.

同一線形系で複数の解が存在する場合,前処理行列やLU分解を格納することで計算時間を短縮することができます.

gmm/gmm_domain_decomp.h の(とても)簡単なプログラムは,一定のオーバーラップ率で,規則的なドメイン分解を構築することができます.加算Schwarz法に対する行列 vB のベクトルを直接生成します.

範囲基底関数

gmm/gmm\_range\_basis.h で定義されている関数 gmm\_range\_basis(B, columns, EPS=1e-12) は,この行列の範囲の基礎となる疎行列 B の列から選択することができます.結果は columns に返されます.これは std::set<size_type> 型で,選択された列のインデックスを含みます.

このアルゴリズムは,線形依存列を持つ大きな行列から独立制約を選択するように特別に設計されています.

実装されたアルゴリズムには4つのステップがあります.

  • NULL列の削除.

  • すでに直交している列の集合の選択.

  • ブロック単位のGram-Schmidtアルゴリズムによる局所依存列の排除.

  • 大域的リスタート付きLanczosアルゴリズムによる残りの零空間のベクトルの計算と除去されるべきいくつかの列の推論.

このアルゴリズムは,局所 Gram-Schmidt アルゴリズムの後,それが低次元null空間のままであれば効率的です.実装されたリスタート付きLanczosアルゴリズムは,null空間ベクトルを1つずつ見つけます.

大域的リスタート付きLanczosアルゴリズムは,ブロックLanczos法( [ca-re-so1994] の例を参照),(並列化するための)ブロックWiedelann法,または単に各イテレーションにおけるnull空間の複数のベクトルの計算によって改良または置換することができます.