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次元目を指す二つの配列添字の値を入れ替えることで実現します．

```Fortran
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`文や剰余を使わずに，非常に単純な演算で添字を切り替えられます．

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

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

```Fortran
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次元移流方程式

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

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

```math
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コードに落とし込むことができます．

```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`を設けて
```Fortran
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$の値が入っていると見なすことにします．


```Fortran
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$への時間積分も同じ式で記述することができたので，これを繰り返しを用いて記述します．

```Fortran
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`に代入してもいいはずです．

```Fortran
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(:)
```
↓
```Fortran
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)
```

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

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

ダブルバッファリングをしていないプログラムの肝の部分を示します．パラメータは表にまとめています．
```Fortran
    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$が取得できるというわけです．

```Fortran
    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`を入れ替える手続です．
```Fortran
    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`とします．

```Fortran
    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`はポインタ変数が指す配列を入れ替える手続です．
```Fortran
    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とは逆になり，添字で切り替える方が高速になりました．

# まとめ
素直が一番．