Fortran の基礎2021年4月4日 | |||||||||||||||||||||||||||||||
Fortran とはいまさら Fortran？Fortran は、死に損ないの古い (ancient) 言語という印象があるが、数値計算の分野ではいまだ現役である。多少知っている人が思い浮かべる典型的な Fortran コードは、Fortran がまだ「生きていた」ころのコードで、それはおそらく FORTRAN 77 であろう。最近のコードはそれとは趣が異なる。 FORTRAN 77 以降には、Fortran 90、95、2003、2008、2018 がある。ふつう、Fortran を学ぶ動機といえば、古いコードを読まざるを得ない場合のために渋々といったところだが、ここでは新しい Fortran にターゲットを絞り、90 から 2008 +α あたりを想定する。 Fortran の利点？数値計算用の言語といったら、現在では Fortran ではなくて MATLAB や NumPy/SciPy モジュールを用いた Python などになる。ただ、この人たちはインタープリタ言語なので、ガチの計算をするには少し遅いという欠点がある。その欠点を克服するものとして出てきたのが Julia であった。 そういう状況で、Fortran をまだかろうじて数値計算用の言語だと呼べなくもないのは、配列が数値計算に特化した仕様になっているところくらいであろうか。言ってしまえば MATLAB っぽい書き方ができる。とはいえ、線形代数ライブラリが標準仕様に含まれるているわけではなく、このあたりは MATLAB などと住み分けられている。また、独特の並列計算の仕組みが言語仕様に追加されている (あまり使われてなさそうだが)。 大抵の言語は数値計算を特に意識していないため、言語仕様が許す範囲でなんとか自然な書き方を試みても結構複雑な事態になりうる。速度と書きやすさを同時に求めるなら、Fortran が一番書きやすいということも (万に一つくらいは) あるかもしれない。 学習コストは、一般的に使われているであろう範囲 (Fortran 90/95 の上澄みくらいか) にとどめておけば、低いほうだと思う。それ以上に頭を突っ込むと、結構厳しくなってくる。Fortran は、利点として、他の言語と比べて単純で余計な機能がないと言われたりする。これは事実とは異なり、実際には奇妙で興味深く不要な機能が色々とある (あった)。このことから察するに、Fortran の秘密は、ユーザーが他の言語どころか Fortran 自体もよく知らないというところにありそうである。すなわち、無知により表現力が単純なレベルに制限されているがゆえの利点がありそうである。したがって、今の Fortran の表現力は C++ にだって劣らないとかいって、下手に Modern Fortran だなんだとやってしまうと、かえって Fortran の利点を損なってしまうかもしれない。 コンパイラフリーで使えるコンパイラとして gfortran がある。ここでは macOS 版 gfortran 8.2.0 で動作を確認し...ていたのだが、Coarray を試そうと思って、OpenCoarrays を Homebrew で入れようとしたら gfortran も入れ替えざるを得なくなって、結局は gfortran 10.2.0 となった。Coarray を試したい場合は、先に OpenCoarrays のセットアップを考えたほうがよいかもしれない。 Ubuntu 18.04 LTS で OpenCoarrays をインストールする場合は、libcoarrays-dev と open-coarrays-bin を入れる。 $ apt-get install libcoarrays-dev open-coarrays-bin プログラムの構造プログラムは program ブロック (Fortran 的にはユニットという) からなる。 program hello print *, "Hello" end program hello program の後には任意の名前をつける。実は program の行自体は必須ではない (プログラムユニットは他から呼べないので)。end program は必要である。end program の後にプログラムの名前を書けるが、program で指定したものと同じものにする。end program は単に end と書いてもよい。仕様にかかわらず、終わりを名前まできっちり書くというお作法があるようだが、それだけ見れば何の終わりかがわかるためということらしい。ただ、そうしないとわかりにくくなっている時点で、すでに何か失敗しているような気がする。 ちなみに、Fortran には大文字と小文字の区別がない。これは、C 言語などと比べると不当な制限に感じるが、接尾辞とかの工夫でなんとかなる。 上記の内容を含むファイル "hello.f90" を用意する。gfortran でコンパイルするには次のようにする。 $ gfortran -o hello hello.f90 実行。 $ ./hello Hello 拡張子の ".f90" は、必ずしも Fortran 90 の意味ではなく、フリーフォーマットの意味のようである。gfortran の場合は、オプション -ffree-form で拡張子とは無関係にフリーフォーマットを指定することができる。 参考 コメント! の後はコメントになる。 ! hello program program hello print *, "Hello" ! say "Hello" end program hello 文文は改行で区切られるが、複数行に渡って文を書きたい場合は、行末に "&" をつければよい。 print *, & "Hello" ! say "Hello" 変数変数の定義変数は、"型 :: 変数1, 変数2, ..." という形でブロック (ユニット) のはじめに定義する。 program example integer :: a, b, c a = 1 b = 2 c = 3 print *, a, b, c end program example この場合、integer のあとの "::" はなくてもよいが、属性というものを設定する場合があり、そのときには必要になる。 初期値を設定することもできる。 integer :: a = 1, b = 2, c = 3 この場合、integer のあとの "::" は必要である。Fortran の場合、C 言語など他の言語と異なり、変数定義時の初期化には特別な (余計な) 意味が付随するので、注意が必要である。初期化は別途行うようにしたほうがよいかもしれない。 変数名などの長さは、昔は最大 6 文字とかだったと思うが、今はそんなことはない。 暗黙の型宣言実のところ、上のコードは変数の宣言なしでも動作する (ただし、出力は少し変わる)。Fortran には暗黙の型宣言といって、変数名の頭文字のアルファベットに応じて型を決める。a, b, c は実数になる。 変数名をうっかり間違っても新しい変数ができるだけなので (Python もそうだが)、結構見つけにくいバグを作り込みかねないという問題があり、"implicit none" と書いて暗黙の型宣言を無効化できるようになっている。 program example implicit none integer :: a, b, c a = 1 b = 2 c = 3 print *, a, b, c end program example 過去のコードを利用するのでなければ、基本的には implicit none を用いた方がよい。JavaScript の "use strict" みたいなものである。gfortran の場合、オプションで -fimplicit-none を指定すれば、デフォルトが implicit none になる。 定数定数 (C 言語でいうところの const) にするには、parameter という属性をつける。 real, parameter :: pi = 3.141592653589793 データ型整数型整数型は integer である。 integer :: a a = 1 これは 4 バイト長である。"integer(kind=8)" のようにバイト長を指定することができる。2、4、8、16 などが指定できる。この kind は種別型パラメータ (kind type parameter) と呼ばれている。ある変数の kind の値が知りたい場合は、組み込み関数の kind を用いる。 print *, kind(a) 使える kind の値は、次のコードで調べられる。 program kinds use iso_fortran_env implicit none print *, "integer ", integer_kinds print *, "real ", real_kinds print *, "character", character_kinds print *, "logical ", logical_kinds end program kinds 出力は以下のようになる。 integer 1 2 4 8 16 real 4 8 10 16 character 1 4 logical 1 2 4 8 16 実数型実数型は real である。 real :: a a = 1. これは単精度である。"real(kind=8)" のようにして精度を指定することができ、4 だと単精度、8 だと倍精度、16 だと 4 倍精度になる。倍精度については、昔ながらの double precision という表現も一応使える。かつて "real*8" といった書き方があった気がするが (実際これでも通るが)、標準的な書き方ではないらしい。 数字の表し方 (リテラル) としては、"12.3" という表記や、"1.23e1" といった指数表記があるが、これらは単精度である。倍精度の場合は "12.3_8"、"1.23e1_8"、"1.23d1"、4 倍精度では "12.3_16"、"1.23e1_16"、"1.23q1" などとする。 整数と実数の変換整数を実数に変換するには、組み込み関数の real を用いる。 print *, real(1)/3 逆に、実数を整数に変換するには、組み込み関数の int を用いる。 print *, int(1.)/3 複素数型複素数型は complex である。 complex :: a a = (1., 2.) これは単精度である。real 同様、kind で精度を指定できる。"(1., 2.)" などの形で実部と虚部をそれぞれ指定する。 文字型と文字列文字列は文字型 character で表す。 character(len=10) :: s s = "1234567890" len で文字列の長さを指定する。昔は "character*10" といった指定ができたと思うが、Fortran 2003 で廃止予定事項になっている。文字列を表す引用符は、シングルクウォート '...' でもよい。 定数の場合は、長さの指定を "*" にできる。 character(len=*), parameter :: s = "hello" 各文字には s(1) などと添字でアクセスできる。また、s(2:3) のような部分文字列へのアクセスもできる。各文字へのアクセスも、s(1) などよりは s(1:1) などとアクセスすることが多いと思う。 文字の長さは len(s) で得られるが、これは変数の長さを返す。上の例では、実際の値に関係なく 10 を返す。本当の文字の長さというか、後ろの空白を除いた長さは、len_trim(s) で得られる。len_trim("omochi") は 6、len_trim(" omochi") は 7、len_trim(" omochi ") は 7 を返す。後ろの空白を除いた文字列自体は trim(s) で得られる。 数値から文字列への変換出力には print や write が使えて、それらは書式を要求する。たとえば、実数 a を出力するには、次のようにする。
print '("a = ", f7.3)', a
上述の例で "*" としていたところが、書式の設定箇所である。上の例では文字列 "a = " の後に実数の出力指定 "f" が指定してあり、幅が 7、小数点以下の桁が 3 としている。整数の場合は "f" の代わりに "i" とする。 write も同様だが、こちらはファイルなどの出力先を指定できる。出力先は文字列変数にもできるので、文字列変数に書き出すことで実数などを文字列に変換できる。 write (s, "(f7.3)") a 文字列から数値への変換入力には read が使える。これはファイルから文字列や数値を読み込むのに使えるが、文字列から数値を読み込むこともできる。 read (s, *) a 論理値型論理値型は logical である。値 (リテラル) は .true. と .false. である。 列挙型C 言語同様の (というかそのものの) 列挙型が使える。 enum, bind(c) enumerator :: sabatora, chatora, mike end enum sabatora、chatora、mike の値は、それぞれ 0、1、2 になる。値を指定することもできる。 enum, bind(c) enumerator :: sabatora, chatora=2, mike end enum この場合、それぞれの値は 0、2、3 となる。 列挙型を型として使いたい場合は、組み込み関数の kind を使う。 integer(kind(sabatora)) type 配列配列は、次のように定義できる。 real :: a(3), b(4) a = [1., 2., 3.] b = [1., 2., 3., 4.] 同じサイズの配列を複数定義するなら、次のようにも書ける。 real, dimension(3) :: a, b 配列の値の表現は、次のような古風な書き方もある。 a = (/1., 2., 3./) 内包表記も可能である。 a = [(2.*i, i = 1, 3)] 多次元配列は、次のように定義できる。 real :: a(2, 2) a = reshape([1., 2., 3., 4.], [2, 2]) C 言語とメモリ上の配置の仕方が異なる。行列で表すと、次のような配置になる。 1.000 3.000 2.000 4.000 次のようにすれば、reshape の配置を変えられる。 a = reshape([1., 2., 3., 4.], [2, 2], order=[2, 1]) これだと、配置は次のようになる。 1.000 2.000 3.000 4.000 ちなみに、上の出力は以下のコードから得られる。 do i = 1, 2 print "(2f7.3)", (a(i, j), j = 1, 2) end do 配列の要素は a(i, j) のように参照するが、添字のはじまりは C 言語などと違って 1 からである。ただし、添字の範囲を変えることもできる。C 言語同様にするには、次のように定義するとよい。 real :: a(0:1, 0:1) ただ、混乱しそうなので、あまりやらないほうがよい気がする。 Fortran の配列の要素は、メモリ上では a(1, 1)、a(2, 1)、... と並んでいることになるので、計算時のループの回し方などには注意が必要なときがある。 単純にすべての要素を 0 で初期化したい場合は、次のようにすればよい。 a = 0. reshape による初期化もそうだが、C 言語の配列などとはかなり趣が異なる。 サイズの取得サイズが必要な時は size 関数を用いる。 print *, size(a) 多次元配列の各次元のサイズが必要であれば、size(a, 1)、size(a, 2) などとする。添字の上限・下限が欲しい場合は、lbound(a, 1)、ubound(a, 1) などを用いる。第 2 引数は次元である。 動的配列配列のサイズを動的に決めることができる。まず、空の変数を用意する。 real, allocatable :: a(:, :) allocate 文で領域を割り当てる。 allocate (a(2, 2), stat=st) ここで、st は整数の変数で、状態を表す数値が入る。0 でなかったら何かしら問題が発生している。 使い終わったら、deallocate 文で領域を解放する。 deallocate (a) 領域が割り当てられているかどうかは、allocated 関数で調べることができる。 すでに割り当てられたもののサイズを変更することはできない (C 言語の realloc 関数に当たるものはない)。代わりに、次のようにする。 real, allocatable :: a_new(:, :) allocate (a_new(3, 3)) a_new(1:2, 1:2) = a(1:2, 1:2) call move_alloc(a_new, a) 別の配列に拡大したサイズを割り当て、それに中身をコピーし、組み込みの move_alloc 関数でメモリを移動する。上の例の場合、a_new の領域は a に渡され、a_new は解放状態になる。また、a がもともと持っていた領域は勝手に解放される。 派生型いわゆるユーザー定義型は、派生型 (derived type) と呼ばれるが、要するに構造体である (だたし、C 言語のよりは C++ のものに近い)。 type :: profile character(len=10) :: name character(len=10) :: type integer :: age end type profile type のあとの "::" は、この場合はなくてもよい (属性を指定するときには必要)。これの変数は、次のように定義する。 type(profile) :: prof 各フィールドにアクセスするには、% を使う。 prof%name = "omochi" prof%type = "sabatora" prof%age = 4 フィールドへの値の設定は、次のようにもできる。
prof = profile("omochi", "sabatora", 4)
キーワードを指定することもできる。 prof = profile(name="omochi", type="sabatora", age=4) デフォルトの値を設定することもできる。 type :: profile character(len=10) :: name = "omochi" character(len=10) :: type = "sabatora" integer :: age = 4 end type profile len や kind などのパラメータを持たせることもできる。 type profile(l, k) integer, len :: l = 10 integer, kind :: k = 4 character(len=l) :: name character(len=l) :: type integer(kind=k) :: age end type profile この場合の変数の定義は、次のようにする。 type(profile(len=20, kind=8)) :: prof パラメータを指定しなければ、デフォルト値が用いられる。ここでの l, k もフィールドなので、値の設定は次のようにする。 prof = profile(10, 4, "omochi", "sabatora", 4) l, k の値は代入先と同じである必要がある。l, k の値がデフォルトでよい場合は、キーワードを指定する方式を使えばよい。 メモリの動的確保上で動的配列の例を見たが、一般の型についてもメモリを動的に確保することができる。 type(profile), allocatable :: prof allocate(prof) ... deallocate(prof) また、Fortran にも C 言語のようなポインタがある。 type(profile), pointer :: prof allocate(prof) ... deallocate(prof) allocatable と pointer の違いは、pointer はポインタなので、別のものを指せる。ただし、指せるのはポインタか、target 指定されたものである。 type(profile), target :: prof type(profile), pointer :: p p => prof ポインタが何かを指しているかどうかは、associated 関数で調べることができる。 メモリまわりは、C 言語に比べると複雑に見えることがある。上述の動的配列のリサイズの例で示した move_alloc 関数の場合、移動元の変数がポインタと結びついていると、ポインタの参照先も移動する (C 言語的には参照先が「移動していない」と言いたくなるが、Fortran の場合はアドレスではなくて変数を参照している)。ただ、移動先が target 属性を持たない場合は未定義のようである (gfortran の場合、問題なく移動される。つまり、target 属性を持たない変数を指すポインタになる)。もともと移動先を指しているポインタがあった場合は、領域は残るようである (これが標準かどうかはわからない)。 Fortran は、領域の割り当て状態やポインタとの関連状態を言語側が把握しているため、メモリ関連はどうもうまいこと処理をしてくれるようである。参照カウンタを持っているのかもしれない (スマートポインタのイメージ？)。後の話になるが、モジュールの final 手続きが常に呼ばれるわけではないのも、そういうことかもしれない。正直よくわからない。 参考 式と演算子文は式からなるが、式はデータと演算子で構成される。 数値については、以下のような演算子が使える。
余りを求める演算子はなくて、関数 mod(a, b) を用いる。modulo(a, b) という似たような関数もあるが、これは引数が負のときに値が異なる。C 言語の演算子 "%" と同じ挙動をするのは mod() のほうである。 文字列の連結用に "//" という演算子が使える。 制御構造block 構文block...end block でブロックを作れる。 block real :: a = 1. print *, a end block print *, a ! エラー: a は見えない ブロック内で定義された変数は、ブロックの外からは参照できない。ブロックは入れ子にすることもできる。 ブロックには名前をつけられる。 b: block real :: a = 1. print *, a end block b 名前をつけておくと、exit 文でブロックから抜けることができる。
b: block
real :: a = 1.
block
exit b
end block
print *, a ! 実行されない
end block b
associate 構文associate 構文は、たとえば、以下のように書く。
program example
implicit none
real :: x = 1., y = 2.
associate (z => sqrt(x**2 + y**2))
print *, z, z + 2., z - 2.
end associate
end program example
変数 z は定義されていないことに注意。sqrt(x**2 + y**2) を z として名前付けしている。マクロのような使い方と言えるかもしれない。ただ、associate 構文の中で x、y の値を書き換えても z の値は変わらない。だから、この場合は変数 z を作って代入したのと結果は同じになる。 associate 構文の外では z の値は参照できないので、z の値が欲しい場合は、他の変数にも関連付けしてそれに値を入れる。
program example
implicit none
real :: x = 1., y = 2., v = 0.
associate (z => sqrt(x**2 + y**2), r => v)
print *, z, z + 2., z - 2.
r = z
end associate
print *, v
end program example
r の値を変えると v の値が変わる (ややこしいが、逆もしかり)。これはポインタというか参照っぽい。式と関連付けるか変数と関連付けるかでイメージが少し違うのである。 派生型の各フィールドに関連付けることもできる。
program example
implicit none
type :: profile
character(len=10) :: name
character(len=10) :: type
integer :: age
end type profile
type(profile) :: prof = profile("omochi", "sabatora", 4)
associate (name => prof%name, type => prof%type, age => prof%age)
print *, name, type, age
age = age + 1
end associate
print *, prof
end program example
値も変更できる。Pascal の With 文 (JavaScript にもあるけど、ないことにされてるやつ) みたいに使えなくもない。 比較演算子<、>、<=、>= などは C 言語と同じである。等しいを表すのは ==、等しくないは /= である (これだけ C 言語と異なる)。昔は .lt.、.gt.、.le.、.ge.、.eq.、.ne. を使ったが (今でも使えるが)、今更使う必要はないだろう。 論理演算子.and. (AND)、.or. (OR)、.not. (NOT) が使える。等しい、等しくないはそれぞれ .eqv.、.neqv. を用いる。==、/= でもよさそうなものだが、論理型には使えない。 if 構文if...end if で条件分岐を作れる。 if (v == 1) then print *, "sabatora" else if (v == 2) then print *, "chatora" else print *, "mike" end if if 構文には名前を付けられる。 select_type: if (v == 1) then print *, "sabatora" else if (v == 2) then print *, "chatora" else select_type print *, "mike" end if select_type 名前を付ける理由は、たぶん構造を明確にするためだと思うが、Fortran のコードは昔からだらだらと長く書かれがちだから、その習慣に対処するものかもしれない。 if 文条件が 1 つなら、if 文が使える。 if (v == 1) print *, "sabatora" あと、算術 if 文という変な機能があったと思うが、Fortran 2003 では廃止予定事項になっている。 case 構文select case...end select で case 構造を作れる。要するに C 言語の switch 文である。 select case (v) case (1) print *, "sabatora" case (2) print *, "chatora" case default print *, "mike" end select 値は文字列でもよい。また、整数の場合、合致する値を複数指定できる。 select case (v) case (1,3,5) print *, "sabatora" case (2,4,6) print *, "chatora" end select 範囲にもできる。 select case (v) case (1:3) print *, "sabatora" case (4:6) print *, "chatora" end select case 構造には名前を付けられる。 select_type: select case(v) case (1) print *, "sabatora" case (2) print *, "chatora" case default select_type print *, "mike" end select select_type do 構文do...end do でループを作ることができる。以下のようにすると、無限ループを作れる。 do print *, "xxx" end do exit 文でループを抜けることができる (C 言語の break にあたる)。 i = 1 do print *, i if (i == 10) exit i = i + 1 end do 1 から 10 まで変ぞえるなら、次のように書ける。 do i = 1, 10 print *, i end do 増分を指定できるので、逆さまに数えることもできる。 do i = 10, 1, -1 print *, i end do do while でいわゆる while 文を書くこともできる。 i = 1 do while (i <= 10) print *, i i = i + 1 end do C 言語の do...while のようにループ条件を後置にはできないようで、そういうときは do...end do を使って if 文で抜ける。 cycle 文でループの次のステップに移行できる (C 言語の continue にあたる)。 do i = 1, 10 if (mod(i, 2) == 0) cycle print *, i end do do 構文には名前を付けられる。
loop: do i = 1, 10
do j = 1, 10
print *, i, j
if (i == 2 .and. j == 2) exit loop
end do
end do loop
do 構文の名前付けには特別な意味がある。上の例のように名前を指定して exit すると、入れ子になったループを一気に抜けられる。これがなかったら goto 文を使う言い訳をひとつ与えてしまうので、いにしえの Fortran ユーザーを成仏させるには重要な事項である。 goto 文考古学的興味のため、goto 文も一応見ておく。 goto 12345 ... 12345 continue 頭にラベルのついた文に飛ぶ。ラベルは最大 5 桁の数字である。ラベルにアルファベットが使えないのが C 言語よりも不便だが、そもそも goto 文自体使うべきではない。 といいつつ、C 言語的発想でいくと、例外処理に使い道が見出せるかもしれない。 if (st /= 0) goto 99999 ... 99999 print *, "allocation error" stop ただ、ラベルが数字しかダメなので、読みにくくなるだけかもしれない。 goto は、もっと変な飛び方をするやつがあったが (計算型 goto)、Fortran 2003 で廃止予定事項に入れられている。 プログラムの停止プログラムを止めるには、stop 文を使う (C 言語の exit 関数にあたる)。 if (st /= 0) then print *, "allocation error" stop end if error stop 文で、エラーメッセージを表示させて止めることもできる。 if (st /= 0) then error stop "allocation error" end if この場合、終了コードは 1 になる。 メッセージの代わりに終了コードを設定することもできる。 if (st /= 0) then error stop 10 end if 手続き手続きとして関数とサブルーチンがある。値を返す必要がある場合 (あるいは返す値が 1 つで済む場合) は関数、値を返す必要がない場合 (あるいは返す値が 1 つでは済まない場合) はサブルーチンが使われる。 関数関数は、function...end function で定義できる。たとえば、引数を 2 倍した値を返す関数 double は、以下のようになる。 function double(x) implicit none real :: x real :: double double = 2.*x end function double end function の仕様は end program と同様である。関数名と同名の変数を定義しており、それに値を入れて返り値を設定している。返り値用の変数を指定することもできる。 function double(x) result(r) implicit none real :: x real :: r r = 2.*x end function double 基本的には関数名で返り値を指定すればよいと思うが、たまに関数の名前を変えようとしたときにめんどくさかったりするので、常に名前を付けるようにするとか、返り値の指定は必ず関数の最後にするとかしておいたほうがよい気がする。 返り値の型は、function の前に指定してもよい。 real function double(x) implicit none real :: x double = 2.*x end function double return 文で関数を終了することができる。ただし、C 言語のようにそれで値を返したりはできない。昔は return 文で変な戻り方ができたと思うが (選択 return 文)、Fortran 2003 で廃止予定事項に入れられている。 呼び出し側は次のようになる。 program example implicit none real :: double print *, double(2.) end program example 呼び出すには関数の型を指定する必要がある。次のように interface 文で関数のインターフェイスをフルセットで書く方法もある。
program example
implicit none
interface
function double(x)
real :: x
real :: double
end function double
end interface
print *, double(2.)
end program example
ところで、昔は 1 行で書く関数定義の方法があったと思うが (文関数)、ふつうの代入文に見えかねないということで、Fortran 2003 で廃止予定事項に入れられている。 関数の定義の場所関数の定義は、program...end program と並列するように定義できるが、program...end program の中に定義することもできる。
program example
implicit none
print *, double(2.)
contains
function double(x)
real :: x
real :: double
double = 2.*x
end function double
end program example
contains 文の後で定義する。この場合、呼び出し側の設定は必要ない。 また、同様にして関数内に関数を定義することもできる。 引数の種類の指定引数の定義に "intent(in)" を付けると、入力値であることを明確にできる。 function double(x) implicit none real, intent(in) :: x real :: double double = 2.*x end function double Fortran の関数の引数は (C 言語などと違って) 参照渡しなので、関数の中で引数 (仮引数) の値を変更すると、呼び出し元の引数 (実引数) の値も変わってしまう。それを意図しない場合、intent(in) という属性を指定する。意図してそうしたい場合は、intent(out) や intent(inout) を指定する。 引数の値渡し引数を値渡しにすることもできる。C 言語とやりとりするためのものだろう。
program example
implicit none
real :: x = 2.
print *, double(x)
print *, x
contains
function double(x)
real, value :: x
real :: double
double = 2.*x
x = 1. ! 影響なし
end function double
end program example
double 関数の引数 x に "value" という属性が指定されている。関数内で x に値を代入しているが、値渡しになっているため、呼び出し側の値は変わらない。ここで、program 内部で関数を定義しているのは、呼び出し側で interface を要求されるためである。 キーワード引数引数は、キーワードを指定して値を指定することもできる。たとえば、2 つの値を掛け合わせる mul(x, y) という関数があったとすると、次のように呼び出すことができる。 print *, mul(2., 3.) print *, mul(x=2., y=3.) print *, mul(y=3., x=2.) print *, mul(2., y=3.) print *, mul(2., x=3.) ! エラー: x はすでに指定してある オプション引数C++ のデフォルト引数のように、引数を指定してもしなくてもよいオプション引数というものが使える。
program example
implicit none
print *, mul(2., 3.)
print *, mul(2.)
contains
function mul(x, y)
real, intent(in) :: x
real, optional, intent(in) :: y
real :: y_
real :: mul
if (present(y)) then
y_ = y
else
y_ = 1.
end if
mul = x*y_
end function mul
end program example
interface を要求されるため、program 内関数として定義している。引数 y について optional が指定してある。y のデフォルト値というのはなくて、present(y) で値が指定してあるかどうかを調べる。指定してあったら、別の変数 y_ に y の値を入れる (y に値を代入できないため)。値が指定してなければデフォルトの値を y_ に入れる。そして y の代わりに y_ を使って計算を行う。 変数の初期化と値の保持関数内でデータを保持するために、C 言語でいう static 変数のようなものを、Fortran でも作ることができる。 function count() implicit none integer, save :: x = 0 ! save integer :: count x = x + 1 count = x end function count "save" で変数 x の値を保持するように指定している。これにより、関数 count を呼び出すたびに、返り値が 1、2、... と大きくなっていく。 ややこしいことに、save を付けなくても同じようになる。 function count() implicit none integer :: x = 0 ! save integer :: count x = x + 1 count = x end function count 変数の定義時に初期化すると、その変数は save 変数だと見なされるようになっている。program の中ならよいが、手続きの中だと意図しない副作用を作り込みかねない。 そういえば、初期化はもともと data 文でやってた気がする。一見ただの初期化に見えるが、data 文による初期化の省略形と考えたほうがよさそうである。 再帰関数再帰的な関数を作る場合は、"recursive" を付ける必要がある。
program example
implicit none
integer :: fibonacci
print *, fibonacci(10)
end program example
recursive function fibonacci(n) result(r)
implicit none
integer, intent(in) :: n
integer :: r
select case (n)
case (0)
r = 0
case (1)
r = 1
case default
r = fibonacci(n - 2) + fibonacci(n - 1)
end select
end function fibonacci
純粋関数副作用がないことを示すために、"pure" を指定できる。 pure function double(x) implicit none real, intent(in) :: x real :: double double = 2.*x end function double これは、副作用のない手続きが必要な場合に使うものだが、副作用のない手続きを作っているつもりのときは、常に指定した方がよいかもしれない。意図しない副作用をうっかり作り込んだときに気が付くので。 pure 手続きの中では、print 文などを実行したり、他の pure でない手続きを呼んだりできない。 要素別処理関数"elemental" を指定すると、配列に対して処理を要素ごとに適用する関数を定義できる。これは引数がスカラーである (配列でない) 必要がある。
program example
implicit none
real :: a(3) = [1., 3., 5.]
a = double(a)
print *, a
contains
elemental function double(x)
real, intent(in) :: x
real :: double
double = 2.*x
end function double
end program example
interface を要求されるため、program の内部関数としている。引数は配列で定義されていないのに、呼び出し側では配列を指定して、さらに返り値を配列で受けている。この場合、[1., 3., 5.] の要素がそれぞれ 2 倍されて [2., 6., 10.] が格納される。 サブルーチンサブルーチンは subroutine...end subroutine で定義する。値を返さないという点以外は、関数とほぼ同じである。ただし、呼び出し時には call を付ける。たとえば、サブルーチン sub を呼び出すには、次のようにする。 call sub() 実は引数がないときは括弧はなくてもよい。 文字列を引数にする文字列を引数で渡す場合は、長さを指定する必要はない。
program example
implicit none
call print_string("hello")
contains
subroutine print_string(s)
character(len=*), intent(in) :: s
print *, s
end subroutine print_string
end program example
サブルーチン print_string の引数 s では、文字の長さを "len=*" としている。 配列を引数にする配列を引数にする際、昔は配列のサイズも一緒に渡していたが、最近の Fortran では必ずしもそうする必要はない。
program example
implicit none
real :: a(5) = [1., 2., 3., 4., 5.]
print *, total(a)
contains
function total(a)
real, intent(in) :: a(:)
integer :: i
real :: total
total = 0.
do i = lbound(a, 1), ubound(a, 1)
total = total + a(i)
end do
end function total
end program example
配列の要素の総和を返す関数 total を定義している (実際には組み込み関数の sum を使うと思うが)。引数は a(:) となっていて、サイズが指定されていない。ループを回す際に、lbound と ubound で添字の上限と下限を取得している (第 2 引数は添字を取る次元の指定)。下限が 1 だと断定できるのなら、size 関数を用いることもできる。 do i = 1, size(a) 2 次元配列の場合は、size の第 2 引数に ubound などと同様に次元を指定できる。 手続きを引数にする手続きの引数に関数などを渡して呼び出してもらうのは、interface を設定してやればできる (前からできた気もするが)。
program example
implicit none
real :: a(3) = [1., 3., 5.]
call map(a, double)
print *, a
contains
function double(x)
real, intent(in) :: x
real :: double
double = 2.*x
end function double
subroutine map(a, f)
real, intent(inout) :: a(:)
interface
function f(x)
real, intent(in) :: x
real :: f
end function f
end interface
integer :: i
do i = lbound(a, 1), ubound(a, 1)
a(i) = f(a(i))
end do
end subroutine map
end program example
上の例では、配列の要素に関数を適用して書き換えるサブルーチン map を定義している。interface で渡せる関数のインターフェイスを定義している。一見めんどくさそうだが、C 言語よりはわかりやすくてよいと思う。 次のような書き方もできる。
program example
implicit none
real :: a(3) = [1., 3., 5.]
interface
function f_int(x)
real, intent(in) :: x
real :: f_int
end function f_int
end interface
call map(a, double)
print *, a
contains
function double(x)
real, intent(in) :: x
real :: double
double = 2.*x
end function double
subroutine map(a, f)
real, intent(inout) :: a(:)
procedure(f_int) :: f
integer :: i
do i = lbound(a, 1), ubound(a, 1)
a(i) = f(a(i))
end do
end subroutine map
end program example
引数 f の型として procedure(f_int) でインターフェイスを指定してる。 common ブロックとかcommon ブロックは、いわゆるグローバル変数を実現する仕組みである。equivalence とかいう変な機能もあった気がする。こういうのはもういらないんじゃないかと思う。グローバル変数が必要な場合は、Fortran 90 以降ではモジュールを使えばよい。 entry とかいう気違いじみた機能もあったが、Fortran 2008 でようやく廃止予定になったらしい。 モジュールmodule...end module でモジュールを作ることができる。モジュールというのは、データや手続きのセットである。C 言語でいう分割コンパイルの高度なやつである。 モジュールのコンパイルたとえば、次のようなファイルを用意する。 profile.f90
module profile_m
implicit none
type :: profile_t
character(len=10) :: name
character(len=10) :: type
integer :: age
end type profile_t
end module profile_m
example.f90
program example
use profile_m
implicit none
type(profile_t) :: prof
prof = profile_t("omochi", "sabatora", 4)
print *, prof
end program example
別々にコンパイルする。 $ gfortran -c profile.f90 $ gfortran -c example.f90 $ gfortran -o example example.o profile.o モジュールをコンパイルすると .mod ファイルが生成されるが、それは利用する側が必要とする。もし参照するモジュールのディレクトリを別にしたい場合は、gfortran では -I オプション (gcc のインクルードパスの指定と同じ) でパスを指定すればよい。モジュールのコンパイル時に .mod ファイルを出力するディレクトリを変更したい場合は、-J オプションで指定する。 モジュールの定義モジュールは次のように定義する。
module profile_m
implicit none
private
type :: profile_t
character(len=10) :: name
character(len=10) :: type
integer :: age
end type profile_t
public :: profile_t, set_name, set_type, set_age
contains
subroutine set_name(prof, name)
type(profile_t), intent(inout) :: prof
character(len=*), intent(in) :: name
prof%name = name
end subroutine set_name
subroutine set_type(prof, type)
type(profile_t), intent(inout) :: prof
character(len=*), intent(in) :: type
prof%type = type
end subroutine set_type
subroutine set_age(prof, age)
type(profile_t), intent(inout) :: prof
integer, intent(in) :: age
prof%age = age
end subroutine set_age
end module profile_m
3 行目の private は、デフォルトの参照許可属性 (要するにアクセス制御) の設定である。private だとデフォルトで非公開、public だとデフォルトで公開になる。何も指定しない場合は public 扱いになる。目的にもよるが、すべて公開するつもりでない限り、デフォルトは private にして、公開したい変数や手続きを個別に public 設定にしていくのがよい。上の例では、contains の前の行の public 文により、派生型と、contains 以下で定義されているサブルーチンについて、個別に公開設定にしている。 変数を定義して使うには次のようにする。
module profile_m
implicit none
private
...
type(profile_t), public :: omochi = profile_t("omochi", "sabatora", 4)
...
contains
...
上の例では omochi という変数を定義しており、public 属性が設定してある。これは、モジュール内手続きから触ることができるだけでなく、モジュールを参照しているすべての手続きからも触ることができるので、実質グローバル変数である。値は参照させたいが、値を変更させたくない場合は、protected 属性を付ける。
type(profile_t), public, protected :: &
omochi = profile_t("omochi", "sabatora", 4)
parameter とどう違うのだろう？ と思ったが、protected はモジュール外で変更できないだけで、モジュール内手続きでは変更できる。定数を提供したい場合ではなく、覗き窓を提供したい場合に使うものだろう。 モジュールの使用モジュールは "use" でモジュール名を指定して使うことができる。 use profile_m JavaScript や Python の import のようなものと言えばよいか (ちなみに、Fortran にも import があって、別の用途に割り当てられている)。 すべての項目を使うわけではない場合は、使うものだけを指定できる。 use profile_m, only : profile_t 取り込む項目の名前を変えることもできる。 use profile_m, only : pdata => profile_t 上の例では、profile_t を pdata という名前で参照するようにしている。 すべての項目を取り込むが、一部だけ名前を変える、ということもできる。 use profile_m, pdata => profile_t サブモジュールサブモジュールは、たぶんインターフェイスと実装を分けるためのものである。上の profile_m をサブモジュールを使って書き直すと、次のようになる。
module profile_m
implicit none
private
type :: profile_t
character(len=10) :: name
character(len=10) :: type
integer :: age
end type profile_t
interface
module subroutine set_name(p, name)
type(profile_t), intent(inout) :: p
character(*), intent(in) :: name
end subroutine set_name
module subroutine set_type(p, type)
type(profile_t), intent(inout) :: p
character(*), intent(in) :: type
end subroutine set_type
module subroutine set_age(p, age)
type(profile_t), intent(inout) :: p
integer, intent(in) :: age
end subroutine set_age
end interface
public :: profile_t, set_name, set_type, set_age
end module profile_m
submodule (profile_m) profile_sm
contains
module procedure set_name
p%name = name
end procedure set_name
module procedure set_type
p%type = type
end procedure set_type
module procedure set_age
p%age = age
end procedure set_age
end submodule profile_sm
かえってわかりにくい気がする。 配列基本操作Fortran の配列は、添字の始まりがデフォルトでは 1 から始まるとか、メモリ上の配列が前の添字のほうから動くとかはあるが、基本的な部分は C 言語と同様である。
program example
implicit none
real :: a(3, 3), v
integer :: i, j
v = 1.
do i = 1, ubound(a, 1)
do j = 1, ubound(a, 2)
a(i, j) = v
v = v + 1.
end do
end do
do i = 1, ubound(a, 1)
do j = 1, ubound(a, 2)
write (*, "(f7.3)", advance="no") a(i, j)
end do
print *
end do
end program example
出力 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000 9.000 ただし、Fortran の配列は、MATLAB に似た操作ができるようになっている。値の初期化部分は、次のように書ける。 a = reshape([(real(i), i = 1, size(a))], shape(a), order=[2, 1]) 出力部分は、次のように書ける。
character(len=10) :: fmt
...
write (fmt, '("(", i1, "f7.3", ")")') ubound(a, 2) ! 書式の作成
do i = 1, ubound(a, 1)
print fmt, (a(i, j), j = 1, ubound(a, 2))
end do
あるいは、もっと単純に、次のようにも書ける。 print fmt, (a(i, :), i = 1, ubound(a, 1)) a(i, :) で a(i, 1)、a(i, 2)、a(i, 3) を表している。 転置を使う方法もある。 print fmt, transpose(a) 演算2 倍したいだけなら、そのまま 2 をかければよい。 a = 2.*a 同じ形状の配列であれば、足したりできる。 real, dimension(3, 3) :: a, b ... a = a + b 同じ形状の配列どうしをかけた場合、要素同士がそれぞれかけ合わされる。 a = a*b 行列の積を実行したい場合は、matmul という関数を用いる。 a = matmul(a, b) ベクトルの内積が欲しい場合は、dot_product という関数がある。 real, dimension(3) :: x, y ... print *, dot_product(x, y) ベクトルのノルムは norm2 で計算できる。 print *, norm2(x), sqrt(dot_product(x, x)) 行列の転置には transpose 関数がある。 a = transpose(a) 要素の参照配列の全要素に同じ値を入れるには、そのまま代入すればよい。 a = 0. これは次のように書くのと同じである。 a(:, :) = 0. a(1, :) とすると、1 行目にアクセスできる。 a(1, :) = 0. a(2:3, :) とすると、2 行目から 3 行目にアクセスできる。 a(2:3, :) = 0. a(:, 1:3:2) と書くと、1 から 3 列目に増分 2 で (つまり 1 列目と 3 列目に)アクセスする。 a(:, 1:3:2) = 0. 特殊な構文forall 構文たとえば、行列の非対角要素を 0 にしたい場合は、次のようにできる。
do i = 1, ubound(a, 1)
do j = 1, ubound(a, 2)
if(i /= j) a(i, j) = 0.
end do
end do
これは、forall を使うと 1 行で書ける。 forall (i = 1:ubound(a, 1), j = 1:ubound(a, 2), i /= j) a(i, j) = 0. ブロックとして書くこともできる。 forall (i = 1:ubound(a, 1), j = 1:ubound(a, 2), i /= j) a(i, j) = 0. end forall 入れ子にしてもよい。
forall (i = 1:ubound(a, 1))
forall (j = 1:ubound(a, 2), i /= j)
a(i, j) = 0.
end forall
end forall
この処理については、どうも入れ子にしないと処理が遅いようである。 forall は、上の例のように、左辺に forall で指定している添字を使用した配列があるパターンで書く。ブロック内で do ループは使えない。 do concurrent 構文forall に似た感じのもので、do concurrent というものがある。
do concurrent (i = 1:ubound(a, 1))
do concurrent (j = 1:ubound(a, 2), i /= j)
a(i, j) = 0.
end do
end do
これは並列計算、というよりベクトル化を想定したもののようで、処理がそれぞれ独立で、どんな順番で処理されてもよいものと想定されているようである。gfortran については、do concurrent で書ける処理をふつうの do ループで書いたとしてもたぶん最適化が効くので、たいして差はなさそうである。do concurrent を使うようにしたほうが、ふつうの do で書くよりも最適化の機会が増えるというのはあるかもしれない。コンパイラやオプションによって、forall と性能が異なると考えた方がよさそうである。 forall と違い、ブロック内で do ループを使うことができる。一見、これがあれば forall は不要なように思えるが、処理が独立で順番関係ないという条件が厳しいときに forall が必要になるのかもしれない。 where 構文たとえば、行列の負の要素を 0 にしたい場合は、次のようにできる。
do i = 1, ubound(a, 1)
do j = 1, ubound(a, 2)
if(a(i, j) < 0.) a(i, j) = 0.
end do
end do
これは、forall で次のように書ける。 forall (i = 1:ubound(a, 1), j = 1:ubound(a, 2), a(i, j) < 0.) a(i, j) = 0. この場合、where 構文を使うと、次のように簡単になる。 where (a < 0.) a = 0. end where あるいは、負でない要素には 1 を入れたいといった場合は、次のようにできる。 where (a < 0.) a = 0. else where a = 1. end where さらに、要素が負、0、正のときにそれぞれ 1、2、3 を入れたいといった場合は、次のように書ける。 where (a < 0.) a = 1. else where (a == 0.) a = 2. else where a = 3. end where 関数の適用forall で関数を適用する場合は、pure 関数でないといけない。また、where で関数を適用する場合は、elemental 関数である必要がある。forall、where ではサブルーチンは使えない。do concurrent は pure な手続きであれば大丈夫なようである。 その他の組み込み関数sum 関数で要素の和をとれる。 a = [(real(i), i = 1, 5)] print *, sum(a) product 関数で要素の積をとれる。 print *, product(a) sum や product はマスクを設定できる。 print *, sum(a, mask=a > 2.) maxval、minval で要素の最大値、最小値を得られる。また、maxloc、minloc で最大値、最小値それぞれの添字を得られる。 print *, maxval(a), maxloc(a) count で条件にあう要素の個数を得られる。 print *, count(a > 2.) all 関数は、すべての要素が条件を満たす場合に真を返す。any 関数は、少なくとも 1 つは条件を満たす要素がある場合に真を返す。 print *, all(a > 2.), any(a > 2.) 計算時間の計測cpu_time 関数で計算時間を計測できる。 real :: start, finish call cpu_time(start) ... call cpu_time(finish) print '(f7.3, " s")', finish - start 単純な場合はこれでよいのだが、マルチスレッドだと正しく計測できない。たとえば、ループのある処理を計測するとして、gfortran にオプション "-O3 -floop-parallelize-all -ftree-parallelize-loops=2" などを渡してみるとよい。 代わりに system_clock 関数を用いる。 integer :: start, finish, time_rate call system_clock(start, time_rate) ... call system_clock(finish, time_rate) print '(f7.3, " s")', real(finish - start)/time_rate 以下の例は、行列の積の実装の比較である。
program example
implicit none
integer, parameter :: n = 1000, ntries = 5
integer :: start, finish, time_rate
real :: time
real, dimension(n, n) :: a, b, c
integer :: i, j, k, t
call random_number(a)
call random_number(b)
time = 0.
do t = 1, ntries
call system_clock(start, time_rate)
do i = 1, ubound(a, 1)
do j = 1, ubound(b, 2)
c(i, j) = 0.
do k = 1, ubound(a, 2)
c(i, j) = c(i, j) + a(i, k)*b(k, j)
end do
end do
end do
call system_clock(finish, time_rate)
time = time + real(finish - start)/time_rate
end do
print '(f7.3, " s")', time/ntries
time = 0.
do t = 1, ntries
call system_clock(start, time_rate)
c = matmul(a, b)
call system_clock(finish, time_rate)
time = time + real(finish - start)/time_rate
end do
print '(f7.3, " s")', time/ntries
end program example
gfortan で -O2 オプション付きでコンパイルして、以下の結果になった。 1.133 s 0.074 s forall や do concurrent も試してみたが、最初、入れ子にしないで書いたら遅かった。入れ子にすれば do と変わらない。matmul はごりっごりにチューニングしてあるようで、ループではまったく歯が立たない。 ファイルファイルの新規作成ファイルを開いて書き込むには、次のように書く。 integer :: uni, sta open (newunit=uni, file="test", status="new", iostat=sta) if (sta /= 0) then error stop "file open error" end if write (uni, *) "test" close (uni) open 文でファイルを開く。ファイルを開くにはユニークなユニット番号を指定する必要があるが、さすがに最近は自動で付番してくれる。引数 file でファイル名を指定する。status="new" だと新規作成になる。ただし、ファイルが存在するとエラーになる。上書きでよければ status="replace" にする。iostat で指定した変数には、問題がなければ 0 が、問題があれば正の数が入る。 write 文では、最初にユニット番号を指定すれば、ファイルに書き出しができる。close 文では、ユニット番号を指定してファイルを閉じる。 ファイルの読み込みファイルを開いて読み込むには、次のように書く。 integer :: uni, sta character(len=128) :: line open (newunit=uni, file="test", status="old", iostat=sta) if (sta /= 0) then error stop "file open error" end if read (uni, *) line print *, trim(line) close (uni) 既存のファイルを開く場合は、open 文で status="old" とする。読み込みは、read 文でユニット番号を指定すればよい。ファイルの内容がなくなるまで 1 行ずつ読み込むには、次のようにする。 do read (uni, "(a)", iostat=sta) line if (sta /= 0) exit print *, trim(line) end do ファイルの追記既存のファイルに追記するには、status="old" と一緒に position="append" を指定する。 open(newunit=uni, file="test", status="old", position="append", iostat=sta) 空白区切やカンマ区切のデータの読み込み空白区切やカンマ区切のデータの並びを読み込むのであれば、変数を並べるだけでよい。 real :: a, b ... read (uni, *) a, b ただ、これだと改行もお構いなしにデータを次々と読み込んでいく。できれば 1 行単位で処理したいとなると、文字列で一度読み込む必要がある。
character(len=128) :: line
character(len=128), allocatable :: items(:)
integer :: cnt
integer :: i
...
read (uni, "(a)") line
cnt = 0
do i = 1, len_trim(line)
if (line(i:i) == ",") then
cnt = cnt + 1
end if
end do
allocate (items(cnt))
read (line, *) items
...
ここで、line(i:i) は line(i) としてしまいそうになるが、これだとうまく行かない。 固定フォーマットの読み込み昔ながらの Fortran プログラムだと、入力ファイルが固定フォーマットになっていたりする。たとえば、次のような内容のファイルがあったとする。 12345678901234567890 これが 10 文字ごとの 2 つの値を表しているとすると、次のように読み込む。 integer :: a, b ... read (uni, "(2i10)") a, b あるいは、実数だと次のようになる。 real(kind=8) :: a, b ... read (uni, "(2f10.0)") a, b namelistFortran には namelist という変なファイルフォーマットがある。たとえば、次のように書く。 &material young = 2.1e9, nu = 0.3 / 1 行目はグループの名前、その後が変数とその値、最後はグループの終わりである。コードのほうでは、次のように書く。 real(kind=8) :: young, nu namelist /material/ young, nu ... read (uni, nml=material) これで、変数に指定した値がそのまま設定される。namelist ファイルには、"!" によるコメントも書ける。 write 文で namelist を書き出すこともできる (大文字になるが)。namelist は、派生型も使えるが、private なメンバーを持つものは使えない。また、allocatable や pointer 属性を持つものも使えない。 標準エラー出力write 文でユニット番号指定を "*" にすると標準出力になるが、標準エラー出力で書き出したいときがある。その場合は、次のようにする。 use iso_fortran_env, only : error_unit ... write(error_unit, '(a)') msg iso_fortran_env モジュールで定義されている error_unit が標準エラー出力のユニット番号である (昔から特定の番号だった気がするが、忘れた)。 書式formatprint 文の * や、write/read 文の二つ目の * には、書式 (format) を設定できる。C 言語の printf 文の "%7.3f" みたいなやつである。昔は次のように書いた。
write (*, 10) a
10 format (f7.3)
今となっては、ラベル自体が悪しき Fortran の象徴のように見えるし、そもそも分離されてわかりにくいので、直接書く。 write (*, "(f7.3)") a "(...)" の形で書けばよい。 データ整数は i、実数は f、文字列は a で表す。
print '("i = ", i5)', i
print '("a = ", f7.3)', a
print '("s = """, a, """")', trim(s)
i には幅を表す数値を指定する。f にはドットを挟んで幅と小数点以下の桁数を指定する。a でも幅を指定できる (しなくてもよい)。書式で直接文字列を指定することができるが、引用符の中で引用符を使う必要があるので、'...' で挟んだ上で "..." と指定すればよい。区切りはカンマである。"..." の中で二重引用符自体を表現したい場合は、"" と二つ書けばよい。二重引用符 1 つなら """" と 4 つ書くことになる。上の例の場合は、次のように "..." を外側にしたほうがまだわかりやすい？
print "('s = ""', a, '""')", trim(s)
実数の指数表示には、e や es を用いる。 print "(e15.3)", a print "(es15.3)", a a = 12.3 とすると、e は "0.123E+02" になり、es は "1.230E+01" になる。ふつうは es を使うだろう。 空白x で空白を表す。 print "(a, x, a)" trim(s1), trim(s2) 空白の数を指定することもできる。 print "(a, 2x, a)" trim(s1), trim(s2) "(a, x, a)" は "(axa)" としてもよいが、"(a, 2x, a)" は "(a2xa)" としてはならない。2 が a の幅と解釈されてしまうので、カンマが必要である。 改行"/" で改行を表す。 print "(/)" この場合、print 文自身の改行もあるので 2 回改行される。 逆に、改行したくない場合は、書式ではどうにもできないようで、write 文のオプションを使う。 write (*, "(/)", advance="no") これだと 1 回しか改行されない。 繰り返し配列を次のように出力する。 print "(f7.3)" a これだと、各要素が 1 行ずつ表示される。1 行に 2 個ずつ表示させたい場合は、次のようにする。 print "(2f7.3)" a 隙間を調整したい場合は、次のようにする。 print "(2(f7.3, 2x))" a 括弧で囲った部分の前に数字をつけると、その部分が繰り返される。 コマンドライン引数コマンドライン引数を処理するには、たとえば、次のようにする。 character(len=32) :: filename if (command_argument_count() < 1) then print "(a)", "usage: prog file" stop end if call get_command_argument(1, filename) print "(3a)", "'", trim(filename), "'" print で書式を指定しているのは、こうしないと変な空白が入るため (FORTRAN 77 からある仕様)。 オブジェクト指向プログラミング派生型の拡張と class ポインタ最近の Fortran にはオブジェクト指向プログラミング用の機能がある (使うべきかどうかは置いといて)。モジュールの例で示した profile モジュールはデータ指向的である。それでも十分ではあるが、もう一歩進めることもできる。 たとえば、profile_t 型に新たなフィールドを追加したい場合は、拡張した型を定義できる。
type, extends(profile_t) :: profile2_t
character(len=10) :: tail_type
end type profile2_t
type(profile2_t) :: prof
prof = profile2_t("omochi", "sabatora", 4, "long")
extends で profile_t 型を拡張 (継承) している。 当然だが、profile_t のポインタで profile2_t を指すことはできない。
type(profile2_t), target :: prof
type(profile_t), pointer :: p
prof = profile2_t("omochi", "sabatora", 2, "long")
p => prof ! エラー: 型が違う
type の代わりに class を用いると、profile_t のポインタで profile2_t を指すことができる。
type(profile2_t), target :: prof
class(profile_t), pointer :: p
prof = profile2_t("omochi", "sabatora", 2, "long")
p => prof
これにより、いわゆる多態性 (あるいは多相性、polymorphism) を実現できる。 データアクセス用の手続きもポインタで呼べる。 call set_name(p, "annko") ただ、これはまだオブジェクトっぽくはない。モジュールを次のように書き換える。
module profile_m
implicit none
private
type, public :: profile_t
character(len=10) :: name
character(len=10) :: type
integer :: age
contains
procedure :: set_name, set_type, set_age
end type profile_t
contains
subroutine set_name(p, name)
class(profile_t), intent(inout) :: p
character(*), intent(in) :: name
p%name = name
end subroutine set_name
subroutine set_type(p, type)
class(profile_t), intent(inout) :: p
character(*), intent(in) :: type
p%type = type
end subroutine set_type
subroutine set_age(p, age)
class(profile_t), intent(inout) :: p
integer, intent(in) :: age
p%age = age
end subroutine set_age
end module profile_m
procedure 指定が type の contains に含まれていることと、サブルーチンの引数の型が type から class に変わっていることに注意。これにより、データアクセス用手続きを次のように呼び出せる。
call p%set_name("annko")
これを型束縛手続き (type-bound procedure) という。手続きの第 1 引数には、呼び出し元のデータ (上の例では p 自身) が渡されるようになっている (Python の self みたいな感じか)。 profile2_t に手続きを追加することもできるが、モジュール化する必要がある。
module profile2_m
use profile_m
implicit none
private
type, extends(profile_t), public :: profile2_t
character(len=10) :: tail_type
contains
procedure :: set_name, set_tail_type
end type profile2_t
contains
subroutine set_name(p, name)
class(profile2_t), intent(inout) :: p
character(len=*), intent(in) :: name
character(len=10) :: temp
write (temp, '(a, "''")') name
call p%profile_t%set_name(temp)
end subroutine set_name
subroutine set_tail_type(p, tail_type)
class(profile2_t), intent(inout) :: p
character(len=*), intent(in) :: tail_type
p%tail_type = tail_type
end subroutine set_tail_type
end module profile2_m
ここでは、サブルーチン set_name をオーバーライドしている。p%profile_t%set_name で、親クラス (にあたるもの) の手続きを呼び出している (ここではそうする必要はないが)。ちなみに、型束縛手続きに non_overridable 属性を指定すると、オーバーライドできなくなる (Java/C++ の final みたいなやつ)。 型束縛手続きは、公開する名前と実装の名前を分けることができる。 procedure :: set_name => profile2_t_set_name 上の例では、公開名を set_name、実装名を profile2_t_set_name としている。 上の例では、派生型 profile_t、profile2_t に public を指定しているが、これは派生型を見えるようにすると同時に、そのメンバーも公開してしまう (デフォルトが public なので)。メンバーを隠蔽したい場合は、改めて private 指定する必要がある。
type, public :: profile_t
private
character(len=10) :: name
character(len=10) :: type
integer :: age
contains
procedure :: set_name, set_type, set_age, get_name, get_type, gat_age
end type profile_t
この場合は setter だけでなく getter もいるだろう。private な手続きが必要なら、contains のほうでも private/public 指定が可能である。 初期化と final 手続きC++ などのクラスだと、初期化と後始末にそれぞれコンストラクタとデストラクタが呼ばれるが、Fortran にはそういうものは明確にはなさそうである。ただ、派生型の初期化の作法は一応あるので、コンストラクタはそれを上書きすればよさそうである。デストラクタについては、後始末手続き (final procedure) という機能があるので、それを設定する。profile モジュールで例を示す。
module profile_m
implicit none
private
type, public :: profile_t
character(len=10) :: name
character(len=10) :: type
integer :: age
contains
procedure :: set_name, set_type, set_age
final :: profile_t_final
end type profile_t
interface profile_t
procedure :: profile_t_init
end interface
contains
function profile_t_init(name, type, age) result(r)
character(len=*), intent(in) :: name
character(len=*), intent(in) :: type
integer, intent(in) :: age
type(profile_t) :: r
print *, "init"
r%name = name
r%type = type
r%age = age
end function profile_t_init
subroutine profile_t_final(p)
type(profile_t), intent(inout) :: p
print *, "final"
end subroutine profile_t_final
...
end module profile_m
モジュール手続きとして profile_t_init と profile_t_final を用意している。interface profile_t で profile_t_init を設定している。ここで profile_t の public 設定がいるかと思ったが、派生型と同名のせいか、不要なようである。profile_t_final は、派生型のほうで final 属性の手続きとして設定している。 たとえば、次のようにブロックの中で profile_t の変数を定義してみる。
block
type(profile_t) :: prof
prof = profile_t("omochi", "sabatora", 2)
print *, prof
end block
これを実行すると、次のような出力になる。 init omochi sabatora 2 final 初期化は、ふつうの派生型の初期化に見えるが、ちゃんと profile_t_init が呼ばれている。また、スコープを抜けると、ちゃんと profile_t_final が呼ばれているのがわかる。 変数がスコープから抜けるタイミングで final 手続きが呼ばれるものだと思うだろう。概ねそうなのだが、実はスコープから抜けていなくても呼ばれることがあるし、スコープから抜けても呼ばれないことがある。 final 手続きが呼ばれるタイミング上の例では、初期化を関数で行なっているが、返り値がコピーされるであろうところが気になる。そのため、サブルーチンで初期化を行いたくなるかもしれない。
subroutine profile_t_init(p, name, type, age)
type(profile_t), intent(out) :: p
character(len=*), intent(in) :: name
character(len=*), intent(in) :: type
integer, intent(in) :: age
print *, "init"
p%name = name
p%type = type
p%age = age
end subroutine profile_t_init
実は、上のサブルーチンでは問題がある。というのも、init の前に final が呼ばれるのである。どういうことかというと、final 手続きの呼び出しは、intent(out) 指定された引数に入れられた場合と、"=" による代入がなされるときらしい。中身をそっくり入れ替えるなら、今あるデータを破棄しておこうということだろう。ただ、初期化時にはデータは持たないので、final 手続きが呼ばれるのは余計である。上記の場合は、intent(out) の代わりに intent(inout) にすればよい。
subroutine profile_t_init(p, name, type, age)
type(profile_t), intent(inout) :: p
...
この例の他にも注意が必要な場合がある。たとえば、初期化で領域確保、final 手続きで領域解放するような派生タイプがあるとする。初期化で領域を確保したあとに、ある手続きにデータを入れてもらおうと intent(out) の引数にその変数を入れてしまうと、final 手続きが呼ばれて領域が解放されてしまう。こちらとしてはデータを入れてもらうだけなので "intent(out)" のつもりなのだが、確保した領域も含めて持ち込み ("intent(inout)") と考えるべきなのだろう。 また、たとえば、関数の返り値を直接別の関数の引数に与えて処理を行う、といったことを考えた場合、返り値を引き受ける明示的な変数がないので、一時的な変数が現れては消える様をイメージするが、Fortran の場合はそのイメージとは異なり、この場合は final 手続きが呼ばれないようである。領域解放を final 手続きに頼っていた場合、領域が宙ぶらりんになる...というのも変なので、自動で解放されてる？ (わからない) final 手続きは、プログラムの最後に達した場合は呼ばれない (上述の例においてブロックで囲っているのはそのため)。C++ などのデストラクタとは動作が結構異なるようである。あまりあてにしないほうがよいかもしれない。 抽象型profile、profile2 の例は、具体的な型を拡張したものだったが、C++ の抽象クラスとか、Java のインターフェイスのように、ガワだけ用意したい場合がある。Fortran では抽象型でそれが実現できる。 たとえば、形状のクラスを定義したいが、とりあえず面積を計算するメソッドだけは欲しい、という場合、次のようになるだろうか。
module shape_m
implicit none
private
type, abstract, public :: shape_t
contains
procedure(get_area_int), deferred :: get_area
end type shape_t
interface
function get_area_int(s) result(r)
import shape_t
class(shape_t), intent(in) :: s
real(kind=8) :: r
end function
end interface
end module shape_m
type に abstract が指定されている。型束縛手続きとして get_area が deferred 指定されているが、実装はあとでするよ、という意味であろう。"procedure(get_area_int)" の "get_area_int" はインターフェイス名である。インターフェイスは interface で定義する。ここで突然出てきた import はなにかというと、どうも、インターフェイスを定義するときにはまだ shape_t の定義が完了していないので、循環を断ち切るため？に必要らしい。 当然、これ自身のインスタンス (にあたるもの) は定義できない。 type(shape_t) :: s ! エラー: 抽象型なので定義できない class ポインタならいける。 class(shape_t), pointer :: s たとえば、rectangle 型を実装すると、次のようになる。
module rectangle_m
use shape_m
implicit none
private
type, extends(shape_t), public :: rectangle_t
real(kind=8) :: width
real(kind=8) :: height
contains
procedure :: get_area
end type rectangle_t
contains
function get_area(s) result(r)
class(rectangle_t), intent(in) :: s
real(kind=8) :: r
r = s%width*s%height
end function get_area
end module_m
ここで get_area を実装しなかったら怒られる。 クラスメソッドC++ の静的メンバ関数 (クラスメソッド) のようなものは Fortran にはないが (クラスがないし)、似たようなことはできなくはない。 クラスメソッドは、要するに特定のインスタンスに結びついていないメソッドである。Fortran の場合、"obj%proc" の形で手続きを呼び出した場合、obj が手続きの第 1 引数として渡される。これを渡さない方法がある。これは Python の @staticmethod デコレータに似ているように思えるが、モジュールの変数や手続きにアクセスできるので、どちらかというと @classmethod のほうである。派生型の手続きの指定で、次のようにする。 procedure, nopass :: proc nopass を指定すると、呼び出し元のデータが手続きに渡されなくなるので、特定のインスタンスに結びつかないようになる。代わりに pass を指定すると、渡されるようになる (デフォルトは pass)。静的変数にあたるものが必要であれば、モジュールの変数を使えばよい。 モジュールの手続きとの違いがよくわからないが、気分の問題？ 型の調査と変換型を調べるには、select type 構文を使う。 select type (v => p) type is (profile_t) print *, "profile_t" print *, v type is (profile2_t) print *, "profile2_t" print *, v class default print *, "???" stop end select v で p を指している形になっているが、この v の型宣言は不要である (associate のような感じ)。class ポインタを print 文に渡そうとすると怒られるのだが (polymorphic だからダメとか)、ここで v はちゃっかり type ポインタに変換される。"type is ..." の代わりに "class is ..." というのもあって、こちらだと v は class ポインタになるようで、print 文で怒られる。つまり、select type は型を調べると同時に型変換 (いわゆるキャスト) もやっているわけである。 型がわかっているのであれば、transfer 関数で直接的に型変換できる。 print *, transfer(p, prof) transfer 関数の第 2 引数は、値ではなくて型を参照している。これと同じ型に変換せよということである。これが Fortran のキャストにあたる (第 2 引数の名前は mold で、だから casting ではなくて molding だとかなんとか)。 class ポインタと type ポインタ (と便宜上呼ぶが) の 2 つがあるのはなぜだろうか。おそらく、C++ だとクラスの定義でメンバ関数を仮想関数にするかどうかで多相性を制御するのに対して、Fortran では、ポインタの側で class ポインタにするか type ポインタにするかで多相性を切り替えることにしたのだろう。 演算子オーバーロードC++ 同様、Fortran でもある型に対して演算子の動作を定義することができる。また、新たな演算子を定義することもできる。 たとえば、"a" という文字式と "b" という文字式を足したら "a + b" という文字式になるような型 expr_t を考える。ここでは "+" と "="、そして出力について実装する。
module expr_m
implicit none
private
type, public :: expr_t
character(len=128) :: expr
end type expr_t
public :: assignment(=)
interface assignment(=)
procedure :: assign
end interface
public :: operator(+)
interface operator(+)
procedure :: add
end interface
public :: write(formatted)
interface write(formatted)
procedure :: print_expr_t
end interface
contains
subroutine assign(s1, s2)
type(expr_t), intent(out) :: s1
type(expr_t), intent(in) :: s2
s1%expr = s2%expr
end subroutine assign
function add(s1, s2) result(r)
type(expr_t), intent(in) :: s1, s2
type(expr_t) :: r
write (r%expr, '(a, " + ", a)') trim(s1%expr), trim(s2%expr)
end function add
subroutine print_expr_t(e, unit, iotype, v_list, iostat, iomsg)
class(expr_t), intent(in) :: e
integer, intent(in) :: unit
character(len=*), intent(in) :: iotype
integer, intent(in) :: v_list(:)
integer, intent(out) :: iostat
character(len=*), intent(inout) :: iomsg
write (*,*) trim(e%expr)
iostat = 0
end subroutine
end module expr_m
expr_t は、文字式を保存する文字列 expr だけを持つ。interface assignment(=) で代入の手続き assign を指定している。interface operator(+) で和の手続き add を指定している。interface write(formatted) は、print や write における派生型の書式 dt のための設定を行う手続きを指定している。 assign は、s1 の文字式を s2 の文字式にコピーしているだけである。add は、s1 の文字列と s2 の文字列から和の文字式を作成している。print_expr_t は出力設定をしている。引数部分がやたら長いが、そういうテンプレートのようである。ちなみに、write の代わりに read の設定を行う場合は、最初の引数を intent(inout) にする。 operator の引数は intent(in) にする必要がある (引数の内容をこっそり書き換えるということはできない)。また、assign の引数は intent(out) (あるいは intent(inout)) と intent(in) にする必要がある。 お試し。
program example
use expr_m
implicit none
type(expr_t) :: a = expr_t("a")
type(expr_t) :: b = expr_t("b")
type(expr_t) :: c = expr_t("c")
type(expr_t) :: d = expr_t("d")
c = a + b
print *, c
d = c
print *, d
end program example
"+" と "=" を試している。どちらも "a + b" と出ればよい。 operator(+) の代わりに、operator(.add.) のようにまったく新しい演算子を定義することもできる。その場合、次のように使えるようになる。 c = a .add. b 総称関数Fortran には昔から総称関数 (generic function) という、C++ でいうオーバーロードのような仕組みがあるにはあったが、組み込み関数に限られていた。最近では、これを自分で作れるようになっている。 たとえば、変数の値を入れ替えるサブルーチン swap を考える。実装は単純で、データ型によらないが、サブルーチンの定義上は実数用、整数用、などと分けて作る必要があるし、呼び出しも別になる。処理がかぶるのはしょうがないが、せめて呼び出しぐらいは統一しようというのが総称手続きである。
module swap_m
implicit none
private
public :: swap
interface swap
procedure :: swap_r, swap_i
end interface
contains
subroutine swap_r(a, b)
real, intent(inout) :: a, b
real :: temp
temp = a
a = b
b = temp
end subroutine swap_r
subroutine swap_i(a, b)
integer, intent(inout) :: a, b
integer :: temp
temp = a
a = b
b = temp
end subroutine swap_i
end module swap_m
interface で、実数用の swap_r と swap_i を swap に統一している。使用例は以下の通りである。 program example use swap_m implicit none real :: a = 1., b = 2. integer :: i = 1, j = 2 call swap(a, b) print *, a, b call swap(i, j) print *, i, j end program example 同じ swap で、実数と整数両方に対応できている。ちなみに、swap_r に整数変数を入れようとすると怒られる (この場合は変換も使えない)。 この程度ならまだよいが、swap_r と swap_i で処理自体が共通している (二重になっている) のが気になる。これは、include 文で対応できる。swap.inc に、次のように書く。 temp = a a = b b = temp これを include する。
module swap_m
implicit none
private
public :: swap
interface swap
procedure :: swap_r, swap_i
end interface
contains
subroutine swap_r(a, b)
real, intent(inout) :: a, b
real :: temp
include "swap.inc"
end subroutine swap_r
subroutine swap_i(a, b)
integer, intent(inout) :: a, b
integer :: temp
include "swap.inc"
end subroutine swap_i
end module swap_m
これで、たとえば swap.inc の内容に変更があった場合、自動的にすべての手続きに適用されることになる。 型に依存しないデータ構造C++ のテンプレートのように、型に依存しないプログラミングについては、Fortran の機能としては上述の総称関数くらいしかない。C 言語でも機能自体はないが、型に依存しないデータ構造については無理やり作れなくもない。Fortran でもそれに近いことができる。 なんでも入る配列 array を定義してみる。
module array_m
implicit none
private
type array_item
class(*), pointer :: value => null()
contains
final :: array_item_final
end type array_item
public :: array_t
public :: set, value_at
type array_t
type(array_item), allocatable :: items(:)
contains
procedure :: set, value_at
final :: array_t_final
end type array_t
interface array_t
procedure :: array_t_init
end interface
contains
subroutine array_item_final(this)
type(array_item), intent(inout) :: this
if(associated(this%value)) deallocate (this%value)
end subroutine array_item_final
function array_t_init(n) result(a)
integer, intent(in) :: n
type(array_t) :: a
allocate (a%items(n))
end function array_t_init
subroutine array_t_final(this)
type(array_t), intent(inout) :: this
if(allocated(this%items)) deallocate (this%items)
end subroutine array_t_final
subroutine set(this, i, value)
class(array_t), intent(inout) :: this
integer, intent(in) :: i
class(*), intent(in) :: value
if(associated(this%items(i)%value)) deallocate(this%items(i)%value)
allocate (this%items(i)%value, source=value)
end subroutine set
function value_at(this, i) result(r)
class(array_t), intent(inout) :: this
integer, intent(in) :: i
class(*), pointer :: r
r => this%items(i)%value
end function value_at
end module array_m
array_t が配列の型で、データは array_item の形でもつ。array_item はデータのポインタとして class(*) のポインタをもつ。class(*) のポインタは、無制限多相的 (unlimited polymorphic) ポインタという。たぶん C 言語の void ポインタのようなものだろう。最初は null() を指している。null() はヌルポインタ (C/C++ でいう NULL、nullptr) である。 コンストラクタ array_t_init で指定した個数の要素を確保する (array_item としては空っぽ)。set で指定した位置にデータを入れるのだが、もし array_item に領域が確保されていれば削除し、その後領域を確保する。associated はポインタが何かを指しているかチェックする (領域が確保されているかどうかチェックするのに allocated というのもあるが、allocatable にしか使えない)。value_at でデータを参照できる (ただし、class(*) のポインタ として返す)。デストラクタ (final procedure) が設定してあるので、後始末は自動で行われる (はず)。 お試し。
program example
use array_m
implicit none
type :: profile
character(len=10) :: name
character(len=10) :: type
integer :: age
end type profile
type(profile) :: prof = profile("omochi", "sabatora", 4)
type(array_t) :: a
integer :: i
a = array_t(4)
call a%set(1, 1)
call a%set(2, 2.)
call a%set(3, "three")
call a%set(4, prof)
do i = 1, 4
select type (v => a%value_at(i))
type is (integer)
print *, v
type is (real)
print *, v
type is (character(len=*))
print *, v
type is (profile)
print *, v
end select
end do
end program example
出力
1
2.00000000
three
omochi sabatora 4
整数、実数、文字列、派生型を入れてみているが、ちゃんと動いている。value_at が class(*) のポインタを返すため、データの取り出し方がちょっと独特になる。ただ、同様の処理を C 言語で書くのに比べたら、おとなしめというか、あまり無理やり感はない (キャストを直接は使っていないし)。 値を transfer 関数で直接取り出すこともできる。 print *, transfer(a%value_at(1), 1) ただ、この方法だと、文字列の取り出し方がわからない (文字列を指定しても 1 文字しか取れないし)。 参考 CoarrayCoarray のコンパイルCoarray は、並列計算を行うためのものである。Open MP よりは MPI に似ている。gfortran で Coarray を用いるには特別なオプションがいる。一番簡単な Coarray のプログラムは、次のようなものになるだろう。 program example implicit none print *, num_images() end program example コンパイルは、次のようにする。 $ gfortran -fcoarray=single example.f90 オプションを付けないと怒られる。実行すると、1 が出力されるはずである。
$ ./a.out
1
Coarray は、image という単位でプログラミングする。image はそれぞれ別の人とみなされる。プロセスと言ったほうがわかりやすいかもしれない。ただ、プロセスとは限らないのだろう。オプションの -fcoarray は、Coarray を使うよという意味で、single は image が 1 つだよという意味である。image が 1 つしかないなら並列計算もなんもないのだが、こうしないと Coarray を用いたコードがコンパイルされないのでしかたがない。 複数の image を使うには、別途環境整備が必要である。OpenCoarrays というものを用いる。こちら を参考にしてインストールすればよい。 インストールできたら、次のようにコンパイルする。 $ caf example.f90 実行。
$ cafrun -n 1 ./a.out
1
$ cafrun -n 2 ./a.out
2
2
実行できる image の数は環境による。中身は結局のところ MPI のようである。
$ caf --show
/usr/local/bin/gfortran -I/usr/local/Cellar/opencoarrays/2.9.2/include/OpenCoarrays-2.9.2_GNU-10.2.0 -fcoarray=lib -Wl,-flat_namespace -Wl,-commons,use_dylibs -L/usr/local/Cellar/libevent/2.1.12/lib -L/usr/local/Cellar/open-mpi/4.0.5/lib ${@} /usr/local/Cellar/opencoarrays/2.9.2/lib/libcaf_mpi.a /usr/local/lib/libmpi_usempif08.dylib /usr/local/lib/libmpi_usempi_ignore_tkr.dylib /usr/local/lib/libmpi_mpifh.dylib /usr/local/lib/libmpi.dylib
$ cafrun --show
/usr/local/bin/mpiexec -n
指定したプロセス数が多すぎると文句を言われる場合は、次のような内容の hostfile を用意する。 localhost slots=10 slots の数はお好きな数。cafrun ではなく、mpiexec を使って実行する。 $ mpiexec --hostfile=hostfile -n 10 ./a.out Coarray の使い方ここでは 2 個の image があるものとして例を示す。Coarray というからには、配列っぽいデータを扱う。たとえば、整数の変数 (スカラー) を定義するには、次のようにする。 integer :: a[*] あるいは、次のように書いてもよい。 integer, codimension[*] :: a たとえば、自分の image の番号を入れる。 a = this_image() print *, a この出力は、次のようになる。
1
2
つまり、image 間で変数 a の内容を共有していない。よその image の変数には割と簡単にアクセスできる。たとえば、2 番目の image の変数 a にアクセスしたければ、a[2] とすればよい。次のようなことができる。 if (this_image() == 1) then a[2] = 1 else a[1] = 2 end if print *, this_image(), a お隣さんの変数に互いに値を入れている。出力は次のようになる。
1 2
2 1
いちいちデータをどこそこに転送するとか、受け取りを待ちますとか、そういう処理を書く必要がない。 同期それぞれの image で配列の和をとって、それらをさらに合計するには、次のようにする。
integer :: a(10)[*]
integer :: s[*]
integer :: i
a = [(i, i = 1, 10)]
s = 0
if (this_image() == 1) then
do i = 1, 5
s = s + a(i)
end do
else
do i = 6, 10
s = s + a(i)
end do
end if
print *, this_image(), s
sync all
if (this_image() == 1) then
s = s + s[2]
print *, s
end if
sync all はすべての image で同期をとるためのものである。もっと細かく image どうしで同期をとりあうこともできる。
if (this_image() == 1) then
sync images (2)
else if (this_image() == 2) then
sync images (1)
end if
たぶん、image の中のいくつかだけが同期をとりあうというのはあまりないと思うが、どれか 1 つのイメージが特別な立場にいて、その他のイメージと同期を取り合うということはありえる。たとえば、image 1 がその他全部と同期を取り合うとしたら、次のようになる。
if (this_image() == 1) then
sync images (*)
else
sync images (1)
end if
データの同期上の配列の和の例では、image 1 ですべての image の値の和をとったが、他の image はそのことを知らない。それを知らせるためには、co_broadcast を用いる。 call co_broadcast(s, 1) 上の例は、image 1 の s の値を各 image に配信するよという意味である。これにより全 image の値が同期する。 データの簡約上の例で見た、全 image の値の和をとって各 image に配信する、というのは co_sum を呼ぶだけでできる。 call co_sum(s) 配信が余計であれば、image 1 だけに値をセットすることもできる。 call co_sum(s, result_image=1) 和ではなくて、全 image の値の中で最大のものとか最小のものが欲しければ、co_max や co_min を用いる。 call co_max(smax) call co_min(smin) たとえば、全 image のデータの最大値が欲しいという場合は、各 image でそれぞれ最大値を出した上で、それらの中での最大値を co_max でとればよい。 より一般的な簡約 (reduction) も co_reduce できる。 call co_reduce(s, f) たとえば、co_sum と同じ結果にしたい場合は、関数 f を次のように定義すればよい。 pure function f(a, b) integer, intent(in) :: a, b integer :: f f = a + b end function f 動的配列配列を動的に確保したい場合は、次のように宣言する。 integer, allocatable :: a(:)[:] Coarray のサイズもコロンにしてあることに注意。領域確保は次のようにする。 allocate (a(10)[*]) deallocate はふつうにすればよい。 C 言語とのやりとりFortran から C 言語の関数を呼ぶFortran には、C 言語の関数を呼び出す仕組みがある。たとえば、文字列を表示するだけの関数 (値を返さないので、Fortran から見たらサブルーチン) hello があるとする。 hello.c #include Fortran では次のように呼ぶ。
program example
use, intrinsic :: iso_c_binding
implicit none
interface
subroutine hello() bind(c)
end subroutine hello
end interface
call hello()
end program example
特殊なモジュール iso_c_binding を読み込んで、C 言語の関数のインターフェイスを指定している。その際、bind(c) というものを指定する。コンパイルは次のようにする。 $ gcc -c hello.c $ gfortran -c example.f90 $ gfortran example.o hello.o 引数がある関数の場合は、次のようになる。 dbl.c
double dbl(double x)
{
return 2.*x;
}
Fortran では次のようになる。
program example
use, intrinsic :: iso_c_binding
implicit none
interface
function dbl(x) bind(c)
import c_double
real(kind=c_double), value :: x
real(kind=c_double) :: dbl
end function dbl
end interface
print *, dbl(2d0)
end program example
ここで、実数の型として c_double というものを用いている。これは C 言語の double 型を意味する。また、引数 x には value 属性がついている。これは、値渡しを想定している。 C 言語の型に対応した型には次のようなものがある。
構造体を扱うには、派生型で bind(c) を指定すればよい。 文字列引数の場合は、次のようにする。 say.c #include Fortran
program example
use, intrinsic :: iso_c_binding
implicit none
interface
subroutine say(msg) bind(c)
import c_char
character(len=1, kind=c_char), dimension(*) :: msg
end subroutine say
end interface
call say("hello")
end program example
C 言語は文字列型がないので、文字型の配列にしている。dimension は dimension(*) とする。 C 言語から Fortran の手続きを呼ぶ逆に、C 言語から Fortran の手続きを呼ぶ場合は、次のようにする。 dbl.f90 function dbl(x) bind(c) use, intrinsic :: iso_c_binding real(kind=c_double), value, intent(in) :: x real(kind=c_double) :: dbl dbl = 2.*x end function dbl この場合のキモは、bind(c) と引数の value 指定である。bind(c) によって C 言語で名前解決できるようになる。引数の型については real(kind=8) でも動きはする。 C 言語 #include この場合は、はじめから C 言語に呼ばれることを想定している。それが想定されていない手続きを呼びたい場合は、ラッパーを作ればよいのかもしれない。 シグナルCtrl-C などの入力を検出したいときがある。標準では対応できないが、gfortan (など) ではシグナルが扱える。たとえば、次の例は、無限ループ状態を Ctrl-C で抜けようとしたことを検出する。
program example
implicit none
integer, parameter :: signum = 2 ! SIGINT
call signal(signum, s)
do
! do nothing
end do
contains
subroutine s()
print *
print *, "Ctrl-C"
stop
end
end
サブルーチ signal でシグナル番号と呼び出されるサブルーチンを設定する。実行すると無限ループに入るので、Ctrl-C を入力すると、設定したサブルーチンが呼び出されて、プログラムが止まる (呼び出されるサブルーチンで stop を呼ばないと止まらない)。 プリプロセッサC/C++ で使われているプリプロセッサは、Fortran の標準にはないが、コンパイラが対応していることがある。gfortran では、"-cpp" オプションをつければよい。ただし、行の先頭に空白は入れられないようである。また、プリプロセッサ上は大文字小文字を区別する。 pure の無効化pure 関数をデバッグしたいときに、print 文を使いたいが、いちいち pure を外さないといけない。めんどくさいときは、gfortran にオプション "-cpp -Dpure" を渡せば、pure をすべて無効化できる。 参考文献
| |||||||||||||||||||||||||||||||
| PENGUINITIS | |||||||||||||||||||||||||||||||