Fles fft
とりあえずドキュメントを読んでみると,ライブラリは拡散項の部分を計算してくれるらしい。FitzHugh-Nagumoのサンプルがあったのでとりあえず走らせてみよう。
コンパイルは,makefile のfort77をg77にしたらあっさり通る。しかし、そこで生成される Type1.out というのを実行しようとしても
>open: No such file or directory apparent state: unit 9 named (null).param lately writing direct unformatted external IO Abort
というエラーが出て止まる。何だ?おちついてSample_Program の中のREADMEをみたら、後ろでパラメータファイルを指定しろ,ということだったので,
>./Type1.out 00
としたらあっさり動いた。ほう。
動いたところで、自分の使いやすいように改造してみる。グラフィックスのライブラリはUNIX式よりも、テキストで吐き出してImageJで読んだ方が楽なので,適当に変えてみる。
メインルーチンがfortmain.fなので開けてみる。まず
> call fname (filename)
というところについてfname.cを開いてみると、どうやらfilenameという配列にパラメータや初期値,グラフィックス等のファイル名を入れているようだ。
>filename[0] = 00.param filename[1] = term00.data filename[2] = init00.data filename[3] = Graph00
ということらしい。ここだけc。fortranだとファイル周りがよくないのか?次の
> open (unit = 9, file = filename(1), status = 'old')
というのもよくわからなかったので説明をあさる。
http://homepage.mac.com/galois21/fortran/tips/tips3.html
unit=…というのは多分ファイル番号であろう。次のparamというサブルーチンを見ると,00.paramからパラメータを読み取るルーチンのようだ。次の preset というサブルーチンはどうやら計算に使う固有ベクトルを計算しているようだ。次のinit.fでデータを読み込んでいる。icaseという変数が1ならデフォルトと一緒,負の場合はinit00.dataから読み取るようになっている。
ImageJで512x512の初期配列を作って,init04.dataというのに書き出して、04.param のicaseを-1にしてみたら
fmt: read unexpected character apparent state: unit 20 named init04.data last format: (10e15.6) lately reading sequential formatted external IO Abort trap
というエラーが出た。フォーマットがおかしいようなので,read文の中の (10e15.6)を*に変えたら通った!結構簡単。次は出力。GLSCでリアルタイムに出すよりは,テキストで吐き出してImageJに読ませる方が楽なので,グラフィックスのルーチンのところをwrite 文に変えてみる。
write (12,*) u c call graph (imax, jmax, u, v, time)
そうすると、term04.dataというでかいファイルが出来てはくれるのだが(簡単だ!)80文字で行を変えてしまうためImageJで簡単に読み込めない。formatに一行あたりの文字数制限があるためらしい。
探していたら
-ffixed-line-length-none
というオプションを入れたらいいらしい,というのがあるので、Makefile のfflags に入れて試してみるが、全然変わらない。オプションを探すのが馬鹿馬鹿しくなって来たのと、uとvを並べて出力するため、
FortranであるWRITE文で改行せず, その次のWRITE文で同じ行に書き続けたいときには 書式 (FORMAT) の最後に$を付ける. こんなかんじ:
program suppress_new_line implicit none write(6,'(a,$)') 'abc' write(6,'(a)') 'def' write(6,'(a)') 'ghi' end program suppress_new_line
ということだったので、
do j = 1, jmax do i = 1, imax write(12, '(f8.4,$)') u(i,j) end do do i = 1, imax write(12, '(f8.4,$)') v(i,j) end do write(12, *) ' ' end do
としてやってまあ普通に出力できるようになった。
ここからは支配方程式のチェック。どうやらscheme.fでやっているらしい。で,陰解法というのがなんだかよくわからなかったのでチェック。…どうやら
(u(x, t+dt) - u(x,t))-dt= d (u(x+dx,t)+u(x-dt,t) -2u(x,t))/dx^2
のかわりに
(u(x, t+dt) - u(x,t))-dt= d (u(x+dx,t+dt)+u(x-dt,t+dt) -2u(x,t+dt))/dx^2
とすると、離散化による誤差がdtがでかくても出なくなる,というもの。原理は線形安定性解析をするとわかる。求めるためには1ステップごとに結構でかい配列の逆行列を求めないといけないが,その部分をフーリエ変換を使うとただのかけ算に出来る,ということらしい(畳み込み積分がただのかけ算になるというあれか?)それとようわからんのが、scheme.fで、拡散している因子が一つしかねえ。uやvの配列をそのまま使っているらしいが,どうせちょっとしかかわらんからそのまま使っているということか? …で、一次元のファイルを少し見てみると,どうやら固有値のマトリックス(eig_u, eig_v)をpreset.f で最初に作ってやるらしい。ということは、domain growth の効果を入れるのはそんなに簡単ではないのか?固有値マトリックスを二つ作って,二度呼び出してやることでいちおうiBookG4では動くようになった。
これをG5マシンに移植しようとするといくつかエラーが出る。 1. glsc.hが読めないというエラーが出てくる>>glsc.h、glsc_ftn.h をソースコードと同じディレクトリにコピーする,もしくはMakefile のオプションに -I/Users /miura/include 等を入れる 2. シンタックスエラーが出る>一行の文字数が72文字になっていると,それより長い行は文法エラーが出てしまう。-ffixed-line-length-none を入れる。
この二つで大概動いた。あとは “ld: table of contents for archive … is out of date;”
というメッセージが出たので,
>sudo ranlib ~/lib/libfles_fft.a
をかけたら通った。
Intel Mac
Intel Mac にかわったので、fles_fft を再コンパイルしてみる。とりあえず、GNU gccではgfortran しか使わなくなるらしいので、まずMakefile のコンパイラをgfortranにかえてみる。すると
>mt19937.f:38.35: > > parameter(UMASK = -2147483648) > 1 >Error: Integer too big for its kind at (1)
というのが出てくる。でかい整数でエラーが出るようなので、
>FFLAGS = -O2 -fno-range-check
というオプションを入れてみたら、コンパイラは通った。
次にglsc をコンパイル。こちらもあっさり通る。sudo make install で/usr/local/lib にインストールし、glsc_ftn.h を/usr/local/include に入れる。
しかし、
Error: Non-numeric character in statement label at (1) /usr/local/include/glsc_ftn.h:6.2: Included at graph.h:12: Included at gprep.f:15: integer G_BLACK, G_RED, G_GREEN, G_BLUE 1
という感じのエラーが大量に出て止まる。何だ??
glsc.hを使うことによって出ているエラー臭いので、このライブラリを全然使っていないSutureLineModel (データをテキストで吐き出す。)のところで試してみる。単純にinclude分をコメントアウトしたら動いた。拍子抜け。