2008年2月17日日曜日

fortran90関連:変数・関数の精度など

しょっちゅう忘れるわりにあんまりまとめが見つからないfortran90の精度関連のメモ。数値計算するときゃ精度が命だ。とはいえ自分以外に参照する人間がいるとは思えない。

・変数の精度について
倍精度は十進で15桁なので大抵の数値計算はこれで良いのだが、倍精度の変数を宣言するときに注意しないといけないことがある。変数の精度は一般的に種別パラメータにより決定される。普通、倍精度の種別パラメーターは8であるので、
倍精度実数 real(8)
倍精度複素数 complex(8)
で大抵は何とかなるが、コンパイラがフツーのじゃ無い場合、種別パラメーター8が違うものをさす場合がある。一番シンプルな解決方法は組み込み関数kindを使って、
倍精度実数 real(kind(0d0))
倍精度複素数 complex(kind(0d0))
てな風にすることだが、見た目があんまり麗しくない&計算の全体の精度を変更したときに、大量に修正する箇所が出てくるという嬉しくない事態が起こりがち。そこでselected_real_kind関数を使うと言う方法がある。selected_real_kind(p,r) のpは有効桁数(10進数)であり、rは指数範囲である。
この組み込み関数の戻り値は整数で、これが、p,rで決められた必要な桁数をみたす種別パラメーターを返している。p,rは省略可能。
module precision_dp
integer,parameter ::dp=selected_real_kind(12)
end module precision_dp
みたいなmoduleを作っておいて、手続き側でuse文でこのmoduleを呼び出して
real(dp) :: xみたいに変数を宣言しておけば、精度の変更がらくちん。

指数形式は
単精度なら 1.345E-03
倍精度なら 0.24577650003D-03
のような表現で可能。

変数の桁の保障は(要はこっちが決めてないところはその桁までちゃんと0を入れておくっていう指示になるぽい)
1.234560000_10
のように最後に_と桁数を書くことで可能になる。桁数は上で出てきたselected_real_kind()の戻り値を使って表現することも可能。しかし、指数形式の倍精度版とこの表現の併用(つまり1.3479847D-04_dpとか0.006D-02_8みたいな)はうまくいかなかった。指数形式をEで書くと、そういう表記も使えるし、桁数を大きくとっておくと値はその桁にあわせた精度で返ってくる。よくわからんなぁ。指数形式をDで書くとそっちのほうで倍精度てのが決定されるからだろうか?

・関数の精度について
fortranの算術関数exp,cos,sinなどは基本的に引数と同じ精度の戻り値を返す。倍精度専用のdexpとかもあるが、非標準になってるコンパイラとかもあるので総称関数exp等を使うほうが無難だ。

気をつけないといけないのが、変数の型の変更とか複素数がらみの組み込み関数である。実数型に変換するreal,複素数の虚部をとるaimag、複素数に代入するcmplxあたりが代表格。aimagは総称関数であるので、とくに注意せずこれを使ってれば問題は起こりにくいが、realとcmplxがやばい。この2つはデフォルトで戻り値が単精度なので、倍精度変数にこれを使ってしまうと、精度ががた落ちという悲しい事態が起こりうる。realの倍精度版dbleは標準になっているので、これを使えばよい。しかし、悲しいかなdcmplxは非標準。じゃあどうするかというと、実はcmplxの引数は第3まであって、第3の引数によって戻り値の精度を制御できる。用はcmplx(1.0d0,2.0d0,kind(0d0))みたくすれば倍精度の戻り値を得られる。

0 コメント: