HOME >  技術情報 >  Fles fft

Fles fft 

反応拡散に関して早いライブラリがあると言うことなので,広島大学の小林亮先生からもらう。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をみたら、後ろでパラメータファイルを指定しろ,ということだったので,

>./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分をコメントアウトしたら動いた。拍子抜け。

ページトップへ戻る
-->
Copyright(c) 2013 九州大学大学院医学研究院 系統解剖学分野 All Rights Reserved.