Orz日記 by Akio Morita

_ [SAD]INTRA EMITのバグ探し

誤差付きOpticsをINTRA EMITすると BINTが NaNを出して まともに動かない件を昨日から真面目に調べていて分かったこと

• bint()に渡される被非積分関数が、積分領域の端でNaNを返す
  • 改良前の台形公式による積分値がNaNなので、bint()内の収束判定が必ず失敗する
  • 被積分関数のパラメータに問題が有り、関数の値が純虚数になる領域がある

というわけで、呼出側(tintrb)での前後関係を洗ってゆくと、NaNが発生する 前の計算で、3x3行列の固有値がおかしいことが発覚。 周辺のコードから、固有値は正定値であることが期待されているのに、 零や負の値が帰ってきている。

固有値分解関数eigs33()へ与えている行列と帰ってきた固有値を 正常時と異常時で比較してゆくと、一番小さい固有値が失われているのが 直接の原因と判明。

内部の2x2部分行列の固有値分解で桁落ちしている部分を修正したら、 多少マシになったが、完治していない。

多分、対角化操作で行列を回転する際に桁落ち or 情報落ちしているものと 思われる。対処法としては、

1. ピポッド選択等のアルゴリズムを改良する
2. 4倍 or 8倍精度実数を使う
3. Error Free Transformationで行列操作を実装する

3x3行列なので、(1)はあまり余地は無いかもしれない。 (2)は、Fortran標準では実装必須にはなっていないので、 動作がコンパイラ依存になるのが難。 (3)は、実装が面倒くさいのが難。C++とかだと class化すれば 表面上のコードは綺麗にかけるけど、FortranからC++を呼び出すのは 木に竹を継ぐような話だから避けたい所である。