配列を返す関数

FORTRAN77 までは、関数の返す値はスカラーのみで配列は返せませんでした。Fortran 90 以降は配列を返せるようになりましたが、返す値を入れる配列として stack に確保される自動変数とすることもできるし、heap に確保される pointer としても返せます。さらに Fortran 2003（正確には Fortran 95 拡張の Technical Report）からは、heap に確保される allocatable な変数として返すこともできるようになりました。allocatable で heap に確保した方が、解放処理を自動でやってくれるので安全性が増します。

ではこれら複数の選択肢のうち、どれを選べばよいのか、少し考えてみました。

stack と heap の別

Fortran では基本的に、サブルーチンや関数内で確保されるスカラー変数や配列、いわゆる自動変数 (automatic variable) は stack に確保されます。Fortran 2008 で導入された block 構文内の局所変数も stack に確保されます。stack はリンク時に静的に確保される一時領域で、普通は単精度 10⁶ 個に満たない容量しかありません。本来、小さい大きさの一時変数を素早く確保・解放するためのものです。

一方 heap は、pointer や allocatable 属性を持つ変数や配列が、必要に応じて OS などに依ってメモリを allocate するもので、大きな領域をもっさりと手間をかけて確保・解放するものです。

しかしながら、一般的に Fortran ユーザーは、stack/heap の使い分けなどを理解しないので、自動変数に非常に大きな配列を取ろうとして stack overflow 実行時エラーが出ることをサポートに苦情を入れるため、最近はコンパイラのオプションでサブルーチンや関数の配列型の自動変数は無条件で heap に確保するスイッチが導入されていたりします。この場合、小さな短寿命の配列を利用する場合余計なオーバーヘッドがかかることになります。

私はとりあえず、stack は手元のメモ帳（素早く取り出せるがいっぱい書けない）、heap はわら半紙（いっぱい書けるがいちいち立ちあがって押し入れに取りに行かねばならぬ）とイメージしてます。

pointer と allocatable の別

allocatable 型の変数は、最近の他言語では unique pointer というものに近い気がします。基本的には単一の unique なメモリー領域を指していて、右辺にあって他の変数に代入する時は deep copy されます。また move_alloc で移動させるときには pointer のように番地がコピーされ元の名前は解放されます。さらに smart pointer 的に、サブルーチンや関数、block 構文内で局所的に宣言された場合は、これらの領域を出る時に自動的に解放処理が行われます。

これに対して pointer 型の変数の場合は古典的なポインタに近く、代入は shallow copy で実体は単一のまま番地だけがコピーされます。また解放処理は明示的に指示せねばならずメモリーリークの危険性が増します。さらに（他言語に比べて制約が非常に強いけれども）自由度の高さがコンパイラの最適化を妨げてしまいます。

しかしながら、関数の返り値として考えた場合、allocatable を用いると heap に確保された結果が deep copy されて余計な時間がかかるのに対し、pointer を用いた場合はアドレスのコピーで済むので無駄が無くなります。

なお stack に返り値を確保した場合も、結果の deep copy が行われます。

実行時間の比較

以下では、stack、heap (allocatable)、heap (pointer) の実行時間の違いを見てみることにします。

結果に有意な差を持たせるために、配列サイズを 10⁸ にしています。これに伴いリンク時に stack size を 1Gbyte ほど取っていますｗ

関数からの結果を受け取る配列は、静的に確保するもの、allocatable で動的に確保するものを用意して比較しています。 pointer を返すものは pointer 番地を受け取っています。

プログラム

    module test_m
    contains
        function fstack(a) result(r)
            real, intent(in) :: a(:)
            real :: r(size(a))
            r = a*a
        end function fstack
        
        function fheap(a) result(r)
            real, intent(in)  :: a(:)
            real, allocatable :: r(:)
            r = a*a
        end function fheap
        
        function fpointer(a) result(r)
            real, intent(in) :: a(:)
            real, pointer :: r(:)
            allocate(r(size(a)))
            r = a*a
        end function fpointer 
        
        subroutine stamp(text)
            character(len=*), intent(in) :: text
            integer :: ic, icrate
            integer, save :: ic0, ic1
            logical, save :: first = .true.
            call system_clock(ic, icrate)
            if (first) then 
                first = .false.
                ic0 = ic
                ic1 = ic
            end if
            print '(2f9.3, 2g0)', real(ic - ic1) / icrate, real(ic - ic0) / icrate, ': ', trim(text)
            ic1 = ic
        end subroutine stamp    
    end module test_m
    
    
    program Console4
        use test_m
        implicit none
        integer, parameter :: n = 10**8
        real :: a(n), b(n), c(n)
        real, allocatable :: e(:), f(:)
        real, pointer :: p(:)
        integer :: i

        
        call stamp('assign to static')
        do i = 1, 3
            print *, '======', i
            call random_number(a)
            call stamp('before fstack')
            b = fstack(a)
            call stamp('after  fstack')
            call stamp('before fheap')
            c = fheap(a)
            call stamp('after  fheap')
        end do

        print *
        call stamp('assign to allocatble')
        do i = 1, 3
            print *, '======', i
            call random_number(a)
            call stamp('before fstack')
            e = fstack(a)
            call stamp('after  fstack')
            call stamp('before fheap')
            f = fheap(a)
            call stamp('after  fheap') 
        end do 

        print *
        call stamp('assign to pointer')
        do i = 1, 3
            print *, '======', i
            call random_number(a)
            call stamp('before fpointer')
            p => fpointer(a)
            call stamp('after  fpointer')
            deallocate(p)
            nullify(p)
        end do 

    
    end program Console4

実行結果

実行結果を見てみますと、結果として返す配列を stack に確保して静的変数に有する場合が最も速くて 0.05 秒くらい、同じく allocatable に返す場合が 0.09 秒くらい、heap に確保して pointer で返す場合が 0.1 秒くらい、heap に確保して静的変数に返す場合が 0.2 秒くらい、同じく allocatable に返す場合が 0.25 秒くらいになっています。これは大体予想通りの順序ですが、思った以上に stack への確保が軽快で、heap への確保がもっさりしているようです。

ただし、よく理解できない点もあって、複数回実行するうち１回目だけは stack/heap に確保した場合は余分な時間がかかっています。何故かよく分かりません。最適化のせいかもしれません。

    0.000    0.000: assign to static
 ======           1
    1.450    1.450: before fstack
    0.195    1.645: after  fstack
    0.000    1.645: before fheap
    0.372    2.017: after  fheap
 ======           2
    1.307    3.324: before fstack
    0.049    3.373: after  fstack
    0.000    3.373: before fheap
    0.214    3.587: after  fheap
 ======           3
    1.318    4.905: before fstack
    0.043    4.948: after  fstack
    0.000    4.948: before fheap
    0.202    5.150: after  fheap

    0.001    5.151: assign to allocatble
 ======           1
    1.288    6.439: before fstack
    0.431    6.870: after  fstack
    0.000    6.870: before fheap
    0.278    7.148: after  fheap
 ======           2
    1.305    8.453: before fstack
    0.088    8.541: after  fstack
    0.001    8.542: before fheap
    0.226    8.768: after  fheap
 ======           3
    1.283   10.051: before fstack
    0.092   10.143: after  fstack
    0.000   10.143: before fheap
    0.240   10.383: after  fheap

    0.001   10.384: assign to pointer
 ======           1
    1.281   11.665: before fpointer
    0.092   11.757: after  fpointer
 ======           2
    1.338   13.095: before fpointer
    0.090   13.185: after  fpointer
 ======           3
    1.341   14.526: before fpointer
    0.096   14.622: after  fpointer