Use BLAS-like subroutines on iterative methods
Closes #702 (closed)
-
hecmw_axpby_R(Y = a * X + b * Y) の作成(MKLに同様のサブルーチンが存在する) - BiCGSTAB, GPBiCG, GMRES, GMRESR, GMRESRENの各サブルーチンにおいて、ベクトル演算のdo loopを、全てBLAS様サブルーチンを用いるように変更した
- axpby, scaleについては、係数が0かの場合分け(※)を行った(
contact_slag_iterなどのテストでの失敗を修正)
※場合分けの理由(GitHub Copilot出力)
原因
hecmw_axpby_R と hecmw_scale_R に IEEE 754 の 0 × NaN = NaN に起因するバグがあります。
hecmw_axpby_R(n, coef, 0.d0, X, Y)(GMRES/GMRESR/GMRESRENで計8箇所使用)
- 元のインラインコード: Y(i) = X(i) * coef → 旧Yを読まない
- 変更後: Y(i) = coef * X(i) + 0 * Y(i) → 旧Yを読む → Yが未初期化/NaN の場合 0 * NaN = NaN で汚染
hecmw_scale_R(n, 0.d0, X)(GMRESで4箇所使用)
- 元のインラインコード: X(i) = 0.d0 → 旧Xを読まない
- 変更後: X(i) = 0 * X(i) → 旧Xを読む → Xが NaN の場合 0 * NaN = NaN で汚染
Fortran の allocate はメモリを初期化しないため、WW 配列の未使用領域に NaN ビットパターンが含まれる可能性があります。環境依存で発生する理由は、メモリアロケータの動作が環境(OS、コンパイラ、最適化レベル)によって異なるためです。
Edited by Kanato Kono