更新日:2008-04-09
作成者:坂口修一 (SAKAGUTI Syuiti)
連絡先:sakaguti@yamaguchi-u.ac.jp

Fortranの玉手箱


このページを読む上での注意

このページの記事は Fortran 90 規格を前提とする。筆者は通常 Fortran 90 規格対応のコンパイラを UNIX 系の環境(Absoft Pro Fortan - MaxOSX か Intel Fortran Compiler - Vine Linux)上のコマンドライン(シェル,bash)から使用している。Mac とある場合は MacOSX (10.x.x) のことで MacOS 9.x 以前ではない。このページの記述が実行環境に依存したものならないように注意しているつもりではあるが,すべての環境について確認することは不可能であるし,実際のところ上記環境の範囲でしか確認していない。

ところどころに記載されている小物(◆)・習作(◇)ソースは,かつて筆者が自分の練習用に書いたものを修正した程度のものである。大したものは無いが,同様の事を初めて行うときには少しは参考になるかも知れないし,備忘用もかねてここに記載する。なお,このページの記事はできるだけ正確で新鮮である様に努めてはいるが,一切の保証はなされないし,このページの情報の品質および使用したことによる損害等について筆者は一切の責任を負わない。


コンパイルエラーやバグの対処法のメモ

1. コンパイルが完了しない場合

ちょっとした書き間違いや宣言忘れなどコーディング上の単純なミスながほとんどのはず。使用するフォントにもよるが,1(数字のイチ)と l (英小文字のエル)のような変数名の打ち間違いは,発見するのは結構厄介。書き捨てのものやよほど小規模なものでない限り,implicit none を指定し,使用する変数をいちいち宣言するようにした方が無難。この段階のエラーはコンパイラが出すものの他にリンカが出すものもある。

1.1 はじめから自分で書いたものや小規模のもの

  1. コンパイラの出すエラーメッセージをよく読み,バグの場所と原因を把握する。大量のエラーが報告される場合も実質のエラーは少ない場合が多いので,ひとつ目のエラーから順番に潰す。大抵のコンパイラには報告するエラーの数とレベルを指定するオプションがあるだろうから,ひとつエラーがあった時点でコンパイルが停止するように設定するのもよい。
  2. implicit none を指定している場合,必要な変数がすべて正しく宣言されているか,型や長さなどを確認する。宣言した変数が使用されていない,という類いのワーニングは潰すに越したことはないが無視してもよい。
  3. do 文や if 文が閉じていなかったなどのミスは,end do や end if を必ず書くようにした上で,Emacs や Mule のような Fortran ソースコードをインデント整形してくれるエディタを使用すると楽に発見・修正ができる。
  4. 文法的なミスがないのにコンパイルが完了しないときは,リンク段階でエラーになっていることがある。このときはコンパイラではなくリンカがエラーを出すので,エラーメッセージ中には f90: などコンパイラ名ではなく ld: などのリンカ名が入る。ソースファイルがひとつのみの場合は,ソース中で使用している関数やサブルーチンが Fortran 規格の組み込み関数なのか,あるいは現在使用している処理系に組み込まれているかそうでないのかを確認する。コンパイル時にその関数等を含むライブラリを呼ぶためのオプション指定が必要な場合もある。雑誌やネットから用例を拾ってきて,とりあえず試しに書いてみたが動かない,という時は見慣れない関数等がないか探す。ソースファイルが複数あり分割コンパイルをしているとき,個々の実行オブジェクトやモジュール(*.o, *.mod)はコンパイルできるのに,すべてをリンクした実行用のモジュールができない場合は,必要なオブジェクト,モジュールやリンクするべきライブラリが揃っているか,コンパイルのし忘れがないか確認する。このミスは make の使用で防げる。
  5. 分割コンパイル時は,コンパイルの順序を確認するのも有効。Makefile を書いておくとよい。
  6. make を使用している場合,Makefile 中の tab であるべき部分を space にしていないか。久しぶりに書き直したりするとやりがち。

1.2 処理系を変えた場合

機種更新や間借りのため,あるいはできるだけ多くのマシンを使うためなど,処理系を変えてコンパイルする場合は1.1 で挙げた点に加えて次の点を確認する。

  1. 漢字コード,行末コード。UNIX 系,Mac 系,Win 系で異なる。
  2. 自作のものも含め,リンクさせているものが揃っているか。
  3. コンパイラのインストール状況(種類,バージョン,コンパイラのコマンド名)。
  4. IMSL, SSL,MKL など使用しているライブラリのインストール状況。
  5. コマンドサーチパス,ライブラリサーチパスやシェル変数の設定。ついでにスタックサイズとコアサイズも確認しておくとよい。
  6. コード中でコンパイラ依存の方言が使用されていないか。またコンパイラ独自のコンパイルオプションが使用されていないか。
  7. コード中で処理系依存ライブラリが使用されていないか。特に UNIX/MacOSX と Win の間で移動する場合。

1.3 自分で書いたものではない場合

他人から引き渡されたものが一発でコンパイラを通らないことはよくある。1.1, 1.2 で挙げた点に加えて次の点を確認する。

  1. Makefile があるならばその内容。
  2. common 文まわり。
  3. マクロ指定などでプリプロセッサが必要か。
  4. そもそも動くコードなのか。動いたものなのであれば,開発時に使用していたコンパイラの種類と,指定してしていたオプション。
  5. Fortran の仕様。FORTRAN IV, FORTRAN 66, FORTRAN 77, Fortran 90, Fortran 95 のいずれか。普通は FORTRAN 77 かFortran 90 規格に準拠したものが多いはずだが,古くからあるサブルーチンを含むファイルなど一部のソースコードが古い規格で書かれていたりしないか。

2. 実行中に異常終了する場合

2.1 実行中に異常終了する場合一般

  1. コンパイラの出すエラーメッセージをよく読み,バグの場所と原因を把握する。
  2. エラーメッセージからだけでは場所と原因がわからない場合は,デバッグライトやデバッガでの変数値トレースや進行状況トレースが効果的。
  3. 実行ファイルを直接実行するのではなく,スクリプトから実行するように設計されていないか,確認する。実行時に必要なパラメータを実行ファイルに渡したり,入出力するファイルを整えたり,いくつかの実行ファイルを続けて起動するなどの目的でドライバスクリプトを書くことがある。UNIX 系 OS ならば普通 sh か bash で書かれている。

デバッグライトとは,デバッグ目的にコード中に書いておく write 文や print 文のこと。デバッグの基本的な方法として非常に有名。基本的には do 文 や call 文 の前後など処理の塊ごとに入れ,変数値の変化(変数値トレース)や実行時にどこまで処理が進んでいるか(進行状況トレース)を調べることが目的。標準エラー出力(装置番号は0のことが多い)に書き出すように書くことが多いが,データの検証が目的の場合はファイルに書き出すこともある。実行速度を落とす原因になるので,デバッグ完了後は必要に応じて削除する。

デバッガとは,デバッグを支援してくれる開発ツールで,コードを1行ずつ実行できたり,実行中の変数値トレースするなどができるもの。gdb (GNUデバッガ)や dbx (BSD系)などが有名。デバッグライトとは異なり,ソース中に write 文や print 文を埋め込む必要はないので,コードはきれいなままにしておける。デバッガを使用するためには専用のデータファイルが必要で,コンパイル時にコンパイラオプションを指定して生成する。一般にこのオプションはコードの実行速度を落とすので,コードが完成したらこのオプションは落とす。

2.2 Segmentation Fault や Access Violation への対処

配列まわりのミスを疑う。配列参照,配列境界,配列名,配列長,配列の次元など。

  1. 宣言されていないような配列を参照・指定していないか確認する。宣言忘れミスは implicit none 方式にすれば防げる。
  2. 配列の次元を間違えていないか確認する。例えば,a(i) と a(i,j) を間違えていないか。
  3. 配列の下限や上限を超えて参照したり,サイズを超えた使用をしていないか確認する。例えば,a(1:n) の a(0) や a(n+1) を参照したりしていないか。特に各次元ごとのサイズが異なる様な多次元配列で do 文を回すときは注意する。また配列を宣言したり割り付けるときには上限のみならず共に下限も明示的に指定するようにする。size() で配列の大きさを確認するのも有効。
  4. 宣言した配列の境界を超えて参照していないかを確認するためのコンパイルオプションを付けてコンパイルしてみる。
  5. シェルのスタックサイズを確認する。大きな配列を宣言した場合,スタックがあふれて止まることがある。bash の場合は ulimit -a でスタックサイズの確認ができる。スタックサイズが原因の場合は,ulimit で割当を変更する。または,スタックを使用しないで記憶領域にすべて静的に割り当てるコンパイラオプション(-s や -static など。コンパイラの説明書を読んで探す)を指定してコンパイルしてみる。
  6. allocatable 配列を使っている場合は,deallocate した後にその配列を使用していないかを確認する。

2.3 入出力関連

  1. 入力が正しいか確認するには read したものをそのまま write してみる。
  2. 実行時に入出力を行う場合,必要なファイルとディレクトリがすべて揃っているか。
  3. ファイルからデータを読み込んでいる場合は,そのファイルに十分なデータが揃っているか,行数などを確認する。
  4. 読み込み方が正しいか。ひとつの read 文でいくつのデータを読み込んでいるのか,データ読み込みの繰り返し方に注意する。
  5. 入出力が書式付きで行われているのかどうか。書式付きならその書式を確認し,期待した通りの入出力が行われているか。

2.4 subroutine, function 関連

  1. subroutine や function の引数の個数は正しいか。Bus error が出た場合など。
  2. 親ルーチンで function の型宣言を忘れていないか。

3. 実行終了するが計算結果や出力がおかしい場合

探し当てるのに苦労する事が多い。すでに結果が分かっている計算を行い,期待通りの結果を返すかなどして絞り込みを行う。必要なら全変数書き出しも行う。

3.1 コーディング上のもの

原因例

  1. 配列演算を間違えていた。順次代入される書き方と一度に代入する書き方の勘違いなど。
  2. 掛け算と割り算を間違えていた。a/b/c, a/(b*c),a/b*c などの勘違いなど。
  3. 整数型変数の割り算をしていた。割り切れない場合に予想外の結果が出る。
  4. 実数型変数の計算に整数型変数が混入していた。
  5. 実数型変数に定数を代入する際,精度を間違えていた。
  6. 関数に渡した引数の型が合っていなかった。
  7. 呼び出している組み込み数学関数の型が合っていなかった。
  8. 数式中の括弧が不整合だった。
  9. ファイル等からのデータ読み込みが正しく行われていなかった。ファイル内容が不正,あるいは入出力の対象の間違い,open - close の不一致など。
  10. 変数宣言はしてあったが,型が間違っていた。
  11. 変数のスコープを勘違いしていた。
  12. 文字列変数のトリミングを忘れていた。
  13. 初期化していない変数を使ってしまっていた。
  14. 配列を参照する際に宣言されている境界を超え,たまたまそこ(別な変数名で確保されている領域)にあった別なデータを読んでしまっていた。
  15. 実数変数の精度設定や浮動小数点計算の精度を変更するようなコンパイルオプションが適切でなかった。
  16. コンパイラにバグがあることが明らかだったが,パッチを当てたりアップデートするなどの修正をしていなかった。
  17. function 名が組み込み関数と同じになっていないか。

対処例

  1. 変数を使用する前には初期化か代入を行う。コンパイラオプションを使用しても良いが,可搬性は落ちる。
  2. 結果が明らかなもので検算する。
  3. 他人からもらったコードの場合,コンパイル時に実数変数の精度を変更するオプションの指定が必要な場合がある。コンパイラによっては,コンパイル時に実数変数(整数変数も可)の精度を一括指定し,real を すべて real(8) だとか real(16) として扱う様にできるものがあるので確認する。
  4. 分割コンパイルをしている場合は,すべての更新が反映されているか。
  5. ループカウンタや分岐判定用など,整数変数とすべきところに実数変数を使用していないか確認する。
  6. 入出力において装置番号 0, 5, 6 番の使われ方を確認する。慣習に基づいて,それぞれ標準エラー出力,標準入力,標準出力の装置番号として使われている事が多い。
  7. コンパイラの説明書を読む。
  8. 他人からコードを受け取る時は Makefile の有無やコンパイルの方法を確認する。

3.2 論理的なもの

原因例

  1. 手計算ミスなど,コーディング以前の論理的な段階で間違いがあった。
  2. 参照した文献の数式に間違いがあった。版が古い教科書なども誤植が残っている場合がある。
  3. 数値計算的な手法の選択に問題があり,十分な精度や速さが出ていない箇所があった。

高速化のためのメモ

設計・コーディング

  1. 数値計算に優しい形式に手計算しておく。
  2. (速い)足し算>掛け算>割り算(遅い)。
  3. 配列変数はメモリ上のアドレスの順に参照すると速い。配列 a(i, j) は i から回す。
  4. 速いアルゴリズムを採用する。
  5. 既に先人が知恵を絞って作った,よく出来たライブラリ関数を使う。使えるものなら,既知の手続きはわざわざ自分で書かない。
  6. 入出力を減らすと速い。CPU time は短くても,HDへの書き出し読み込みなど I/O に時間がかかると遅い。 real time が長くなる。
  7. 扱っている問題にもよるが,ベクトル化・パラレル化できると速い。専用のコンパイラが必要なことがある。
  8. ホットスポットを潰せると速い。ホットスポットはプロファイラで探す。

コンパイル・リンク

  1. g77よりは商用のコンパイラの方が速い実行ファイルが出来る。
  2. デバッグやプロファイリングのためのオプションを付けたままだと遅い。行列境界チェックオプションも同様。
  3. 最適化オプション(いろいろある)を付けてコンパイルすると速い。 -O2 のようなものが大抵用意されている。 -On の n が大きくなると,コンパイラが高速化のためにいろいろなことを陰でしてくれる。コンパイルにかかる時間は長くなる。説明書を読んでおくとよい。
  4. 使っている CPU と相性のよいコンパイラを使用すると速い。CPU の種類(PPC だとか Pentium IV だとか)を指定してチューンするオプションを使用する。
  5. 使っている CPU に最適化された(その CPU のマシン語で書かれた)数値計算ライブラリを使うと速い。
  6. 数学関数で精度重視と速度重視のどちらかが選べる様なオプションがある場合がある。

※ここでオプションと言っているのは,コンパイル時にコンパイラ等に渡すコンパイラオプションのこと。

プロファイリング

実行時間を短縮したくなった場合,まずはどの部分(subroutine や function)で時間がかかっているかを次の方法等を用いて把握する。

  1. gprof などのプロファイラを使う。これが一番楽。コンパイル時にオプション(コンパイラによるが大体 -P とか -p とか)を与えて,プロファイラで解析するためのファイルを作成する必要がある。
  2. プログラムに実行時間を吐くコードを入れておく。パラメーターを振りながら計算するときに有効な事もある。
  3. time コマンドを使う。ただし複数の実行ファイルを連続して実行する時くらいしか意味が無い。
  4. 無駄な入出力は削る。

特に時間のかかる部分(ホットスポット)や特に多数回実行されている部分を改善すると効率が良い。逆に,どの部分も均等に働いているようなら,ちょっとした作業で劇的な高速化を期待するのは難しい。小手先の修正で得られる以上の高速化を望む場合はアルゴリズムやデータ構造を含めてプログラムの再検討が必要。

高速化作業は意外と楽しいのでついつい時間をつぎ込みがちだが,何回も走らせないようなコードだと,その作業時間分の時間が回収できないかも知れないので要注意。科学技術計算用コードで使用予定人数が極端に少ない場合,開発時間も実行時間のようなものである。人間が寝てようが遊んでいようがマシンはずっと計算してくれる。実行実時間で数日以上時間短縮できないなら,放っておいてもよいのではないだろうかと思う。初めて高速化方法を勉強している時は当然別問題で大切な作業。


入出力に関するメモ

装置番号

書式(データ編集記述子)に関するメモ

変数値を write するときに形を整える。何回も使うような形式は format 文に定義してもよい。

入出力関係のコーディング例

◆変数の標準入力からの読み込み

変数を標準入力から読み込む。教科書にもよく載っているだろう。要求された分の読み込みが終わらない限り,先には進まないので,どのような変数をいくつ入れる必要があるのか,write しておく方がよい。入力は区切ってあれば一行で与えても複数行で与えてもよい。UNIXマシンで実行する場合で,変数入力の部分をシェルスクリプトに書き込みたいときは,echo "var1 var2 ... " | ./a.out のようにして echo や cat とパイプで入力値を渡してやればよい。

!---------------------------------------------------------------------
!     read.f90 (Fortran90)
!     SAKAGUTI Syuiti 
!     06-Mar-2005
!---------------------------------------------------------------------

program read

  implicit none
  integer :: i, k
  integer, dimension(1:6) :: na

  read *, (na(i), i=1,6)
  write(*,'(1i4,1i4,1i4)') (na(i), i=1,3)
  write(*,'(1i4,1i4,1i4)') (na(i), i=4,6)

end program read

◆何行あるか分からないファイルを最後まで読み込みたい。

iostat - if exit 法を使う。下に例を示す。この例では入力ファイルは 11 とし,open-close は省略した。read 文中にある iostat は read 文の結果(終了状態)を返す整数型変数で,この例では ios をあててある。入力ファイルの最後まで行って read 文がエラーを返して ios=-1 となったとき, exit 文で do ループを抜ける。

   do i=1, ndm
     read(11,*,iostat=ios) a(i), b(i)
     if (ios .lt. 0) exit
   end do

◆コマンドライン引数の読み込み・文字列変数からの数値読み込み。

C言語でよくやる様に,プログラム実行時にコマンドラインから引数を読み込む事ができる。read されるファイルを用意する必要がないので,実行時にパラメータを渡したい場合などに大変便利である。この機能は UNIX 系 OS 上の Intel (旧 Compac, 更にその前は DEC) の Fortran および MacOSX 上の absoft Pro Fortran で実装されていることを確認した。UNIX 系 OS ならば使用できる公算が高い。将来これと同等の機能が Fortran の標準として組み込まれる予定の様だが,現時点ではまだ拡張機能に当たると思われるので,安易に使用すると移植性が落ちるので注意が必要。

この方法でコマンドライン引数の読み込みを行うと,文字列変数である argv に値が格納される。しかしながら,数値計算コードにおいて実行時に渡したいものはものは大抵は数値で,文字列であることはむしろ稀だろう。do文ではできないような方法で,あるいはもっと柔軟にパラメータを振りながら計算を行いたいから,実行時に決定したいような数値が出てくるのだと思う。さて,argv に入った値を integer, real 等の数値変数として使うには,内部入出力の方法を使う。read 文や write 文の引数は普通 file 番号か標準入出力を示す * (とフォーマット指定)だが,ここには内部変数を書いても構わず,これを内部入出力と言う。ここでは入力の例を示すが,後述する様に format 文と write 文で数字の入った定型文字列を作るということもできる便利な方法である。

!---------------------------------------------------------------------
!     getarg.f90 (Fortran90) (UNIX Library)
!     SAKAGUTI Syuiti
!     24-Aug-2006
!
!     compile with absoft pro Fortran 90:
!     f90 -lU77 -o getarg getarg.f90
!
!     argc: total number of argument
!     argv: string for arguments
!     getarg(i,argv): subroutine to get arguments from command line.
!            0th argument is the name of excutable and 1st argument 
!            is the first word follow to excutable.
!---------------------------------------------------------------------

program gtarg
  
  implicit none
  character(32) :: argv
  integer :: argc, iargc, i
  real(8), allocatable :: a(:)
  
  argc=iargc()
  write(*,'(a,1x,1i1)') 'argc:', argc  
  
  allocate(a(0:argc))  
  do i=0, argc
     call getarg(i,argv)
     ! Internal io
     read(argv,*) a(i)
     write(*,'("argv(",1i1,"): ",1a,1es10.4)') i, argv, a(i)
  end do
  deallocate(a)

end program gtarg

このソースをコンパイルして実行した例を示す(コンパイル条件は Absoft Pro Fortran 90 - MacOSX)。コマンドライン引数として, 0, 2, 1.0, 1e3, 1d2, a を与えた。ソース中で a(i) を実数変数として定義してあるので,文字変数である a を読ませた際にエラーを出して止まっている。実用する際にはこの例の様におざなりに配列を do 文で回すのではなく,適切な変数を用意して適切に read 文を書いてやれば問題はない。

$ f90 -O -lU77 -o getarg getarg.f90; ./getarg 0 2 1.0 1e3 1d2 a

argc: 6
argv(0): ./getarg                        0.0000E+00
argv(1): 0                                 0.0000E+00
argv(2): 2                                 2.0000E+00
argv(3): 1.0                              1.0000E+00
argv(4): 1e3                             1.0000E+03
argv(5): 1d2                             1.0000E+02

? FORTRAN Runtime Error:
? Illegal character in numeric input
? READ(UNIT=INTERNAL,...

◇ファイルの開閉・読み込みと書き出しの練習

幾つかのファイルを開いて,読み込み,書き出しをして閉じる。

!---------------------------------------------------------------------
!     file_io_several.f90 (Fortran 90)
!     SAKAGUTI Syuiti
!     24-Aug-2006
!---------------------------------------------------------------------

program file_io_several
  
  implicit none
  integer :: i, n=9
  character(32) :: a(1:2)
  
  open(21,status='unknown',file="aaa.out")
  open(22,status='unknown',file="bbb.out")
  do i=1, n
     write(21,*) i, ': test write for file 21'
     write(22,*) i, ': test write for file 22'
  end do
  close(21)
  close(22)
  
  open(21,status='old',file="aaa.out")
  open(22,status='old',file="bbb.out")
  do i=1, n
     read(21,'(1a)') a(1)
     read(22,'(1a)') a(2)
     write(*,*) i, a(1), a(2)
  end do
  close(21)
  close(22)
  
end program file_io_several

◆連番ファイルの開閉

名前に数字を含んだ,シリアライズされたファイル名のファイルをいくつも open / close する際に便利な方法。数値計算の際の中間結果出力やパラメータに応じた入力用ファイルを扱うときに,数値パラメータをファイル名に含めることを念頭に置いている。 このような場合,例えば a というパラメータの値に応じて, a0010.outfile, a0020.outfile, a0030.outfile, ...., a1000.outfile というような出力ファイルを作りたいとする。このとき,a の値に関係なく同じファイル名で write しておいて後からファイル名を手で変えるようにしてしまうと後で混乱が起きやすい。またパラメータを振る場合,ファイルが大量になるので,いちいち open 文を書き換えて名前を変えるのは面倒だ。ファイル名が書き込まれたファイルを作成して,そこからファイル名を読み込むという手もあるが,それもやはり面倒くさい。そこで,実行時にパラメータの数字を自動的にファイル名に取り込む方法を示す。

内部入出力の方法を使用する。 まずファイル名のための文字配列 filename を定義する。ファイル名を決定する際には,フォーマット付きで write 文を使用して filename に書き込む。最後に open 文でファイル名に filename を指定してファイルを開く。ファイル番号は使い回せるので,close 文は普通に書いておけばよい。具体的なコーディング例を次に示す。

!---------------------------------------------------------------------
!     openfiles.f90 (Fortran 90)
!     SAKAGUTI Syuiti
!     24-Aug-2006
!---------------------------------------------------------------------

program openfiles
  
  implicit none
  integer :: i, n=10
  character(64) :: filename
  
  open(10,file='openfiles_namelist.out')
  do i=1, n
     write(filename,"('openfiles_',1i2.2'.out')") i
     write(10,*) filename
     open(21,status='unknown',file=filename)
     write(21,'(1i2.2,": ",a)') i, "A test write of openfiles."
     close(21)
  end do
  
  close(10)
  
end program openfiles

日付と時間の処理に関するメモ

Fortran 90 以降,実時間時計からの情報を返す手続きが組み込み手続きとなった。手続きは date_and_time と system_clock の2種類がある。これらは function ではなく subroutine として定義されているので call して使用する。2個目以降の引数は省略可能である。従来マシンから日付情報やCPU時間を取得するには IMSL などのライブラリや etime, ctime といった外部関数を用いた事が多いと思われるが,ソースの可搬性を考えるならこの組み込み手続きを使用した方がよい。

character(10) :: date, time, zone
integer :: values(1:8)
integer :: count, count_rate, count_max)

call date_and_time(date, time, zone, values)

call system_clock(count, count_rate, count_max)

date_and_time は日付処理のためのもの。date, time, zone のいずれも 文字型で長さは 10 あれば十分。それぞれ,現在の年月日,時分秒.ミリ秒,時間帯を返す。時間帯はGMTとの差なので,マシンが日本時間(JST)設定なら GMT+0900 なので,+0900 が返ってくる。values は integer 配列,配列長は8で十分。年,月,日,分表示の時差,時,分,秒,ミリ秒が返される。おそらく values が一番使い勝手がよいだろう。計算結果として多数のファイルを作成するような場合,自動生成したファイルに生成日付を write して記録すると計算日時が特定でき,どの時点のコードで計算したか分かるので便利である。

system_clock はCPU時間を取得するためのもの。ctime や etime と同様。count, count_rate, count_max のいずれも integer のスカラーで,それぞれ,処理系の時計の現在値,時計が1秒間に刻む回数,count の取りうる最大値を返す。count は count_max までカウントされると 0 に戻って再カウントされる様なので要注意。system_clock(count) は初めにcall したときに count=0 を,その後は call するたびにそこからの経過時間を返す。

◆日付を文字列で返す・実時間を整数スカラーで返す

次にこの2つの手続きの用例を示す。

!---------------------------------------------------------------------
!     date_time_clock.f90 (Fortran90)
!     SAKAGUTI Syuiti
!     25-Aug-2006
!---------------------------------------------------------------------

program date_time_clock
  implicit none
  integer :: i, j
  integer, parameter :: n=100
  integer, allocatable :: k(:,:)
  ! for date_and_time
  integer :: iv(1:8)
  character(10) :: d,t,z
  ! for system_clock
  integer :: ic, icr, icm
  !----------------------------------
  
  call date_and_time(d,t,z,iv)
  write(*,*)
  write(*,*) 'date:   ', d
  write(*,*) 'time:   ', t
  write(*,*) 'zone:   ', z
  write(*,*) 'values: ', iv
  write(*,*)
  write(*,'(1i4.4,"-",1i2.2,"-",1i2.2)') iv(1),iv(2),iv(3)
  !----------------------------------
  
  call system_clock(ic,icr,icm)
  write(*,*)
  write(*,*) '1st'
  write(*,*) 'count:      ', ic
  write(*,*) 'count rate: ', icr
  write(*,*) 'count max:  ', icm

  allocate(k(1:n,1:n))
  do j=1, n
     do i=1, n
        k(i,j) = 1
     end do
  end do
  deallocate(k)

  call system_clock(ic,icr,icm)
  write(*,*)
  write(*,*) '2nd'
  write(*,*) 'count:      ', ic
  write(*,*) 'count rate: ', icr
  write(*,*) 'count max:  ', icm

end program date_time_clock

!---------------------------------------------------------------------
!     eof
!---------------------------------------------------------------------

実行結果を次に示す。

$ f90 date_time_clock.f90; ./a.out 
 
 date:   20060825  
 time:   143118.443
 zone:   +0900     
 values:   2006  8  25  540  14  31  18  443
 
2006-08-25
 
 1st
 count:        0
 count rate:   1000000
 count max:    2147483647
 
 2nd
 count:        346
 count rate:   1000000
 count max:    2147483647

◆年月から年度を返す関数

整数の年と月を渡すと整数の年度を返す関数。ここで言う年度とは日本の会計年度で, n 年4月から n+1 年3月までを n 年度と呼ぶもののこと。

integer function nendo(iyear,imonth)
  implicit none
  integer, intent(in) :: iyear, imonth
  if(imonth .ge. 4) then
     nendo=iyear
  else if(imonth .le. 3) then
     nendo=(iyear-1)
  end if
end function nendo

数値計算に関するメモ

計算時間

◇多次元配列に対し do 文を使う際の添字の順序による計算時間の差

多次元配列といえども,メモリ上では一次元的に配置されている。do 文で連続的にアクセスする際には,メモリ上を連続的にアクセスした方が当然速い。Fortran の場合,二次元配列 a(i,j) は a(1,1), a(2,1), a(3,1), ..., a(1,2), a(2,2) ... の様にメモリに格納されているので,内側の添字から回す方が速い。なお, Fortran90 以上に限るが,配列に同じ値を代入する際には a = 0 の様な書き方が許されていて,配列 a の全要素への一括代入ができる。こちらの方が do 文による逐次代入よりも速い様である。

!---------------------------------------------------------------------
!     array_i_j.f90 (Fortran 90)
!     SAKAGUTI Syuiti
!     28-Aug-2006
!---------------------------------------------------------------------
      
program array_i_j
  
  implicit none
  integer, parameter :: N=10000
  integer :: i, j, ic, ir,itime0, ia(N,N)
6 format(1a,1i10,1x,"counts")

  ! i -> j
  call system_clock(ic,ir); itime0=ic
  do j=1,N
     do i=1,N
        ia(i,j)=0
     enddo
  enddo
  call system_clock(ic)
  write(*,6) 'i -> j:     ', ic-itime0
  
  ! j -> i
  call system_clock(ic); itime0=ic
  do i=1,N
     do j=1,N
        ia(i,j)=0
     enddo
  enddo
  call system_clock(ic)  
  write(*,6) 'j -> i:     ', ic-itime0
  
  ! vector
  call system_clock(ic);  itime0=ic
  ia=0
  call system_clock(ic)
  write(*,6) 'f90 vector: ', ic-itime0
  
end program array_i_j

上記ソースを Absoft Pro Fortran 90 - MacOSX - PowerMacG5 でコンパイル・実行した結果。このマシンでは 1 s に1,000,000 counts なので,1 counts = 1 μs。明らかな差が分かる。配列記法を使うのが最も速い。

i -> j: 1131083 counts
j -> i: 24387157 counts
f90 vector: 345928 counts

計算精度・数値計算の誤差

整数と実数の有効桁数

◇多数回和を取ったときの誤差の確認 + 定数代入による誤差の確認

◇小さな数の多数回和を取ったときの誤差の確認(積み残し誤差)

非常に小さな数多くの数の和を取る際,大きな方から足した場合と小さな方から足した場合で結果が異なることがある。 これは有効な桁数以上の大きさの差を持つ2 つの実数の加算を行った際,小さい方の値は丸め誤差により反映されないからである。級数を計算する時などは注意。

!---------------------------------------------------------------------
!     numerical_error.f90 (Fortran90)
!     SAKAGUTI Syuiti
!     28-Aug-2006
!---------------------------------------------------------------------

program numerical_error
  implicit none
  integer :: i, j
  integer, parameter :: k=1000, m=1000000
  real(8) :: a, b
  !----------------------------------
  
  ! 10万回以上の和を取った時の丸め誤差

  do j=1,1000
     do i=1,1000000
        a = a + 1.2d-9
     enddo
  enddo
  
  write(*,'(1a,1es24.10)') 'Sum_{k=1}^{1,000,000}, real(8): ', a
  
  
  ! 大きな数に小さい数を加算していった時の積み残し誤差
  
  a = 1d+0
  b = 0d+0
  
  do j=1, m
     a = a + 1.0d-16
     b = b + 1.0d-16
  end do
  
  write(*,*) a
  write(*,*) b + 1d+0
  
end program numerical_error

実行結果

Sum_{k=1}^{1,000,000}, real(8):         1.2000000081E+00
   1.00000000000000
   1.00000000010000

◆階乗計算の精度

倍精度で階乗の自然対数を計算し,配列に格納する。階乗の値は非常に速く大きくなり,20! 程度でも倍精度の有効桁数である 16 桁を超えてしまい,情報の欠落が起こる。これを防ぐためには階乗の自然対数を取ったものを用いて必要な計算を行い,最後に指数を取るのがよい。

!     準備中

複素変数の使用に関するメモ

複素変数の型宣言

基本型複素数 complex は実数部4バイト虚数部4バイトあわせて8バイトである。 では,倍精度型複素数はどうなるのだろうか?倍精度実数型ならば real(8) で,これはどのコンパイラでも共通の様だ。しかし倍精度複素数型に関してはどうもそうではないらしく,コンパイラ毎に定義が違うようで,筆者が今まで見たものには complex(8) の場合と complex(16) の二つの場合があった。確保される記憶領域は実8虚8の16バイトで共通であった。

以下, この記事においては complex(8) は実部・虚部それぞれ real(8) あわせて16バイトの記憶領域を確保するとする。

◇虚数単位の宣言

虚数単位 (imargnal unit) に当たるものを, 次の様に定数として宣言しておくと便利である。

complex(8), parameter :: iu = (0.0d0,1.0d0)

◇複素変数に対する代入操作

定数を代入する際には括弧を用いて書けばよい。しかし,実部と虚部を変数を用いて代入する場合には 実部+虚部*iu の様にする。

complex(8), parameter :: iu = (0.0d0,1.0d0)
complex(8) :: z
real :: a, b

z = (1.0d0,0.0d0)
z = a + b*iu

◇複素配列の割り付け

allocatableで宣言された複素配列を要素数 k 個で割り当てる為には, 2*k 個のサイズで割り当てなくてはならない。

complex(8), allocatable :: z(:)
real(8), allocatable :: a(:)
integer :: k=10
k2=2*k
allocate(a(k),z(k2))

文字列操作に関するメモ

2次元の文字配列は下記の様に宣言する。部分文字列を指定したい時は,a(行の位置)(列の始点位置:列の終点位置)の様に指定する。

character(文字列の長さ) :: a(行数)

文字列を表示するための編集記述子は a 。

!      31-Jan-2007

Program char
  implicit none
  integer :: i
  integer, parameter :: NN=20, MM=10
  character(NN) :: a(MM)='_________0_________0', b='123456789abcdefghijk'

!  write(*,'(10a/)') a
  write(*,'(1a,/)') a
  write(*,'(1a,/)') b

  do i=1, MM
     a(i)=b(i:i+10)
     write(*,'(1i2,": ",1a)') i, a(i)
  end do

  write(*,'(/,1a,1a,/)') 'a(5)(6:8) = ', a(5)(6:8)

  write(*,*) 
  write(*,'(1a)') a

  write(*,*) len(a(1)),  len(a)

end program char

上記コードの出力。format で '1a' を指定しているので,a(i) が表示されている。複数個表示させたい時は,'10a' の様に指定する。

_________0_________0

_________0_________0

_________0_________0

_________0_________0

_________0_________0

_________0_________0

_________0_________0

_________0_________0

_________0_________0

_________0_________0

123456789abcdefghijk

 1: 123456789ab         
 2: 23456789abc         
 3: 3456789abcd         
 4: 456789abcde         
 5: 56789abcdef         
 6: 6789abcdefg         
 7: 789abcdefgh         
 8: 89abcdefghi         
 9: 9abcdefghij         
10: abcdefghijk         

a(5)(6:8) = abc

 
123456789ab         
23456789abc         
3456789abcd         
456789abcde         
56789abcdef         
6789abcdefg         
789abcdefgh         
89abcdefghi         
9abcdefghij         
abcdefghijk         
   20  20

基礎的な習作コード

関数・サブルーチン・一般

◇簡単なプログラムの練習

単純に write 文を使うだけ。C言語でもよくやる,お約束の様なもの。

!---------------------------------------------------------------------
!      hello_fortran.f90 (Fortran90)
!      18-Jan-2003
!---------------------------------------------------------------------

program hello_fortran

  implicit none
  integer :: i

  do i=1, 5
     write(*,*) i, ': Hello FORTRAN'
  end do

end program hello_fortran

◇module 文の練習

module 文を使用して module 変数を宣言し,独立した複数のプロシージャ間で値を共有する練習。

!---------------------------------------------------------------------
!     module_tes_1.f90 (Fortran90)
!     21-Aug-2001
!---------------------------------------------------------------------

module md
  implicit none
  real (8) :: a=3090.0912, b
end module md

program module_tes_1

  use md
  implicit none

  write(*,'(1d24.16)') a
  write(*,'(1d24.16)') b

end program module_tes_1

bashの設定のメモ

ユーザーリソースの設定

bash の場合,ユーザーリソースの設定は ulimit で行う。現在の設定は,ulimit -a で確認することができる。大抵初期設定のままで問題が無い。変更の必要が出てくるのは,core file size と stack size くらいのものだろう。デバッガを使用するなど core が必要なときは ulimit -c (適当な数値などを -c の後に続ける)とすればコードが死ぬときに core を吐く様になる。要らないのであれば ulimit -c 0 とすれば,core は作成されなくなるので,ディスク容量と消去の手間の節約になる。stack size は大きな配列を使用したコードを走らせるときに変更する必要があるかも知れない。

$ ulimit -a
core file size        (blocks, -c) unlimited
data seg size         (kbytes, -d) 6144
file size             (blocks, -f) unlimited
max locked memory     (kbytes, -l) unlimited
max memory size       (kbytes, -m) unlimited
open files                    (-n) 256
pipe size          (512 bytes, -p) 1
stack size            (kbytes, -s) 8192
cpu time             (seconds, -t) unlimited
max user processes            (-u) 100
virtual memory        (kbytes, -v) unlimited

※ MacOSX のコマンドラインでの例。

シェル変数

シェル変数で重要なものは,コマンドサーチパス (PATH) やライブラリ関連の設定である。 また,Fortran のコンパイラ名を各マシンの .bashrc でいちいちシェル変数定義しておくと,マシン毎の Fortran のパスやコマンド名の差異を吸収できる。シェル変数の定義方法はコマンドラインから行う時も,.bashrc などシェル起動時に実行される設定ファイルに書き込む時も同じである。例えば,FC (通常,Fortran コンパイラのコマンド名として使用されるシェル変数。C ならば CC がよく使われる。)を定義する場合,次の様に書く。このようにすると,Makefile 中の ${FC} が /Applications/Absoft/bin/f95 に対応する。

FC=/Applications/Absoft/bin/f95

※ MacOSX のコマンドラインでの例。


コンパイル/ポーティングに関するメモ

コンパイラの種類とオプション

同じことを指定するにも,コンパイラにより形式が異なることがある。

Make の活用

Makefile の例

Makefile_sample.pdf

参考文献

Fortran,数値計算

UNIX 関連の開発用ツール


更新履歴

2008-04-09 コンパイルエラーやバグの対処法のメモを更新(2.4 を追加,その他小修正・小追記数箇所)。日付と時間の処理に関するメモ(年度関数)を追加。
2007-01-31 文字列処理のためのメモを追加。
2006-08-31 章分けを少し変更。入出力に関するメモ,多次元配列に対し do 文を使う際の添字の順序による計算時間の差,計算精度・数値計算の誤差を更新。
2006-08-25 高速化のためのメモ,日付を文字列で返す・実時間を整数スカラーで返すを更新。
2006-08-24 コマンドライン引数の読み込み・文字列変数からの数値読み込み,ファイルの開閉・読み込みと書き出しの練習,連番ファイルの開閉を更新。
2006-08-24 bashの設定のメモ,Makefileの例を追加。
2006-03-20 pv カウンタ設置。
2005-05-31 Fortranの玉手箱を公開。


pv: (scince 2006-03-20)