ペナルティ法の固有値解析で使えない線形解が黙って使われる問題を修正

概要

ペナルティ法（MPCMETHOD=1）で固有値解析を行うと、固有値が 0 や 1.8E+308 のような値になる。原因はペナルティ法そのものではなく、線形ソルバーが要求された許容残差に到達できないまま MAXIT まで反復し、その間に発散した解を Lanczos ループが黙って使うことである。

ペナルティ法は剛性行列に「対角最大値 × penalty」（!MPC の既定値 1.0E+4）を加えるため条件数が悪化し、tests/analysis/eigen/mpc_elim の lever が指定している許容残差 1.0E-12 に届かなくなる。hecmw_solve_CG は許容残差を下回るたびに真の残差を計算し直して判定するので反復が止まらず、再帰更新の残差は 1E-11 から 1E+19 まで発散する。hecmw_solve_iterative はこれを警告 HEC-MW-SOLVER-W-3001 として出すだけである。到達不能な許容残差を指定すれば消去法でも同じことが起きるので、ペナルティ法に固有の問題ではない。

許容残差を到達可能な値にすれば、ペナルティ法でも消去法と一致する固有値が得られる。

モデル	条件	固有値
lever	MPCMETHOD=1, resid=1.0E-12	0.0000E+00 / 7.0854E-20 / 1.7977E+308
lever	MPCMETHOD=1, resid=1.0E-8	3.8208E-02 / 4.0838E-01 / 5.9792E-01
lever	MPCMETHOD=3, resid=1.0E-12	3.8208E-02 / 4.0838E-01 / 5.9792E-01
tie	MPCMETHOD=1, resid=1.0E-8	2.4623E-02 / 2.1799E-01 / 5.8578E-01
tie	MPCMETHOD=3, resid=1.0E-12	2.4623E-02 / 2.1799E-01 / 5.8579E-01

したがってペナルティ法の固有値解析は動作する。対応したのは、使えない線形解が黙って使われることと、この組み合わせのテストが無いことである。

未収束または発散のフラグが立ったときにだけ相対残差を計算し直し、1.0E-2 以上であれば停止する。未収束そのものを致命的にしていないのは、既存の固有値テストのいくつかが、要求した許容残差を満たさないまま正しい固有値を出しているためである。閾値は実測から決めた。判定に掛かったときの相対残差を一時的に出力して固有値解析を含むテスト44件を実行したところ、20回の該当があり、値は 2.47505E-07 から 1.92107E-05 の範囲だった。1.0E-2 はその最大値に対して520倍の余裕があり、かつ解が1%程度の誤差を含む水準なので、そこから計算される固有値はモード解析に使えない。

閾値は要求許容残差に依存しない絶対値1つとした。まずは単純な判定で明らかにおかしな結果を出さないようにすることにとどめ、それ以外は警告を出して計算を継続する方針である。要求許容残差との比を併用する案も検討したが、停止させる範囲を広げることになるため採らなかった。要求を満たしていないこと自体は下記の警告で残差の値とともに伝わる。!EIGEN の収束判定値と比べる案は、上記のテストがその判定値より大きい線形残差のまま正しい固有値を出しているため成立しない。

上限に達しない場合は停止せず、到達した相対残差を警告として出す。既存の HEC-MW-SOLVER-W-3001 は反復回数内に収束しなかったとしか言わず、どれだけ外れているかを示さないためである。この警告が出る状況は HEC-MW-SOLVER-W-3001 が出る状況と一致しており、警告の総量は増えない。

線形ソルバー側の改修（到達限界で停止して最良の解を返す）は行っていない。全解析が使う中核ソルバーに手を入れることになるため、別issueとするのが適当である。

Close #865 (closed)

変更内容

コミット	要約
e08739db eigen: report a linear solve that fell short of its tolerance, and stop on an unusable one	線形解の使用可否を判定し、使えない場合は停止、使える場合は残差を警告するようにした
8bc89fcd tests: add eigen analysis with !EQUATION handled by the penalty method	MPCMETHOD=1 の固有値解析の回帰テストを追加