反復法の残差再計算の周期を !SOLVER から指定できるようにする

概要

CG法・BiCGSTAB法・GPBiCG法は、残差ベクトルを漸化式で更新し、一定の反復ごとに b - Ax から計算しなおしている。漸化式の残差が丸め誤差の蓄積で真の残差から乖離し、収束していない解を収束と判定するのを防ぐためである。この周期は解法ごとにソースコード内の定数で固 定されていた。CGは50、BiCGSTABは100、GPBiCGは20である。適切な値は行列の条件数と前処理の効き方に依存するため、利用者が調整できるようにした。

!SOLVER のヘッダ行に RECOMPUTE_RESIDUAL を追加した。正の整数はその周期、0または未指定は解法ごとの従来値、負の値は周期的な再計算を行わないことを表す。未指定時の反復履歴は変更前と完全に一致するため、既存の入力ファイルの挙動は変わらない。

PipeCG も同じ漸化式と固定周期50を持つため、あわせて対象にした。issue #61 (closed) が挙げる3解法には含まれないが、残差を再計算する解法で周期の出どころが分かれるのを避けたい。

従来の固定値が最適でない例は両方向にある。tests/solver/CG の入力を1スレッドの逐次実行で測定した結果である。

解法	収束判定値	未指定（従来値）	1	-1
BiCGSTAB	1.0e-11	521反復で破綻（NaN）	194反復で収束	441反復で収束
CG	1.0e-11	102反復で収束	10000反復で収束せず	104反復で収束
PipeCG	1.0e-8	91反復で収束	84反復で収束	86反復で発散

BiCGSTAB法は従来の周期100ではこの精度の解が得られないが、周期1なら収束する。CG法は逆に、周期1にすると収束しなくなる。PipeCG法は周期を長くする側で、再計算を行わないと不定値の判定にかかって発散する。周期は問題と収束判定値に応じて調整する価値がある。

無効化は周期へ MAXIT + 1 を代入することで表現した。反復ループはいずれも do iter = 1, MAXIT であり剰余が0になることがないため、再計算の要否を判定する2か所をそのまま残せる。判定式に周期が正であることのガードを置く案もあるが、Fortranは論理式の短絡評価を保証 しないため、0による剰余演算が実行される可能性が残る。

収束判定が成立したときの確認のための再計算は、本パラメータによらず常に行う。負の値を指定しても、収束と報告された解は真の残差で検証されている。

Close #61 (closed)

変更内容

コミット	要約
3acfdfa8 solver: make the residual recomputation period configurable	!SOLVER から周期を読み、hecMAT%Iarray 経由で3解法へ渡す
420e69ad tests: add a CG case with a non-default recomputation period	solver/CG に RECOMPUTE_RESIDUAL=10 のケースを追加
0e8abb0f solver: take the recomputation period from !SOLVER in PipeCG too	PipeCG の固定周期50を同じ経路に置き換える