fortran 備忘録

確保した配列をオーバーした箇所を調べる

コンパイルオプションに -fbounds-check を入れる。

nan や inf の発生場所を調べる

1) プログラム trapfpe.c を作る 

/*******************    trapfpe.c    **********************/
#define _GNU_SOURCE 1 
#include <fenv.h> 
static void __attribute__ ((constructor))
  trapfpe ()
{ 
  feenableexcept (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW); 
} 
/**********************************************************/
2) trapfpe.c をコンパイル
gcc -c -o libtrapfpe.a trapfpe.c

3) 対象となるプログラムを libtrapfpe をリンクしてコンパイル。
オプションの -g を忘れないこと
g77 -o tmp -g tmp.f -L/ディレクトリ -ltrapfpe

4-1) 実行ファイルを gdb 上で実行
4-2) gdb ./tmp
4-3) gdb 上で run を実行。問題箇所が表示される。

例 1: 不正な計算を行った場合 

------------------- サンプルプログラム --------------------
      program samp1
      real a,b 
      a=0.0 
      b=1.0
      write(*,*) b/a 
      end
-------------------------------------------------------------
b を 0 で割っているので INF となる。上記の方法でデバッグすると gdb に 0x0804873e in MAIN__ () at samp1.f:5 5 write(*,*) b/a と表示される。 例 2: 共有ライブラリの関数に不正な値を渡した場合 
------------------- サンプルプログラム --------------------
      program samp2
      write(*,*) acos(2.0)
      end
-------------------------------------------------------------
gdb 上で実行すると 
(gdb) run
        :
Program received signal SIGFPE, Arithmetic exception.
0x40051a9e in fegetexcept () from /lib/i686/libm.so.6
となる。これだけだとよく分からない(この場合は明らかだが）。 そこで backtrace コマンドで問題にたどり着く前の情報をみる。 
(gdb) backtrace
#0  0x40051a9e in fegetexcept () from /lib/i686/libm.so.6
#1  0x0804877c in MAIN__ () at samp2.f:2
#2  0x08048806 in main ()
#3  0x4008bc1f in __libc_start_main () from /lib/i686/libc.so.6
これなら samp2.f の 2 行目が怪しいことが分かる。この方法はCでも使える。

実行ファイルに引数を使う

iargc,getarg を使う。
------------------- サンプルプログラム --------------------
      program sample
      integer iargc,nargc
      external iargc
      character*16 argv
      integer itmp
      real rtmp
      nargc=iargc()   ! iargcは引数の数を返す
      write(*,*) '引数の数は',nargc,'個'
      do i=0,nargc
        call getarg(i,argv)  ! argv に i番目の引数を代入
        write(*,*) argv
      enddo
      call getarg(1,argv)
      read(argv,*) itmp  !argv を整数:itmpに代入
      call getarg(2,argv)
      read(argv,*) rtmp  !argv を実数:rtmpに代入
      write(*,*) itmp,rtmp,itmp+rtmp
      end
-------------------------------------------------------------
注：gfortran では iargc が無くなったので COMMAND_ARGUMENT_COUNT() を使う。
Ref.: COMMAND_ARGUMENT_COUNT

乱数

Fortan で使える乱数は g77 標準の rand や Mersenne Twister , CERNLIB の rndm, grndm 等がある。
Mersenne Twisswter のFortran77 用コード(バックアップ)。
コンパイラ、オプションによる速度の違い(Pen4 3.0 GHz でテスト)。
Mersenne Twister:grnd() を 2^27 回計算した場合

time ./a.out
ifort 3.3 sec
g77 -O3 3.8 sec
g77 13.3 sec

関数による速度の違い

rand, rndm, grndm を 2^27 回計算した場合(最適化なし)

time ./a.out
rand 5.6 sec
rndm 2.9 sec
grndm 10.2 sec

どの関数も周期性を持っているのでどの関数を使うかは自分で判断。 rndm の周期性は n 回目と n+1 回目 を２次元プロットすることで簡単に確認できる。
乱数の周期性
乱数に周期構造が現れる様子。 横軸は n 回目、縦軸は n+1 回目に call されたときの値。rndm は周期的な構造が見える。
(cernlib とGEANT3 のドキュメント)

可変長のファイルを読む

読みたいファイルの行数が可変の場合は read の引数に end を使うと便利。 配列はファイルの十分な大きさを確保しておく(配列も可変にしても良いけど）。
------------------- サンプルプログラム --------------------
      open(10,FILE='samp.dat')
      do i=1,200
         read(10,*,end=99) a(i),b(i),c(i)  ! end of file で 99 番に飛ぶ
      enddo
  99  continue
      close(10)
---------------------------------------------------------

文字列操作

文字列操作で知っておくと便利なこと
Fortran90 以降では他にも便利な関数があるけど環境が限られるので省略。
------------------- サンプルプログラム --------------------
      program sample
      implicit none
      character*4 ci
      character*6 cnum
      character*16 str1,str2
      integer num

      str1='Othello'
      str2='world?'
      write(*,*) str1//str2
      write(*,*) str1(3:8)//str2(1:5)//'!'

* integer to character
      num=9801
      write(cnum,'(I4.4)') num
      write(*,*) 'PC'//cnum
      write(*,*) 'len=',len(cnum)
      write(*,*) 'len_trim=',len_trim(cnum)

* character to integer
      ci='9'
      read(ci,*) num
      write(*,*) num**2

      end
-------------------------------------------------------------
実行結果
Othello world?
hello world!
PC9801
len= 6
len_trim= 4
81

追記: "/"(スラッシュ)は制御編集記述子なのでそのままでは文字列変数に代入できない。 "/"を文字として使いたい場合は write(cs,'(A)') '/' の様に一度変数に代入すると 大丈夫。もっとスマートな方法は無いのか？

cernlib を使ったサンプルプログラム

CERN Program Library に含まれる乱数やヒストグラム作成機能を使ったサンプルプログラムを置いておきます。詳しいことはcernlib のマニュアルを見て下さい。

プログラムのコード
main.f: モンテカルロでπを計算
hbook.f: ヒストグラムを定義
Makefile: Makefile

上のファイルをどこか適当なところにコピーして make でコンパイル。 (cernlib が無ければ当然コンパイルできない。Makefile 中のライブラリパス は各自の環境に合わせる)



[戻る]