球面調和のライブラリって何が一般的に何が使われるのだろうか。
Fortranではないけれど,データ分析用言語Interactive Data Language(IDL)に球面調和関数が実装されているらしい。が,波数lが大きい所(l>150)での挙動がおかしいらしい。
http://ryukyu.astr.tohoku.ac.jp/pukiwiki/index.php?Members%2Fchinone%2F%B3%D0%BD%F1%2FIDL#s2fdd458
科学技術計算関数のライブラリにGNU Scientific Library(GSL)がある。Fortran用インターフェースとしてFGSLもある。
インストールはsudo port install fgslで行けた。
球面調和関数はgsl_sf_legendre_sphPlmで計算できる。しかしIDLと同様,lが約150よりも大きくなるとオーバーフローするため計算できない(すくなくともver.1.11までは,http://www.cbrc.jp/~tominaga/translations/gsl/gsl-1.11/gsl11.pdf)。
日本語だと地球流体電脳ライブラリ(DCL)がある。名前の通り地球流体計算向けのサブルーチンが充実している。緯度は等間隔緯度のみでガウス緯度には非対応。データ解析用ということらしい。結論から言うとdcl-5.4.8を入れようとしたら挫折したが, dcl5.9betaはどうやら導入できたようだ。
まずはdcl5.4.8の導入から。ソースコードからコンパイルしようとしたら
http://d.hatena.ne.jp/Koshida_Tomoki/20111023/1319342778
のようなエラーが出る。どうやらX11に関するライブラリが無いらしい。
http://www.gfd-dennou.org/arch/gtool/gtool4-tools-library/gtool4-tutorial/trmm/doc/dcl_ifc.html#whats
/usr/X11R6/include 以下に Xlib.h, X.h といったヘッダファイルが存在しているか
と書かれているのだが,そもそも/usr/X11R6/include/X11以下に Xlib.h, X.h といったヘッダファイルが存在しているのでは? XINCPATHを/usr/X11R6/include/X11としてもダメ。仕方ないので諦める。
バグだらけの安定版よりもバグ修正された開発版を使った方がいいというのは吉里吉里で学んだので, dcl-5.9.0 betaで試してみる。
export FC=gfortran
./configure
make
sudo make install
したら,サンプルコードのビルドと実行まで成功した。
球面調和に関するリファレンスはhttps://www.gfd-dennou.org/arch/dcl/pdf/math2.pdfにある。
SHTOOLSというライブラリがある。こちらは主に地球力学・地球電磁気学での使用を想定していることが,サンプルコードや球面調和関数の正規化の方法のデフォルトからわかる。基本的にhttp://d.hatena.ne.jp/hotta_hideyuki/20121217/1355729259を参照した。依存関係のあるFFTWとLAPACK&BLAS。緯度は等間隔緯度と,ガウス緯度の両方に対応。
FFTWはmacportsから普通に行けた。/opt/local/includeにヘッダが,/opt/local/libにライブラリがインストールされる。したがって,コンパイル時には
$ gfortran test.f90 -I/opt/local/include -L/opt/local/lib -lfftw3
LAPACKのインストールはちょっと手間取った。先にBLASをmakeする必要がある。詳しくは上記のページへ。ビルドが終わるとliblapack.a, librefblas.a, libtmglib.aの3つが生成される。この3つのファイルは好きな場所において使う。僕はとりあえず~/lib直下に置くことにした。LAPACKを使う際は
$ gfortran test.f90 -L$HOME/lib -llapack -lblas
とすればよい。"-L~/lib"とするとエラーが出る。Makefileをつくるときは"-L$(HOME)/lib"としない動かなかった。(実行されるシェルが違うのかも)
ここまではhttp://www.rcs.arch.t.u-tokyo.ac.jp/kusuhara/tips/linux/fortran.htmlが参考になる。
SHTOOLSのビルドはまずDL&解凍してできたMakefileをいじる。
gfortranのオプションのうち'-march=native'を削除
そして$ make を実行。これでビルドできた。箇条書き2つめを入れずにやると
/var/folders/p_/3mcbcnfd6f5fp5xclbvb9slh0000gn/T//ccYTqNNk.s:508:no such instruction: `vmovsd (%rdx), %xmm0'
みたいな感じのエラーが吐かれる。このエラーに対してはgfortranのコンパイルオプション'-march=native'を削除すればいいらしい。http://stackoverflow.com/questions/10327939/erroring-on-no-such-instruction-while-assembling-project-on-mac-os-x-lion
Compile your code with the following flags:
gfortran -Imodpath -m64 -O3 -Llibpath -lSHTOOLS2.10 -lfftw3 -lm
where modpath and libpath are replaced with their respective paths.
のように使い方が表示される。こちらもまた好きな場所に置ける為,~/libに置いた。
コンパイルオプションをいちいち書いてられないので,Makefileを作る。Makefileの入門にはhttp://www.jsk.t.u-tokyo.ac.jp/~k-okada/makefile/が分かりやすい。
shtools_test: shtools_test.f90
gfortran $^ -o $@ -I$(HOME)/lib/SHTOOLS-2.9.1/modules -m64 -O3 -L$(HOME)/lib/SHTOOLS-2.9.1/lib -lSHTOOLS2.10 -I/opt/local/include -L/opt/local/lib -lfftw3 -L$(HOME)/lib -llapack -lblas -lm -latlas
Undefined symbols for architecture x86_64:
"_ATL_cGetNB", referenced from:
_ATL_ilaenv in liblapack.a(ATL_ilaenv.o)
"_ATL_dGetNB", referenced from:
_ATL_ilaenv in liblapack.a(ATL_ilaenv.o)
"_ATL_dgemaxnrm", referenced from:
_ATL_dgels in liblapack.a(ATL_dgels.o)
Undefined symbols for architecture x86_64:
"_dfftw_destroy_plan_", referenced from:
石岡圭一氏の作成したFortran用数値計算パッケージにISPACKがある。緯度はガウス緯度のみ対応。