はじめに

以前、C++ parallel algorithmを使ったNVIDIA C++ compilerでのGPU並列化を紹介しました。
同様にNVIDIA Fortran compilerでは、Fortran 2008で定義されたdo concurrent文により、言語仕様レベルでのGPU並列化コードを実装可能です。
C++と異なり標準ライブラリではなく構文としてのサポートなので、より簡単な記述になっています。

NVIDIA technical blogでは説明されていますが、
NVIDIA HPC SDKのドキュメントには詳細な説明がないため、今回はこの機能を検証してご紹介します。
NVIDIA HPC SDKのコードの検証には、NVIDIA HPC SDK 23.7を使用しています。

シンプルな並列化

1次元のループの場合には、do文をdo concurrent文に置き換えるだけです。
do i=1,nとdo concurrent(i=1:n)が同等で、あとはコンパイラオプションによって並列化方法を指定します。

subroutine saxpy(n, a, x, y)
    implicit none
    integer, intent(in) :: n
    real, intent(in) :: a, x(:)
    real, intent(inout) :: y(:)
    integer :: i

    do concurrent(i=1:n)
        y(i) = a * x(i) + y(i)
    end do
end subroutine saxpy

OpenACCと同じように、コンパイル時に-Minfoオプションで並列化情報を取得できます。
do concurrentによる並列化を有効化するには、C++の場合と同じように-stdparを使います。CPUで並列化する場合には-stdpar=multicoreを指定します。

$ nvfortran -stdpar -Minfo -c saxpy.f90
saxpy:
      8, Generating NVIDIA GPU code
          8, Loop parallelized across CUDA thread blocks, CUDA threads(128) blockidx%x threadidx%x
      8, Generating implicit copy(y(:n)) [if not already present]
         Generating implicit copyin(x(:n)) [if not already present]

2つ目の並列化メッセージの通りimplicit copyやcopyin (OpenACCの構文と同等の意味と思われます) が生成されているので、do concurrentによる並列化はUnified memoryが必須ではないようです。ただ、-stdparを有効化するとUnified memoryも同時に有効になっていそうです。

そのため、do concurrentで計算する配列は静的確保したものでも可能です。
実行してみたところ、静的配列で実行した場合にCPU-GPU間のmemcpyが走っていることも確認できました。

$ NVCOMPILER_ACC_NOTIFY=2 ./a.out
upload CUDA data  file=/home/prometech/run_saxpy.f90 function=saxpy line=25 device=0 threadid=1 variable=x(:n) bytes=48
upload CUDA data  file=/home/prometech/run_saxpy.f90 function=saxpy line=25 device=0 threadid=1 variable=y(:n) bytes=48
download CUDA data  file=/home/prometech/run_saxpy.f90 function=saxpy line=27 device=0 threadid=1 variable=y(:n) bytes=48
...

（NVCOMPILER_ACC_NOTIFY環境変数はOpenACC用ですが、バックエンドがOpenACCとdo concurrentで同一実装と予測されます）

より複雑な並列化

次に多重ループの並列化方法として行列ベクトル積 (sgemv) を実装します。以下が実装例です。

subroutine sgemv(m, n, A, x, y)
    implicit none
    integer, intent(in) :: m, n
    real, intent(in) :: A(:, :), x(:)
    real, intent(out) :: y(:)
    integer :: i, j
    real :: ret

    do concurrent(i=1:n)
        ret = 0
        do concurrent(j=1,m)
            ret = ret + A(j, i) * x(j)
        end do
        y(i) = ret
    end do
end subroutine sgemv

行列ベクトル積の場合二重ループとなりますが、どちらもdo concurrentを使用できます。
CUDA GPUで考えると、外ループをblock方向に、内ループをthread方向に並列化するイメージです。

内ループで使っている変数retは畳み込みになっているため、OpenACCやOpenMPなどと同様にreductionの指示が必要ですが、OpenACCと同様にreductionの有無はコンパイラが判定しやすいため、ほとんどの場合は自動的に生成されます。
上記のコードをコンパイルすると、以下の通り意図通りの並列化が実現されていることがわかります。

$ nvfortran -stdpar -Minfo -c sgemv.f90
sgemv:
      9, Generating implicit private(ret)
         Generating NVIDIA GPU code
          9, Loop parallelized across CUDA thread blocks ! blockidx%x
         11, Loop parallelized across CUDA threads(128) ! threadidx%x
             Generating implicit reduction(+:ret)
      9, Generating implicit copyin(a(:m,:n)) [if not already present]
         Generating implicit copyout(y(:n)) [if not already present]
         Generating implicit copyin(x(:m)) [if not already present]

暗黙的に定義されている指定を明示した場合は、以下の書き方になります。
do concurrentのreduceは今月ISOの言語規格としてpublishされる予定のFortran 2023で定義されるとのことで、OpenACCやOpenMPのreduction句と同等です。reduction処理はほとんどの場合で必須ですから、NVIDIA Fortran compilerは先行して実装しているようです。
Fortran 2018で定義されるlocal句が、OpenACCやOpenMPのprivate句に相当する変数のlocality指定です。

do concurrent(i=1:n) local(ret)
    ret = 0
    do concurrent(j=1,m) reduce(+:ret)
        ret = ret + A(j, i) * x(j)
    end do
    y(i) = ret
end do

説明を読む限り、local等は以下の対応関係と思われます。

Jacobi法の並列化

最後にNVIDIA technical blogでも紹介されていますが、2D Jacobi法を並列化した場合について説明します。

do concurrent(i=2 : n-1, j=2 : m-1)
    a(i,j) = w0 * b(i,j) + &
            w1 * (b(i-1,j) + b(i,j-1) + b(i+1,j) + b(i,j+1)) + &
            w2 * (b(i-1,j-1) + b(i-1,j+1) + b(i+1,j-1) + b(i+1,j+1))
enddo
do concurrent(i=2 : n-1, j=2 : m-1)
    b(i,j) = w0 * a(i,j) + &
            w1 * (a(i-1,j) + a(i,j-1) + a(i+1,j) + a(i,j+1)) + &
            w2 * (a(i-1,j-1) + a(i-1,j+1) + a(i+1,j-1) + a(i+1,j+1))
enddo

※NVIDIA technical blogからの引用

上記はdo concurrentの例ですが、OpenACCでは以下の書き方になり、計算コードは完全に同一です。OpenMPも同様で、acc kernelsの代わりにomp parallel doをそれぞれの二重ループの頭に記述します。

!$acc kernels
    do i=2,n-1
    do j=2,m-1
      a(i,j) = ...
    enddo
    enddo

    do i=2,n-1
    do j=2,m-1
      b(i,j) = ...
    enddo
    enddo
!$acc end kernels

コンパイルは以下のオプションで行い、問題サイズは8192×8192、OpenMP並列化したCPUでの計算時間を1として算出しました。ただし、CPU性能はよりフェアな比較となるように強力な最適化を実施できるIntelコンパイラを使って評価しています。

# OpenMP
ifort -O3 -march=native -fopenmp jacobi_openmp.f90

# do concurrent
nvfortran -stdpar -fast jacobi_do_concurrent.f90

# OpenACC
nvfortran -acc -gpu=managed -fast jacobi_openacc.f90

以下の表からもわかる通り、do concurrentによりOpenACCと同等の性能が得られているようです。
より大きな問題や3次元問題などメモリバンド幅が更に効くような計算では、CPUに比べ更に高い性能が期待されます。

まとめ

今回はFortran do concurrentによるGPU並列化を紹介しました。
do concurrentはFortran 2008から規格化された仕様のため、移植性が高い方法と言えます。

またC++ parallel algorithm異なりライブラリ機能ではなく言語機能のため、よりシンプルな書き方であることは強いメリットと言えそうです。
NVIDIA Fortran compilerに限った話ですがCUDA GPUでどのように並列化されるかOpenACCと同じ枠組みで見れますので、最適化はしやすいと感じました。

C++ parallel algorithmと同様、より簡易なGPUの利用方法としてご検討ください。