Fortranにおける配列のダブルバッファリングにdo concurrentを使う

概要

Fortranで配列のダブルバッファリングを実装する，Fortranにおける配列のダブルバッファリングのうまい使い方を調査するの続きです．

ダブルバッファリングにポインタを利用する方法＋配列を用いる処理にdo concurrentを用いる組合せが，手軽に高速化の恩恵を受けられることがわかりました．

タイトルが「ダブルバッファリングに」となっていますが，正確には「ダブルバッファリングの対象となる配列を使った処理に」です．

環境

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

ダブルバッファリングする配列を使った計算にdo concurrentを使ってみる

Fortranで配列のダブルバッファリングを実装するにおいて，配列のダブルバッファリングを行うと実行時間が遅くなることを報告し，Fortranにおける配列のダブルバッファリングのうまい使い方を調査するでその原因を考察しました．

ダブルバッファリングを導入すると，配列を使った処理の中で読み取りのみの変数と変数と書き込まれる変数を区別できず，コンパイラが最適化しなくなると推察しました．そして，それを回避するためには，配列を使った処理をサブルーチンとして切り出すことを示しました．

Fortranにおける配列のダブルバッファリングのうまい使い方を調査するに対して，@cure_honey さんから

do concurrent で左辺が右辺と独立していることが伝えられそうなので、@cometscome_phys さんの記事のコードを intel fortran で試してみたところ、それっぽい結果が出て subroutine を読んだ場合と同程度の時間になりました。

とのコメントが寄せられました．

本記事でも，do concurrentを使って処理を記述することで，ダブルバッファリングを導入したコードにどの程度影響があるかを調査しました．

1次元移流方程式

知見の反映

Fortranにおける配列のダブルバッファリングのうまい使い方を調査するで得た，ダブルバッファリングを導入する際の知見は，事の発端となった1次元移流方程式（Fortranで配列のダブルバッファリングを実装する）には適用していませんでした．

do concurrentを使った場合に，1次元移流方程式の検証が不十分だった（ダブルバッファリングで高速化したのではなくdo cuncurrentの影響で高速化した）と指摘されるのを避けるために，1次元移流方程式でもダブルバッファリングの効果があることを確認します．

1次元移流方程式のプログラムを知見に従って書き換え，高速化することを確認します．
詳しくはこれまでの記事を見ていただくとして，基準となるコードスニペットを示しておきます．このコードでは，時間積分の処理を配列式でベタ書きし，配列の入れ替えに単純なコピーを利用しています．

naive/copy

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

先の記事（Fortranにおける配列のダブルバッファリングのうまい使い方を調査する）では，ダブルバッファリングを導入する際は下記の点に気をつけるよう述べていました．

• 配列を使う処理は手続とせよ
• 手続の引数にはintentを指定せよ
• ポインタを使う場合はcontiguousを指定せよ

それに従って，1次元移流方程式のプログラムを修正します．まずは時間積分の処理を手続に切り出します．

subroutine integrate_convection_equation(f_curr, u, f_next)
    implicit none
    real(real64), intent(in) :: f_curr(:)
    real(real64), intent(in) :: u
    real(real64), intent(inout) :: f_next(:)

    f_next(2:Nx-1) = f_curr(2:Nx-1) - dt*u*(f_curr(3:Nx) - f_curr(1:Nx-2))/(2d0*dx)
end subroutine integrate_convection_equation

ダブルバッファリングの速度を調査するために利用した3個のプログラム（copy, array, poiter）を，時間積分手続を呼ぶように変更します．

proc/copy

do n = 1, Nt
    call integrate_convection_equation(f_curr, U_conv, f_next)
    f_curr(:) = f_next(:)
end do

proc/array

do n = 1, Nt
    call integrate_convection_equation(f(:, curr), U_conv, f(:, next))
    call update(curr, next)
end do

proc/pointer

do n = 1, Nt
    call integrate_convection_equation(f_curr(:), U_conv, f_next(:))
    call update(f_curr, f_next)
end do

実行時間の比較

実行時間を比較するにあたって，時間積分の処理を配列式でベタ書きし，配列の入れ替えに単純なコピーを利用する方法（naive/copy）を基準とします．比較するプログラムでは，手続に切り出した時間積分の処理を呼出しており，配列の入れ替え方法だけが異なります．配列をコピーする方法（proc/copy），添字で切り替える方法（proc/array），ポインタを使う方法（proc/pointer）を用いています．測定方法はこれまでの記事と同様に各プログラムを10回実行し，最大値と最小値を除いた平均を算出しています．コンパイルオプションには，fpmの--profile releaseオプションを利用した際に与えられるオプションを記載しています．

	proc/copy [s]	proc/array [s]	proc/pointer [s]	naive/copy [s]	コンパイルオプション
gfortran	1.631	1.268	1.262	1.635	-O3 -Wimplicit-interface -fPIC -fmax-errors=1 -funroll-loops -fcoarray=single
ifort	1.572	1.287	1.289	1.582	/fp:precise /align:all /error-limit:1 /reentrancy:threaded /nogen-interfaces /assume:byterecl
nagfor	1.625	1.287	1.33	1.654	-O4 -coarray=single -PIC

配列をコピーする方法であっても，配列を使う時間積分の処理を手続に切り出すことで，僅かですが高速化する傾向が見られました．添字で切り替える方法，ポインタを使う方法は，先の記事の通り高速化しました．ダブルバッファリングが1次元移流方程式（1次元配列）に対しても有効であることが確認できました．

do concurrentを使う

ダブルバッファリングの有効性が確認できたので，1次元移流方程式の計算にdo concurrentを導入します．do concurrentはFortran 2008から導入された構文です．その名前から，並列に実行されるdo文だと思われるかもしれませんが，正確には構文内に書ける内容に制限のあるdo文です．do concurrentを用いると，構文内のコードがその制限に沿っているため，自動並列化も含めてより高い最適化を行ってもよいとコンパイラに知らせることができます．つまり，do concurrentは，並列に実行しても問題無いdo文だと捉えることができます．

do concurrentを使うと，構文内でのcycyleやexit，stopなどフロー制御が許されません．また，各反復における変数の更新などに，順序の依存性がないことが前提になります¹．

以前の記事では，ダブルバッファリングを導入するには，配列を使う処理を手続に切り出す必要があることを見てきました．この理由は，intentを指定して，更新される配列と更新されない配列をコンパイラに教えるためでした．do concurrentを使うと，少なくともある反復において，右辺と左辺が独立であることを伝えられるので，処理を手続に切り出すという作業を回避できる可能性があります．

ところで，これまで1次元移流方程式の時間積分の処理には，配列式を用いてきました．do concurrentを用いると，配列式ではなくいわゆる繰り返し処理を用いることになります．do concurrentの影響ではなく，配列式を繰り返し処理に変更した影響の方が大きいのでは？という指摘を受けそうなので，do concurrentだけでなく，通常のdo文を用いるプログラムも作成し，実行時間を調査します．

1次元移流方程式のソースは，それぞれ下記のように変更しました．

do/copy

do n = 1, Nt
    do i = 2, Nx - 1
        f_next(i) = f_curr(i) - dt*U_conv*(f_curr(i+1) - f_curr(i-1))/(2d0*dx)
    end do
    f_curr(:) = f_next(:)
end do

do/array

do n = 1, Nt
    do i = 2, Nx - 1
        f(i,next) = f(i,curr) - dt*U_conv*(f(i+1,curr) - f(i-1,curr))/(2d0*dx)
    end do
    call update(curr, next)
end do

do/pointer

do n = 1, Nt
    do i = 2, Nx - 1
        f_next(i) = f_curr(i) - dt*U_conv*(f_curr(i+1) - f_curr(i-1))/(2d0*dx)
    end do
    call update(f_curr, f_next)
end do

do concurrent/copy

do n = 1, Nt
    do concurrent(i=2:Nx-1)
        f_next(i) = f_curr(i) - dt*U_conv*(f_curr(i+1) - f_curr(i-1))/(2d0*dx)
    end do
    f_curr(:) = f_next(:)
end do

do concurrent/array

do n = 1, Nt
    do concurrent(i=2:Nx-1)
        f(i,next) = f(i,curr) - dt*U_conv*(f(i+1,curr) - f(i-1,curr))/(2d0*dx)
    end do
    call update(curr, next)
end do

do concurrent/pointer

do n = 1, Nt
    do concurrent(i=2:Nx-1)
        f_next(i) = f_curr(i) - dt*U_conv*(f_curr(i+1) - f_curr(i-1))/(2d0*dx)
    end do
    call update(f_curr, f_next)
end do

実行時間の比較

実行時間の測定では，Fortranで配列のダブルバッファリングを実装するで測定した3個のプログラム（copy, array, pointer，時間積分の処理には配列式を用い，手続ではなくベタ書き）を基準として，配列式からdo文に変更した影響（手続），do文をdo concurrentに変更した影響を確認します．do concurrentを用いるにあたって，コンパイラが並列化したコードを出力しても1スレッドで実行されるよう，環境変数OMP_NUM_THREADSを1に設定しました（これはあくまで念のためで，OpenMPを有効にするコンパイルオプションを付けないと，並列化はされません）．

	copy [s]	array [s]	pointer [s]	do/copy [s]	do/array [s]	do/pointer [s]	do concurrent/ copy [s]	do concurrent/ array [s]	do concurrent/ pointer [s]	コンパイルオプション
gfortran	1.635	1.674	1.674	1.646	1.277	1.270	1.641	1.273	1.271	-O3 -Wimplicit-interface -fPIC -fmax-errors=1 -funroll-loops -fcoarray=single
ifort	1.582	2.516	1.592	1.574	2.471	1.275	1.578	1.273	1.270	/fp:precise /align:all /error-limit:1 /reentrancy:threaded /nogen-interfaces /assume:byterecl
nagfor	1.654	1.719	1.787	1.764	1.632	1.336	1.633	1.398	1.334	-O4 -coarray=single -PIC

まず，配列式をdo文に変更した影響（copy, array, pointerとdo/copy, do/array, do/pointer）を見てみます．

配列をコピーする方法を用いた場合，gfortran, Intel Fortranは配列式をdo文に変更した影響は僅かです．一方で，NAG Fortranは，do文に変更すると実行時間が長くなりました．NAG Fortranは，配列式の最適化が上手いのかもしれません．

添字で切り替える方法を用いると，gfortranは大幅に高速化しました．この実行時間は，前小節の結果（proc/array）とほとんど同じ（むしろそれよりも高速）であることから，ダブルバッファリングの効果が得られていると判断できます．Intel FortranおよびNAG Fortranは，do文に変更する事で僅かに高速化しました．

ポインタを利用する方法を用いると，いずれのコンパイラも大幅に高速化し，ダブルバッファリングの効果が得られるようになりました．

ここまでの比較によって，まずは以下のことが判りました．

• 配列をコピーする方法を用いる場合，配列式とdo文のどちらを使っても実行時間はほとんど変わらない．
• 添字で切り替える方法を用いる場合，配列式よりはdo文を用いる方が有利である．コンパイラによってはダブルバッファリング導入による高速化が達成できる．
• ポインタを利用する方法を用いる場合，配列式を用いるとダブルバッファリング導入の影響は全くないが，do文を用いるとダブルバッファリング導入によって高速化される．

次に，do文をdo concurrentに変更した効果（do/copy, do/array, do/pointerとdo concurrent/copy, do concurrent/array, do concurrent/pointer）を確認します．

配列をコピーする方法を用いた場合，gfortran, Intel Fortranの実行時間に変化はほとんどありません．NAG Fortranは若干高速化していますが，これは配列式をdo文に置き換えた際の実行時間の増加分が短くなっただけです．配列式を用いた実装（copy）と実行時間はほぼ同じです．

添字で切り替える方法では，gfortranは（do文の導入によって高速化していたので）実行時間に変化はありませんでした．少なくとも，do concurrentにしたから遅くなるということはないようです．Intel Fortranは，do concurrentにすることで大幅に高速化し，ダブルバッファリング導入の効果を得られるようになりました．@cure_honey さんの指摘通り，do concurrentを用いる事で左辺と右辺の独立性をコンパイラに示すことに成功していると言えるでしょう．NAG Fortranも高速化はしていますが，他の2コンパイラよりは時間がかかっています．先の記事でもそうでしたが，NAG Fortranは最適化の仕方が異なるのでしょうか．

ポインタを利用する方法の場合は，全てのコンパイラで実行時間はほぼ同じでした．そもそも配列式からdo文に切り替えた時点で大幅に高速化していたので，do concurrentに切り替えるペナルティが無いことは示せたでしょう．

1次元のまとめ

1次元移流方程式の計算にダブルバッファリングを導入するにあたって，先の記事の知見（配列を使う処理を手続に切り出す，ポインタを使う場合はcontiguous属性を付与する）を反映しました．その結果，これらの知見は1次元配列に対しても有効でした．

一方，配列を使う処理を，配列式ではなくdo文で記述する場合，先の記事の知見とは異なる結果になりました．

• 配列を使う処理をdo文で記述する場合，ポインタを利用する方法であれば，処理を手続に切り出さなくてもダブルバッファリングによって高速化される．これは今回利用した全てのコンパイラで有効である．
• do文をdo concurrent文に変更することで，配列を使う処理を手続に切り出さなくてもダブルバッファリング導入の効果を得られる場合がある．
  • Intel Fortranの場合，添字で切り替える方法に対しては極めて有効である．
• do concurrent文を用いるペナルティはない．

2次元移流方程式

1次元移流方程式を用いた検証によって，do concurrentはダブルバッファリングとの親和性が高いことが示されました．2次元移流方程式でも有効と期待されるので，do concurrentを導入してみます．

1次元移流方程式と同様，2次元移流方程式も配列式を用いて実装していたので，まずは配列式をdo文を用いた繰り返し処理に変更し，その後do concurrentに変更します．

先の記事で作成・測定した結果（手続に切り出した時間積分の処理を呼出し，配列をコピーする方法（copy），添字で切り替える方法（array），ポインタを使う方法（pointer）を用いたプログラムの実行時間）を比較対象として，配列を使う処理（時間積分処理）を手続とはせず，また配列式ではなくdo文，do concurrent文による繰り返し処理で記述したプログラムの実行時間を測定します．

比較対象のプログラムの肝の部分は下記のようになっています．

copy

do n = 1, Nt
    call integrate_convection_equation(f_curr, u, v, f_next)
    f_curr(:, :) = f_next(:, :)
end do

array

do n = 1, Nt
    call integrate_convection_equation(f(:, :, curr), u, v, f(:, :, next))
    call update(curr, next)
end do

pointer

do n = 1, Nt
    call integrate_convection_equation(f_curr, u, v, f_next)
    call update(f_curr, f_next)
end do

配列を使う処理（時間積分処理）を手続とはせず，do文による繰り返し処理で記述したプログラムをdo/copy（配列をコピーする方法），do/array（添字で切り替える方法），do/pointer（ポインタを使う方法）と称し，do concurrent文による繰り返し処理で記述したプログラムをそれぞれdo concurrent/copy, do concurrent/array, do concurrent/pointerと称します．

do/copy

do n = 1, Nt
    do j = 2, Ny-1
    do i = 2, Nx-1
        f_next(i,j) = f_curr(i,j) &
                        - dt*u(i,j)*(f_curr(i+1, j) - f_curr(i-1, j))/(2d0*dx) &
                        - dt*v(i,j)*(f_curr(i, j+1) - f_curr(i, j-1))/(2d0*dy) !&
    end do
    end do
    f_curr(:, :) = f_next(:, :)
end do

do/array

do n = 1, Nt
    do j = 2, Ny-1
    do i = 2, Nx-1
        f(i,j,next) = f(i,j,curr) &
                        - dt*u(i,j)*(f(i+1, j,curr) - f(i-1, j,curr))/(2d0*dx) &
                        - dt*v(i,j)*(f(i, j+1,curr) - f(i, j-1,curr))/(2d0*dy) !&
    end do
    end do
    call update(curr, next)
end do

do/pointer

do n = 1, Nt
    do j = 2, Ny-1
    do i = 2, Nx-1
        f_next(i,j) = f_curr(i,j) &
                        - dt*u(i,j)*(f_curr(i+1, j) - f_curr(i-1, j))/(2d0*dx) &
                        - dt*v(i,j)*(f_curr(i, j+1) - f_curr(i, j-1))/(2d0*dy) !&
    end do
    end do
    call update(f_curr, f_next)
end do

do concurrent/copy

do n = 1, Nt
    do concurrent(j=2:Ny-1, i=2:Nx-1)
        f_next(i,j) = f_curr(i,j) &
                        - dt*u(i,j)*(f_curr(i+1, j) - f_curr(i-1, j))/(2d0*dx) &
                        - dt*v(i,j)*(f_curr(i, j+1) - f_curr(i, j-1))/(2d0*dy) !&
    end do
    f_curr(:, :) = f_next(:, :)
end do

do concurrent/array

do n = 1, Nt
    do concurrent(j=2:Ny-1, i=2:Nx-1)
        f(i,j,next) = f(i,j,curr) &
                        - dt*u(i,j)*(f(i+1, j,curr) - f(i-1, j,curr))/(2d0*dx) &
                        - dt*v(i,j)*(f(i, j+1,curr) - f(i, j-1,curr))/(2d0*dy) !&
    end do
    call update(curr, next)
end do

do concurrent/pointer

do n = 1, Nt
    do concurrent(j=2:Ny-1, i=2:Nx-1)
        f_next(i,j) = f_curr(i,j) &
                        - dt*u(i,j)*(f_curr(i+1, j) - f_curr(i-1, j))/(2d0*dx) &
                        - dt*v(i,j)*(f_curr(i, j+1) - f_curr(i, j-1))/(2d0*dy) !&
    end do
    call update(f_curr, f_next)
end do

do concurrentを用いるにあたって，コンパイラが並列化したコードを出力しても1スレッドで実行されるよう，環境変数OMP_NUM_THREADSを1に設定しました．

	copy [s]	array [s]	pointer [s]	do/copy [s]	do/array [s]	do/pointer [s]	do concurrent/ copy [s]	do concurrent/ array [s]	do concurrent/ pointer [s]	コンパイルオプション
gfortran	1.307	1.043	1.035	1.270	1.004	1.037	1.279	1.000	1.029	-O3 -Wimplicit-interface -fPIC -fmax-errors=1 -funroll-loops -fcoarray=single
ifort	1.305	1.041	1.301	1.291	1.584	1.625	1.270	1.008	1.020	/fp:precise /align:all /error-limit:1 /reentrancy:threaded /nogen-interfaces /assume:byterecl /heap-arrays:1000
nagfor	1.473	1.049	1.207	1.881	1.611	1.639	1.273	0.998	1.117	-O4 -coarray=single -PIC

配列式＋手続と，ベタ書きのdo文（copy, array, pointerとdo/copy, do/array, do/pointer）の結果を比較します．

配列をコピーする方法を用いた場合，gfortran, Intel Fortranは，配列式＋手続をベタ書きのdo文に変更しても実行時間はほとんど変わりませんでした．一方で，NAG Fortranは，手続をベタ書きのdo文に変更すると実行時間が長くなりました．これは，do文をベタ書きしたことで読み取りのみの変数と変数と書き込まれる変数が区別できなくなったこと，1次元の時と同様に配列式とdo文の最適化の違いが関係していると推察しています（つまりは，よくわからないということです）．

添字で切り替える方法を用いると，gfortranは大幅に高速化し，配列式＋手続で書くよりも速くなりました．Intel FortranとNAG Fortranは，配列式＋手続よりもベタ書きのdo文の方が遅くなっています．これは，先の記事で配列式をベタ書きしていた時と同じ傾向なので，読み取りのみの変数と変数と書き込まれる変数を区別した最適化ができないことが，速度低下の要因だと考えられます．

ポインタを利用する方法の結果は，添字で切り替える方法の結果よりも若干時間がかかるものの，同じ傾向を示します．

次に，do文をdo concurrentに変更した影響（do/copy, do/array, do/pointerとdo concurrent/copy, do concurrent/array, do concurrent/pointer）を見てみます．

配列をコピーする方法を用いた場合，gfortran, Intel Fortranの実行時間に変化はほとんどありません．NAG Fortranは，顕著に高速化し，他の2コンパイラと同じ程度の時間になりました．

添字で切り替える方法では，全てのコンパイラで高速化が確認され，ポインタを用いる方法（do/pointer）と同じ程度まで高速化しました．do concurrent文を用いると，配列を使う処理をベタ書きをしていたとしても，右辺と左辺の独立性が示されるので高速化することが判りました．また，配列式で書くよりも40 msほど高速です．これは計算に用いる点の参照方法の違い（特に(f_curr(2:Nx-1, 3:Ny) - f_curr(2:Nx-1, 1:Ny-2)と(f_curr(i, j+1) - f_curr(i, j-1)）に起因していると予想しています．

ポインタを利用する方法の場合は，doからdo concurrentに切り替えることで高速化が達成されました．1次元の場合，do concurrent文を用いるペナルティはないと結論していましたが，2次元の場合も，do concurrent文を用いるペナルティは無いと言い切ってもよいでしょう．
添字で切り替える方法と比較すると，少し遅いという結果でした．速度を追求する場合は，ポインタを利用する方法ではなく添字で切り替える方法を利用するべきでしょう．

まとめ

Fortranでダブルバッファリングを使う際のうまいやり方を模索し，do concurrentを使ってみました．do concurrentを使うことで，式の右辺と左辺の独立性をコンパイラに伝えることができるようになり，手続に切り出さなくてもダブルバッファリングによる高速化（処理が十分に最適化されることに加えて，配列入れ替えのためのメモリコピーの削減）が期待できます．

do concurrentを使う前に，配列式で書いていた処理を全てdoを使った繰り返し処理に書き換えました．配列次元を一つ増やして添字で切り替える方法とポインタを利用する方法は，一部のコンパイラでは高速化しましたが，ダブルバッファリング導入の効果は得られないと考えた方が安全です．

doをdo concurrentに変更すると，添字で切り替える方法とポインタを利用する方法の両方で高速化しました．do concurrentを用いたことで，式の右辺と左辺の独立性をコンパイラに伝えることができ，読み取りのみの変数と変数と書き込まれる変数を区別した最適化が行われたと考えられます．その最適化に加え，ダブルバッファリングによってメモリコピーが削減されたことが高速化の理由です．

ポインタを利用する方法は，添字で切り替える方法と比較すると，少し遅いという結果でした．

これらの情報をまとめると，ダブルバッファリングを導入するには下記のように検討すればよいということでしょう．

• 処理を配列式で書いており，あまり大きな変更は入れたくない→処理を手続に切り出し，ポインタを利用してダブルバッファリングを行う．
  • ただし，高速化の効果はあまり大きくない．
• 処理を配列式で書いており，大きな変更を入れてもいい→処理を手続に切り出し，配列次元を一つ増やして添字で切り替えることでダブルバッファリングを行う．
  • 高速化の効果は大きい．
• 処理をdo文で書いており，あまり大きな変更は入れたくない→doをdo concurrentに変更し，ポインタを利用してダブルバッファリングを行う．
  • 手軽で，かつ高速化の効果も大きい．ただし，do concurrentの制約に引っかかるような処理を書いている場合は，対応できない．
• 処理をdo文で書いており，大きな変更を入れてもよい→処理を手続に切り出し，配列次元を一つ増やして添字で切り替える方法を用いる．
  • 検証はしていないが，性能を低下させる要因が見当たらない．
  • do concurrentの制約に引っかかるような処理も書ける．

1. そういう処理を書くとコンパイルエラーがでるのではなく，依存性が無いことを前提として，コンパイラがコードを生成してよいという意味です．これとは別に，構文内で呼ぶことができるのは，pureな手続に限られます． ↩