線形代数ソルバ

構造解析の連立一次方程式は、f64 の疎行列に対する直接法（Cholesky LLᵀ 分解、LU 分解）で解きます。通常の静解析では拘束縮約後の剛性行列が対称正定値であることを前提とし、自由度規模によらず疎 Cholesky を用います。特異・非正定値の行列、または倍精度で信頼できる正定値求解を確認できない行列では解析を停止します。荷重がゼロ、または不安定な変位方向を励起しない荷重でも分解を省略しません。

疎行列の扱い（スカイライン法との対比）

構造解析で古くから使われる疎行列の格納・求解法に、スカイライン法（バンド法・エンベロープ法）があります。スカイライン法は各列の対角から最上端の非ゼロまでの「帯」を連続配列で保持する方式で、帯の内側にあるゼロ要素も格納・演算の対象になるため、効率が節点番号付け（Cuthill–McKee 法などによるバンド幅最小化）に強く依存します。

本プログラムはスカイライン法ではなく、汎用疎行列格納と現代的な疎行列直接法を用います:

1. 非ゼロ要素のみの格納：全体剛性行列は要素剛性の Triplet（行・列・値）として組み立て、CSC（圧縮列格納）形式へ変換します。同一位置の重複は加算マージし、帯の内側か否かに関係なく、非ゼロ要素だけを保持します。
2. フィルイン低減オーダリング：疎 Cholesky（faer）は数値分解に先立つ記号解析（symbolic analysis）で、フィルイン（分解で新たに生じる非ゼロ）を減らす行列の並べ替えを自動で適用し、因子 \( L \) の非ゼロパターンを事前に確定します。節点番号付けの良し悪しに性能が左右されるスカイライン法と異なり、入力モデルの節点順序に依存せずフィルインが抑えられます。
3. 記号解析と数値分解の分離：非ゼロパターンの解析（記号解析）と数値分解を分離しているため、一度 factorize すれば複数右辺（多数の荷重ケースや時刻歴の各ステップ）は前進・後退代入のみで解けます。

計算量の目安として、スカイライン法の分解コストは自由度数 \( n \)・平均バンド幅 \( b \) に対しおおむね \( O(n b^2) \) で帯内のゼロにも演算が及ぶのに対し、オーダリング後の疎 Cholesky は因子の非ゼロ数 \( \mathrm{nnz}(L) \) に応じたコストで済みます。立体骨組のように帯内に多くのゼロを含む系では、この差が求解時間・メモリの双方に効きます。

算定式

対称正定値の \( K \cdot u = F \) には疎 Cholesky を用います。通常の静解析で Cholesky 分解に失敗しても LU への自動切替は行いません。

正の対角成分から対角行列 \( D=\operatorname{diag}(\sqrt{K_{ii}}) \) を作り、合同平衡化した行列

\[ A=D^{-1}KD^{-1},\qquad b=D^{-1}F \]

を用います。並進と回転の単位や対角剛性の大きさによる見かけの差を取り除きます。対角が非正または非有限、平衡化した成分が非有限の場合は停止します。フィルイン低減の並べ替え \( P \) の後、

\[ PAP^\mathsf{T}=LL^\mathsf{T} \]

へ分解し、前進・後退代入で \( Az=b \) を解き、\( u=D^{-1}z \) として元の変位へ戻します。剛性や拘束を付加する操作ではありません。疎 LU は不定値・非対称行列を扱う別の直接法であり、この正定値診断とは区別します。

数値条件の確認と逆1ノルム推定

浮動小数点の丸めにより、不安定な架構の零剛性が微小な正値になって分解が成功する場合があります。そのため、荷重に対する解を公開する前に、平衡化後の行列について推定逆条件数を確認します。対象は全体行列を拘束縮約した後の行列です。

行列の1ノルムは列ごとの絶対値和の最大値です。

\[ \lVert A\rVert_1=\max_j\sum_i|A_{ij}|,\qquad \widehat{r}=\frac{1}{\lVert A\rVert_1,\widehat{\lVert A^{-1}\rVert_1}} \]

逆行列の1ノルムは、逆行列を作らずに標準的な反復推定法で評価します。一定値の初期ベクトル、符号ベクトル、選ばれた座標方向への求解を繰り返し、符号の反復・推定値の停滞・反復回数で終了し、最後に交互符号のベクトルも確認します。対称行列では転置方向への作用も同じ因子で求解できます。これは標準的な逆1ノルム推定に用いられる DLACN2 の方法です。既存の因子による数回の前進・後退代入とベクトル演算を追加し、逆行列の生成や追加の因数分解は行いません。

\( \widehat{r}<\epsilon \) の場合、または推定途中に非有限値が生じた場合は、数値精度上信頼できる正定値求解を確認できないものとして停止します。ここで \( \epsilon=2^{-53}\simeq1.1102230246251565\times10^{-16} \) は倍精度の丸め単位です。条件数推定は DPOCON、判定水準は DPOSVX の作業精度で特異とみなす判定を根拠とします。DPOSVX は警告後も求解を続けますが、本プログラムでは結果を公開せず停止します。

この推定法が与える逆1ノルムは下界であり、推定逆条件数はすべての不安定な架構を検出する数学的保証ではありません。また、数学的に正定値でも作業精度で識別できない行列は停止対象です。診断は物理的な不安定性の証明でも、建物の安全性の判定でもありません。診断通過は要求する変位精度を保証せず、各求解で元の行列に対する残差を常時検査する仕様でもありません。

数値処理の一次資料：LAPACK 3.12.1 DPOCON、DLACN2、DPOSVX、DLAMCH の目的・計算手順・機械精度の定義によります。

記号・単位

分解の再利用と適用条件

• 線形静解析：通常・長期の固定剛性は採用条件ごとに分解し、同じ条件の荷重ケースに再利用します。引張専用ブレースはケースごとの有効集合を反復更新し、剛性が変わるたびに数値分解します。
• 固有値解析・参考固有周期：静解析の準備で分解した基準剛性を全規模で再利用します。長期の軸剛性低減や引張専用のケース別最終剛性には置き換えません（固有値解析）。
• 線形時刻歴応答（Newmark-β）：全ステップ共通の有効剛性を1回分解し、各ステップの求解に再利用します。
• 増分解析・非線形時刻歴：接線剛性が変わるたびに数値分解し、同じ疎パターンの記号解析を再利用します。

上記の正定値診断は、これらの Cholesky 分解で実際に求解する行列へ共通に適用します。時刻歴では有効剛性や初期加速度の質量行列が対象であり、静的な不安定性だけから動的な解の存在を一律に否定するものではありません。初期加速度に用いる質量行列の既存の LU への切替は、この診断とは別の契約です。

行列が更新された際の再分解に失敗した場合は、Cholesky・LU ともに以前の行列の分解結果を使って求解しません。求解を再開するには、対象行列の分解が正常に完了する必要があります。

線形連立方程式を直接法で解く方針と、固有値反復・Newton 反復・引張専用の有効集合反復の収束条件は別です。反復が規定回数内に収束しない場合はエラーで停止します。因子のフィルインによって分解の時間とメモリが増えるため、大規模解析の資源量は構造・拘束・疎パターンに依存します。