Fortranで配列のダブルバッファリングを実装する

概要

Fortranで配列のダブルバッファリングを行うには，2通りの方法があります．

• 配列のrankを一つ増やして，添字で切り替える
• ポインタを利用する

愚直に配列をコピーする方法，配列のrankを増やす方法，ポインタを使う方法の3通りで実行時間を比較すると，愚直に配列をコピーする方法が最も高速でした．

環境

• Windows 10
• gfortran 10.3
• Intel Fortran classic 2021.5.0
• NAG Fortran 7.1.0

配列のダブルバッファリング

ダブルバッファリングは，動画像処理等の分野において，画面表示のチラつきを低減するために二つのバッファを交互に描画する手法として知られています．本記事で対象とするのはそれではなく，時間発展問題の数値計算において，二つの配列の役割を各時刻で入れ替える手法を意味します．

非常に簡単に言うと，二つの配列を入れ替える方法だと考えてください．

配列のrankを増やして添字で切り替える方法

二つの配列f_1(:), f_2(:)があった場合，f_1(:), f_2(:)の代わりに新しく配列f(:,0:1)を新たに設けます．配列の入れ替えは，2次元目を指す二つの配列添字の値を入れ替えることで実現します．

integer(int32) :: x, y, i
integer(int32), allocatable :: f(:, :)

x = 0; y = 1
allocate (f(5, x:y))

f(:, x) = [(i, i=1, ubound(f, 1))]
f(:, y) = [(i, i=ubound(f, 1), 1, -1)]

print *, f(:, x) ! 1 2 3 4 5
print *, f(:, y) ! 5 4 3 2 1

y = x
x = 1 - x

print *, f(:, x) ! 5 4 3 2 1
print *, f(:, y) ! 1 2 3 4 5

Fortranの配列は，無指定であれば1始まりですが，0始まりとすることで添字の値の入れ替えを簡単に実現できます．
添字の切替の肝はx = 1 - xです．添字は0か1しか取らないことを前提にすると，xが0のとき，x = 1 - xを計算するとx=1が得られます．逆に，xが1のときはx=0となります．if文や剰余を使わずに，非常に単純な演算で添字を切り替えられます．

Fortranの配列は，低次元の要素が連続になる（f(1,0), f(2,0), ...の順に連続になる）ので，配列の切替は最高次元で行います．配列のダブルバッファリングの情報を調べると，C言語系統の情報がよく見つかるわけですが，その通りに実装すると配列要素へのアクセスが不連続になり，著しい性能低下を招いてしまいます．
配列添字の下限を0にするのはC言語系統と同じですが，多次元配列のどの次元で切り替えを行うかは，それらとは異なることに気をつけてください．

ポインタを使う方法

二つの配列f_1(:), f_2(:)があった場合，それぞれにtarget属性を付与するとともに，f_1(:), f_2(:)を指すポインタを新たに設けます．配列の入れ替えは，ポインタを付け替えることで実現します．

integer(int32), allocatable, target :: f_1(:), f_2(:)
integer(int32), dimension(:), pointer :: x, y, swap
integer(int32) :: i

allocate (f_1, source=[(i, i=1, 5)])
allocate (f_2, source=[(i, i=5, 1, -1)])

x => f_1
y => f_2

print *, x ! 1 2 3 4 5
print *, y ! 5 4 3 2 1
swap => x
x => y
y => swap

print *, x ! 5 4 3 2 1
print *, y ! 1 2 3 4 5

配列のダブルバッファリングの使いどころ

1次元移流方程式

\frac{\partial f}{\partial t} + c\frac{\partial f}{\partial x}=0

を，時間に前進差分，空間に中心差分を用いるFTCSスキームで離散化すると，次のような離散式が得られます．

f^{n+1}_{i}=f^{n}_{i} - c\Delta t\frac{f^{n}_{i+1}-f^{n}_{i-1}}{2\Delta x}

$f^{n}, f^{n+1}$ を保持するために配列（例えばf_1, f_2）を宣言したとすると，$f$の時間発展は次のようなFortranコードに落とし込むことができます．

f_2(2:Nx-1) = f_1(2:Nx-1) - dt*c*(f_1(3:Nx) - f_1(1:Nx-2))/(2d0*dx) 

ここで，Nxは配列の要素数，cは$c$，dxは$\Delta x$，dtは$\Delta t$を意味しています．f_1(1), f_1(Nx)は，中心差分で計算できないため，ここでは計算していません．

$t$から$t+\Delta t$への時間積分ができました．f_1には時刻$t$における$f$の値が格納されており，f_2には時刻$t+\Delta t$の$f$の値が格納されることになります．では，$t+\Delta t$から$t+2\Delta t$への時間積分をするにはどうすれば良いでしょうか？

新たにf_3を設けて

f_3(2:Nx-1) = f_2(2:Nx-1) - dt*c*(f_2(3:Nx) - f_2(1:Nx-2))/(2d0*dx) 

と書くのは非常に非効率ですし，時間積分を行う回数だけプログラムが肥大化していきます．

$t$から$t+\Delta t$への時間積分も，$t+\Delta t$から$t+2\Delta t$への時間積分も，原理的には全く同じはずです．そこで，f_1にf_2の値を代入し，f_1には時刻$t+\Delta t$の$f$の値が入っていると見なすことにします．

f_2(2:Nx-1) = f_1(2:Nx-1) - dt*c*(f_1(3:Nx) - f_1(1:Nx-2))/(2d0*dx) ! 時刻t -> t+dt
f_1(:) = f_2(:) ! f_1に時刻t+dtの値を代入

f_2(2:Nx-1) = f_1(2:Nx-1) - dt*c*(f_1(3:Nx) - f_1(1:Nx-2))/(2d0*dx) ! 時刻t+dt -> t+2dt
f_1(:) = f_2(:) ! f_1に時刻t+2dtの値を代入

! 必要な回数繰り返す

$t$から$t+\Delta t$への時間積分も，$t+\Delta t$から$t+2\Delta t$への時間積分も同じ式で記述することができたので，これを繰り返しを用いて記述します．

do n = 1, Nt !Ntは時間積分の回数
    f_2(2:Nx-1) = f_1(2:Nx-1) - dt*c*(f_1(3:Nx) - f_1(1:Nx-2))/(2d0*dx) ! 時刻t -> t+dt
    f_1(:) = f_2(:) ! 次の時間積分の準備
end do
! 必要な回数だけ時間積分を行うと，f_1に欲しい結果が格納されている

時間積分の過程では，f_1からf_2を計算した後，f_2の値をf_1に代入しています．計算に使う数値を考えると，f_1からf_2を計算した後，f_2を使って計算した値をf_1に代入してもいいはずです．

f_2(2:Nx-1) = f_1(2:Nx-1) - dt*c*(f_1(3:Nx) - f_1(1:Nx-2))/(2d0*dx)
f_1(:) = f_2(:)
f_2(2:Nx-1) = f_1(2:Nx-1) - dt*c*(f_1(3:Nx) - f_1(1:Nx-2))/(2d0*dx) ! ここのf_1の数値はf_2と全く同じ．f_2は参照されない．
f_1(:) = f_2(:)

↓

f_2(2:Nx-1) = f_1(2:Nx-1) - dt*c*(f_1(3:Nx) - f_1(1:Nx-2))/(2d0*dx)
f_1(2:Nx-1) = f_2(2:Nx-1) - dt*c*(f_2(3:Nx) - f_2(1:Nx-2))/(2d0*dx)

しかし，これでは時間積分回数が偶数回に制限されますし，同じような式を手動で繰り返すのは，良い作法ではありません．

そこで，ダブルバッファリングを用いて配列を入れ替えることで，配列のコピーを回避しつつ，手動での繰り返しを避けることができます．

ダブルバッファリングをしていないプログラムの肝の部分を示します．パラメータは表にまとめています．

    real(real64), allocatable :: f_curr(:) ! f at current time t
    real(real64), allocatable :: f_next(:) ! f at next time t+dt

    init: block
        integer(int32) :: i
        real(real64), allocatable :: x(:)

        x = [(dble(i - 1)*dx, i=1, Nx)]
        f_curr = merge(((1d0 - cos(2d0*PI*x))/2d0)**2, 0d0, x < 1d0)
        allocate (f_next, source=f_curr)
    end block init

    time_marching: block
        integer(int32) :: n
        real(real64) :: start_s, end_s

        call cpu_time(start_s)
        do n = 1, Nt
            f_next(2:Nx-1) = f_curr(2:Nx-1) - dt*U_conv*(f_curr(3:Nx) - f_curr(1:Nx-2))/(2d0*dx) !&
            f_curr(:) = f_next(:)
        end do
        call cpu_time(end_s)
        print *, end_s - start_s
    end block time_marching

    deallocate (f_curr)
    deallocate (f_next)

変数名	意味	値
Lx	計算領域の長さ	$5$
Nx	計算領域に設けた格子点数	$10001$
dx	格子点間隔	Lx/dble(Nx - 1)
U_conv	移流速度	$1$
Co	Courant数	$0.01$
t_end	計算終了時間	$4$
dt	計算時間間隔	Co*dx/U_conv
Nt	時間積分回数	nint(t_end/dt)

配列のrankを増やして添字で切り替える方法

二つの配列f_curr(:), f_next(:)を，配列f(:,0:1)にまとめます．f(:,0),f(:,1)を参照する配列添字の変数には，f_curr(:), f_next(:)の役割に応じて，curr, nextと付けました．f(:,0)とf(:,1)のどちらに現在の時刻の$f$の値が入っているのかを意識する必要は全くなく，f(:,curr)を参照すれば，現在の時刻の$f$が取得できるというわけです．

    real(real64), allocatable :: f(:, :)
    integer(int32) :: curr, next

    curr = 0; next = 1
    allocate (f(Nx, curr:next), source=0d0)

    init: block
        integer(int32) :: i
        real(real64), allocatable :: x(:)

        x = [(dble(i - 1)*dx, i=1, Nx)]
        f(:, curr) = merge(((1d0 - cos(2d0*PI*x))/2d0)**2, 0d0, x < 1d0)
    end block init

    time_marching: block
        integer(int32) :: n
        real(real64) :: start_s, end_s

        call cpu_time(start_s)
        do n = 1, Nt
            f(2:Nx-1,next) = f(2:Nx-1,curr) - dt*U_conv*(f(3:Nx,curr) - f(1:Nx-2,curr))/(2d0*dx)
            call update(curr, next)
        end do
        call cpu_time(end_s)
        print *, end_s - start_s
    end block time_marching

    deallocate (f)

updateは，配列添字currとnextを入れ替える手続です．

    subroutine update(a, b)
        implicit none
        integer(int32), intent(inout) :: a, b

        b = a
        a = 1 - a
    end subroutine update

ポインタを使う方法

二つの配列f_curr(:), f_next(:)にtarget属性を付与するとともに，f_curr(:), f_next(:)を指すポインタを新たに設けます．ここでは，修正を最小限に抑えるために，f_curr(:), f_next(:)をf_1(:), f_2(:)に変更し，ポインタ変数の名前をf_curr, f_nextとします．

    real(real64), allocatable, target :: f_1(:)
    real(real64), allocatable, target :: f_2(:)

    real(real64), dimension(:), pointer, contiguous :: f_curr
    real(real64), dimension(:), pointer, contiguous :: f_next

    allocate (f_1(Nx), source=0d0)
    allocate (f_2(Nx), source=0d0)
    f_curr => f_1
    f_next => f_2

    init: block
        integer(int32) :: i
        real(real64), allocatable :: x(:)

        x = [(dble(i - 1)*dx, i=1, Nx)]
        f_curr = merge(((1d0 - cos(2d0*PI*x))/2d0)**2, 0d0, x < 1d0)
    end block init

    time_marching: block
        integer(int32) :: n
        real(real64) :: start_s, end_s

        call cpu_time(start_s)
        do n = 1, Nt
            f_next(2:Nx-1) = f_curr(2:Nx-1) - dt*U_conv*(f_curr(3:Nx) - f_curr(1:Nx-2))/(2d0*dx) !&
            call update(f_curr, f_next)
        end do
        call cpu_time(end_s)
        print *, end_s - start_s
    end block time_marching

    f_curr => null()
    f_next => null()
    deallocate (f_1)
    deallocate (f_2)

updateはポインタ変数が指す配列を入れ替える手続です．

    subroutine update(a, b)
        implicit none
        real(real64), dimension(:), pointer, intent(inout) :: a, b
        real(real64), dimension(:), pointer :: s

        s => b
        b => a
        a => s
    end subroutine update

実行時間の比較

配列間のメモリコピーを回避できたので，実行時間の短縮が期待できます．配列をコピーする方法（copy），添字で切り替える方法（array），ポインタを使う方法（pointer）の3種類の実装で実行時間を測定しました．それぞれ10回実行し，最大値と最小値を除いた平均を算出しました．コンパイルオプションには，fpmの--profile releaseオプションを利用した際に与えられるオプションを記載しています．単位はすべて秒です．

	copy [s]	array [s]	pointer [s]	コンパイルオプション
gfortran	1.635	1.674	1.674	-O3 -Wimplicit-interface -fPIC -fmax-errors=1 -funroll-loops -fcoarray=single
ifort	1.582	2.516	1.592	/fp:precise /align:all /error-limit:1 /reentrancy:threaded /nogen-interfaces /assume:byterecl
nagfor	1.654	1.719	1.787	-O4 -coarray=single -PIC

いずれのコンパイラも，配列をコピーする方法が最も速いという結果になりました．添字で切り替える方法とポインタを使う方法のどちらが速いかは，コンパイラによって変わります．

• gfortranはどちらも同じ結果になりました．ただし，結果には出ていませんが，実行時間のばらつきはポインタを使う方法の方が大きく，また除外した最大値もポインタを使う方法の方が大きいという結果でした．
• Intel Fortran classicでは，配列をコピーする方法とポインタを使う方法はほとんど同じです．実行時間のばらつきや除外した最大値については，gfortranと同じ傾向を示しました．添字で切り替える方法は，他のどのコンパイラよりも顕著に遅いという結果になりました．
• NAG Fortranによる実行結果はIntel Fortran classicとは逆になり，添字で切り替える方が高速になりました．

まとめ

素直が一番．