1. はじめに

前回の記事では姫野ベンチマークの実装についてあまり説明ができていなかったので、姫野ベンチマークをOpenACCで並列化してみて、一連の流れからメリット・デメリットを感じていただければと思います。

OpenACCの使い方や、コンパイラの使い方などはGDEPソリューションズにて連載されている東京大学の星野先生による入門記事がございますので、こちらも並行して読んでいただくことが望ましいです。旧PGIコンパイラを用いていますが、多くはNVIDIA HPCコンパイラも同様の動作です。GPUの特徴や、OpenACCのメリットとデメリット、具体例によるOpenACCの実装と最適化などを解説いただいております。OpenACCに関する詳細なドキュメントもHPCWORLDで管理しておりますので、こちらも参照いただけると尚良いと考えます。

今回は執筆時の最新版であるNVIDIA HPC SDK 20.11を使用します。前回は20.9でしたが、性能にはほとんど影響しません。

2. 対象のコード

まずOpenACC化を試みるコードですが、姫野ベンチマークのC, static allocate versionの問題サイズL版を使用します。どの問題サイズでもコードは同一です。283行目付近にある関数float jacobi(int nn)が、ベンチマークの測定対象となるJacobi法の計算です。a, pなどはfloatの4次元ないし3次元配列で、Jacobi法の計算と次ステップのために配列wrk2のデータをpにコピー、という2つの3重ループで構成されています。nnは反復回数で、呼び出し元で自動計算されます。

float
jacobi(int nn)
{
  int i,j,k,n;
  float gosa, s0, ss;

  for(n=0 ; n<nn ; ++n){
    gosa = 0.0;

    for(i=1 ; i<imax-1 ; i++)
      for(j=1 ; j<jmax-1 ; j++)
        for(k=1 ; k<kmax-1 ; k++){
          s0 = a[0][i][j][k] * p[i+1][j  ][k  ]
             + a[1][i][j][k] * p[i  ][j+1][k  ]
             + a[2][i][j][k] * p[i  ][j  ][k+1]
             + b[0][i][j][k] * ( p[i+1][j+1][k  ] - p[i+1][j-1][k  ]
                              - p[i-1][j+1][k  ] + p[i-1][j-1][k  ] )
             + b[1][i][j][k] * ( p[i  ][j+1][k+1] - p[i  ][j-1][k+1]
                               - p[i  ][j+1][k-1] + p[i  ][j-1][k-1] )
             + b[2][i][j][k] * ( p[i+1][j  ][k+1] - p[i-1][j  ][k+1]
                               - p[i+1][j  ][k-1] + p[i-1][j  ][k-1] )
             + c[0][i][j][k] * p[i-1][j  ][k  ]
             + c[1][i][j][k] * p[i  ][j-1][k  ]
             + c[2][i][j][k] * p[i  ][j  ][k-1]
             + wrk1[i][j][k];

          ss = ( s0 * a[3][i][j][k] - p[i][j][k] ) * bnd[i][j][k];

          gosa+= ss*ss;
          /* gosa= (gosa > ss*ss) ? a : b; */

          wrk2[i][j][k] = p[i][j][k] + omega * ss;
        }

    for(i=1 ; i<imax-1 ; ++i)
      for(j=1 ; j<jmax-1 ; ++j)
        for(k=1 ; k<kmax-1 ; ++k)
          p[i][j][k] = wrk2[i][j][k];

  } /* end n loop */

  return(gosa);
}

まず、nvcコンパイラでコンパイル可能か確認してください。このとき-Dで問題サイズを指定します。下記のようにプログラムを実行して、MFLOPS measuredなどの実行ログが得られればOKです。このプログラムはシングルコア・シングルスレッドで動作します。

$ which nvc
/opt/nvidia/hpc_sdk/Linux_x86_64/20.11/compilers/bin/nvc
$ nvc -DLARGE -O3 -o xbmt himenoBMTxps.c
$ ./xbmt
...
 Loop executed for 367 times
 Gosa : 6.965881e-04
 MFLOPS measured : 6923.526678  cpu : 59.300072
 Score based on Pentium III 600MHz : 84.433252

3. OpenACCを用いてGPUで動作させる

まずはシンプルにOpenACCのディレクティブを追記し、GPUで動作するようにします。OpenACC化したコードはBitbucketに公開しています。

計算をGPU化するときは、#pragma acc kernelsまたは#pragma acc parallelを使用します。両者の違いですが、kernelsはOpenACCコンパイラに並列化やループ構造などの実装方法を任せるのに対し、parallelはそれらの実装手段をユーザーが明示的に実装する必要があります。parallelの動作はOpenMPの仕様に近く、明示的な並列化を行う場合に使います。多くの場合kernelsを使い、コンパイラに実装を任せるほうが便利と考えます。変数gosaはスレッド間でリダクションが必要なので、reduction句を使います。reduction句はOpenMPと同様の使い方です。ただし、kernelsにはreduction句を指定できないので、直下に#pragma acc loopを追加します。loopは次に続くループブロックの並列化方法を指定するための機能で、たとえばreductionやループの並列化に割り当てるスレッド数などが指定できます。#pragma omp parallel forと同様に、2つのディレクティブを1行で記述することが可能です。

...
#pragma acc kernels loop reduction(+:gosa)
    for(i=1 ; i<imax-1 ; i++)
      for(j=1 ; j<jmax-1 ; j++)
        for(k=1 ; k<kmax-1 ; k++){
          s0 = a[0][i][j][k] * p[i+1][j  ][k  ]
          ...

#pragma acc kernels
    for(i=1 ; i<imax-1 ; ++i)
      for(j=1 ; j<jmax-1 ; ++j)
        for(k=1 ; k<kmax-1 ; ++k)
          p[i][j][k] = wrk2[i][j][k];
...

NVIDIA HPCコンパイラの場合、-accオプションでOpenACCを有効化します。OpenMPにおける-fopenmpや-qopenmpオプションに該当します。NVCOMPILER_ACC_NOTIFY環境変数に1を設定すると、GPUカーネルの起動ログが出力されます。同環境変数はデータコピーやメモリ確保など色々なログが取れますので、ドキュメントを参照してください。

$ nvc -DLARGE -O3 -acc -o xbmt himenoBMTxps.c
$ NVCOMPILER_ACC_NOTIFY=1 ./xbmt
...
 Now, start the actual measurement process.
 The loop will be excuted in 27 times
 This will take about one minute.
 Wait for a while

launch CUDA kernel  file=/home/hirokawa/himenobmt-oacc/himenoBMTxps.c function=jacobi line=203 device=0 threadid=1 num_gangs=65535 num_workers=1 vector_length=128 grid=65535 block=128 shared memory=1024
launch CUDA kernel  file=/home/hirokawa/himenobmt-oacc/himenoBMTxps.c function=jacobi line=203 device=0 threadid=1 num_gangs=1 num_workers=1 vector_length=256 grid=1 block=256 shared memory=1024
launch CUDA kernel  file=/home/hirokawa/himenobmt-oacc/himenoBMTxps.c function=jacobi line=229 device=0 threadid=1 num_gangs=65535 num_workers=1 vector_length=128 grid=65535 block=128
...

4. データ転送を最小化し十分な性能を得る

先程のプログラムを実行していただければわかりますが、CPUでの実行に比べ性能が悪いです。なぜかというと、上記のコードではCPU-GPU間のデータ移動が最適ではないためです。OpenACCで単純に並列化してしまうと、kernelsディレクティブが現れたときに、GPUのメモリ確保や配列のデータコピーなどが行われ、コードブロックの末尾でこれらのリソースの解放処理が行われてしまいます。これらは処理コストが大きいため、なるべくメモリの確保と解放は繰り返したくないし、データコピーも最小限にしたいところです。GPUカーネルの実行コストがメモリ確保と転送の処理コストよりも非常に大きければ気にする必要はありませんが、多くの場合は逆の状態になります。

Jacobi法のループでは計算中CPUでデータにアクセスしないため、計算している間はGPUにデータの所有権があります。したがって、関数jacobiの呼び出し時と戻り時の2回だけCPU-GPU間のデータ転送があればいいことがわかります。

GPUカーネルを立ち上げずデータ転送のみを行う場合は、#pragma acc dataを用います。#pragma acc dataは指定した変数を次のコードブロックの開始時にCPUからGPUにコピーし、コードブロックの終了時にGPUからCPUに書き戻します。前者はcopyin句で、後者はcopyout句で指定し、また両方を指定するcopy句もあります。上記の関数を見ると配列a, b, c, bnd, wrk1は代入文の左辺に現れないことから、関数jacobiの中では読み取り専用の配列であることがわかり、copyin句で良いと判断できます。配列pは計算結果を書き出す変数なので、copy句の指定が必要です。配列wrk2は計算の途中結果を保持するだけの一時領域のようなので、CPU-GPU間のデータ転送は不要です。ただし配列自体はGPUに確保されていないといけないので、create句を使ってメモリの確保は必要になります。またスカラ変数（omega）は自動的にスレッドプライベート変数としてキャプチャ（firstprivate指定）されるので、明示的に指定する必要はありません。

#pragma acc data copyin(a,b,c,wrk1,bnd) copy(p) create(wrk2)
  for(n=0 ; n<nn ; ++n){
    gosa = 0.0;

#pragma acc kernels loop reduction(+:gosa)
    for(i=1 ; i<imax-1 ; i++)
      for(j=1 ; j<jmax-1 ; j++)
        for(k=1 ; k<kmax-1 ; k++){
          s0 = a[0][i][j][k] * p[i+1][j  ][k  ]
          ...
    ...
  }

NVIDIA Quadro GV100で実行した結果が以下です。前回の記事は他に並列化方法を変更しさらに性能を向上していますが、V100と今回使用したQuadro GV100はほぼ同じ性能を持っているため、前回の結果に対し9割に近い性能が得られていると判断できます。データ転送の最適化がいかに重要かおわかりいただけたでしょうか。

$ ./xbmt
...
 Loop executed for 464 times
 Gosa : 6.766306e-04
 MFLOPS measured : 223677.662736        cpu : 2.320662
 Score based on Pentium III 600MHz : 2727.776375

5. CUDA実装と比較した場合

BitbucketにCUDAの実装例も公開しています。OpenACCよりも細かい最適化を行いましたが、コード行数は元のコードに対し約1.25倍とOpenACCに比べて非常に多くのコード修正や追加が必要になっています。CUDA実装の性能が以下の通りで、OpenACC実装はこれに対し約83%の性能を得られています。

$ make
nvcc -DLARGE -O3 -Xcompiler '-fopenmp' -gencode arch=compute_60,code=sm_60 -gencode arch=compute_70,code=sm_70 -c himenoBMTxps.cu
nvcc -DLARGE -O3 -Xcompiler '-fopenmp' -gencode arch=compute_60,code=sm_60 -gencode arch=compute_70,code=sm_70 -o xbmt himenoBMTxps.o
$ ./xbmt
...
 Loop executed for 843 times
 Gosa : 6.151068e-04
 MFLOPS measured : 268823.624361        cpu : 3.508138
 Score based on Pentium III 600MHz : 3278.336882

もちろんCUDA実装をより細かく最適化していけば差は広がる可能性がありますが、OpenACCは非常に少ない労力でGPUの計算能力を扱えるわけですから、かなり魅力的な開発環境ではないでしょうか。

6. おわりに

今回はまずGPUでプログラムを動かすこと、そしてデータ転送を最小化し実用に耐えられる性能を得ました。OpenACCはCUDAに比べ容易にGPUにアクセスでき、非常に強力なツールです。しかし、データ転送を適切に記述しないと性能を得ることは難しいと感じると思います。これを単純化する仕組みとして、Pascal architecture世代から追加されたUnified memoryによるヒープメモリの自動管理が追加されています。

次回は、Unified memoryを使ったC, dynamic allocate versionの実装についてご紹介します。