固有値解析の収束判定で要求された低次モードを見るようにする (反復回数が3〜12倍に膨らみ、maxiter 到達時に擬似ゼロ固有値が出る)

Close #786 (closed) 詳しくは下記のとおりですが、固有値解析についての2017年頃の作業で移植ミスがあったのが原因で、もともと（v4系では）正しく計算できていたようです。 テストでは解はそのままに、反復回数の増加は解消しています。

問題

!SOLUTION, TYPE=EIGEN のランチョス反復において、収束判定が「利用者が要求した低次モード」ではなく、QL法が返した未ソートの Ritz 値の先頭 nget 個という無関係な部分集合の残差を見ている。

シフト・インバート形 A = (K - sigma*M)^-1 * M では、最低次の固有値 lambda が最大の Ritz 値 theta に対応する（lambda = 1/theta - sigma）。しかし fstr_EIG_tridiag.f90 の判定は

do i = 1, min(nget, iter)
  resid = dabs(Tri%beta(iter+1)*L(iter,i))/dabs(alpha(i))

と配列先頭から機械的に nget 個を取っており、この alpha はソートされていない。

その結果:

• 求めたいモードは8〜12反復で収束しているのに、判定が通らず33〜101反復まで回り続ける（3.5〜12.6倍の過剰反復）
• MPI領域分割数を変えると反復回数が 101 → 357 → 195 と大きく振れる
• maxiter に到達すると無警告で「収束」扱いになり、擬似ゼロ固有値 0.0000E+00 や刺激係数 NaN が出力される

報告は積層シェル要素についてだが、要素種別には依存しない。ソリッド要素（341/361）でも3.5〜4.9倍の過剰反復が出る。積層シェルで顕在化しやすいのは、集中質量の対角に並進成分と回転成分で数桁の開きがあり（static_LIB_shell.f90 の mass_Shell）、A のスペクトルが広く密になって theta の小さい側に密集帯ができるためと考えられる。

原因

QL_decomposition の並べ替えが移植ミスで機能していない。移植元 EISPACK tql2.f の

if (d(j) .ge. p) go to 260

のラベル260は内側 do の終端 continue であり cycle 相当だが、固定形式から自由形式への変換（3c7a60ef, 2017-08-12）の際に exit aa と読み替えられた。このため選択ソートが最初の非減少要素で打ち切られる。乱数3重対角200件で検証したところ、199/200で昇順にならず、最大 theta は200/200で末尾に来る。すなわち alpha(1:nget) に最重要モードは一度も入らない。

!GP: Get unordered eigenvalues and eigenvectors というコメントは、2011年版（8627228c）で du / zu に未ソートのコピーを取っていた処理に付いていたもので、b4ee8093 で du / zu が削除された際にソート処理の直上に取り残されたものである。同じ関数のヘッダコメントが現在も "contains the eigenvalues in ascending order" と書いていることが元の意図を示している。

ソートを直すだけでは解決しない

QL_decomposition を EISPACK どおり昇順に直しても解決しない。昇順に並ぶのは theta であり、求めたいモードは theta が最大の側＝配列の末尾に来る。do i = 1, nget のままでは、今度は確実に「theta 最小＝最高周波数側」を判定することになり、現状より悪化する。

求めたいモードの定義は |1/theta - sigma| の小ささであって theta の大小ではないため、theta をどうソートしてもモード選択の問題は解けない。判定側でモードを選択する必要がある。

なお2011年版（V4系）は fstr_solve_eigen.f90 で EVAL(i) = 1/LLDIAG(i) + lczsgm を全モードについて算出し、CALL EVSORT(EVAL,NEW,LTRIAL) で並べ替えたうえで JITER = NEW(JJITER) により KITER = NGET+2 個（!Extra 2 values as a safety feature）の残差を見ていた。2018年の書き換え（37ce461a）でこの選択機能が失われている。本MRはこれを回復するものである。

修正内容

fistr1/src/analysis/dynamic/mode/fstr_EIG_tridiag.f90 の1ファイルのみ（45行追加・63行削除）、4コミット。

コミット	内容
Remove the unused variable max_eigval ...	#775 の修正（a55f69a0）以降まったく参照されていないデッドコードの削除
Check convergence of the requested modes ...	本MRの中核。eigval の算出と evsort を判定より前に移し、in = iparm(i) で低次側 nget+2 個を選ぶ
Drop the eigenvalue ordering from QL_decomposition	判定・出力とも順序に依存しなくなったため、壊れた選択ソートを削除。ヘッダコメントの矛盾も解消
Report non-convergence at maxiter ...	maxiter 到達時に警告を出し、得られた Ritz 対の数だけを出力する

主要な差分:

+    sigma = 0.0d0
+    if(fstrEIG%is_free) sigma = fstrEIG%sigma
+    do i = 1, iter
+      if(alpha(i) /= 0.0d0)then
+        eigval(i) = 1.0d0/alpha(i) - sigma
+      else
+        eigval(i) = huge(0.0d0)
+      endif
+    enddo
+
+    call evsort(eigval, iparm, iter)
+
     is_converge = .true.
-    do i = 1, min(nget, iter)
-      if (dabs(alpha(i)) > 0.0d0) then
-        resid = dabs(Tri%beta(iter+1)*L(iter,i))/dabs(alpha(i))
+    do i = 1, min(nget+2, iter)
+      in = iparm(i)
+      if (dabs(alpha(in)) > 0.0d0) then
+        resid = dabs(Tri%beta(iter+1)*L(iter,in))/dabs(alpha(in))

補足:

• evsort は判定前に1回だけ呼び、収束時の固有ベクトル組み立てでは同じ iparm を再利用する。二度呼ぶと L(j,in) が必要とする QL 側の元添字が失われる
• 判定対象を nget+2 個としたのはV4系に倣ったもので、要求集合の境界にあるモードが並べ替えで判定から外れるのを防ぐ
• alpha(i) == 0 のモードには huge(0.0d0) を与えてソート末尾へ送る。判定前に全モードの eigval を作る構造にすると、代入されずに0のまま残った eigval が先頭に来て求めたいモードを押しのけるため。あわせて、メッセージなしにMPI全体を落としていた hecmw_abort の呼び出しを削除し、未収束として反復を続ける扱いに変えた

効果

K731ORTHOISO（積層シェル731、nget=5、tolerance=1.0e-6）:

	修正前	修正後
反復回数（逐次）	101	16
反復回数（np=1/2/4）	101 / 357 / 195	16 / 16 / 16

固有値は修正前後、および np=1/2/4 で一致する。

判定残差の推移（修正前は96反復目で3桁悪化しており、判定対象が反復ごとに入れ替わっていたことが分かる）:

iter    修正前        修正後
   7  1.14086E+01  1.14086E+01    <- ここまでは全モードを判定するので一致
   8  1.50886E+01  4.12996E-01
   9  2.33247E+01  4.85579E-01
  ...
  13       -       7.34617E-05
  16       -       7.28933E-08    <- 収束
  ...
  95  6.61416E-05       -
  96  9.42594E-02       -         <- 3桁悪化
 101  4.07641E-07       -         <- 収束

8ケースの比較（「必要反復」は厳密解と4桁一致する最小の maxiter）:

ケース	要素	nget/tol	必要反復	修正前	修正後
K731FIBERM45	積層シェル731	5 / 1e-6	8	33	16
K731FIBERP45	積層シェル731	5 / 1e-6	9	79	16
K731ORTHOISO	積層シェル731	5 / 1e-6	8	101	16
K741FIBER180	積層シェル741	5 / 1e-6	9	56	19
K741LAMINATE3	積層シェル741	5 / 1e-6	12	101	16
exK/K741	シェル741 単層	3 / 1e-8	8	30	12
exK/K341	ソリッド341	3 / 1e-8	7	34	12
exK/K361	ソリッド361	3 / 1e-8	8	28	16

全8ケースで固有値は修正前と5桁一致する。過剰反復の比は3.5〜12.6倍から1.3〜2.0倍に下がった。

maxiter 到達時（nget=5、maxiter=3）:

	修正前	修正後
出力モード数	5	3
モード4	0.0000E+00、振動数0、刺激係数 NaN	出力されない
モード5	2.3715-322（デノーマル）	出力されない
警告	無し	標準出力と 0.log に ### WARNING: eigen analysis stopped at maxiter without convergence.

影響範囲

QL_decomposition の呼び出しは tridiag の1箇所のみであり、他の解析種別からは呼ばれていないため、影響は固有値解析に閉じる。

入出力の仕様に変更は無く、frontistr_documents への反映は不要と判断した。ただし以下は利用者から見える挙動変化である。

• maxiter に到達した場合、未収束の警告が出るようになる（従来は無警告で正常終了扱い）
• 未収束のまま maxiter に到達し、かつ得られた Ritz 対が nget に満たない場合、出力されるモード数が nget ではなく実際に得られた数になる