PGIアクセラレータによる姫野ベンチマークの最適化手順(その1)

キーワード GPGPU、CUDA、姫野ベンチマーク、PGIアクセラレータ、Fermi 性能

 前回のコラムでは、行列積の計算をPGIが提供する複数のGPU用のプログラミング方法によって性能がどの程度異なるのかを示しました。今回のコラムは、非圧縮性流体の Poisson 方程式ソルバーの一つである姫野ベンチマークを例題に、PGIアクセラレータ・プログラミングモデルを使用してどの程度まで性能を向上できるかを実際の作業を説明しながらお見せします。この中の説明では、具体的な PGIアクセラレータのディレクティブの説明も行っていますので、実践的にPGIアクセラレータを使おうとする方も見様見真似で自分のプログラムに適用してみてはいかがでしょうか。ここでは、PGI 10.6/10.4 を使用した Fermi の性能を呈示しています。今回はその 1 で、次回のコラムで、PGIアクセラレータのディレクティブを使用して最終的な Kernel ループスケジューリング等のチューニングした結果をお見せします。最終的な性能は、PGI 10.4 コンパイラを使用して、55 GFLOPS の性能を記録しました。
2010年7月7日 Copyright © 株式会社ソフテック 加藤

姫野ベンチマーク・プログラムへの適用

 姫野ベンチマークの性能は、メモリ帯域で律速となるアプリケーション特性を有します。従来のインテルやAMDのマルチコアCPUでは、複数のコア数を備えスレッド並列しようとも、メモリ性能が悪いシステムでは、その(スレッド並列の)実効性能はメモリ帯域の如何で頭打ちとなってしまう特性を有します。従来の CPU に較べ、NVIDIAの GPU が持つ高いメモリ帯域を利用することで、うまくいけば、より高い性能を得ることができることは容易に想像がつきますが、PGIアクセラレータ・プログラミングモデルであるコンパイラの指示行レベルで、どこまで性能が伸びるかは興味のあるところです。

 PGIアクセラレータによる姫野ベンチマークの最適化手順 コラム(全2回)

Performance Summary

単精度 性能加速 13.5x

 最初に、今回 PGIアクセラレータコンパイラを用いて得られた性能を示します。姫野ベンチマークの実行に使用した計算サイズは、M (256 x 128 x 128) です。

表-1 PGI 10.6 における 姫野ベンチマーク(Fortran90) の GPU 性能(その1)
作業(最適化)ステップ GTX 480 GTX 285
STEP 1 反復ループの内側にアクセラレータ計算領域ディレクティブ 1.29 GFLOPS 1.27 GFLOPS
STEP 2 反復ループの外側にアクセラレータ計算領域ディレクティブ 24.20 GFLOPS 18.29 GFLOPS
STEP 3 STEP 2+データ転送Clauseを追加指定 24.41 GFLOPS 23.62 GFLOPS
STEP 4 反復ループの内側にアクセラレータ計算領域を置き、別途、データ領域ディレクティブ使用 24.33 GFLOPS 23.54 GFLOPS
STEP 5 STEP 3+データ領域ディレクティブを別途指定 24.40 GFLOPS 23.62 GFLOPS
STEP 6以降 続編のコラムへ --- ---

性能実測するシステム環境について

 まず始めに、性能実測に使用する私のテストマシン環境について説明します。

ホスト側プロセッサ Intel(R) Core(TM) i7 CPU 920 @ 2.67GHz x 1 プロセッサ
ホスト側メモリ PC3-10600(DDR3-1333)2GB x 2 枚
ホスト側 OS Cent OS 5.4、カーネル 2.6.18-164.el5
PGI コンパイラ PGI 10.6 (PGI 2010)
GPU ボード1 NVIDIA Geforce GTX480(Fermi)、1401.0 MHz clock, 1535.7 MB memory
GPU ボード2 NVIDIA Geforce GTX285、1476.0 MHz clock, 1023.3 MB memory
NVIDIA CUDA バージョン CUDA 3.0ドライバと3.0 Toolkit
表-1.1 NVIDIA GPU 理論性能の比較
ハードウェア仕様 Geforce GTX285 Geforce GTX480
Compute Capability Compute Capability 1.3 2.0
総CUDAコア(SP)数 CUDA cores
(Streaming Processors)
240 cores 480 cores
プロセッサクロック CUDA cores
(Streaming Processors)
1476 MHz 1401 MHz
単精度理論FP演算性能 Single Precision GFLOPS 708 GFLOPS 1344 GFLOPS
倍精度理論FP演算性能 Double Precision GFLOPS 88 GFLOPS 1344x1/8 = 168 GFLOPS

使用する Fortran 90 姫野ベンチマークプログラムと本稿での適用ステップの説明

 使用する姫野ベンチマークプログラムは、Fortran90構文で記述された himenoBMTxp.f90 をベースとして使用します。このプログラムの性能に関しての説明は、前の技術コラムである「自動並列化性能」でも紹介しました。ここでは、このコラムで使用した各最適化ステップ毎のソースファイルも提供することにします。

 本稿のPGIアクセラレータの適用ステップは次の通りです。

  1. STEP1:単純に 並列の対象ループ・ブロックを !$acc region /end region で囲む、悪い例
  2. STEP2:反復ループの外側を !$acc region /end region で囲む
  3. STEP3:STEP2に明示的にホストとGPU間のデータ転送を指定するクローズ(Clause)を使用する
  4. STEP4:STEP 3 と同等な形態であるが、反復ループの内側にアクセラレータ計算領域を置き、別途、データ領域ディレクティブ使用する例
  5. STEP5:STEP 3 と同等な形態であるが、データ領域ディレクティブで別途指定する例
● Makefile の中の FFLAGS (PGIオプション)
# for CUDA 3.0 driver on your system
 FFLAGS = -fastsse -Minfo=accel -ta=nvidia,cuda3.0,time
# for CUDA 2.3 driver on your system
#FFLAGS = -fastsse -Minfo=accel -ta=nvidia,time

● Makefile の方法
(例)step1.F90 をコンパイルし、step1.exe と言う実行モジュールを作成する場合
   ソースファイル名は step*.F90 と言う形式を提供している。* は、数字。
   実行モジュール名は step*.exe と言う名前としている。* は、数字。
   コンパイルする際は、make の引数にビルドしたい実行モジュール名を指定する。
$ make step1.exe 

※ PGIアクセラレータ用コンパイル・オプションについて

※ PGI Windows版の makefile の利用に関して

姫野ベンチマークプログラムの実行条件

 姫野ベンチマークを実行する上で、ソースプログラムの変更は以下の通りです。

  • 実行に使用する計算サイズは、M (256 x 128 x 128) とする。以下の各ステップの実行結果・性能は GTX480(Fermi) を使用した場合の値である。なお。GTX480 の倍精度演算ピーク性能は、単精度の 1/8 に制限されていますのでご注意下さい
  • ヤコビの反復回数を 800 回の固定で実行する。これは、最適化ステップの実行時間の推移を定量的に比較するためにこのような措置をとる。
  • プログラムのとおり、gosa の総和計算は、各反復の度に計算を行う。すなわち、以下で示す実行時間はプログラムの通り、反復毎にリダクション計算を行っている場合の性能表示である

 それでは、PGI アクセラレータの指示行を入れて、順番に最適化作業を行ってみましょう。最初に、単純に「アクセラレータ計算領域」を指示行で指示してコンパイル・実行してみましょう。

ステップ 1 (単純に 並列の対象ループ・ブロックを !$acc region /end region で囲む、悪い例)

 まず、GPU をアクセラレータとして使用するために、GPUでの処理を行うプログラムの「アクセラレータ計算領域」を決めます。姫野ベンチのプログラムは、subroutine jacobi(nn,gosa) のルーチンの処理で、ほぼ100% の時間を使用しているため、ここが GPU 並列処理の対象となります。下記のリストは、そのルーチンの抜粋をしたものです。左端に、行番号が付されています。ソースプログラムの全体を見て、1) ループ構造を見出し2) その中で使用されている配列(変数)が各ループ・インデックスに対して、データ依存があり得るかを判別します。(依存性に関しては、今後のコラムで説明するつもりです。)依存性がないループを有している場合、GPU 並列化の対象ブロックとすることができます。対象とするループは、1次元のループよりも多次元のループ構造の方が、並列化対象とする実行時間量(これを一般に、並列対象の「粒度」と言う。)が大きいため、並列性能の加速性が良くなります(これは、GPUに限らず、一般的な並列化における特徴・概念です)。姫野ベンチマークの場合は、3次元の i, j, k の空間のポイント・ヤコビ計算のため、並列の粒度としては申し分ない特徴を有します。さて、ソースを見てお分かりのとおり、配列化対象のループは、303行目の k ループ以下の三重ループとなります。これは、ループ内で使用されているデータの依存性はありません。正確に言えば、 gosa の総和計算する部分を除き、依存性がないと言えます。一般に、並列化対象ループ内に、データ依存性のある演算が、一つでも含まれている場合、そのループ・ブロックは並列化できません。ただ、今回のような総和計算(sum)や最大値・最小値(max,min) 等の「リダクション演算」は、シミュレーションの計算においては常に付きまとうため、PGIコンパイラソフトウェア側で、並列リダクション用のルーチンを内部的に作成します。コンパイラはこれを組込み、このループも全体として並列化できるようになっています。なお、リダクション演算は、独立の CUDA kernel として実行されます。

 以下のソースには、303行目からのループ・ブロックと、328 行目からのループ・ブロックの二つに分かれています。NVIDIA の CUDA プログラミングの概念で考えると、これらの二つのループ・ブロックの間には、「同期処理」が必要です。詳細に言えば、一時配列 wrk2(I,J,K) と言う配列に値をストアすることが前者のブロックの目的であり、後者は、この配列データを正規な p(I,J,K) 配列にコピーすることが目的です。前者のループ・ブロックで,独立で行われる全ての並列計算処理が全て終了しない限り、後者の計算ブロックへの進めません。並列処理では、こうした状況において「同期をとる」ことが行われます。この同期処理は、CUDAの概念で言う全ての「スレッド・ブロック」間で行う必要があるため、こうした同期は、一旦、kernel での処理を終えて、結果的にこれにて同期処理をする必要があります。従って、PGIコンパイラは、一つの PGI 並列リージョン内で、二つの明示的な kernel プログラムを生成します。(正確には、総和のリダクション処理も一つの独立した kernel 処理となっていますので、三つの CUDA kernel として、PGIコンパイラは生成します。)姫野ベンチマークは、CUDA kernel 的には大体、こう言った構造となります。

 もう一つ、一番外側の 299行目のループは、一般には収束状況を見るための反復計算となりますので、ここの大きなループはシリアル処理をしなければなりません。また、他のシミュレーションでは、一番外側に時間積分のためのループがある場合もあります。外形的なループ構造の形は、多かれ少なかれ同じ形となります。従って、299行目のループでは、並列化を行うことはできません。

 最初に、PGIアクセラレータ・プログラミングモデルのコンパイラ指示行(ディレクティブ)を並列化対象となるループ・ブロックの前(302行)と後(335行)に挿入してコンパイルしてみましょう。一般には、並列化ができる領域と言うことですから、こう言った挿入方法となりますね。この領域を「アクセラレータ計算領域」と言います。(参考、この領域は、「構造化ブロック(Structured block)」でなければならない。)GPU 並列化対象となる領域の前後に、!$acc region と!$acc end region (これを「アクセラレータ・リージョン・ディレクティブ」と言う。)を挿入することで、コンパイラはそのアクセラレータ計算領域を分析し、CUDA kernel を生成します。ここで、!$acc region と !$acc end region の重要な役割の一つを理解しておいて下さい。これは、!$acc region の時点(領域のエントリ)で、GPU側で必要な変数・配列データの転送を行うと言うことと、!$acc end region の時点(領域の出口)で、GPU側で計算された(ストアされた)変数と配列をホスト側のメモリへ転送すると言う操作を行うと言うことです。もし、この「アクセラレータ計算領域」の外側で、ループを回すような処理が行われる場合は、このループ1回毎に、ホストとGPU間のデータ転送が発生します。以下のソース・プログラムの例は、まさにこのタイプです。299 行目の loop インデックスのループ内で「アクセラレータ計算領域」の処理が行われます。以下に述べるように、この実行性能は、CPUで演算するよりも遅い性能を記録しています。結果的に、データ転送の繰り返しで多くの時間を消費していると言うことです。
 しかし、このヤコビ計算の本質を考えると、299行目の loop ループ の前に、1回限り ホスト から GPU へ必要なデータを送り込み、そのデータを使用し302行以降の3次元のヤコビ計算を収束するまでGPU内部で繰り返し、収束したら最後に1回で ホスト へ計算されたデータを返すと言う計算をすれば良い訳です。この方法を次の STEP 2 で説明します。

INFORMATION
PGIアクセラレータ・プログラミング・モデルの基本コンパイラ指示行

PGI Accelerator Compute Region Directive(アクセラレータ計算領域)

$acc region ~ !$acc end region

アクセラレータ計算領域ディレクティブは、GPU上で行う計算領域を明示的に指定するために使用します。その領域の対象となる実行文の前に !$acc region を指定して、領域の最後の実行文の後に、!$acc end region を指定して、「アクセラレータ計算領域」を指定します。一般には、Fortran/C の do/forループ構造を持つ構造化ブロックをその領域とすることが多いです。この領域の中では、call文は使用できません。もし、subroutine を call するような領域となる場合は、そのルーチンをインライン展開するか、あるいは、明示的に、「アクセラレータ計算領域」を複数に分断して、GPU並列化の阻害要因を避ける必要があります。

  • コンパイラは、この領域内のプログラムを分析し、GPU用の並列化ができるものであれば、自動的にGPU CUDA のコードの生成とGPUへのデータ転送の命令をコード化し実行ファイルを作成します。
  • $acc region は、領域のエントリであり、このポイントでGPU側で必要な変数・配列データの転送を行う。
  • $acc end region は、領域の出口であり、このポイントでGPUからホスト側へ変数・配列データの転送(戻し)を行う。
PGF90 (Version     10.6)          07/03/2010  22:53:07
(  283) !*************************************************************
(  284) subroutine jacobi(nn,gosa)
(  285) !*************************************************************
(  286)   use pres
(  287)   use mtrx
(  288)   use bound
(  289)   use work
(  290)   use others
(  291) !
(  292)   implicit none
(  293) !
(  294)   integer,intent(in) :: nn
(  295)   real(4),intent(inout) :: gosa
(  296)   integer :: i,j,k,loop
(  297)   real(4) :: s0,ss
(  298)
(  299)  do loop=1,nn  ! 本来、収束反復のループ(ここはシリアル実行しなければならない)
(  300)    gosa = 0.0
(  301) !アクセラレータ計算領域の開始宣言。
(  302) !$acc region       ! 以下の三重ループは並列化可能。アクセラレータ領域の開始宣言。
(  303)      do k=2,kmax-1       ! 左辺のs0は、一時待避変数。
(  304)         do j=2,jmax-1    ! s0は、一時待避変数ss スカラ値に反映され、
(  305)            do i=2,imax-1 ! wrk2(i,j,k)にストアされる。wrk2は iループ内で依存性無し。
(  306)               s0=a(I,J,K,1)*p(I+1,J,K) &
(  307)                    +a(I,J,K,2)*p(I,J+1,K) &
(  308)                    +a(I,J,K,3)*p(I,J,K+1) &
(  309)                    +b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K) &
(  310)                                -p(I-1,J+1,K)+p(I-1,J-1,K)) &
(  311)                    +b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1) &
(  312)                                -p(I,J+1,K-1)+p(I,J-1,K-1)) &
(  313)                    +b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1) &
(  314)                                -p(I+1,J,K-1)+p(I-1,J,K-1)) &
(  315)                    +c(I,J,K,1)*p(I-1,J,K) &
(  316)                    +c(I,J,K,2)*p(I,J-1,K) &
(  317)                    +c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
(  318)               ss=(s0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
(  319)               gosa= gosa + ss*ss              ! gosa 総和をとるリダクション演算
(  320)               wrk2(I,J,K)=p(I,J,K)+OMEGA *SS  ! gosa の演算は本来、並列化できないが、
(  321)            enddo                               ! コンパイラが、並列reductionに置き換える
(  322)         enddo
(  323)      enddo
(  324)   ! F90の「配列構文」からDOループ構造へ変更する
(  325) !     p(2:imax-1,2:jmax-1,2:kmax-1)= &
(  326) !         wrk2(2:imax-1,2:jmax-1,2:kmax-1)
(  327)   ! これにより今後、kernel tuning の directive を入れることが可能
(  328)      do k=2,kmax-1
(  329)         do j=2,jmax-1
(  330)           do i=2,imax-1
(  331)               p(I,J,K)=wrk2(I,J,K)
(  332)            enddo
(  333)         enddo
(  334)      enddo
(  335)  !$acc end region
(  336)  !アクセラレータ計算領域の終了宣言。
(  337)  enddo
(  338)
(  339) !! End of iteration
(  340)
(  341)   return
(  342) end subroutine jacobi

 PGIアクセラレータ・コンパイラで、コンパイルしてみましょう。以下は、用意した Makefile でコンパイラが出力した「PGIアクセラレータ翻訳」に関するメッセージです。以下のメッセージのリスト内に、コンパイラが行った処理に対して、その意味を説明しています。大雑把に説明すると、まず、302行目のアクセラレータ計算領域開始のディレクティブ(C言語ではプラグマ)の指定ポイント(これは領域、あるいはリージョンの entry pointと言う。)で、並列領域内で使用するデータ(配列や変数)を ホスト のメモリから GPU のデバイス・メモリへコピーする操作を行います。これは、コンパイラが自動的に使用する配列・変数を分析し、最も少ないデータ転送量となるように、データ転送命令を生成します。「最も少ないデータ転送量」と言う処理がコンパイラのデフォルト操作となります。以下の例で言うと、例えば、(wrk1(2:imax-1,2:jmax-1,2:kmax-1))と言うように、計算で必要な配列の添字の部分を転送する(2からimax-1と言った部分配列の転送)と言う意味です。すなわち、部分配列の転送ですのでデバイス・メモリ上での配列の割付は、ホスト上での割付レイアウトとは異なる形になります。このことは続編で話題にしますが、デバイス・メモリのアクセス性能に係わってくる場合があります。また、これについては、以前のコラムでも説明しています。

 次に生成された CUDA kernel に関する情報が記述されています。もし、コンパイラが自動で並列性を抽出できない(並列化を阻害する要因がある)場合は、kernel を生成できなかったと言うメッセージ(No parallel kernels found. accelerator region ignored.) を出します。以下の例では、大きく二つの kernel が生成されたことを告げ、それぞれの Kernelのスレッド・ブロックやグリッドの構成を記しています。また、PGI 10.5 以降は、生成したスレッド・ブロックの場合の使用するスレッド当たりのレジスタ量、shared memory の使用量やCUDA Occupancy(Warp の占有率)に関する情報を表示します。これは、CUDA Kernel パラメータ(スレッド・ブロックサイズ等)の調整に活かすことができます。

●コンパイル時の GPU翻訳に関するメッセージ -Minfo=accel 
[kato@photon29 STEPS]$ make step1.exe
pgfortran -o step1.exe step1.F90 -fastsse -Minfo=accel -ta=nvidia,cuda3.0,cc13,cc20,time
jacobi: ! ホストのメモリからGPUのメモリへ部分配列のデータ転送
    302, Generating copyin(wrk1(2:imax-1,2:jmax-1,2:kmax-1))
         Generating copyin(c(2:imax-1,2:jmax-1,2:kmax-1,1:3))
         Generating copyin(b(2:imax-1,2:jmax-1,2:kmax-1,1:3))
         Generating copyin(a(2:imax-1,2:jmax-1,2:kmax-1,1:4))
         Generating copyin(bnd(2:imax-1,2:jmax-1,2:kmax-1))
         Generating copyin(p(1:imax,1:jmax,1:kmax))
         Generating copyout(p(2:imax-1,2:jmax-1,2:kmax-1))
         Generating copyout(wrk2(2:imax-1,2:jmax-1,2:kmax-1))
         Generating compute capability 1.0 binary
         Generating compute capability 1.3 binary
         Generating compute capability 2.0 binary
    303, Loop is parallelizable ! 行番号303~305ループまで並列化可能
    304, Loop is parallelizable
    305, Loop is parallelizable
         Accelerator kernel generated ! CUDA kernel プログラムを生成した
        303, !$acc do parallel, vector(4) ! kernel scheduling parameter を自動設定した
        304, !$acc do parallel, vector(4)  ! Block=(4x4x16)のサイズ
        305, !$acc do vector(16)
             Cached references to size [18x6x6] block of 'p'
             ! Shared Memoryに[18x6x6]のブロックのキャッシュを設定
             ! 以下は、このスレッド・ブロック特性の場合のOccupancyを支配するパラメータを表示
             CC 1.0 : 32 registers; 2616 shared, 204 constant, 0 local memory bytes; 33 occupancy
             CC 1.3 : 32 registers; 2616 shared, 204 constant, 0 local memory bytes; 50 occupancy
             CC 2.0 : 34 registers; 2600 shared, 200 constant, 0 local memory bytes; 50 occupancy
    328, Loop is parallelizable
    329, Loop is parallelizable
    330, Loop is parallelizable ! もう一つの CUDA kernel プログラムを生成した
         Accelerator kernel generated
        319, Sum reduction generated for gosa
        328, !$acc do parallel, vector(4)
        329, !$acc do parallel, vector(4)
        330, !$acc do vector(16)
             CC 1.3 : 13 registers; 24 shared, 188 constant, 0 local memory bytes; 100 occupancy
             CC 2.0 : 13 registers; 8 shared, 196 constant, 0 local memory bytes; 100 occupancy

 以下に、STEP 1 の実行結果を示しました。この場合の実行性能は、1288 MFLOPS です。一般的な CPU の性能よりも悪い性能です。この理由は、前述の通り、ホストとGPU間のデータの転送を外側のループ1回ごとに行われているため、その時間のほとんどをデータ転送に費やされたことによるものです。従って、PGIアクセラレータ計算領域は、ループ繰り返しの内部に置いてはなりません。さて、以下の出力のうち、「Accelerator Kernel Timing data」と言うタイトルで、GPU 処理のプロファイル情報が記されています。302行開始のアクセラレータ計算領域の総合時間情報が最初に記されています。ここで、kernels=4406962 と表記された時間は、純粋に GPU の演算器を利用した計算時間を示します。また、その隣の data=77023164 は、データ転送の総消費時間を表し 77秒ほど転送に費やしたと言うことになります。それでは、このデータ転送をアクセラレータ計算領域の前後で1回限り行って本来の GPU 活用を行う方法を次の STEP 2 で説明します。

●STEP1 の実行結果
[kato@photon29 STEPS]$ step1.exe
 For example:
 Grid-size=
            XS (64x32x32)
            S  (128x64x64)
            M  (256x128x128)
            L  (512x256x256)
            XL (1024x512x512)
  Grid-size =
  mimax=          257  mjmax=          129  mkmax=          129
  imax=          256  jmax=          128  kmax=          128
  Time measurement accuracy : .10000E-05
  Start rehearsal measurement process.
  Measure the performance in 3 times. ! 3回の実行でトライアル・ラン
 initialize nvidia GPU
   MFLOPS:    1280.981       time(s):   0.3210940000000000        1.6939555E-03
 Now, start the actual measurement process.
 The loop will be excuted in          800  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for           800  times  800回の反復回数に固定しています。
  Gosa :   8.3822059E-04
  MFLOPS:    1288.662       time(s):    85.11471800000000
  Score based on Pentium III 600MHz :    15.55604

Accelerator Kernel Timing data
/home/kato/GPGPU/Himeno/STEPS/step1.F90 
  jacobi ! GPUの実行プロファイル情報(時間単位はμ秒)
    302: region entered 803 times 
        time(us): total=85435123 init=232 region=85434891
                  kernels=4406962 data=77023164 ! ホストとGPU間のデータ転送に 77秒掛かっている
                  ! kernel部分の総実行時間は 4.4秒
        w/o init: total=85434891 max=108401 min=106206 avg=106394
        305: kernel launched 803 times
            grid: [32x32]  block: [16x4x4]! グリッドとブロックのサイズ
            time(us): total=4147047 max=5211 min=5125 avg=5164 ! kernel の実行時間特性、平均5.1msec/1反復
        319: kernel launched 803 times
            grid: [1]  block: [256]
            time(us): total=8048 max=14 min=10 avg=10
        330: kernel launched 803 times
            grid: [32x32]  block: [16x4x4]
            time(us): total=251867 max=322 min=310 avg=313
acc_init.c
  acc_init
    41: region entered 1 time ! GPUボードの初期化オーバーヘッドは 0.57秒
        time(us): init=573803

ステップ 2 (反復ループの外側を !$acc region /end region で囲む)

 STEP 1 では、GPUを使用する上で悪い例を示しました。非圧縮性流体の Poisson 方程式ソルバーの一つである姫野ベンチマーク本来のロジックでは、GPU内部で使用できるデータは、1回GPU側へ転送して、収束した結果を最後にホスト側に戻すことを行えば良い訳で、収束のための計算本体は GPU 内部のデバイス・メモリ上のデータを活用できることになります。(注意:姫野ベンチはベンチマーク用途コードですので、収束判定を行わず単なる繰り返しをしている構造となっています。実際は、もちろん収束判定のロジックが含まれても並列化は可能です。)こうした動作を指示するには、アクセラレータ計算領域の開始(!$acc region)を 301行目の do loop=1,nn の反復ループの前に挿入します。「アクセラレータ計算領域の開始」は前述したとおり、GPU側へオフロードする領域のエントリを意味するもので、主な目的の一つは、GPUで必要なデータの転送を行うポイントであると言うことです。「真に並列化できる対象の開始ポイント」を指すことではありません。301 行目の上位で、!$acc regionを指定すると、望みとおりのデータ転送が 1 回限り行われます。しかし、待てよ、 do loop=1,nn と言う反復計算の外側ループは、並列実行してもらうと困るループなので、コンパイラは並列化を施さないのではないかと言う疑問が生まれるかもしれません。PGIアクセラレータ用の指示行(ディレクティブ)には、さらに、アクセラレータ・ループ・マッピング・ディレクティブと言う、ループの並列化を細かに指示するための指示行が存在します。その中の一つに、!$acc do hostと言う指示があります。アクセラレータ計算領域内にはあるが、この指定されたループは、ホスト側の制御を行い、シーケンシャルに実行しなさいと言う意味となります。これにてコンパイルするとどうでしょうか。(実は、ユーザがこうしたディレクティブを入れなくても、コンパイラは自動的に分析して以下と同じ並列化を行います)

INFORMATION
PGIアクセラレータ・プログラミング・モデルの基本コンパイラ指示行

PGI Accelerator Loop Mapping Directives (ループマッピング)

!$acc do [clause [,clause]...]

  • $acc do ... は、ループ構文の直前に挿入し、そのループがどのような形で並列化を行うかを指示するもの。並列のベクトル長(SIMD 幅)や並列分割数(グリッド、ブロックのサイズ)を指定できる。
  • 以下に使用する $acc do host は、GPU上での並列処理はせず、ホスト側制御によるシーケンシャル実行をすることを指示するものです。

ループマッピング・ディレクティブの詳細説明はこちらにあります

【注意】
PGI 10.5 では、!$acc do host ディレクティブの機能に関してバグがあります。
PGI 10.6 で解決されています。このディレクティブを使用する場合は、
PGI 10.4以前 か 10.6 以降のバージョンをご利用下さい。
(  299) !$acc region ! ここのポイントでGPUへのデータの転送が行われる。
(  300) !$acc do host
(  301)  do loop=1,nn   ! ホスト側の制御のループとして指示する。シーケンシャル実行。
(  302)    gosa = 0.0
(  303)         ! 以下の3重ループは、GPU並列処理される。
(  304)      do k=2,kmax-1
(  305)         do j=2,jmax-1
(  306)            do i=2,imax-1
(  307)               s0=a(I,J,K,1)*p(I+1,J,K) &
(  308)                    +a(I,J,K,2)*p(I,J+1,K) &
(  309)                    +a(I,J,K,3)*p(I,J,K+1) &
(  310)                    +b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K) &
(  311)                                -p(I-1,J+1,K)+p(I-1,J-1,K)) &
(  312)                    +b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1) &
(  313)                                -p(I,J+1,K-1)+p(I,J-1,K-1)) &
(  314)                    +b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1) &
(  315)                                -p(I+1,J,K-1)+p(I-1,J,K-1)) &
(  316)                    +c(I,J,K,1)*p(I-1,J,K) &
(  317)                    +c(I,J,K,2)*p(I,J-1,K) &
(  318)                    +c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
(  319)               ss=(s0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
(  320)               gosa= gosa + ss*ss
(  321)               wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
(  322)            enddo
(  323)         enddo
(  324)      enddo
(  325)
(  326) !     p(2:imax-1,2:jmax-1,2:kmax-1)= &
(  327) !         wrk2(2:imax-1,2:jmax-1,2:kmax-1)
(  328)
(  329)      do k=2,kmax-1
(  330)         do j=2,jmax-1
(  331)           do i=2,imax-1
(  332)               p(I,J,K)=wrk2(I,J,K)
(  333)            enddo
(  334)         enddo
(  335)      enddo
(  336)
(  337)  enddo
(  338) !$acc end region
!$acc end region  ! ここのポイントでGPU から ホストへのデータの転送が行われる。

 以下のリストは、アクセラレータ関係のコンパイル・メッセージです。詳細は、以下にインラインで説明していますが、アクセラレータ計算領域内に、本来並列化できない loop (301行)がありますので、これに関する警告メッセージが表示されています。コンパイラはもちろん、このループをシーケンシャルに実行すると言う形で、コードを生成します。304行目以降は、GPU上でスレッド並列実行を行うためのコードを生成した際の情報を表示しています(これは上述 STEP 1 の場合と同じです)。

● コンパイル時の GPU翻訳に関するメッセージ -Minfo=accel 
[kato@photon29 STEPS]$ make step2.exe
pgfortran -c -fastsse  -Minfo=accel -ta=nvidia,cuda3.0,cc13,cc20,time himenoBMTxp.F90
jacobi:  ! ホストのメモリからGPUのメモリへ部分配列のデータ転送
    299, Generating copyout(wrk2(2:imax-1,2:jmax-1,2:kmax-1))
         Generating copyin(p(1:imax,1:jmax,1:kmax))
         Generating copyout(p(2:imax-1,2:jmax-1,2:kmax-1))
         Generating copyin(bnd(2:imax-1,2:jmax-1,2:kmax-1))
         Generating copyin(a(2:imax-1,2:jmax-1,2:kmax-1,1:4))
         Generating copyin(b(2:imax-1,2:jmax-1,2:kmax-1,1:3))
         Generating copyin(c(2:imax-1,2:jmax-1,2:kmax-1,1:3))
         Generating copyin(wrk1(2:imax-1,2:jmax-1,2:kmax-1))
         Generating compute capability 1.0 binary
         Generating compute capability 1.3 binary
         Generating compute capability 2.0 binary
          ! 以下のメッセージは 301 の反復ループでの並列化は依存性があるためできないと言うことの警告
           ! 並列化するつもりはないため、以下のメッセージは無視する
    301, Loop carried dependence due to exposed use of 'p(2:imax-1,2:jmax-1,2:kmax-1)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(1:imax-1,3:jmax,2:kmax-1)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(3:imax,1:jmax,2:kmax-1)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(1:imax-2,1:jmax-2,2:kmax-1)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(2:imax-1,3:jmax,1:kmax-2)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(2:imax-1,2:jmax,3:kmax)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(3:imax,2:jmax-1,3:kmax)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(1:imax-2,2:jmax-1,1:kmax)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(2:imax-1,1:jmax-2,1:kmax)' prevents parallelization
         Loop carried dependence due to exposed use of 'p(2:imax,2:jmax-1,1:kmax-2)' prevents parallelization
          ! 以下のメッセージも並列化の対象と考えていないため、無視する。
         Accelerator restriction: scalar variable live-out from loop: gosa
         Parallelization would require privatization of array 'wrk2(2:imax-1,i3+2,i2+2)'
          ! コンパイラは以上の並列化阻害要因の分析から並列化を行わず、
          ! !$acc do hostの指示に従い、host 側でシーケンシャル実行する。
         Sequential loop scheduled on host
    304, Loop is parallelizable ! 行番号304~306ループまで並列化可能
    305, Loop is parallelizable
    306, Loop is parallelizable
         Accelerator kernel generated ! CUDA kernel プログラムを生成した
        304, !$acc do parallel, vector(4)
        305, !$acc do parallel, vector(4)
        306, !$acc do vector(16)
             Cached references to size [18x6x6] block of 'p'
             CC 1.0 : 32 registers; 2624 shared, 204 constant, 0 local memory bytes; 33 occupancy
             CC 1.3 : 32 registers; 2624 shared, 204 constant, 0 local memory bytes; 50 occupancy
             CC 2.0 : 34 registers; 2600 shared, 204 constant, 0 local memory bytes; 50 occupancy
    329, Loop is parallelizable
    330, Loop is parallelizable
    331, Loop is parallelizable
         Accelerator kernel generated ! CUDA kernel プログラムを生成した
        320, Sum reduction generated for gosa
        329, !$acc do parallel, vector(4)
        330, !$acc do parallel, vector(4)
        331, !$acc do vector(16)
        	 CC 1.0 : 13 registers; 28 shared, 188 constant, 0 local memory bytes; 66 occupancy
             CC 1.3 : 13 registers; 28 shared, 188 constant, 0 local memory bytes; 100 occupancy
             CC 2.0 : 13 registers; 8 shared, 200 constant, 0 local memory bytes; 100 occupancy

 以下に、実行結果を示しました。実行性能は、24.2GFLOPS に一気に向上しました。この理由は、ホストとGPU間のデータの転送が反復ループの外側で1回のみ行われたことによるデータ転送時間の大幅な削減によるものです。以下の出力のうち、「Accelerator Kernel Timing data」と言う、GPU 処理のプロファイル情報を見て下さい。299行開始のアクセラレータ計算領域の総合時間情報が最初に記されています。ここで、kernels=4402983(μ秒) と表記された時間は、純粋に GPU の演算器を利用した計算時間を示します。また、その隣の data=193225 は、ホストとGPU間のデータ転送の総消費時間を表し 0.19秒のみ転送に費やしたと言うことになります。前述の STEP 1 の時は、77秒 であったことから考えると大幅な削減となります。これによって、実効性能が大幅に向上したことになります。次のステップ 3 では、ホストとGPU間のデータ転送時間の最小化を行います。

●STEP 2 の実行結果
 Now, start the actual measurement process.
 The loop will be excuted in          800  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for           800  times
  Gosa :   8.3822059E-04
  MFLOPS:    24197.57       time(s):    4.532856000000000
  Score based on Pentium III 600MHz :    292.1001

Accelerator Kernel Timing data
/home/kato/GPGPU/Himeno/STEPS/step2.F90
  jacobi
    299: region entered 2 times  ! トライアル実行と本番実行の2回のみ
        time(us): total=4651760 init=1 region=4651759
                  kernels=4402983 data=193225 ! ホストとGPU間の双方向データ転送は 0.192秒のみ
                  ! kernel部分の総実行時間は 4.4秒
        w/o init: total=4651759 max=4532854 min=118905 avg=2325879
        306: kernel launched 803 times  ! 並列ブロックは 803回実行(トライアル実行3回含む)
            grid: [32x32]  block: [16x4x4]
            time(us): total=4142545 max=5200 min=5121 avg=5158
        320: kernel launched 803 times
            grid: [1]  block: [256]
            time(us): total=8034 max=13 min=10 avg=10 
        331: kernel launched 803 times
            grid: [32x32]  block: [16x4x4]
            time(us): total=252404 max=318 min=311 avg=314
acc_init.c
  acc_init
    41: region entered 1 time
        time(us): init=571633

ステップ 3 (明示的にホストとGPU間のデータ転送を指定するディレクティブ)

 STEP 2の場合は、コンパイラがプログラムを分析し、ホストとGPU間のデータ移動すべき配列、変数とその添字範囲を特定して自動でデータ転送のコードを作成しました。コンパイラのデフォルトの動作は、プログラムのロジックのとおり、例えば、GPU上で実行に必要な添字範囲(部分配列)だけを抜き出して、データの転送を行うための命令を作成します。また、GPUからホスト側にデータを戻す場合も、GPU 上で値が更新された変数あるいは、値が更新された配列(部分的に更新された場合はその部分配列のみ)をホスト側にコピーする命令を作成します。上述の STEP 1,2のコンパイルメッセージを見ると分かるとおり、Generating copyout(a(2:imax-1,2:jmax-1,2:kmax-1)) と言った形で、配列の全ての要素ではなく、ループ内で使用している添字範囲の部分配列のみを転送する命令をコード化します。ホストとGPU間のデータの転送は、いわゆるPCIバスを利用した I/O 処理であるため低速なものです。データ転送のオーバーヘッドを小さくするには、できるだけ I/O の DMAモードが有効に機能する転送を行うことが望ましいのですが、これはメモリのデータを連続で転送する場合に有効になります。部分配列の転送の場合は、データの跳びが発生するため、多少、転送性能が劣ります。但し、データ転送が無視できる程のものであれば、あまり気にする必要はありません。
 PGIアクセラレータのディレクティブ・プラグマにはホスト~GPU間のデータの転送を明示的に指定できるディレクティブがあります。このディレクティブを使用して、明示的に配列の添字範囲を指定して、コンパイラにその転送の形態を指示できます。一般的には、配列の一部ではなく全部を転送するように指定すれば、連続アクセス転送となりますので転送時間が小さくなります。  

 以下の例は、!$acc region アクセラレータ計算領域ディレクティブの中に備わる、データ転送を定義する 節(clause)を使用して転送定義した例です。姫野ベンチマークの配列の定義は、ベクトル型スパコン時代のメモリ・バンクコンフリクトを避けるために、各配列の添字を1ずつ増やして配列割付(mimax,mjmax,mkmax) を行っています。実際に使用する範囲は、1:imax,1:jmax,1:kmax と言う変数で定義されています。こうした場合、連続アクセスデータの転送を行う場合は、実際の配列割付を行っている添字の最大値 mimax,mjmax,mkmax を指定して転送指示を行う必要があります。

INFORMATION
PGIアクセラレータ・プログラミング・モデルの基本コンパイラ指示行

PGI Accelerator Compute Region Directive(アクセラレータ計算領域)のクローズ

!$acc region [clause [,clause]...]

アクセラレータ計算領域ディレクティブのデータ転送定義のためのクローズ(節)として、主に以下のものが用意されています。!$acc region の後に続いて指定します。複数のクローズ(節)をブランクを区切って指定できます。

  • "copy (list)" は、アクセラレータ計算領域の入口でホストからGPUへlistで指定された配列を転送し、アクセラレータ計算領域の出口で、その配列をホスト側に戻すと言う処理を指示するものです。
  • "copyin (list)"は、アクセラレータ計算領域の入口でホストからGPUへlistで指定された配列を転送することのみを指示するものです。
  • "copyout (list)"は、アクセラレータ計算領域の出口で、指定された list の配列をホスト側に戻すと言う処理を指示するものです。
  • "local (list)"は、アクセラレータ計算領域内で、ローカルに使用することを宣言し、これをコンパイラに指示します。GPU側の配列の割付を行いますが、データの転送は行いません。
    ! 反復ループの前でGPUへのデータの転送が行われる。
    ! ホスト側で割付けた各配列をその添字 1 から最後まで全て転送する。
    ! 転送範囲が、1:imax,1:jmax,1:kmax ではないことに注意。
    ! 連続データとして DMA モードで転送できる(速い)。
    ! wrk2配列はループ内のみで使用する一時変数のため、!$acc localを指定する。
(  299) !$acc region &
(  300) !$acc copyin (a(1:mimax,1:mjmax,1:mkmax,1:4)) &
(  301) !$acc copyin (b(1:mimax,1:mjmax,1:mkmax,1:3), c(1:mimax,1:mjmax,1:mkmax,1:3))  &
(  302) !$acc copy   (p(1:mimax,1:mjmax,1:mkmax))  &
(  303) !$acc copyin (wrk1(1:mimax,1:mjmax,1:mkmax), bnd(1:mimax,1:mjmax,1:mkmax)) &
(  304) !$acc local  (wrk2(1:mimax,1:mjmax,1:mkmax)) &
(  305) !$acc copyin(omega),copy(gosa),local(ss,s0)
(  306)
(  307) !$acc do host
(  308)  do loop=1,nn ! 反復ループ
(  309)    gosa = 0.0
(  310)
(  311)      do k=2,kmax-1
(  312)         do j=2,jmax-1
(  313)            do i=2,imax-1

          (省略)
  
(  344)  enddo
(  345) !$acc end region
    ! 反復ループが終了したこのポイントでGPU から ホストへのデータの転送が行われる。

 以下のリストは、データ転送の部分のコンパイル・メッセージを抜き出したものです。ディレクティブの指示通り、配列の下限~上限まで全ての要素が転送されています。

● コンパイル時の GPU翻訳に関するメッセージ -Minfo=accel 
[kato@photon29 STEPS]$ make step3.exe
pgfortran -o step3.exe step3.F90 -fastsse -Minfo=accel -ta=nvidia,cuda3.0,time
jacobi:  ! ホストのメモリからGPUのメモリへ全体配列のデータ転送
    299, Generating local(wrk2(1:mimax,1:mjmax,1:mkmax))
         Generating copy(p(1:mimax,1:mjmax,1:mkmax))
         Generating copyin(bnd(1:mimax,1:mjmax,1:mkmax))
         Generating copyin(a(1:mimax,1:mjmax,1:mkmax,1:4))
         Generating copyin(b(1:mimax,1:mjmax,1:mkmax,1:3))
         Generating copyin(c(1:mimax,1:mjmax,1:mkmax,1:3))
         Generating copyin(wrk1(1:mimax,1:mjmax,1:mkmax))
         Generating compute capability 1.0 binary
         Generating compute capability 1.3 binary
         Generating compute capability 2.0 binary
         (以下、省略) STEP2の場合と同じ
     

 以下に、実行結果を示しました。実行性能は、24.41GFLOPS と STEP 2 に較べて若干ではありますが性能が向上しています。ただ、誤差範囲の向上とも言える程度ですが、ここで着目して欲しいのは、data=86256 の時間です。これは、ホストとGPU間のデータ転送の総消費時間を表します。STEP 2 の時は data=193225 であったので、転送時間が半分以下になったと言うことです。実際の転送バイト数は、全配列要素を転送する STEP 3 の方が多いのですが、DMA 転送が効いて転送時間の向上がなされたようです。どうしてもホストとGPU間の転送を数多く行わなければいけないロジックでは、これは GPGPU 処理のオーバーヘッドとなるため、こうした細かなチューニングが必要となります。

●STEP 3 の実行結果
 Now, start the actual measurement process.
 The loop will be excuted in          800  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for           800  times
  Gosa :   8.3822059E-04
  MFLOPS:    24418.85       time(s):    4.491779999999999
  Score based on Pentium III 600MHz :    294.7713

Accelerator Kernel Timing data
/home/kato/GPGPU/Himeno/STEPS/step3.F90
  jacobi
    299: region entered 2 times  ! トライアル実行と本番実行の2回のみ
        time(us): total=4557715
                  kernels=4415100 data=86256 ! ホストとGPU間の双方向データ転送は 0.086秒
                  ! kernel部分の総実行時間は 4.4秒
        306: kernel launched 803 times  ! 並列ブロックは 803回実行(トライアル実行3回含む)
            grid: [32x32]  block: [16x4x4]
            time(us): total=4142807 max=5205 min=5127 avg=5159
        327: kernel launched 803 times
            grid: [1]  block: [256]
            time(us): total=8043 max=13 min=10 avg=10
        338: kernel launched 803 times
            grid: [32x32]  block: [16x4x4]
            time(us): total=264250 max=334 min=325 avg=329
acc_init.c
  acc_init
    41: region entered 1 time
        time(us): init=562479

ステップ 4 (STEP 3 と同等な形態であるが、ディレクティブの使い方が異なる例)

 PGIアクセラレータ・プログラミング・モデルにおいて、STEP 3 の場合と同等でありながら、別の方式のディレクティブの使い方をした場合の例をお見せします。一般的には、こちらの方法の方がプログラム・フローにマッチしたコンパイラの指示行の使い方と思いますが、STEP 3 よりも若干だけ性能が劣るようです。  

 前のステップ 3 では、!$acc region 「アクセラレータ計算領域」ディレクティブの中に備わるデータ転送用のクローズ(Cluase) を使用して、データ転送の定義を行いました。このステップ 4 は、別の「アクセラレータデータ領域」ディレクティブを使用して、「データの転送」の指示と「GPU並列計算」の指示を明確に分けて行う方法です。ホストとGPU間でデータの移動を行う配列の宣言を「アクセラレータ計算領域」の前で、「アクセラレータデータ領域」ディレクティブを使用して行います。これによって、「アクセラレータ計算領域」内では、ホストとGPU間の配列データの転送を行わず、GPUのデバイス・メモリ内にデータを留めたまま、繰り返しの処理を行うことができます。一般のシミュレーションのロジック・フローはこうした形が多いですので、このような方法でホストとGPU間でデータ転送時間を極力少なくすることが可能です。後は、下のリストに示すように、309行目の反復ループの内側に、「アクセラレータ計算領域」が存在するような形になっています。最も自然なディレクティブの入れ方です。

INFORMATION
PGIアクセラレータ・プログラミング・モデルの基本コンパイラ指示行

PGI Accelerator Data Region Directive(アクセラレータデータ領域)

!$acc data region [clause [,clause]...] ~ !$acc end data region

アクセラレータデータ領域ディレクティブは、ホストとGPU間のデータ転送を行う場所を明示的に指定するために使用します。一般に、!$acc region (アクセラレータ計算領域の定義)よりも外側で指定します。「アクセラレータデータ領域」の内側には、複数の「アクセラレータ計算領域」を有していても構いません。このディレクティブを指定すると、その内側にある「アクセラレータ計算領域」でのデータの転送は行われません。データ転送定義のためのクローズ(節)として、主に以下のものが用意されています。!$acc data region の後に続いて指定します。複数のクローズ(節)をブランクを区切って指定できます。

  • "copy (list)" は、アクセラレータデータ領域の入口でホストからGPUへlistで指定された配列を転送し、アクセラレータデータ領域の出口で、その配列をホスト側に戻すと言う処理を指示するものです。
  • "copyin (list)"は、アクセラレータデータ領域の入口でホストからGPUへlistで指定された配列を転送することのみを指示するものです。
  • "copyout (list)"は、アクセラレータデータ領域の出口で、指定された list の配列をホスト側に戻すと言う処理を指示するものです。
  • "local (list)"は、アクセラレータデータ領域内で、ローカルに使用することを宣言し、これをコンパイラに指示します。GPU側の配列の割付を行いますが、データの転送は伴いません。
    ! 「アクセラレータデータ領域」の指示行で配列データの転送を指示。
    ! 反復ループの前でGPUへのデータの転送が行われる。
(  299) !$acc data region &
(  300) !$acc copyin (a(1:mimax,1:mjmax,1:mkmax,1:4)) &
(  301) !$acc copyin (b(1:mimax,1:mjmax,1:mkmax,1:3), c(1:mimax,1:mjmax,1:mkmax,1:3))  &
(  302) !$acc copy   (p(1:mimax,1:mjmax,1:mkmax))  &
(  303) !$acc copyin (wrk1(1:mimax,1:mjmax,1:mkmax), bnd(1:mimax,1:mjmax,1:mkmax)) &
(  304) !$acc local  (wrk2(1:mimax,1:mjmax,1:mkmax))
(  305)
(  306)  do loop=1,nn  ! 反復ループ(シーケンシャル実行)
(  307)    gosa = 0.0
(  308)        ! ここが「アクセラレータ計算領域」の入口で、その指定を行う。
              ホストから GPU へのスカラ変数等の転送の記述のみ行う。
(  309) !$acc region copyin(omega),copy(gosa),local(ss,s0)
(  310)      do k=2,kmax-1
(  311)         do j=2,jmax-1
(  312)            do i=2,imax-1
(  313)               s0=a(I,J,K,1)*p(I+1,J,K) &
(  314)                    +a(I,J,K,2)*p(I,J+1,K) &
(  315)                    +a(I,J,K,3)*p(I,J,K+1) &
(  316)                    +b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K) &
(  317)                                -p(I-1,J+1,K)+p(I-1,J-1,K)) &
(  318)                    +b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1) &
(  319)                                -p(I,J+1,K-1)+p(I,J-1,K-1)) &
(  320)                    +b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1) &
(  321)                                -p(I+1,J,K-1)+p(I-1,J,K-1)) &
(  322)                    +c(I,J,K,1)*p(I-1,J,K) &
(  323)                    +c(I,J,K,2)*p(I,J-1,K) &
(  324)                    +c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
(  325)               ss=(s0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
(  326)               gosa= gosa + ss*ss
(  327)               wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
(  328)            enddo
(  329)         enddo
(  330)      enddo
(  331)
(  332) !     p(2:imax-1,2:jmax-1,2:kmax-1)= &
(  333) !         wrk2(2:imax-1,2:jmax-1,2:kmax-1)
(  334)
(  335)      do k=2,kmax-1
(  336)         do j=2,jmax-1
(  337)           do i=2,imax-1
(  338)               p(I,J,K)=wrk2(I,J,K)
(  339)            enddo
(  340)         enddo
(  341)      enddo
(  342) !$acc end region ! 「アクセラレータ計算領域」の終了ポイントだが,
              GPUからホストへのデータの転送が行われない。
(  343)
(  344)  enddo
(  345) !$acc end data region !「アクセラレータデータ領域」の出口
    ! end data regionを指定したこのポイントで、GPUからホストへのデータの転送が行われる。

 コンパイラのメッセージも、STEP2/STEP 3 のような余計な警告メッセージが出ず、素直な GPU 並列化のためのメッセージとなっています。

● コンパイル時の GPU翻訳に関するメッセージ -Minfo=accel 
[[kato@photon29 STEPS]$ make step4.exe
pgfortran -o step4.exe step4.F90 -fastsse -Minfo=accel -ta=nvidia,cuda3.0,time
jacobi:             ! 306行ループの前にデータの転送を予め行う。
    299, Generating local(wrk2(1:mimax,1:mjmax,1:mkmax))
         Generating copyin(bnd(1:mimax,1:mjmax,1:mkmax))
         Generating copyin(wrk1(1:mimax,1:mjmax,1:mkmax))
         Generating copyin(c(1:mimax,1:mjmax,1:mkmax,1:3))
         Generating copyin(b(1:mimax,1:mjmax,1:mkmax,1:3))
         Generating copyin(a(1:mimax,1:mjmax,1:mkmax,1:4))
         Generating copy(p(1:mimax,1:mjmax,1:mkmax))
    309, Generating compute capability 1.0 binary
         Generating compute capability 1.3 binary
         Generating compute capability 2.0 binary
              ! 310行ループ以降、アクセラレータ計算領域の指定がされている。
    310, Loop is parallelizable
    311, Loop is parallelizable
    312, Loop is parallelizable
         ! コンパイラが自動的にGPU 並列化を行った際の CUDA kernel情報。
         Accelerator kernel generated
        310, !$acc do parallel, vector(4)
        311, !$acc do parallel, vector(4)
        312, !$acc do vector(16)
             Cached references to size [18x6x6] block of 'p'
             CC 1.0 : 32 registers; 2616 shared, 196 constant, 4 local memory bytes; 33 occupancy
             CC 1.3 : 33 registers; 2616 shared, 196 constant, 0 local memory bytes; 25 occupancy
             CC 2.0 : 34 registers; 2600 shared, 192 constant, 0 local memory bytes; 50 occupancy
    335, Loop is parallelizable
    336, Loop is parallelizable
    337, Loop is parallelizable
         Accelerator kernel generated
        326, Sum reduction generated for gosa
        335, !$acc do parallel, vector(4)
        336, !$acc do parallel, vector(4)
        337, !$acc do vector(16)
             CC 1.0 : 10 registers; 24 shared, 180 constant, 0 local memory bytes; 100 occupancy
             CC 1.3 : 10 registers; 24 shared, 180 constant, 0 local memory bytes; 100 occupancy
     

 以下に、実行結果を示しました。実行性能は、24.33GFLOPS と STEP 3 に較べて若干ではありますが性能が低下しています。STEP 3 と STEP 4 では、ディレクティブの指示の方法が違いますので、コンパイラは異なるコードを生成します。これに伴う、若干の性能差と見ることができます。

●STEP 4 の実行結果
 Now, start the actual measurement process.
 The loop will be excuted in          800  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for           800  times
  Gosa :   8.3822059E-04
  MFLOPS:    24334.20       time(s):    4.507404999999999
  Score based on Pentium III 600MHz :    293.7495

Accelerator Kernel Timing data
/home/kato/GPGPU/Himeno/STEPS/step4.F90
  jacobi
    309: region entered 803 times
        time(us): total=4476331 init=61 region=4476270
                  kernels=4417192 data=0
        w/o init: total=4476270 max=5770 min=5537 avg=5574
        312: kernel launched 803 times   ! jacobi計算の核心の部分の時間情報。
            grid: [32x32]  block: [16x4x4]  ! CUDA Grid/Grid の形
           time(us): total=4145260 max=5206 min=5123 avg=5162
        326: kernel launched 803 times
            grid: [1]  block: [256]
            time(us): total=8042 max=13 min=10 avg=10
        337: kernel launched 803 times
            grid: [32x32]  block: [16x4x4]
            time(us): total=263890 max=334 min=325 avg=328
/home/kato/GPGPU/Himeno/STEPS/step4.F90
  jacobi    ! アクセラレータデータ領域領域の情報。
    299: region entered 2 times
        time(us): total=4572383 init=2 region=4572381
                   data=85644 ! データ転送時間
        w/o init: total=4572381 max=4507402 min=64979 avg=2286190
acc_init.c
  acc_init
    41: region entered 1 time
        time(us): init=562202

ステップ 5 (STEP 3の形態、但しデータ領域ディレクティブを使用する)

 STEP 3 は、!$acc region アクセラレータ計算領域ディレクティブの中に備わるデータ転送用クローズを使用して、配列のデータ転送定義しました。この方法とは別に、「アクセラレータデータ領域」ディレクティブで転送命令を記述し、このデータ領域の中に、「アクセラレータ計算領域」を定義する方法を STEP 4 に説明しました。STEP 4 の変形として、STEP 3 で述べた、シーケンシャル実行を行わなければならない反復ループの前に、!$acc regionを指定する場合の例も示しておきましょう。STEP 4の形式でも、STEP 5 の形式でもコンパイラは正しく、並列領域を解釈します。この STEP 3 の性能は、ほぼ STEP 3 と同じです。  

    ! 「アクセラレータデータ領域」の指示行で配列データの転送を指示。
    ! 反復ループの前でGPUへのデータの転送が行われる。
(  299) !$acc data region &
(  300) !$acc copyin (a(1:mimax,1:mjmax,1:mkmax,1:4)) &
(  301) !$acc copyin (b(1:mimax,1:mjmax,1:mkmax,1:3), c(1:mimax,1:mjmax,1:mkmax,1:3))  &
(  302) !$acc copy   (p(1:mimax,1:mjmax,1:mkmax))  &
(  303) !$acc copyin (wrk1(1:mimax,1:mjmax,1:mkmax), bnd(1:mimax,1:mjmax,1:mkmax)) &
(  304) !$acc local  (wrk2(1:mimax,1:mjmax,1:mkmax))
(  305)
    ! ここが「アクセラレータ計算領域」の入口で、その指定を行う。
    ホストから GPU へのスカラ変数等の転送の記述のみ行う。
(  306) !$acc region copyin(omega),copy(gosa),local(ss,s0)
(  307) !$acc do host
(  308)  do loop=1,nn ! 反復ループ(シーケンシャル実行)
(  309)    gosa = 0.0
(  310)
(  311)      do k=2,kmax-1   ! GPU並列領域
(  312)         do j=2,jmax-1
(  313)            do i=2,imax-1
(  314)               s0=a(I,J,K,1)*p(I+1,J,K) &
(  315)                    +a(I,J,K,2)*p(I,J+1,K) &
(  316)                    +a(I,J,K,3)*p(I,J,K+1) &
(  317)                    +b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K) &
(  318)                                -p(I-1,J+1,K)+p(I-1,J-1,K)) &
(  319)                    +b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1) &
(  320)                                -p(I,J+1,K-1)+p(I,J-1,K-1)) &
(  321)                    +b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1) &
(  322)                                -p(I+1,J,K-1)+p(I-1,J,K-1)) &
(  323)                    +c(I,J,K,1)*p(I-1,J,K) &
(  324)                    +c(I,J,K,2)*p(I,J-1,K) &
(  325)                    +c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
(  326)               ss=(s0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
(  327)               gosa= gosa + ss*ss
(  328)               wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
(  329)            enddo
(  330)         enddo
(  331)      enddo
(  332)
(  333) !     p(2:imax-1,2:jmax-1,2:kmax-1)= &
(  334) !         wrk2(2:imax-1,2:jmax-1,2:kmax-1)
(  335)
(  336)      do k=2,kmax-1
(  337)         do j=2,jmax-1
(  338)           do i=2,imax-1
(  339)               p(I,J,K)=wrk2(I,J,K)
(  340)            enddo
(  341)         enddo
(  342)      enddo
(  343)
(  344)  enddo
(  345) !$acc end region ! 「アクセラレータ計算領域」の終了ポイント
(  346) !$acc end data region !「アクセラレータデータ領域」の出口
    ! end data regionを指定したこのポイントで、GPUからホストへのデータの転送が行われる。

姫野ベンチマークの各作業ステップの性能(その1)

 ここまでの各ステップの性能を整理しておきましょう。ここまで行った作業は、単なる「アクセラレータ計算領域」の設定と「アクセラレータデータ領域」を指定した場合の、単純にGPUで計算すべきプログラム領域を指定しただけの作業です。GPUの CUDA 用の kernel のチューニングは一切していません。次のコラムでこの続編を書きますが、そこでは kernel のチューニング方法について説明します。PGIアクセラレータプログラミングモデルによる最初の最適化作業フェーズにおける性能は以下の通りです。GTX480(Fermi) と GTX 285 の性能を示しました。なお、以下の表は PGI 10.6 バージョンを使用した際の性能です。

表-2 PGI 10.6 における 姫野ベンチマーク(Fortran90) の GPU 性能(その1)
作業(最適化)ステップ GTX 480 GTX 285
STEP 1 反復ループの内側にアクセラレータ計算領域ディレクティブ 1.29 GFLOPS 1.27 GFLOPS
STEP 2 反復ループの外側にアクセラレータ計算領域ディレクティブ 24.20 GFLOPS 18.29 GFLOPS
STEP 3 STEP 2+データ転送Clauseを追加指定 24.41 GFLOPS 23.62 GFLOPS
STEP 4 反復ループの内側にアクセラレータ計算領域を置き、別途、データ領域ディレクティブ使用 24.33 GFLOPS 23.54 GFLOPS
STEP 5 STEP 3+データ領域ディレクティブを別途指定 24.40 GFLOPS 23.62 GFLOPS

 次回へ続きます。
 PGIアクセラレータによる姫野ベンチマークの最適化手順 コラム(全2回)