概要
Fortranで配列のダブルバッファリングを実装するの続きです．Fortranでダブルバッファリングを使う場合，いくつか重要な点がありました．
- ダブルバッファリングの対象となる配列を使う処理は，手続にして仮引数の
intent属性を明記する． - 配列のrankを一つ増やして添字で切り替える方法を用いる場合，当該配列を実引数として手続に渡す際に，切替用の添字を明記して二つの配列に分ける．
- ポインタを利用する場合は，ポインタ変数には
contiguous属性を付与する． - ポインタを用いると，どの程度高速化するかがコンパイラに依存する．配列のrankを一つ増やして添字で切り替える方法を用いた方がよい．
詳しい考察は，コンパイラによって生成されたアセンブリコードなどを読み解く必要があります．
環境
- Windows 10
- gfortran 10.3
- Intel Fortran classic 2021.5.0
- NAG Fortran 7.1.0
ダブルバッファリングのうまい使い方を考える．
先の記事（Fortranで配列のダブルバッファリングを実装する）において，Fortranで配列のダブルバッファリングを実現する方法を示しました．ダブルバッファリングを用いると，配列のコピーが不要になるので，実行時間を短縮できると予想していました．しかし，1次元移流方程式の数値計算を対象にダブルバッファリングを実装してみたところ，配列をコピーする方法が最も高速であるという結果が得られました．
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
curr = 0; next = 1
allocate (f(Nx, curr:next), source=0d0)
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
f_curr => f_1
f_next => f_2
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
原因はよくわかりませんが，いくつか性能に関わりそうな要因を考えてみました．
- ポインタ変数に付与している
contiguous属性が何らかの影響を与えているのではないか． - 問題サイズが小さすぎるのではないか．配列サイズがもっと大きい問題であれば，メモリコピーの負荷が相対的に大きくなり，ダブルバッファリングの効果が現れるのではないか．
- 配列次元が影響するのではないか．（そんな訳ないんですが）
- ダブルバッファリングを使うと，コンパイラが最適化をうまくできなくなるのではないか．配列をコピーする方法でも，時間積分を手続に切り出せば，何らかの変化が現れるのではないか．
コンパイラ最適化については，ビルド過程で生成されるアセンブリコードなどから確認できるのかも知れませんが，議論できるだけの十分な知見を有していません．
contiguousの影響
Fortranにおいて，配列を指すポインタは，Fortranのリッチな配列アクセスの表現に対応して，不連続なメモリを指すことができます．contiguous属性を付与すると，当該のポインタの指す配列が連続なメモリを占有していることをコンパイラに知らせることができるので，コンパイラはそれを前提とした最適化を行うことが可能になります．
real(real64), target :: f(1:10)
real(real64), dimension(:), pointer :: f_no_contiguous
! [f(2), f(4), f(6), f(8)]を持つ配列f_no_contiguous(1:4)と同じように振る舞うが，実態は不連続メモリアクセスになる．
f_no_contiguous => f(2:8:2)
先の記事で実装したポインタを利用する方法において，ポインタ変数の宣言の際に付いていたcontiguousを削除しました．
- real(real64), dimension(:), pointer, contiguous :: f_curr
- real(real64), dimension(:), pointer, contiguous :: f_next
+ real(real64), dimension(:), pointer :: f_curr
+ real(real64), dimension(:), pointer :: f_next
contiguousを付与した方が高速になると考える方が自然でしょうが，先の記事で検討していなかったので，どの程度性能に影響があるかを確認するという目的もあります．
それぞれ10回実行し，最大値と最小値を除いた平均の実行時間を算出しました．コンパイルオプションには，fpmの--profile releaseオプションを利用した際に与えられるオプションを記載しています．
contiguousあり [s] |
contiguousなし [s] |
コンパイルオプション | |
|---|---|---|---|
| gfortran | 1.674 | 1.676 | -O3 -Wimplicit-interface -fPIC -fmax-errors=1 -funroll-loops -fcoarray=single |
| ifort | 1.592 | 2.516 | /fp:precise /align:all /error-limit:1 /reentrancy:threaded /nogen-interfaces /assume:byterecl |
| nagfor | 1.787 | 2.143 | -O4 -coarray=single -PIC |
gfortranでは，影響は全くないと言ってよいでしょう．Intel Fortran classicとNAG Fortranではcontiguousを削除した影響が大きく，実行時間は顕著に遅くなりました．
ダブルバッファリングにポインタを使う場合は，contiguous属性は必須と考えるべきでしょう．
問題サイズの影響
先の記事では，10000要素の配列を用いていました．これを大きくしていってもよいのですが，1次元問題ではそこまで大きな配列サイズを使うことはあまりありません．
問題サイズを大きくしても不自然でないようにするために，問題を2次元に変更します．
2次元移流方程式
\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x}+ v\frac{\partial f}{\partial y} =0
を，時間に前進差分，空間に中心差分を用いるFTCSスキームで離散化すると，次のような離散式が得られます．
f^{n+1}_{i,j}=f^{n}_{i,j} - \Delta t\left( u_{i,j}\frac{f^{n}_{i+1,j}-f^{n}_{i-1,j}}{2\Delta x} +v_{i,j}\frac{f^{n}_{i,j+1}-f^{n}_{i,j-1}}{2\Delta y}\right)
2次元移流方程式に対する例題として，Rotating Coneを用います．これは主に移流項の離散化のチェックに用いられる問題です．
$x \in [-1, 1], y\in [-1,1]$の計算領域中に，パッシブスカラ$f=\frac{1}{2}\left[1+\cos(2\pi r)\right], r=\sqrt{x^2+y^2}, r\le 0.5$が分布しています．計算領域には$(u, v) = (-2\pi y, 2\pi x)$で表現される速度が存在しており，その速度によって$f$がメリーゴーラウンドのように移動していくという問題です．
計算条件は下記のように定めました．格子点数は，おおよそ$500\times 500$ということで，先の記事の25倍になっています．Rotating Coneでは$f$が計算領域を1周する様子を計算するので，計算終了時間は通常$1$としますが，計算に要する時間が膨大になるため，$0.01$に縮めました．このように条件を設定すると，1回のプログラム実行が，おおよそ1~2秒で終わります．
| 変数名 | 意味 | 値 |
|---|---|---|
x_min |
計算領域の$x$座標の最小値 | $-1$ |
y_min |
計算領域の$y$座標の最小値 | $-1$ |
Lx$\times$Ly
|
計算領域の長さ | $2\times 2$ |
Nx$\times$Ny
|
計算領域に設けた格子点数 | $501\times 501$ |
dx, dy
|
格子点間隔 |
Lx/dble(Nx-1), Ly/dble(Ny-1)
|
U_conv_max |
移流速度の最大値 | $2\pi$ |
Co |
Courant数 | $0.01$ |
t_end |
計算終了時間 |
|
dt |
計算時間間隔 | Co*dx/U_conv |
Nt |
時間積分回数 | nint(t_end/dt) |
cone_center_x |
$f$の中心位置の$x$座標 | $0$ |
cone_center_y |
$f$の中心位置の$y$座標 | $0.5$ |
初期値の設定，移流方程式の計算部分を示します．
real(real64), allocatable :: u(:, :)
real(real64), allocatable :: v(:, :)
real(real64), allocatable :: f_curr(:, :)
real(real64), allocatable :: f_next(:, :)
init: block
integer(int32) :: i, j, unit
real(real64), allocatable :: x(:, :), y(:, :), r(:, :)
x = reshape([((dble(i - 1)*dx + x_min, i=1, Nx), j=1, Ny)], [Nx, Ny])
y = reshape([((dble(j - 1)*dy + y_min, i=1, Nx), j=1, Ny)], [Nx, Ny])
r = sqrt((x - cone_center_x)**2 + (y - cone_center_y)**2)
f_curr = merge((1d0 + cos(2d0*PI*r))/2d0, 0d0, r <= 0.5d0)
allocate (f_next, source=f_curr)
u = -2d0*PI*y
v = 2d0*PI*x
end block init
time_marching: block
real(real64) :: start_s, end_s
integer(int32) :: n
call cpu_time(start_s)
do n = 1, Nt
f_next(2:Nx-1, 2:Ny-1) = f_curr(2:Nx-1, 2:Ny-1) &
- dt*u(2:Nx-1, 2:Ny-1)*(f_curr(3:Nx, 2:Ny-1) - f_curr(1:Nx-2, 2:Ny-1))/(2d0*dx) &
- dt*v(2:Nx-1, 2:Ny-1)*(f_curr(2:Nx-1, 3:Ny) - f_curr(2:Nx-1, 1:Ny-2))/(2d0*dy)
f_curr(:, :) = f_next(:, :)
end do
call cpu_time(end_s)
print *, end_s - start_s
end block time_marching
deallocate (u)
deallocate (v)
deallocate (f_curr)
deallocate (f_next)
配列のrankを一つ増やして添字で切り替える方法は，下記のように変更されます．
real(real64), allocatable :: u(:, :)
real(real64), allocatable :: v(:, :)
real(real64), allocatable :: f(:, :, :)
integer(int32) :: curr, next
curr = 0; next = 1
allocate (f(Nx, Ny, curr:next), source=0d0)
init: block
integer(int32) :: i, j, unit
real(real64), allocatable :: x(:, :), y(:, :), r(:, :)
x = reshape([((dble(i - 1)*dx + x_min, i=1, Nx), j=1, Ny)], [Nx, Ny])
y = reshape([((dble(j - 1)*dy + y_min, i=1, Nx), j=1, Ny)], [Nx, Ny])
r = sqrt((x - cone_center_x)**2 + (y - cone_center_y)**2)
f(:, :, curr) = merge((1d0 + cos(2d0*PI*r))/2d0, 0d0, r <= 0.5d0)
u = -2d0*PI*y
v = 2d0*PI*x
end block init
time_marching: block
real(real64) :: start_s, end_s
integer(int32) :: n
call cpu_time(start_s)
do n = 1, Nt
f(2:Nx-1, 2:Ny-1, next) = f(2:Nx-1, 2:Ny-1, curr) &
- dt*u(2:Nx-1, 2:Ny-1)*(f(3:Nx, 2:Ny-1, curr) - f(1:Nx-2, 2:Ny-1, curr))/(2d0*dx) &
- dt*v(2:Nx-1, 2:Ny-1)*(f(2:Nx-1, 3:Ny, curr) - f(2:Nx-1, 1:Ny-2, curr))/(2d0*dy)
call update(curr, next)
end do
call cpu_time(end_s)
print *, end_s - start_s
end block time_marching
deallocate (u)
deallocate (v)
deallocate (f)
updateは先の記事と全く同じです．
subroutine update(a, b)
implicit none
integer(int32), intent(inout) :: a, b
b = a
a = 1 - a
end subroutine update
ポインタを用いる方法を用いる場合，下記のようになります．
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, Ny), source=0d0)
allocate (f_2(Nx, Ny), source=0d0)
f_curr => f_1
f_next => f_2
init: block
integer(int32) :: i, j, unit
real(real64), allocatable :: x(:, :), y(:, :), r(:, :)
x = reshape([((dble(i - 1)*dx + x_min, i=1, Nx), j=1, Ny)], [Nx, Ny])
y = reshape([((dble(j - 1)*dy + y_min, i=1, Nx), j=1, Ny)], [Nx, Ny])
r = sqrt((x - cone_center_x)**2 + (y - cone_center_y)**2)
f_curr = merge((1d0 + cos(2d0*PI*r))/2d0, 0d0, r <= 0.5d0)
u = -2d0*PI*y
v = 2d0*PI*x
end block init
time_marching: block
real(real64) :: start_s, end_s
integer(int32) :: n
call cpu_time(start_s)
do n = 1, Nt
f_next(2:Nx-1, 2:Ny-1) = f_curr(2:Nx-1, 2:Ny-1) &
- dt*u(2:Nx-1, 2:Ny-1)*(f_curr(3:Nx, 2:Ny-1) - f_curr(1:Nx-2, 2:Ny-1))/(2d0*dx) &
- dt*v(2:Nx-1, 2:Ny-1)*(f_curr(2:Nx-1, 3:Ny) - f_curr(2:Nx-1, 1:Ny-2))/(2d0*dy)
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)
deallocate (u)
deallocate (v)
updateでは，ポインタの付け替えを行います．
subroutine update(a, b)
implicit none
real(real64), dimension(:, :), pointer, intent(inout) :: a
real(real64), dimension(:, :), pointer, intent(inout) :: b
real(real64), dimension(:, :), pointer :: s
s => b
b => a
a => s
end subroutine update
先の記事と同様に，配列をコピーする方法（copy），添字で切り替える方法（array），ポインタを使う方法（pointer）の3種類の実装で実行時間を測定しました．それぞれ10回実行し，最大値と最小値を除いた平均を算出しました．コンパイルオプションには，fpmの--profile releaseオプションを利用した際に与えられるオプションを記載しています．Intel Fortranについては，stack overflowが出たので，1 MBを越えるメモリはheapに置くようにオプションを追加しました．ポインタを使う方法では，contiguous属性を付与しています．
| copy [s] | array [s] | pointer [s] | コンパイルオプション | |
|---|---|---|---|---|
| gfortran | 1.305 | 2.174 | 2.197 | -O3 -Wimplicit-interface -fPIC -fmax-errors=1 -funroll-loops -fcoarray=single |
| ifort | 1.359 | 1.596 | 2.203 | /fp:precise /align:all /error-limit:1 /reentrancy:threaded /nogen-interfaces /assume:byterecl /heap-arrays:1000 |
| nagfor | 1.301 | 1.350 | 1.416 | -O4 -coarray=single -PIC |
先の記事と結果の傾向は類似しています．問題サイズも配列の次元も影響は無さそうです．1次元問題の時との違いとしては，Intel Fortranでポインタを用いる方法が顕著に遅くなったことでしょう．また，1次元の時に最も高速だったIntel Fortranが一番遅くなっています．heapを使うように指定したことが影響している可能性があるものの，これ以外にstack overflowを回避する手立てがないので，仕方ありません．もしかしたら，速度が定数から配列に変化したことも影響があるかもしれません．
時間積分を手続に切り出す
時間積分を行っている箇所
f_next(2:Nx-1, 2:Ny-1) = f_curr(2:Nx-1, 2:Ny-1) &
- dt*u(2:Nx-1, 2:Ny-1)*(f_curr(3:Nx, 2:Ny-1) - f_curr(1:Nx-2, 2:Ny-1))/(2d0*dx) &
- dt*v(2:Nx-1, 2:Ny-1)*(f_curr(2:Nx-1, 3:Ny) - f_curr(2:Nx-1, 1:Ny-2))/(2d0*dy)
を手続integrate_convection_equationに切り出します．この手続は，配列をコピーする方法，添字で切り替える方法，ポインタを使う方法の全てで使い回すことができます．
subroutine integrate_convection_equation(f_curr, u, v, f_next)
implicit none
real(real64), intent(in) :: f_curr(:, :)
real(real64), intent(in) :: u(:, :), v(:, :)
real(real64), intent(inout) :: f_next(:, :)
f_next(2:Nx-1, 2:Ny-1) = f_curr(2:Nx-1, 2:Ny-1) &
- dt*u(2:Nx-1, 2:Ny-1)*(f_curr(3:Nx, 2:Ny-1) - f_curr(1:Nx-2, 2:Ny-1))/(2d0*dx) &
- dt*v(2:Nx-1, 2:Ny-1)*(f_curr(2:Nx-1, 3:Ny) - f_curr(2:Nx-1, 1:Ny-2))/(2d0*dy)
end subroutine integrate_convection_equation
do n = 1, Nt
call integrate_convection_equation(f_curr, u, v, f_next)
f_curr(:, :) = f_next(:, :)
end do
do n = 1, Nt
call integrate_convection_equation(f(:, :, curr), u, v, f(:, :, next))
call update(curr, next)
end do
do n = 1, Nt
call integrate_convection_equation(f_curr, u, v, f_next)
call update(f_curr, f_next)
end do
添字で切り替える方法は，f(:,:,curr)とf(:,:,next)のように添字を指定することで，配列を仮想的に二つに分けています．
ところで，添字で切り替える方法に対しては，integrate_convection_equationの別の実装を考えることができます．配列fを添字無しで手続に渡し，同時に添字curr, nextを渡して手続の中で現在時刻と次の時刻を指定する方法です．
subroutine integrate_convection_equation(f, u, v, curr, next)
implicit none
real(real64), intent(inout) :: f(1:Nx, 1:Ny, 0:1)
real(real64), intent(in) :: u(:, :), v(:, :)
integer(int32), intent(in) :: curr, next
f(2:Nx-1, 2:Ny-1,next) = f(2:Nx-1, 2:Ny-1,curr) &
- dt*u(2:Nx-1, 2:Ny-1)*(f(3:Nx, 2:Ny-1,curr) - f(1:Nx-2, 2:Ny-1,curr))/(2d0*dx) &
- dt*v(2:Nx-1, 2:Ny-1)*(f(2:Nx-1, 3:Ny,curr) - f(2:Nx-1, 1:Ny-2,curr))/(2d0*dy)
end subroutine integrate_convection_equation
do n = 1, Nt
call integrate_convection_equation(f, u, v, curr, next)
call update(curr, next)
end do
上の説明と同じ条件（コンパイルオプション，平均の算出）で実行時間を測定しました．配列をコピーする方法（copy），添字で切り替える方法（array），ポインタを使う方法（pointer）に加え，添字で切り替える方法において，添字なしで手続に渡す実装（array wo index）の実行時間を測定しました．また，参考のため，問題サイズの影響で示した表内の，配列をコピーする方法の（時間積分の処理を手続にしていない）結果（naive copy）も掲載します．
| copy [s] | array [s] | pointer [s] | array wo index [s] | naive copy [s] | |
|---|---|---|---|---|---|
| gfortran | 1.307 | 1.043 | 1.035 | 2.170 | 1.305 |
| ifort | 1.305 | 1.041 | 1.301 | 1.596 | 1.359 |
| nagfor | 1.473 | 1.049 | 1.207 | 1.367 | 1.301 |
コンパイラによって明確な違いが現れました．
以降の考察は，あくまでテスト問題を実行した結果からの推察であり，正確でない可能性があります．
gfortran
gfortranでは，配列をコピーする方法に対しては，時間積分の処理を手続に切り出しても実行時間は変わりませんでした．一方で，添字で切り替える方法，ポインタを使う方法共に高速化しており，ダブルバッファリングの導入の効果が見られました．
時間積分の処理を手続に切り出しても実行時間が変わらなかったのは，gfortranが読み取りのみの変数と変数と書き込まれる変数をうまく区別して最適化をしていたのだろうと考えられます．ダブルバッファリングを導入した場合，時間積分の処理を手続に切り出す前は，読み取りのみの変数と変数と書き込まれる変数を区別できずに最適化したコードを生成できなかったが，手続に切り出すことでintentの指定によって最適化が可能になり，加えてメモリコピーを回避したことによって高速化したと考えられます．
添字で切り替える方法において，添字なしで手続に渡す実装は遅くなりました．上の推察（高速化の要因は，intentの指定によって読み取りのみの変数と変数と書き込まれる変数を区別した最適化が行われた）が正しいとすると，添字なしで手続に渡す実装では，f(:,:,curr)とf(:,:,next)のどちらもintent(inout)に指定せざるを得ず，読み取りのみの変数と変数と書き込まれる変数を区別した最適化が行えなかったため，遅くなったと考えられます．
Intel Fortran
Intel Fortranでは，全ての場合に実行時間が短縮されました．
配列をコピーする方法に対しては，時間積分の処理を手続に切り出すと実行時間が短縮しており，仮引数のintent指定が有効に作用したものと考えられます．問題サイズの影響の考察で，
速度が定数から配列に変化したことも影響があるかもしれません．
と書きましたが，あながち的外れではなかったのかもしれません．
ポインタを用いる方法の実行時間は配列をコピーする方法とほぼ同じであり，ダブルバッファリングの効果が得られていません．2.203 sから1.301 sまで短縮されただけでもヨシとする必要があるでしょうか．Intel Fortranのポインタは，高速化目的で使うのはお薦めできないかもしれません．
添字で切り替える方法は無事に高速化されており，メモリコピーを回避した効果が現れています．高速化の理由は，gfortranと同じくintentを指定したことで最適化されたコードを生成できたためと考えられます．
添字で切り替える方法において，添字なしで手続に渡す実装の実行時間は，時間積分の処理を手続に切り出さない実装（問題サイズの影響）と同じです．つまり，添字で切り替える方法において，時間積分の処理を手続に切り出さない実装と添字なしで手続に渡す実装は，共に読み取りのみの変数と書き込まれる変数を区別した最適化が行われていないコードを出力しているのだと考えられます．
NAG Fortran
NAG Fortranでは，速くなる場合と遅くなる場合が見られました．
配列をコピーする方法に対しては，時間積分の処理を手続に切り出すと実行時間が伸びてしまいました．時間積分の処理と配列のコピーを連続して書くことで，特別な最適化を行っていたのでしょうか？これについては理由がわかりません．
ダブルバッファリングを導入すると高速化されます．特に，添字で切り替える方法は，他のコンパイラとほぼ同じ実行時間になっています．高速化の要因は，他のコンパイラと同じくメモリコピーの回避とintentの指定による最適化だと考えられます．
ポインタを用いる方法では，配列をコピーする方法よりは高速であるものの，添字で切り替える方法ほど高速ではありません．この時間の差がどこから来るのでしょうか？
他のコンパイラの結果において，配列をコピーする方法と添字で切り替える方法の実行時間の差は0.264 sです．この時間が配列のコピーに要する時間であると考え，NAG Fortranの結果に当てはめると，1.473-0.264=1.209となって，ポインタを用いる方法の実行時間とほぼ一致します．
つまり，配列をコピーする方法とポインタを用いる方法では，時間積分の処理に要する時間は同じで，メモリコピーを回避した分だけポインタを用いる方法が高速化したと結論付けられます．そうすると，添字で切り替える方法だけで特別に高速な（最適化された）時間積分手続が生成されていることになります．ますます判らなくなってきました･･･
添字で切り替える方法において，添字なしで手続に渡す実装の実行時間は，Intel Fortranと同じく時間積分の処理を手続に切り出さない実装（問題サイズの影響）と同じです．
まとめ
Fortranでダブルバッファリングが効果を発揮する書き方を模索しました．その結果，ダブルバッファリングをうまく使うためには，以下のことを意識する必要があるとわかりました．
- ダブルバッファリングの対象となる配列を使う場合，処理は手続として実装し，仮引数の
intent属性を明記する． - 添字で切り替える方法を使っている場合，実引数として渡す時は添字を明記する．（
intent(in)の配列と，intent(inout)の配列に仮想的に分ける）
f(:, :, next) = f(:, :, curr) - ...
↓
call do_something(f(:, :, curr), ..., f(:, :, next))
subroutine do_something(f, ..., g)
real(real64), intent(in) :: f(:, :)
real(real64), intent(inout) :: g(:, :)
end subroutine do_something
- ポインタを利用する場合は，ポインタ変数には必ず
contiguous属性を付与する．
上記を考慮すると，
- 添字で切り替える方法が最も高速である．
- ポインタを利用する方法は，コンパイラによってどの程度高速化されるかが変化する．最良で添字で切り替える方法と同じ，最悪で配列をコピーする方法と同じ．