反応拡散に関して早いライブラリがあると言うことなので，広島大学の小林亮先生からもらう。implicit でdtを大きくとれるので速い，とのこと。インストールまで小林先生がやってくれてちょっと恐縮。
とりあえずドキュメントを読んでみると，ライブラリは拡散項の部分を計算してくれるらしい。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をみたら、後ろでパラメータファイルを指定しろ，ということだったので，
としたらあっさり動いた。ほう。
動いたところで、自分の使いやすいように改造してみる。グラフィックスのライブラリはUNIX式よりも、テキストで吐き出してImageJで読んだ方が楽なので，適当に変えてみる。
メインルーチンがfortmain.fなので開けてみる。まず
というところについて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')
というのもよくわからなかったので説明をあさる。
unit=...というのは多分ファイル番号であろう。次のparamというサブルーチンを見ると，00.paramからパラメータを読み取るルーチンのようだ。
次の preset というサブルーチンはどうやら計算に使う固有ベクトルを計算しているようだ。
次のinit.fでデータを読み込んでいる。icaseという変数が１ならデフォルトと一緒，負の場合はinit00.dataから読み取るようになっている。
ImageJで５１２x５１２の初期配列を作って，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というでかいファイルが出来てはくれるのだが（簡単だ！）８０文字で行を変えてしまうためImageJで簡単に読み込めない。formatに一行あたりの文字数制限があるためらしい。
探していたら
というオプションを入れたらいいらしい，というのがあるので、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がでかくても出なくなる，というもの。原理は線形安定性解析をするとわかる。求めるためには１ステップごとに結構でかい配列の逆行列を求めないといけないが，その部分をフーリエ変換を使うとただのかけ算に出来る，ということらしい（畳み込み積分がただのかけ算になるというあれか？）
それとようわからんのが、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. シンタックスエラーが出る＞一行の文字数が７２文字になっていると，それより長い行は文法エラーが出てしまう。-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分をコメントアウトしたら動いた。拍子抜け。
最終更新：2008年05月07日 02:16
