総称名ではヤバイ場合
https://www.nag-j.co.jp/fortran/FI_10.html
には「内部関数について総称名を使うことを強く推奨する」とある。ぼくの意見では
「安全であると知ってるものだけ総称名を使う。個別名で対応するほうが明瞭で良い」
が正しいです。型変換に関わる部分が危険です。
強く推奨する
というので、どういうことかと思うとあまり良くない例が載せてある(間違ってるともいえなくはない)。それは以下のコードです。危険です。
program generic_vs_specific
implicit none
double precision :: a = 2.0
complex(kind(0d0)) :: b = (2.0,1.0)
print *, sqrt(a)
print *, sqrt(b)
end program generic_vs_specific
complex(kind(0d0)) :: b = (2d0,1d0)としないとまずいです。complex(kind(0d0)) :: b = (2,1)でもいいですけどいつもd0つけると覚えとくのが良い。できるだけ脳の無駄遣いをしなくていいコードを書きましょう。上のコードでa=2.1,b=1.1として
print *, a
print *, b
を追加して実行するとわかるけど、a,bは単精度しか出ません。a=2.0,b=1.0だと一見精度が出ているわけでとても危険な例です。
上の例は比較的トリビアルだけど、もっと危険な例があります。以下のコードを実行してみてください。
complex(kind(0d0)) :: ccc = (5.1,3.2),c,cc,cc1,cc2,cc0,cc3
real(8)::a,b
write(6,*) 'ccc=',ccc
a=5.1d0
b=3.2d0
cc1=cmplx(a,b)
write(6,*) 'cc1=',cc1
cc2=dcmplx(a,b)
write(6,*) 'cc2=',cc2
cc3=complex(a,b) ! complexはgfortranでしか動かないー使うべきでない。
write(6,*) 'cc3=',cc3
end
倍精度が保てるのはcc2,cc3です。cc1も単精度ですー「総称名を使いましょう」の罠です。
real(a)と書くと当然単精度への変換となるわけだけど、cmplxだと複素化する機能だからね、と錯覚しちゃうわけです。
とにかく「cmplxは注意して使う」ということは記憶しておかないといけない。これはkindを与えない場合、精度を単精度に落としてしまう(real,imag,conjgも同様)。
ではどうする
ぼくの考えでは「dreal,dimag,dconjg,dcmplx」を使ったほうがいい。fortran2008ではA%REなども使えるらしいが使ったことはない。混合精度なコードだとkind入れて、cmplxを使うのは致し方ない。
ただconjgはkindを取れない.
gnufortranだとcomplexもあるが使おうとは思わない。
一般関数なら総称名でよい。
整数を倍精度に変換するときは、a=3d0で良い(simple is best).
整数化はnint,ceiling,floorを使う。
丸め誤差と離散化誤差:
誤差について考える必要もある。
数値計算は「丸め誤差と離散化誤差の狭間での戦い」なので状況によりヤバイことにもなりえます。普通に微分を取るとまあ3桁ぐらいは桁落ちするー2階微分なら6桁ぐらいって感じかもしれない。そうすると単精度だと精度が2桁しか残らないということになる。あるいは100万回ぐらい積算するとその最悪の値が悪さをして全体の精度に響くということもあり得る、負符号も取りうる量を足すときには悩まされます。正定値の計算では平均化でいいように働くかもしれない。
系統誤差と偶然誤差:
系統誤差(打ち切りパラメータ、離散化パラメータ)はパラメータを変えて観察すれば良さそうなものだが、実際はなかなかに制御しずらく測りきれない。入力に対する出力のスムーズネスをみれば偶然誤差(ランダム誤差)は計測できるはずー倍精度だと偶然誤差などないかもしれないが桁落ちすることもあるのでなんとも言えない。収束計算においては収束安定性において問題だったりもする。実際のところ、なかなかに制御しづらい。
台形公式とシンプソン公式:
シンプソン公式は2次関数によるフィッティングに基づいている。差分も見ていることになり、倍精度は必須だと思える。丸め誤差がそれなりに積算してる状況、とか被積分関数のスムーズネスが期待できない場合なら台形公式のほうが計算的な安定性がよい。