OpenACCとcuBLASを用いてGPGPU

はじめに

GPUをグラフィックスの処理だけでなくより一般的な用途にも活用する，いわゆるGeneral-Purpose GPU (GPGPU) は10年以上前からありました．GPGPUを行うためにはいくつかの方法があり，（うまくやれば）最も性能向上が期待できる方法がCUDAなどの専用言語を用いて開発する方法です．この方法は特に既存のプログラムのGPGPU化に適用するにはハードルが高いと考えられます．もう一つの方法としては，ディレクティブベースのOpenACCとベンダーが提供する最適化されたライブラリー用いる方法があり，こちらの方がより簡便にGPGPU化が行えると期待できます．この記事では，CUDAを使わずにGPGPU化を行う手法を，著者が実アプリケーションで遭遇したことのあるパターンの処理を行うサンプルコードを用いて説明したいと思います．

基本的な方針

数値計算では汎用的な行列演算や高速フーリエ変換が計算時間の多くを占めることがあり，このような演算は最適化されたライブラリーを用いることによって高速に処理することができます．

GPGPU環境にもcuBLASというBLASやcuFFTというFFTのライブラリーが存在するため，これらのライブラリーへのライブラリーコールを差し替えるだけで高速化が実現できるかもしれません．しかし，よく言われているようにCPUとGPUの間のデータ転送は低速なため，ライブラリーコールを差し替えるだけではここがボトルネックになり高速化が実現できないかもしれません．そこで，CPUとGPUの通信をなるべく少なくすることが望ましいです．CPUとGPUの通信を少なくするということは，すなわち一度GPUにデータを渡したらCPUに処理が帰ってくるまでになるべく多くの処理を行うことです．この部分をOpenACCによって記述することを基本的な考えとします．OpenACCによる並列化で対象箇所の高速化が実現できれば御の字ですが，仮にこの部分はさほど高速化されなかったとしても性能を阻害するCPUとGPU間のデータのやり取りを抑制できるだけでも有用なのではないかと考えています．

OpenACCについて

OpenACCってなに？

並列計算を行うための規格の一つで，別にGPGPUに特化した仕様ではありませんがGPGPUでよく使われるものです．ディレクティブベースなのでCUDAなどの専用言語にくらべ簡単に導入することができます．特に既存のプログラムをGPGPU化する場合コードの可搬性はできるだけ維持したいので，ディレクティブベースのOpenACCを採用するというのは（性能が十分にでるのであれば）有力な選択肢だと思います．

使える環境

Wikipediaによると「PGIとClayのコンパイラーで使える」とされていますが，PGIのコンパイラーは現在NVIDIAが所有し，NVIDIA HPC SDKの一部となっています．OpenACC自体は名前が示すようにオープンな規格のはずですが，NVIDIA以外のベンダーはあまり力を入れていないようです．

OpenACCのディレクティブ

今回の試作プログラムで用いるOpenACCのディレクティブについて説明します．コード例はいずれもFotranの場合です．

基本の書き方

Fortranの場合，!$accのあとにディレクティブを記述します．各ディレクティブはクローズ(clause)を持つ場合があり，それらはクローズ名(値)のような記法で指定します．複数のクローズを指定したい場合スペース区切りで記述します．長くなりそうな場合，Fortranの継続行のルールと同じ方法で継続させることができます．

data

CPUとGPUの間でデータをやり取りするためのディレクティブです．たとえば以下のように記述します．

!$acc data copyin(a) copyout(b) copy(c) create(work) present(d)
...
!$acc end data

!$acc dataにおいてGPU上でメモリーが確保され，!$acc end dataまで保持されます．dataは以下のクローズを持ちます

copyin

CPUからGPUへデータを送ります．入力配列をCPUで作り，GPUへ送る場合などに使います．

copyout

end dataのタイミングでGPUからCPUへデータを送られる配列を指定します．GPUの計算結果をCPUで用いたい場合に使います．

copy

開始時にCPUからGPUへデータを送り，終了時GPUからCPUへデータを送ります．入力と出力両方を担う配列に使います．

create

GPU上にメモリーを確保します．転送は行われません．一時配列などに用います．

present>

GPU上にすでに確保されている配列であることを表します．必須ではありませんが，たとえば実際には確保していない配列を利用しようとする際に分かりやすいエラーで終了するので，なるべく指定した方が間違いを減らせそうです．

C言語の場合は渡す配列の大きさも同時に指定しないと高い確率で異常終了してしまいます．Fotranの場合は配列自身が大きさの情報を持っているので，大きさの指定は不要です．部分配列をやり取りしたい場合，Fortranの通常の方法で指定することもできます(例：A(2:n-3) など)

kernels

自動並列のためのディレクティブです．たとえば以下のように記述します．

!$acc kernels
...
!acc loop independent
do i=1, n
...
end do
!$acc end kernels

kernelsディレクティブはコンパイラーによる自動並列を行わせるためのものであり，対象コードが必ず並列実行されるわけではありません．データの独立性があるにも関わらずコンパイラーにはわからない場合もあり，そのようなループに対して!acc loop independent を挿入すると並列実行してくれるようになります

parallel

並列実行を指定するためのディレクティブです．たとえば以下のように記述します．

!$acc parallel
...
!acc loop gang
do i=1, n
!acc loop vector
  do j=1, m
    ...
  enddo
enddo
...
!$acc end parallel

kernelsとよく似ていますが，parallelディレクティブの場合明示的に並列実行対象のループを指定する必要があります．kernelsディレクティブはコンパイラーに対してヒントを与えるものであるのに対し，parallelディレクティブは命令を与えるものです．したがって，並列実行できないループに対して指定すると正しくない計算が行われることになってしまいます．

host_data

cuBLASなどのCUDAで記述されたサブルーチンへのサブルーチンコールを行う際に，対象の配列がGPU上のメモリーに確保されていることを伝えるためのディレクティブです．このディレクティブを使うことによってCUDAで記述されたライブラリーを非CUDAコードから呼ぶことができます．たとえば以下のように記述します．

!$acc host_data use_device(a, b, c)
call cublasDgemm('N', 'N', n, m, k, 1.d0, a, lda, b, ldb, 1.d0, c, ldc)
!$acc end host_data

CUDAライブラリーの別の呼び方としては(Fortranの場合) device属性のついた配列にデータをコピーし，cuBLASのルーチンに渡すという方法もあります．しかし，OpenACCにおいては配列を「ホスト用」と「デバイス用」とで分けて用いるのは推奨されないやり方になっています．

cuBLASについて

cuBLASってなに？

cuBLASはCUDAで記述されたBLASの実装で，GPU上で高速に動作することが期待できます．BLASは代表的な行列演算を一通りカバーするライブラリーです．ベクトル―ベクトル演算を担うLevel 1, ベクトルー行列演算を担うLevel 2, 行列―行列演算を担うLevel 3に分けられますが，特にLevel 3 BLASはプラットフォームに最適化されたそれを用いることによってコードを加速させることができます．

使える環境

cuBLASはNVIDIAのGPUにおいてしか用いることができませんが，BLAS自体はあらゆる計算機環境において最適化されたものが用意されています．そのため，BLASのルーチンに置き換えることによってあらゆるプラットフォームで高速化が期待できます．

今回用いるcuBLASの関数

今回のサンプルプログラムでは行列―行列積の計算を行うので，Level 3 BLASのルーチンであるdgemmを用います．cuBLASでのルーチン名はcublasDgemmのようです．cuBLASのバージョン2というものやcuBLASXtというものもあるそうですが，今回は通常のBLASに最も近いレガシーcuBLASを用いることにします．

試作コードとそのGPGPU化

試作コードの処理内容

試作コードは何か意味のある計算を行うものではありませんが，実アプリケーションにおいて実際に現れる(著者が見たことがある)パターンを簡略にしたものです．想定としては行列―行列積を行うものなのでLevel 3 BLASの関数dgemmを用いるのですが，対象とする行列のサイズのうち一つは非常に大きくなり得るため場合によってはメモリーが枯渇してしまいます．そこで，小行列を小分けに作っていきながら行列行列積を実行し，結果の配列に格納していきます．

入力行列$A$および$B$の大きさはそれぞれ$(n, k)$および$(m, k)$，結果を格納する配列$C$の大きさは$(n, m)$とします．$nmk$は典型的には$n > k > m$という関係があるものとします．また，行列$A$はその全要素をメモリ上に確保することはなく，小分けにつくっていけるものであるとします．行列$A$のを構成する要素は問題規模に依存しない行列 $A_p$と依存するベクトル$a_p$であるとし，$A$の行列要素は$A_p$と$a_p$の行列要素を用いて$A(i, k+(j-1) \times N_{A_p} ) = A_p(i, j) a_p(k)$ のような感じで計算できるとします(ここで$N_{A_p}$は$A_p$のカラム数)．つまりメモリ上に保持する配列は$A_p$と$a_p$のみであり，行列$A$の一部を構築しながらdgemmを適用していくことを行います．

言葉では分かりづらいので，具体的なコードを以下に示します．

cpu.f90

program dgemm_test
  implicit none
  integer :: k_partial_var, k_partial_fixed, k, n, m, nblock, niter_inner, niter_outer
  real(kind=8), allocatable, dimension(:,:) :: mat_c, mat_b, mat_a_partial, work
  real(kind=8), allocatable, dimension(:)   :: mat_a_partial_vec
  integer :: i, j
  integer :: stime, etime, rate

  niter_outer     = 10
  niter_inner     = 5
  k_partial_fixed = 5
  k_partial_var   = 500
  m               = 2000
  n               = 50000
  nblock          = 1000
  k               = k_partial_fixed*k_partial_var

  allocate(mat_c(n, m))
  allocate(mat_b(m, k))
  allocate(mat_a_partial(n, k_partial_fixed))
  allocate(mat_a_partial_vec(k_partial_var))
  allocate(work(nblock, k))
  call set_vector(0.2d0, 1, k_partial_var, mat_a_partial_vec)
  do i=1, niter_outer
    call set_matrix(0.015d0, i, 1, n, k_partial_fixed, mat_a_partial)
    do j=1, niter_inner
      call set_matrix(0.045d0, i, j, m, k, mat_b)
      call system_clock(stime)
      call dowork(n, m, k_partial_var, k_partial_fixed, nblock, mat_b, mat_a_partial, mat_a_partial_vec, work, mat_c)
      call system_clock(etime, rate)
      write(0,'(2i5, f25.5, f8.3)') i, j, sum(mat_c), dble(etime-stime)/dble(rate)
    enddo
  enddo
  deallocate(mat_c)
  deallocate(mat_b)
  deallocate(mat_a_partial)
  deallocate(mat_a_partial_vec)
  deallocate(work)

  contains

  subroutine set_vector(factor, i, n1, mat1d)
    real(kind=8), intent(in) :: factor
    integer, intent(in) :: i, n1
    real(kind=8), dimension(n1), intent(out) :: mat1d
    integer :: jj
    do jj=1, n1
      mat1d(jj) = (jj-n1/2)*factor+i
    enddo
  end subroutine set_vector

  subroutine set_matrix(factor, i, j, n1, n2, mat2d)
    real(kind=8), intent(in) :: factor
    integer, intent(in) :: i, j, n1, n2
    real(kind=8), dimension(n1,n2), intent(out) :: mat2d
    integer :: ii, jj
    do ii=1, n2
      do jj=1, n1
        mat2d(jj, ii) = (ii-n2/2)*factor+(jj-n1/2)+i-j
      enddo
    enddo
  end subroutine set_matrix

  subroutine dowork(n, m, k_partial_var, k_partial_fixed, nblock, mat_b, mat_a_partial, mat_a_partial_vec, work, mat_c)
    integer, intent(in) :: n, m, k_partial_var, k_partial_fixed, nblock
    real(kind=8), intent(in),  dimension(:, :) :: mat_b, mat_a_partial
    real(kind=8), intent(in),  dimension(:) :: mat_a_partial_vec
    real(kind=8), intent(out), dimension(:, :) :: mat_c
    real(kind=8), allocatable, dimension(:,:) :: work
    integer :: i, j, l, k, iend, isize, ksize
    mat_c = 0.d0
    ksize = k_partial_var * k_partial_fixed
    do i=1, n, nblock
      iend = min(n, i+nblock-1)
      isize = iend-i+1
      do j=1, k_partial_var
        do k=1, k_partial_fixed
          do l = 1, isize
            work(l, k+(j-1)*k_partial_fixed) = mat_a_partial(i+l-1, k)*mat_a_partial_vec(j)
          enddo
        enddo
      enddo
      call dgemm('N', 'N', isize, m, ksize, 1.d0, work, nblock, mat_b, ksize, 1.d0, mat_c(i,1), n)
    enddo

  end subroutine dowork

end program dgemm_test

• 配列mat_Aが$A$, mat_Bが$B$, mat_cが$C$に相当します
• 大きさ$nmk$はそのままスカラー変数n, m, kに対応します
• 配列mat_a_partialが$A_p$に対応します。また，ベクトルmat_a_partial_vecが$a_p$に対応します．
• k_partial_fixedが$N_p$に対応します．k_partial_varは$a_p$のサイズです．
• サブルーチンset_vectorおよびset_matrix は配列に適当なデータを設定するものです．実アプリケーションなら何かしら意味のある計算を行い，配列に格納する役割を果たすものでしょう．
• 処理は二重ループで行われるようになっています．mat_a_partial_vecはプログラム実行中一度設定されると変化されることはありません．mat_a_partialは外側のループで値が設定され，内側のループでは変化しません．
• サブルーチン dowork が行列を小分けに作成し，dgemm をコールするサブルーチンです．小分けの単位がnblockで，do文をnblockの単位でスキップしながら回しています．その下の三重ループにおいて配列work に行列$A$の部分行列をコピーし，dgemm に渡しています．

試作コードのGPGPU化

その1. サブルーチンで完結するケース

最も簡単なGPGPU化として，サブルーチンdoworkで完結するケースを紹介します．これで十分な高速化が見込めるのであればサブルーチン単位で少しずつ対処していけばよいので，実アプリケーションのGPGPU化を少しずつ進めることができ魅力的です．

サブルーチンdowork の変数宣言のあとに，以下のようなディレクティブを挿入しました．

gpu1.f90

    !$acc data copyin(mat_b, mat_a_partial, mat_a_partial_vec)  &
    !$acc      copyout(mat_c) &
    !$acc      create(work)
    !$acc host_data use_device(work, mat_b, mat_c)

dataディレクティブでGPUにデータの扱い方を指示しています．copyinによってGPUにデータを転送しています．copyout でdataディレクティブのスコープが終了した際GPUから転送するデータを指定しています．createによって一時配列をGPU上で確保しています．その次は結果を格納する配列mat_cの初期化をkernelsディレクティブを用いてGPUで行うようにしています．

gpu1.f90

    !$acc kernels
    mat_c = 0.d0
    !$acc end kernels

さらに，dgemmに渡す小行列を作る部分をparallelディレクティブでGPUにおいて並列処理を行うようにしています．

gpu1.f90

      !$acc parallel
      !$acc loop gang collapse(2)
      do j=1, k_partial_var
        do k=1, k_partial_fixed
          !$acc loop vector
          do l = 1, isize
            work(l, k+(j-1)*k_partial_fixed) = mat_a_partial(i+l-1, k)*mat_a_partial_vec(j)
          enddo
        enddo
      enddo
      !$acc end parallel

最上位のjループはgang単位の並列化にするため!$acc loop gangを入れています．ここではその次のループも展開するようcollapse(2)も追加しています．最内のlループはvector単位の並列化にするため !$acc loop vectorを入れています．dgemmはcublasDgemmに変更しています．

gpu1.f90

call cublasDgemm('N', 'N', isize, m, ksize, 1.d0, work, nblock, mat_b, ksize, 1.d0, mat_c(i,1), n)

cublasDgemmに渡す配列はdevice属性が付与されている必要がありますが，このことはhost_data ディレクティブによってうまい具合に処理されるようになっています．

最後に，dataとhost_dataを終了します．

    !$acc end host_data
    !$acc end data

その2. 二重ループ中のデータ転送を結果配列のみにする

試作プログラムではサブルーチンの呼び出しを二重ループ内で行います．内側のループで結果の配列をCPUが用いますが，それ以外のデータはあらかじめGPUに置いておくことができます．特に一時配列は実際にはCPUとGPUの転送は不要です．そこで，結果配列mat_c以外はあらかじめGPUに確保するようにしてみます．

まず，二重ループの外側で一時配列をGPUに確保します．

gpu2.f90

  !$acc data create(work, mat_a_partial_vec, mat_a_partial, mat_b, mat_a_partial_vec)
  call set_vector(0.2d0, 1, k_partial_var, mat_a_partial_vec)
  do i=1, niter_outer
    call set_matrix(0.015d0, i, 1, n, k_partial_fixed, mat_a_partial)
    ...
    ...
  enddo

結果配列mat_cのcopyout指定はサブルーチン呼び出しの前後に移動します．

gpu2.f90

      !$acc data copyout(mat_c)
      call dowork(n, m, k_partial_var, k_partial_fixed, nblock, mat_b, mat_a_partial, mat_a_partial_vec, work, mat_c)
      !$acc end data

サブルーチンset_vector, set_matrixもGPGPU化しないと通信が発生してしまいます．以下のように，dataディレクティブで用いる配列がGPUにあるものと宣言しておき，サブルーチン本体はkernelsディレクティブで囲います．

  subroutine set_vector(factor, i, n1, mat1d)
    real(kind=8), intent(in) :: factor
    integer, intent(in) :: i, n1
    real(kind=8), dimension(n1), intent(out) :: mat1d
    integer :: jj
    !$acc data present(mat1d)
    !$acc kernels
    do jj=1, n1
      mat1d(jj) = (jj-n1/2)*factor+i
    enddo
    !$acc end kernels
    !$acc end data
  end subroutine set_vector

  subroutine set_matrix(factor, i, j, n1, n2, mat2d)
    real(kind=8), intent(in) :: factor
    integer, intent(in) :: i, j, n1, n2
    real(kind=8), dimension(n1,n2), intent(out) :: mat2d
    integer :: ii, jj
    !$acc data present(mat2d)
    !$acc kernels
    do ii=1, n2
      do jj=1, n1
        mat2d(jj, ii) = (ii-n2/2)*factor+(jj-n1/2)+i-j
      enddo
    enddo
    !$acc end kernels
    !$acc end data
  end subroutine set_matrix

サブルーチンdoworkにおけるdataディレクティブは以下のように変更します．

gpu2.f90

    !$acc data present(mat_a_partial, mat_a_partial_vec, work, mat_b, mat_c)

presentクローズを用いることによって，GPUにおいて確保されている配列を宣言しています．

以上の変更によって，二重ループにおけるデータ転送は結果を格納する配列mat_cのみにすることができました．

その3. CPUとGPUの通信をなくす

結果配列mat_cは二重ループ中でCPUが用いるのでサブルーチンコールの前後で転送が必要でした．仮に二重ループ後に用いるのであればこの転送も省くことができます．まず配列mat_cをdata createに含めるようにします．

!$acc data create(work, mat_a_partial_vec, mat_a_partial, mat_b, mat_a_partial_vec, mat_c)

ついで，二重ループ中のmat_cのcopyoutを削除します．また，試作プログラムではサブルーチンdworkをコールしたあと配列mat_cの要素の和をエラー出力に出力するようにしていましたが，これはCPUの処理なので削除します．以上の変更によってCPUとGPUの通信をなくしてしまうことができます．

コンパイル・リンク

コンパイルとリンクはNVIDIA HPC SDKのバージョン22.5を用いました．以下の要領でコンパイルとリンクを行います．

nvfortran -Minfo=accel -acc -cuda -cudalib=cublas -o gpu1 gpu1.f90

-cudalib=cublasをつければcuBLASとリンクすることができます．-Minfo=accelを指定するとOpenACCによってどのように並列化がなされたかの詳細が以下の要領で報告されます．

dowork:
     72, Generating copyin(mat_a_partial_vec(:),mat_a_partial(:,:),mat_b(:,:)) [if not already present]
         Generating create(work(:,:)) [if not already present]
         Generating copyout(mat_c(:,:)) [if not already present]
     77, Loop is parallelizable
         Generating NVIDIA GPU code
         77,   ! blockidx%x threadidx%x auto-collapsed
             !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
     83, Generating NVIDIA GPU code
         85, !$acc loop gang collapse(2) ! blockidx%x
         86,   ! blockidx%x collapsed
         88, !$acc loop vector(128) ! threadidx%x
     88, Loop is parallelizable

このような情報を参考に，仮にコンパイラーが並列化できなかったループがあった場合に対応を検討したりすることができます．

検証

CPU版と3通りのGPGPU版を用いてプログラムを実行してみます．CPU版は非並列およびOpenMP 16並列，GPUは1枚利用します．用いたCPUはAMD EPYC 7742 64コア，GPUはNVIDIA A100です．

検証 1.

次のような配列サイズで試作プログラムを実行してみました．

  k_partial_fixed = 5
  k_partial_var   = 500
  m               = 2000
  n               = 50000
  nblock          = 1000

何度か述べたように今回の試作プログラムは著者が実際に遭遇した実アプリケーションのパターンの処理をもとに作ったものですが，そのアプリケーションにおいて実際に現れ得る程度の大きさを選んでいます．外側ループ一回あたりにかかった計算時間は下記の通りです．単位は秒です．

CPU 非並列	CPU16並列	GPGPU その1.	GPGPU その2.	GPGPU その3.
10.2	1.53	0.131	0.128	0.03

GPUによってCPU 16並列の場合と比較して結果配列を転送する場合でも一桁程度の高速化を実現することができました．しかしCPU-GPU間転送がない GPGPU その3.と比較すると転送のあるケースは相当程遅く，計算時間の大部分が通信に用いられていることが分かります．

検証 2.

上の検証ではある程度の高速化が実現できましたが，一方でほとんどの時間がCPU-GPU間の転送で使われているということも分かりました．この影響が小さくなるような配列サイズを採用し，どの程度結果が変化するか検証してみます．具体的には先のパラメーターのk_partial_var を500から5000にすることによって転送量はそれほど増やさず，演算量のみを増やしてみます．

CPU 非並列	CPU16並列	GPGPU その1.	GPGPU その2.	GPGPU その3.
102	13.66	0.404	0.379	0.289

CPU非並列の場合は計算時間がちょうど10倍になりました．CPU 16並列の場合は，演算量が増えた分並列化効率が向上したため10倍にはなっておらず9倍程度計算時間が増えました．GPUその1.とその2. は通信量はほぼ変わらず演算量のみが増えたので3倍程度の計算時間の増加にとどまっています．GPGPUその3.はもともと通信はなかったのですが，演算量が増えた分効率が上がったためか計算時間の増加は10倍に少し満たないくらいになりました．

考察

GPUを用いることによって高速化が実現できることが確認できました．ただ，実は今回の検証はCPUに対してフェアなものではない部分があります．HPC用のCPUとGPUとではお値段が違いますし，CPUは64コア積んでいるにも関わらず16並列までしか実施していません．これは，これ以上並列数を増やすとOpenMPでは逆に遅くなってしまったからなのですが，OpenMPだけでなくMPIも併用すれば64コアをより効率よく使えたかもしれません．他方GPU側に不利な点としては，プログラムとしてあまりにシンプルで演算を行っている箇所が少ないため，CPU-GPU間の転送が隠蔽されにくいという点が挙げられるかもしれません．実アプリケーションでは転送後もっといろいろな処理が行われるはずなので，この例ほどはCPU-GPU間の転送が目立つことはないことも多いと思います．

終わりに

OpenACCとcuBLASを用いてサンプルプログラムをGPGPU化してみました．ディレクティブベースの改変とライブラリーの差し変えだけで（配列の大きさの組み合わせによっては）相当程度の高速化が実現できました．チューニングの余地はいろいろとあるかもしれませんが，この程度の手間で高速化が実現できることが確認できたのはよかったと思います．