2009年6月7日日曜日

f77とf90の混在

さて、最近良くも悪くも忙しい。論文なぞ書いていて、ふと、「この数値計算ほんとに大丈夫なんだろうか」という不安が頭をよぎりだしたりして、D論を終わらせる前に計算の妥当性をちゃんと調べておこうという気になった。まーすでに何本かジャーナルに出てるんで妥当でないと困るんだが。

なんでそんなことが気になるかというと、NRGにおけるちょとしたテクがいまの系ではうまく働かないからだ。NRG(数値くりこみ群)は伝導電子バンドを対数的に離散化することで、問題を有限サイズ系の対角化にはめこんでいる。離散化の影響を最低限に抑えるために、いろんなメッシュとろうぜってのを試すとちょっとアレ?ってな感じだったりする。いまの自分のコードでは一般的な場合は影響がほとんどないんだが、伝導バンドとの混成にエネルギー依存をつけると、あやしい。一般的な場合に影響がでてないので、NRGの本体ではなくそこにほりこむパラメータ類の計算が怪しい。

パラメータの類はMathematicaで計算していた。Mathematicaは高性能かもしれんが、ブラックボックスなので不安である。なんでMathematicaかというと、計算に特殊関数のBessel関数とその数値積分がいるからだ。ニューメリカルレシピは単精度だったのでつかえねぇ…。というわけでMathematicaつかってたんだが、つい最近、高精度のBessel関数のフリーなライブラリを見つけた。
http://www.kurims.kyoto-u.ac.jp/~ooura/index-j.html
オオウラ先生ありがとぅぅぅ。Bessel関数だけじゃなく、高精度数値積分のパッケージまで作っておられる。まさにネ申。というかここのパッケージ使えば大半の物理計算ができてしまうと思われ。

恥ずかしながら、DE(二重指数関数型)公式という数値積分法は知らなかった。もう一個のIMT公式もだけど、日本人が発見した方法らしい。日本ってやっぱ数学つよいんだな。変数変換をもちいて区分求積することで猛烈に精度がよくなるらしい。なんでうちの専攻の授業ではガウス積分とかばっかでもっとつかえそうなこっち教えないんだ?

さて、ありがたくもらってきたんだが、ソースはf77で書かれている。私の普段使ってるのはf90。でも大丈夫。基本的にはおなじコンパイラ(g95,ifort)でコンパイルするので、-cオプションでオブジェクトファイルにしておいてリンクすれば.f90のファイルから呼び出すことができる。

ただし、問題が一つ。f90は結構型宣言について厳しい。私は、先頭にimplicit noneつけないと気持ち悪いんだが、これをするとfunctionを使うときにやっかい。たとえば、外部のファイルにfunction hogeがあるとして


program main

implicit none

integer :: i
double precision :: x,res
external :: hoge

res=hoge(x)

end program main

は、hogeが型宣言されてないとしてはじかれる。関数を double precision :: hoge みたいに型宣言してもあいてがf77なので仕様が違う→結構ヘンなことになりそうで怖い。

program main

implicit none

integer :: i
double precision :: x,res

interface
double precision function hoge(x)
double precision , intent (in) :: x
end function
end interface


res=hoge(x)

end program main

とinterface文でhogeの属性を指定してやるのがおそらくベスト。うぅ。いちいち使うたびにinterfaceかくのめんどくさいので、


subroutine ldbesj0(x,res)
implicit none

double precision :: x,res
interface
double precision function dbesj0(x)
double precision , intent (in) :: x
end function
end interface

res=dbesj0(x)

end subroutine ldbesj0

というのを各ベッセル関数について作った。こうすると、ライブラリにしたときに、ldbesj0とかを呼べばinterfaceについて考えなくてもよいしなぁ。というわけで、明日はライブラリ作ろうと思います。

0 コメント: