5.1 線形代数ソルバ
構造解析の基礎となる連立一次方程式を、数値線形代数の標準手法(Cholesky LLᵀ 分解、LU 分解、前処理付き共役勾配法 PCG)で解きます。対称正定値系のソルバ選択は既定で自動(AUTO)とし、自由度規模に応じて直接法(疎 Cholesky)と反復法(PCG)を使い分けます。PCG が収束しない場合は直接法へ自動でフォールバックするため、最終的に解が得られないことはありません。
疎行列の扱い(スカイライン法との対比)
構造解析で古くから使われる疎行列の格納・求解法に、スカイライン法(バンド法・エンベロープ法)があります。スカイライン法は各列の対角から最上端の非ゼロまでの「帯」を連続配列で保持する方式で、帯の内側にあるゼロ要素も格納・演算の対象になるため、効率が節点番号付け(Cuthill–McKee 法などによるバンド幅最小化)に強く依存します。
本プログラムはスカイライン法ではなく、汎用疎行列格納と現代的な疎行列直接法を用います:
- 非ゼロ要素のみの格納:全体剛性行列は要素剛性の Triplet(行・列・値)として組み立て、CSC(圧縮列格納)形式へ変換します。同一位置の重複は加算マージし、帯の内側か否かに関係なく、非ゼロ要素だけを保持します。
- フィルイン低減オーダリング:疎 Cholesky(faer)は数値分解に先立つ記号解析(symbolic analysis)で、フィルイン(分解で新たに生じる非ゼロ)を減らす行列の並べ替えを自動で適用し、因子 \( L \) の非ゼロパターンを事前に確定します。節点番号付けの良し悪しに性能が左右されるスカイライン法と異なり、入力モデルの節点順序に依存せずフィルインが抑えられます。
- 記号解析と数値分解の分離:非ゼロパターンの解析(記号解析)と数値分解を分離しているため、一度
factorizeすれば複数右辺(多数の荷重ケースや時刻歴の各ステップ)は前進・後退代入のみで解けます。 - 反復法という選択肢:PCG は行列を CSR 形式(非ゼロのみ)で保持し、因子分解を一切行わず、1 反復あたり非ゼロ数に比例するコストの疎行列ベクトル積だけで解きます。フィルインが原理的に発生しないため、直接法の分解メモリが支配的になる大規模系で有利になります。AUTO 選択(後述)では、大規模系でこの PCG を自動的に用います。
計算量の目安として、スカイライン法の分解コストは自由度数 \( n \)・平均バンド幅 \( b \) に対しおおむね \( O(n b^2) \) で帯内のゼロにも演算が及ぶのに対し、オーダリング後の疎 Cholesky は因子の非ゼロ数 \( \mathrm{nnz}(L) \) に応じたコストで済みます。立体骨組のように帯内に多くのゼロを含む系では、この差が求解時間・メモリの双方に効きます。
算定式
対称正定値の \( K \cdot u = F \) には疎 Cholesky を用い、非対称やラグランジュ乗数付き拘束では疎 LU にフォールバックします。
直接法(疎 Cholesky)では、フィルイン低減の並べ替え \( P \) を施した剛性行列を下三角行列 \( L \) へ分解し:
\[ P K P^\mathsf{T} = L L^\mathsf{T} \]
前進代入 \( L y = P F \)、後退代入 \( L^\mathsf{T} \tilde{u} = y \) の 2 段の代入計算で解 \( u = P^\mathsf{T} \tilde{u} \) を得ます。分解は 1 回で済み、右辺(荷重ケース)ごとの求解は代入計算のみを繰り返します。疎 LU も同様に \( K = L U \)(行・列の置換付き)へ分解し、前進・後退代入で解きます。
Jacobi 前処理付き PCG(内部 f32 演算)では対角前処理を次とします:
\[ d[i] = 1/a_{ii} \]
収束は \( \|r\|/\|b\| < \text{tol} \) で判定します。 反復中の内積が underflow 域(\( 10^{-30} \) 未満)に落ちた場合は停滞とみなして打ち切ります。 以降の除算が NaN・Inf になり反復が破綻するためです。 対角成分が同域の場合は Jacobi 前処理の当該成分を 1(前処理なし)とします。 右辺ノルムが同域の場合はゼロ解を返します。
直接法/反復法の自動選択(AUTO)
対称正定値系の線形静解析は、既定で AUTO ソルバを用います。しきい値は自由度数 50,000 で、1 節点あたり 6 自由度(並進 3 + 回転 3)なので、節点数に換算するとおよそ 8,300 節点(支持条件による拘束で有効自由度は減るため、実際はもう少し大きいモデル)に相当します。選択規則は次のとおりです:
- 自由度数 \( n < 50{,}000 \)(目安: 約 8,300 節点未満):疎 Cholesky 直接法(f64 厳密解)。中小規模では分解コストが十分小さく、精度面でも直接法が有利なため。
- 自由度数 \( n \geq 50{,}000 \)(目安: 約 8,300 節点以上):Jacobi 前処理付き PCG(tol \( 10^{-6} \)、最大 10,000 反復)を先に試みる。収束しない場合は疎 Cholesky へ自動フォールバックし、フォールバック分解は以降の右辺で再利用する。
ただし、分解の再利用が効く解析は直接法(疎 Cholesky)に固定します:
- 線形時刻歴応答(Newmark-β): 有効剛性が全ステップ共通で、1 回の分解を全ステップの求解で再利用するため。
- 固有値解析(部分空間反復): 同一分解を(部分空間サイズ×反復回数)回の求解で再利用するため(求解のたびに反復をやり直す PCG では分解再利用の利得が失われる。5.3 固有値解析 参照)。
- 増分解析(プッシュオーバー)の Newton 反復・非線形時刻歴の有効剛性: 反復のたびに数値は変わるものの剛性行列の非ゼロパターンは同一のため、直接法の記号解析(フィルイン低減の並べ替えと因子パターンの確定)を初回の 1 回だけ行い、以降の反復は数値分解のみを繰り返します。反復回数の多い非線形解析では、この記号解析の再利用が支配的な効率化になります。
また、非対称やラグランジュ乗数付き拘束の系は AUTO の対象外で、疎 LU を用います。
実装:squid_n_math(sparse.rs(Triplet→CSC 組立)/solver.rs/cholesky.rs/lu.rs/pcg.rs/auto.rs(AUTO 選択)、faer ライブラリ利用)が算定します。