FortranのPureなFunction

はじめに

pure

　モダンなFortran90になって、言語仕様としてpure属性が定義されました。

pure functionは日本語で「純粋関数」です。

   f(x) = y

という計算をするときに、純粋関数$f$は、

• 副作用がない：関係ない変数 $z$ に影響されないし、影響しない
• 参照透過性：$x$ の値が同じなら $y$ の値は必ず同じ

といった性質があります。ちなみに、プログラミング言語での「関数」はサブルーチン的な役割があるので純粋でないことも多いですが、数学で出てくる「関数」は普通は純粋関数じゃないかなと思います。

　副作用がない関数というのはプログラミングにおいて利点がいろいろあります。これについては関数型言語（LISPやHaskellなど）での議論が詳しく、あまり深入りするとFortranから離れてしまいますので、ここではFortranにおいての利点に注目します。
　Fortrannerにとっての純粋関数の利点は次の２つかと思います。

　* 関数の外の変数を意識しなくて良いのでデバッグがやりやすい
　* 関数の実行が独立なので並列化しやすい

　２つ目の利点は「pureは並列計算の時にさらに活きる属性」ということです。どういうことかという説明の前に、elementalについて触れておきます。

elemental

Fortranは配列を扱うのが得意な言語です。特化しすぎてリストとか辞書型とか木構造を使うと辛い
このとき、配列の各成分に関数を適用して別の配列を作るというのはよくあると思いますが、そのたびにdoループを回すのは冗長です。そこでfortranではelementalという属性が定義されています。
例えば三角関数のsinがありますが、これに配列を渡すことを考えます。

double precision, allocatable :: v(:), w(:)
double precision, allocatable :: A(:, :), B(:, :)
! 中略
w = sin(v)
B = sin(A)

結果のwとBは、vやAと同じ形状の配列で、各成分ごとにsinを適用したものだと想像できます。Fortranのelemental関数は、上記の例のようにその直感をそのまま実現しています。
　なんか「W = v.map(sin)としたほうが直感的だ」という声が聞こえた気がしましたが、Fortranタグの前には無力です。¹

ここで注目したいこととして、各成分にsin関数を適用する際、sinは成分を超えて情報をやりとりしないということがあります。sin(v(1))の計算時にsin(v(2))を考えなくても大丈夫です。このため、sin(v(1))、sin(v(2))、・・・、sin(v(i))を独立に計算できます。これは並列計算にもってこいな性質です。

　逆に言うとある成分での計算が他の成分に依存する場合、色々な問題が発生します。なので、このsinのようなelemental関数は引数だけに依存するように、つまりpure関数としなければなりません。(追記：最近のFortran2008ではimpureなelementalが可能になりました。記事の最後にまとめておきます。)
例えば、配列の要素をどのような順番で適用するかとかが問題になりますが、順番というのは関数の外にあって、配列ループを回すときのインデックスと同値です。これでは「引数が決まれば結果が決まる」という純粋関数の性質を満たしません。また、コンパイラにどうやって順番を教えるかという問題もあります。

　別の利点として、上記の例のように引数の配列の形状を限定しません。もしelementalでない関数で上記のように複数の形状の配列に対応しようとすると、スカラー値用のsin、1次元配列用のsin1d、２次元配列用のsin2d、・・・と関数の定義がとんでもなく面倒になります。

　まとめますと、pureやelementalというのは配列の取り扱いや並列計算において重要なキーワードです。

この記事では主にpureな関数の振る舞いを紹介したいと思います。

環境

以下の環境で確認しました。

$ gfortran --version
GNU Fortran (Ubuntu 5.4.0-6ubuntu1~16.04.10) 5.4.0 20160609
Copyright (C) 2015 Free Software Foundation, Inc.

$ uname -rosi
Linux 4.4.0-139-generic x86_64 GNU/Linux

申し訳ありませんがintel fortranは現在所持していないので、確認できません。

基本

PUREの条件

ある関数がpureであるための条件は、導入の話を反映させると以下のようになります。

• 副作用がない
  • 引数がintent(in)である²
  • 外部の変数への代入がない
  • ファイル・端末への出力がない
• 同じ引数なら必ず同じ結果を返す
  • 中でrandom_numberのような呼ぶたびに結果が違う関数は呼べない
  • ファイル・端末からのreadはできない

端末への出力が副作用だというのは、副作用についての議論に馴染みがないとわかりにくいかもしれませんが、ここでは詳しくは割愛します。

具体例

早速pureやelementalな関数の定義例を見てみます。

module func
    implicit none
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        r = n*n
    end function square

    pure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        r = n*10 + m
    end function addmul
end module func

square(n)、addmul(n, m)は引数のnやmが決まれば結果が一意に決まり、また関数外の変数にも手を出しておらず、副作用もありません。また、nやmもintent(in)の効果で変更されません。pureとしての資格は十分です。

これを以下のプログラムのように呼び出します。

program test
    use func
    implicit none
    integer :: n, a(10), b(10), i

    n = 10
    a = (/(i, i=0, 9)/)
    b = (/(i, i=9, 0, -1)/)
    print *, square(n)
    print *, addmul(a, b)
end program test

これはコンパイルして出力すると以下のようになります。

出力

         100
           9          18          27          36          45          54         63          72          81          90

　ひとつ目のprintではsquare(10)の結果、つまり100が出力されます。挙動自体は普通の関数ですね。
　ふたつ目のprintでは、09, 18, 27, ..., 90 が出力されます。というのはaddmulはひとつ目の引数を10倍して、ふたつ目の引数を足した整数を出力しますが、aは0, 1, 2, ..., 9で、bは9, 8, 7, ..., 0が代入されていますので、addmulの結果は9, 18, 27, ..., 90と９の倍数が並んだものになります。elementalなので配列を直接引数においてやると、同じ形の配列が返って来ました。

バリエーション例

elementalではpure宣言は無くてもよい

    elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        r = n*10 + m
    end function addmul

pureがなくてもpureとしての性質は満たさないといけません。

文字列への変換のためのwrite

addmulについて、結果を整数型ではなく文字列で返してみます。

    elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        character(len=4) :: r
        write(r, '(i4)') n*10 + m
    end function addmul

このwrite文は副作用のない使用法なので、pure関数内で使えます。ただし、文字列長は固定です。

「値を変えない」 共有変数

たとえばモジュールの変数やCOMMONブロックでも、左辺に現れなければ使用可能です。

module func
    implicit none
    integer :: i_local
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        r = i_local*n*n
    end function square

    pure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        integer :: c1, c2
        common /hoge/ c1, c2
        r = c1*n*10 + c2*m
    end function addmul
end module func
 
program test
 ! 中略
    common /hoge/ c1, c2
    i_local = 10
    c1 = 2
    c2 = 4
    n = 10
    a = (/(i, i=0, 9)/)
    b = (/(i, i=9, 0, -1)/)
    print *, square(n)
    print *, addmul(a, b)
end program test

出力

        1000
          36          52          68          84         100         116         132         148         164         180

これは注意してやらないと参照透過性が保てないかもしれません。

NG例

左辺に現れた共有変数

pure関数の中で外部の変数を変更することはできません。

module func
    implicit none
    integer :: i_local
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        i_local = i_local + 1
        r = n*n
    end function square

コンパイルエラー

         i_local = i_local + 1
        1
Error: Variable ‘i_local’ can not appear in a variable definition context (assignment) at (1) in PURE procedure

COMMONブロックの変数についても、同様のエラーが出ます。

SAVE

    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        integer, save :: j = 0
        r = n*n
    end function square

コンパイルエラー

 integer, save :: j = 0
                     1
Error: SAVE attribute at (1) cannot be specified in a PURE procedure

使っていようといまいと、現れた時点でOUTです。

Write文

IOは副作用の代表格です。

    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        write(*, *) n
        r = n*n
    end function square

コンパイルエラー

         write(*, *) n
                     1
Error: IO UNIT in WRITE statement at (1) must be an internal file in a PURE procedure

INTENT忘れ

INTENT(IN)もしくはVALUEを厳密に要求します。

    pure function square(n) result(r)
        integer :: n
        integer :: r
        r = n*n
    end function square

コンパイルエラー

     pure function square(n) result(r)
                          1
Error: Argument ‘n’ of pure function ‘square’ at (1) must be INTENT(IN) or VALUE

引数をターゲットとしたポインタの割当

基本的にpure関数の中でもallocatableな配列を使ったりpointerを割り当てたりが可能ですが、使えないケースというのもあり、例えば引数が関係するようなポインタ割当はエラーとなります。

    pure function square(n) result(r)
        integer, intent(in), target :: n
        integer :: r
        integer,pointer :: pn
        pn => n
        r = pn*pn
    end function square

コンパイルエラー

         pn => n
              1
Error: Bad target in pointer assignment in PURE procedure at (1)

乱数

pureの中で乱数は使えません。

    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        double precision :: x
        call random_number(x)
        r = int(x*n*n)
    end function square

コンパイルエラー

call random_number(x)
                    1
Error: Subroutine call to intrinsic ‘random_number’ at (1) is not PURE

elementalの形状間違い

任意の配列形状を受け取れるelementalですが、複数の引数をとるときにその次元が一致しない場合はどう計算してよいかわかりませんので、エラーとなります。

    a = (/(i, i=0, 4)/)
    b = (/(i, i=9, 0, -1)/)
    print *, addmul(a, b)

コンパイルエラー

print *, addmul(a, b)
                       1
Error: Different shape for elemental procedure at (1) on dimension 1 (10 and 5)

他にも

• 中で呼び出される関数・サブルーチンはすべてpureでなければならない。
• プログラムのSTOP文を入れられない。
• ファイルの読み書きなどができない。
• externalによる関数渡しはできない。

などといった決まりがあります。まあ、純粋関数の性質を考えたら妥当な仕様かなと思います。

もし何らかのコンパイラで、こういったことができるものがあったとしても、pureの目的を考えると使うべきではないでしょう。

デバッグ

pureはメリットが多い属性ですが、IO文が使えないため、printデバッグができなくなってしまうデメリットがあります。
すべての関数にPUREをつけるという遊びをしている時にいかにデバッグするかを考えてみます。

プリプロセッサ

main.f90

module func
    implicit none
contains 
#ifdef DEBUG
#define PUREF function
#else
#define PUREF pure function
#endif

    PUREF square(n) result(r)
        integer, intent(in) :: n
        integer :: r
#ifdef DEBUG
        write(*, *) "DEBUG"
        write(*, *) "input: ", n
#endif
        r = n*n
    end function square
end module func

program test
    use func
    implicit none
    integer :: n

    n = 10
    print *, square(n)
end program test

要はプリプロセッサのマクロでDEBUGが定義されているとデバッグモードとして適切にソースコードが切り替わるようにしています。

• DEBUGあり：PUREFはfunctionと定義、関数内のwrite文を有効に
• DEBUGなし：PUREFはpure functionと定義、関数内のwrite文を取り除く

デバッグありコンパイル

$ gfortran main.f90 -cpp -DDEBUG
$ ./a.out
 DEBUG
 input:           10
         100

デバッグなしコンパイル

$ gfortran main.f90 -cpp
$ ./a.out
         100

なお、この方法ではelementalを外すとインターフェースもしくは呼び出し側も大幅に変えなければならない可能性があり、工夫が必要です。

デバッガ

正直慣れるとこれが一番かと思います。
基本のプログラムについて、デバッグのデモ目的で中間変数を追加します。

main.f90

module func
    implicit none
contains 
    pure function square(n) result(r)
        integer, intent(in) :: n
        integer :: r
        integer :: s
        s = n*n
        r = s
    end function square

    pure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        integer :: t10
        integer :: t1
        t10 = n*10
        t1 = m
        r = t10 + t1
    end function addmul
end module func

program test
    use func
    implicit none
    integer :: n, a(10), b(10), i

    n = 10
    a = (/(i, i=0, 9)/)
    b = (/(i, i=9, 0, -1)/)
    print *, square(n)
    print *, addmul(a, b)

end program test

デバッガ情報を持たせてコンパイルします。

gfortran -g main.f90

デバッガから起動します。

% gdb ./a.out 
GNU gdb (Ubuntu 7.11.1-0ubuntu1~16.5) 7.11.1

(中略)

For help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from ./a.out...done.
(gdb) run
Starting program: /path/to/a.out 
         100
           9          18          27          36          45          54          63          72          81          90
[Inferior 1 (process 29438) exited normally]
(gdb) q

qでデバッガが終了します。

さて、中間の変数sにちゃんと平方数が入っているか確認したいとなったとき。

        s = n*n

これは8行目にありました。

(gdb) break 8
Breakpoint 1 at 0x4007c5: file main.f90, line 8.

これで8行目にブレークポイントを設定できました。実行してみます。

(gdb) run
Starting program: /path/to/a.out 

Breakpoint 1, func::square (n=10) at main.f90:8
8               s = n*n

8行目のソースを表示しつつ止まりました。ここで変数の内容を出力します。

(gdb) print n
$1 = 10
(gdb) print s
$2 = 0

引数nには正しく10が入っていますね。ですが、まだsには代入されていないようです。1行進めてみます。それにはstepを用います。

(gdb) step
9               r = s
(gdb) print s
$3 = 100

8行目をすぎて9行目に到達し、sにちゃんと代入できていることが確認できました。見たいものが終わったのでプログラムを次のブレークポイントもしくはプログラム最後まで実行します。

(gdb) continue
Continuing.
         100
           9          18          27          36          45          54          63          72          81          90
[Inferior 1 (process 30346) exited normally]

次はt10にちゃんと値が入っているか確認しようと思います。

17行目

        t10 = n*10

　ここでもうひとつ、gdbのコマンドはgnuplotのように１文字まで省略できます。それを示しながらやってみます。

　１７行目のブレークポイントの設定と実行。ブレークポイントより前にあるので、square(10)の出力が出ていますね。

(gdb) b 17
Breakpoint 1 at 0x400792: file main.f90, line 17.
(gdb) r
Starting program: /path/to/a.out 
         100

まずは一回目のブレークポイント到着。addmulの関数内で、引数もちゃんと表示されています。

Breakpoint 1, func::addmul (n=0, m=9) at main.f90:17
17              t10 = n*10
(gdb) p n
$1 = 0
(gdb) s
18              t1 = m
(gdb) s
19              r = t10 + t1
(gdb) p t10
$2 = 0
(gdb) c
Continuing.

pはprint、sはstep、cはcontinueと、最初の1字で大体OKです。
t10には0が入っていますが、これは初期化されたままなのか、正しく0が代入されたのかわかりません。なので別の引数で確認したいところ。

今回addmul(n, m)はelementalで、引数に長さ10の配列a、bを渡しており、内部的にはループとなりますので、17行目のブレークポイントは何度も通ります。ということで2回めのブレークポイントです。

Breakpoint 1, func::addmul (n=1, m=8) at main.f90:17
17              t10 = n*10
(gdb) p n
$5 = 1
(gdb) p t10
$6 = 0
(gdb) s
18              t1 = m
(gdb) p t10
$7 = 10
(gdb) p m
$8 = 8
(gdb) s
19              r = t10 + t1
(gdb) c
Continuing.

引数も配列の次の成分のものとなっていること、またt10にも正しく10が入っていることが確認できました。

3回目のブレークポイント。

Breakpoint 1, func::addmul (n=2, m=7) at main.f90:17
17              t10 = n*10
(gdb) p n
$9 = 2
(gdb) s
18              t1 = m
(gdb) s
19              r = t10 + t1
(gdb) s
20          end function addmul
(gdb) p r
$10 = 27

end functionの行ではまだ関数の中のようですので、rを出力することでこの関数がどんな値を返しているかも見られます。

continueコマンドに整数を渡してやるとその数だけブレークポイントをスキップします。

(gdb) c 10
Will ignore next 9 crossings of breakpoint 1.  Continuing.
           9          18          27          36          45          54          63          72          81          90
[Inferior 1 (process 30419) exited normally]

これで最低限の変数の確認ができるかと思います。

ほかにも関数名で指定してステップ実行など、便利な機能がいろいろあります。たとえば今いる行の前後のプログラムをその場で確認できるコマンドlistがあります。

(gdb) list
12          pure elemental function addmul(n, m) result(r)
13              integer, intent(in) :: n, m
14              integer :: r
15              integer :: t10
16              integer :: t1
17              t10 = n*10
18              t1 = m
19              r = t10 + t1
20          end function addmul
21      end module func

ブレークポイントの設定についても、現在位置からの相対行数、関数名などでの設定などもあります。
http://flex.phys.tohoku.ac.jp/texi/gdb-j/gdb-j_18.html

FortranのみならずCやC++でも使えるコマンドですので、Linuxで生きている人なら覚えて損はないと思います。
計算化学の人みてるとprintデバッグばっかでデバッガなんて使っちゃいない

Cとの連携

Fortran2003からISO_C_BINDINGが導入され、C言語との連携が強化されました。C言語で定義された関数に対しPUREが使えるか試してみます。

なお、先に言っておきますとelemental属性はBIND(C)と衝突しますので、使えません。

     pure elemental function addmul(n, m) result(r) bind(C, name="addmul")
                                  1
Error: BIND(C) attribute conflicts with ELEMENTAL attribute at (1)

普通に

func.c

int square(int n){
    return n*n;
}

main.f90

module func
    use iso_c_binding
    implicit none
    interface
        pure function square(n) result(r) bind(C, name="square")
            import :: c_int
            integer(c_int), value :: n
            integer(c_int) :: r
        end function square
    end interface
end module func

program test
    use func
    implicit none
    integer :: n

    n = 10
    print *, square(n)
end program test

コンパイルと実行

gcc -c func.c
gfortran main.f90 func.o
./a.out

出力

         100

普通のC関数の呼び出しと同じように使えますね。

ルール違反

さて、ここでFortranの目が届かないところで良くないことをしてみます。例えばPURE関数内で端末出力。

func.c

#include <stdio.h>

int square(int n){
    printf("inside C: %d\n", n);
    return n*n;
}

先と同様にコンパイルをしてみますと、何の問題もなくリンクできます。
実行してみます。

出力

inside C: 10
         100

あれ、pure関数が画面出力してる（棒）。

うまく使えばpure関数用の画面出力関数を作ってprintデバッグできるかもしれません。

まとめ

PURE関数をうまく使うといろいろ良いことがあります。

その他1

組み込みの数学関数は大体 pure です。

その他2

Visual StudioのFortran対応きぼん。ソリューション管理とProfilerマジ便利。

追記：impureなelemental

Fortran2008になって、impureキーワードが導入されました。

これを用いると、一部のimpureな操作がelementalでできるようになります。

    impure elemental function addmul(n, m) result(r)
        integer, intent(in) :: n, m
        integer :: r
        write(*, *) n  ! コンパイルエラーが出ない！実行できる！
        r = n*10 + m
    end function addmul

考えてみれば、elementalとpureは関連が深いとはいえ同値ではないですね。

pure	elemental	例
○	○	各種数学関数
○	×	sumなどのreduce操作関数
×	○	（理想的な）乱数
×	×	状態を保持する擬似乱数、I/O

なおこの表では数学的な定義に基づいて考えていますので、Fortranのimpure elementalではまた別の判断基準（要は仕様書）に基づいて判断したらいいかなと思います。

@cure_honey さん、情報ありがとうございました。

1. v.map(sin)において、mapは高階関数で引数に関数をとっていて、sin関数を引数に渡していますが、Fortranは関数が第一級オブジェクトではないので、mapのような高階関数は使えません。まあ関数ポインタとか使えばなんとかできるんですが、正直めんどい。しかも実数型配列vにメンバ関数を追加することもできないので、まあFortranではできない書きかたかなと。 ↩

2. subroutineではintent(inout)とかintent(out)が可能なようです。確かに、これができないとpure subroutineの存在意義が・・・ ↩