PGIアクセラレータによるヤコビ反復法のGPU性能高速化
GPGPU、CUDA、ヤコビ反復法、PGIアクセラレータ、Fermi 倍精度性能
前回は、PGIアクセラレータのディレクティブの使い方や Kernel ループ・スケジューリングの方法などを説明し、姫野ベンチマークを題材にその性能最適化手順を説明しました。本稿では、PGIアクセラレータを使って、一般的な二次元の問題をヤコビ反復法で解く「アプリケーション」を使用して、実際に GPU 用に最適化してゆくステップを説明します。具体的には、簡単にプロファイリングを行い、CPU時間のホットスポットを特定します。次に、そのホットスポットの並列化対象部分に「データの依存性」がないかを確認する方法を説明します。そして、PGIアクセラレータのディレクティブを入れて性能を最適化します。この例題は、倍精度プログラムであるため、Fermi（GTX480) における実際の倍精度性能が分かります。
2010年7月23日 Copyright © 株式会社ソフテック 加藤
ヘルムホルツ式を差分で離散化したものをヤコビ反復法で解く倍精度プログラムへの適用
本稿で使用するプログラムは、ヘルムホルツ式を差分で離散化したものをヤコビ反復法で解くものです。普通のユーザが、特殊なソースレベル最適化を施さないでコーディングした一般的なプログラムです。本プログラムは、以前のコラム「マルチコアCPU上の並列化手法、その並列性能と問題点」の中の性能実験で使用したものと同じです。同じプログラムを使用することで一般的なインテル/AMD の x64 プロセッサの性能（シングルコア性能、マルチスレッド性能）と比較し、CPUとGPUの性能の違いを理解する一助にもなります。（但し、今回の計算は、以前のコラムで説明した解析サイズの半分で計算しています。）
【参考コラム】
「マルチコアCPU上の並列化手法、その並列性能と問題点」
現在のマルチコア・プロセッサ(IntelやAMDのMPU）に焦点を当て、そのプロセッサ上で実現できている並列化手法とその並列性能に関する知見、現在のプロセッサ技術の問題点等を簡単に整理
Performance Summary
最初に、今回 PGIアクセラレータコンパイラを用いて得られた性能を示します。
| 作業（最適化）ステップ | GTX 480 (Fermi) | |
|---|---|---|
| STEP 1 | 単純に、ループをアクセラレータ計算領域ディレクティブで囲む | 18.30 秒 |
| STEP 2 | ホスト～GPU間のデータ転送の最小化 | 1.78 秒 |
| STEP 3 | STEP 2＋ループスケジューリングの変更 | 1.75 秒 |
| STEP 4 | STEP 3 + その他のチューニング（ループ内の除算の排除） | 1.37 秒 |
| Nehalemホスト側の1コアでの性能 | 11.64 秒 |
|---|
性能実測するシステム環境について
まず始めに、性能実測に使用する私のテストマシン環境について説明します。
| ホスト側プロセッサ | 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 |
（注意）GTX480 の倍精度性能スペックは、同じ Fermi アーキテクチャのTesla C2050 のものとは異なります。
使用するプログラムと本稿での最適化適用ステップの説明
ここで使用するプログラムは、以前のコラム「マルチコアCPU上の並列化手法、その並列性能と問題点」の中の性能実験で使用したものと同じです。プログラムの概要は、そのページをご参照下さい。今回は、「倍精度性能評価」を目的としますので、倍精度版 jacobi.f をベースとして使用しますが、GPUの搭載メモリ制約により、解析する系の大きさを半分に落としています（以前のコラムで呈示した計算時間の半分です）。ここで使用した各ステップのソースファイルと Makefile を纏めたアーカイブを提供しますので、お使い下さい。
本稿のPGIアクセラレータの適用ステップは次の通りです。
- STEP1：単純に 並列の対象ループ・ブロックを !$acc region /end region で囲む、悪い例
- STEP2：ホスト～GPU間のデータ転送の最小化 + ループスケジューリング調整
- STEP3： ループスケジューリングの変更
- STEP4：その他のチューニング（ループ内の除算の排除）
● 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 の方法 （例）jacobi1.F をコンパイルし、jacobi1.exe と言う実行モジュールを作成する場合 ソースファイル名は jacobi*.F と言う形式を提供している。* は、数字。 実行モジュール名は jacobi*.exe と言う名前としている。* は、数字。 コンパイルする際は、make の引数にビルドしたい実行モジュール名を指定する。 $ make jacobi1.exe
※ PGI Windows版の makefile の利用に関して
プログラムの実行条件
本ベンチマークを実行する上で、ソースプログラムの変更は以下の通りです。
- 実行に使用する計算サイズは、n = 5120、m = 5000 とする。
- 時間計測では、ヤコビ反復は収束させるまで計算せず、100回反復の固定で計算時間を評価する。
- 実行は、NVIDIA Geforce GTX480(Fermi)のみを使用して行う。
ホスト側のCPUの 1コアの性能
まずホスト側の実行時間を計測しておきます。make ユーティリティで、最大の最適化オプションでコンパイルします。その後、実行し Nehalem 1コアでの実行時間は、11.6 秒となりました。
● make でコンパイル [kato@photon29 DOUBLE]$ make host.exe pgfortran -o host.exe jacobi4.F -fastsse -Minfo 最終ステップ 4 のソースを使用する initialize: 150, Invariant if transformation 152, Invariant assignments hoisted out of loop Loop not vectorized: mixed data types Unrolled inner loop 4 times jacobi: 224, Memory copy idiom, loop replaced by call to __c_mcopy8 232, Generated 4 alternate loops for the loop Generated vector sse code for the loop Generated 4 prefetch instructions for the loop error_check: 275, Invariant if transformation 277, Invariant assignments hoisted out of loop Loop not vectorized: mixed data types Unrolled inner loop 4 times Generated a prefetch instruction for the loop ● 実行 [kato@photon29 DOUBLE]$ host.exe Input n,m - grid real*8 in x,y direction N= 5120 M= 5000 Input alpha - Helmholts constant Alpha= 1.000000000000000 Input relax - Successive over-relaxation parameter relax= 0.5000000000000000 Input tol - error tolerance for iterative solver error tolerance= 1.0000000000000000E-013 Input mits - Maximum iterations for solver Max iterations= 100 Time measurement accuracy : .10000E-05 Total Number of Iterations 101 Residual 3.8512793897632485E-011 Solution Error : 1.0538681005932186E-004 Elpased Time (Initialize + Jacobi solver + Chack) : 11.641 FORTRAN STOP
| Nehalemホスト側の1コアでの性能 | 11.64 秒 |
|---|
ホスト側のCPUでプロファイル情報を得る
まず、アプリケーションの実行時間の中で、どのルーチンがどの程度の時間を消費しているかを判断します。この中で消費時間が多いものから、その計算部分のGPUへのオフロード化を進めます。なお、このプログラム題材に関しては、別のコラム「プログラムのプロファイルを取得する（時間の掛かる部分を特定する）」で、プロファイル情報の取得方法を説明しております。このときは、プロファイルのサンプリングを別の方法を使用した場合の説明でした。今回のプロファイルのサンプリングは、「実時間計測ベース」のプロファイル方法で情報を取得します。これで、時間の掛かっているルーチンを実時間ベースで特定します。以下に示すように、今回は、三つのルーチン（jacobi, initialize, error_chack）の GPU 並列化を行えるかと言う方針で臨みます。但し、この中で、jacobiルーチンが最も時間を使っているルーチンであることは、念頭に入れておく必要があります。このルーチンを速くすれば、大きな性能向上をもたらすと言うことでもあります。
● 関数分布表示レベルのプロファイル用にコンパイル [kato@photon29 DOUBLE]$ pgfortran -fastsse -Mprof=func -o jacobi1.exe jacobi1.F ● 実行して、プロファイル情報ファイル pgprof.out を生成 [kato@photon29 DOUBLE]$ jacobi1.exe Input n,m - grid real*8 in x,y direction N= 5120 M= 5000 （中略） .... Elpased Time (Initialize + Jacobi solver + Chack) : 11.641 FORTRAN STOP [kato@photon29 DOUBLE]$ ls -l pgprof.out （pgprof.out ファイルが生成された） -rw-r--r-- 1 kato softek 382 Jul 15 15:13 pgprof.out ● PGPROFユーティリティでプロファイル結果を表示 [kato@photon29 DOUBLE]$ pgprof -text pgprof.out GUI を使わず、コマンドベースで PGPROF を使用する pgprof 10.6-0 Copyright 1989-2000, The Portland Group, Inc. All Rights Reserved. Copyright 2000-2010, STMicroelectronics, Inc. All Rights Reserved. Datafile : pgprof.out Processes : 1 Threads : 1 pgprof> p (サマリー情報をプリント、ルーチン毎に消費時間を％表示で） Profile output - Thu Jul 15 15:13:14 JST 2010 Program : /home/kato/GPGPU/OpenMP/DOUBLE/a.out Datafile : pgprof.out Process : 0 Total Time for Process : 11.641867 secs Sort by max time Select all Routine Source Line Calls Time(%) Name File No. 1 94 jacobi jacobi1.F 162 一番時間が掛かっている 1 4 initialize jacobi1.F 127 1 2 error_check jacobi1.F 243 1 0 main jacobi1.F 1 1 0 driver jacobi1.F 69 2 0 second jacobi1.F 277 pgprof> t raw （絶対時間を単位とした表示に切り替え） pgprof> p Routine Source Line Calls Time(s) Name File No. 1 10.9721 jacobi jacobi1.F 162 10秒を縮めることができるか？ 1 0.478491 initialize jacobi1.F 127 1 0.19044 error_check jacobi1.F 243 1 0.000738 main jacobi1.F 1 1 0.000071 driver jacobi1.F 69 2 0.000028 second jacobi1.F 277
GPU で計算を行わせたいループが並列化可能であるか？
そもそも、GPUで処理をさせたいループ領域は、その本質として「並列化」が可能なループでなければなりません。さらに、スレッド・ブロックの対象となるループ範囲では、「ベクトル化」可能でなければなりません。ここで、並列化とかベクトル化とか言っても、ピンと来ないと思います（この概念の説明は、こちらで説明しています。）。厳密論で言えば、並列化とベクトル化は異なる概念であり、その処理体系における「データの依存性が起こりうる状況」も若干異なるのですが、依存性があるかどうかの判断においては、この二つを別物として考えなくても良いです。依存性を考える上では、「ベクトル依存性」がないループは、「並列の依存性」もない、あるいは、「並列の依存性」がないループは、「ベクトル化」もできると考えても結構です（但し、厳密には例外もありますが、細かいことは気にしない）。
「データの依存性がある」とか「ない」とか言う概念は、どう言う意味でしょうか？例えば、並列実行における依存性を考える場合は、まず、あるループ内の処理において、それぞれのインデックスに係わるループ内計算を順序バラバラに並列に処理するとした場合を考えます。そして、その並列計算中にストアされる（配列）データ値が、並列処理の順番の変化によってその値が変わり得るかどうかと言うことで判断します。順番次第で値が変化するような場合は、「並列依存性」があるとされます。この場合、並列計算は不可能です。プログラムのループ内にある、右辺と左辺に置かれている同一配列の添字関係を調べて、ループ・インデックスの順番を変えて実行した場合、どう言う風になるかと言うことを考えれば、依存性あるかないかは直ぐに分かると思います。一方、ベクトル依存性を考える場合は、次のように考えます。ベクトル処理とは、一つの演算命令で同時に実行できる複数要素（これを「ベクトル長」と言う。）の計算処理を言いますが、ベクトル単位の一塊の計算を順番に行った時に、ストアされる左辺に位置する配列データの値が変わり得るかと言うことを考えれば良いのです。ベクトル処理とシリアル処理の違いは、シリアル処理では1要素ずつ順番に添字を変化させて計算するのに対して、ベクトル処理は、ベクトルの長さの複数の要素（添字）を一度に処理することです。こうした「依存性」を考える場合は、常に、シーケンシャルの場合とベクトル・並列処理の場合とで、実行時に左辺側の配列の値が変化する可能性があるかと言うことを対比して考えれば、考え易いです。
ともかく、依存性の存在については、同一配列の各次元の右辺と左辺の添字関係を追っていけば良いのですが、まずは、コンパイラのメッセージから、ループの中の依存性を判断したいと言う場合も多いと思います。簡易的に調べるときは、コンパイラのメッセージで確認します。PGIコンパイラは、従来の x64 プロセッサに対するマルチコア並列化対応のコンパイラでもあります。自動でプログラム内の並列化可能なループを抽出し、自動並列化用の実行モジュールを作成できます。さらに、今のインテル、AMDのプロセッサの中には、 SSE と言うベクトル演算機構があり、ベクトル処理できるようになっています。PGIコンパイラは x64 プロセッサ上で動かすプログラムに対しても、「自動ベクトル化」と「自動並列化」を行う機能を持ち備えます。PGIコンパイラが持つこの分析機能を利用して、プログラムのどの部分がベクトル化、並列化できるのかを知るのも一つの方法です。「並列化、ベクトル化ができる」と言ったメッセージであれば、GPU上でも並列化実行対象となると判断できます。もし、並列化できないと言ったメッセージの場合は、並列化が可能となるようにソースプログラムを変更する必要があります。
pgfortran / pgf77 / pgcc / pgCC コマンド
-Minfo[=all|accel|ccff|ftn|hpf|inline|
intensity|ipa|loop|lre|mp|opt|par|pfo|stat|time|unified|vect]
この中で、ベクトル化と並列化に関するメッセージを出力するサブオプション
-Minfo=vec
-Minfo=par
-Minfo オプションは、コンパイラが行う各種最適化に関するメッセージを出力するためのオプションです。単純に、-Minfo だけを指定すると -Minfo=all と同じ意味で、全てのメッセージを出力します。もし、ある最適化事象だけに絞ってメッセージを出したい場合は、-Minfo=vec と言ったサブフラグを指定します。複数のフラグを指定する場合は、-Minfo=vec,par,accel と言った形でカンマで区切ります。
- -Minfo=vec は、ベクトライザを有効にするオプション -fastsse を指定してコンパイルしたとき、最内側のループに対してその SSE ベクトル化の可能性を確認します。そして、このベクトライザの処理に関してのメッセージが出力されます。このメッセージを見ると、ベクトル化可能なループかどうかが分かります。必ず、-fastsse を付けてコンパイルして下さい。
-
-Minfo=par は、パラレライザを有効にするオプション -Mconcur を指定してコンパイルしたとき、プログラム内の多重ループに対して並列化が適用できるかその可能性を確認します。さらに、
-Mconcur=innermost を指定すると一番内側のループも並列化しようとします。こうしたパラレライザが行った並列化処理に関してのメッセージが出力されます。このメッセージを見ると、並列化可能なループかどうかが分かります。必ず、-Mconcur あるいは -Mconcur=innermost を付けてコンパイルして下さい。
● PGIコンパイラの「ベクトライザ」と「パラレライザ」を利用してメッセージを出す
[kato@photon29 DOUBLE]$ pgfortran -fastsse -Mconcur=innermost -Minfo=vec,par jacobi1.F
initialize:
148, Parallel code generated with block distribution if trip count is greater than or equal to 50
149, Loop not vectorized: mixed data types
jacobi:
210, Loop not vectorized/parallelized: multiple exits
214, Parallel code generated with block distribution for inner loop if trip count is greater ...
220, Parallel code generated with block distribution if trip count is greater than or equal to 50
221, Generated 4 alternate loops for the loop
Generated vector sse code for the loop
error_check:
261, Parallel code generated with block distribution if trip count is greater than or equal to 50
262, Loop not vectorized: mixed data types
実際に、プログラムのベクトル依存性あるいは並列依存性が存在するかを判断するために、実際に -Minfo=vec,par オプションを付けてコンパイルしてみた結果が、上記のメッセージ・リストとなります。この場合は、ベクトライザと自動パラレライザを明示的に起動させるために、必ず、-fastsse -Mconcur=innermost の二つのオプションも同時に指定しています。
さて、どうでしょう。以下にその対象となったプログラムの部分を抜粋しました。左端の番号は、ソース行番号です。上記のコンパイル・メッセージは、この番号に沿ってそのガイダンスを明示しています。
以下のサブルーチン initialize では、行番号148 の外側ループで並列化できること、149 行の内側ループでは、ベクトル化できないことを指摘しています。一般にコンパイラは、多重ループの存在があれば、内側ループで「SSEベクトル化」を外側ループで「並列化」を行おうとします。これがデフォルトの形です。149 行のループはなぜベクトル化できないかと言う原因についてですが、149 行ループの外形的な見え方では、何らデータの依存性は存在しません。左辺の f(i,j) と言う配列は、右辺式の中に存在しないため、添字関係によるデータのストアの干渉がないからです。実際の原因は、xx, yy と言う変数が、実は「整数宣言」されていることにあります。153 行の演算で左辺は「整数演算」であり、この結果をストアするものは、f(i,j) と言う real 配列です。インテル/AMD のプロセッサの SSE 演算機構（ベクトル演算機構）は、整数・浮動小数点の混合演算ができないという制約があるため、x86/x64 用の PGI コンパイラは、そのベクトル化を抑止しています。ここで使用しているコンパイルオプション -fastsse は、従来の x86/x64プロセッサ用のベクトライザですので、そのメッセージもそのプロセッサ特有のメッセージです。この149 行ループは、実際は、演算式におけるデータの依存性はないため、GPU上ではベクトル化可能です。従って、この計算ブロックは、性能効率は兎も角として、GPU並列化対象となります。
● サブルーチン initialize のリスト抜粋 ( 138) integer i,j, xx,yy ( 147) !$acc region ( 148) do j = 1,m ( 149) do i = 1,n ( 150) xx = -1.0 + dx * real(i-1) ! -1 < x < 1 ( 151) yy = -1.0 + dy * real(j-1) ! -1 < y < 1 ( 152) u(i,j) = 0.0 ( 153) f(i,j) = -alpha *(1.0-xx*xx)*(1.0-yy*yy) ( 154) & - 2.0*(1.0-xx*xx)-2.0*(1.0-yy*yy) ( 155) enddo ( 156) enddo ( 157) !$acc end region initialize: 148, Parallel code generated with block distribution if trip count is greater than or equal to 50 149, Loop not vectorized: mixed data types
下記の jacobi ルーチンは、GPUによって時間を短縮できる最大の候補です。行番号 210 の do while ～ endo ループで、収束判定ループを形成しています。このループ内に実際のヤコビ式の本体ループ本体が置かれています。ここでは、コンパイラのメッセージに頼ることなく、まずアルゴリズムとして並列の対象を考えてみることをお勧めします。210 行の収束判定ループは、本来並列化はできません。必ずシリアル処理を行わなければなりません。従って、その内部のループ・ブロックがベクトル化・並列化できるかどうか判断します。コンパイル・メッセージを見てみましょう。
コンパイラは、確かに 210 行の do while ループはベクトル化・並列化できないと認識しています。これは、プログラムが外形的に問題があるからです。並列化対象ブロックは、そもそもアルゴリズム以前に、「構造化ブロック(Structured block)」でなければならないと言う大原則があります。この詳細は、以前のコラムに説明しています。この原因は、 do while ループの条件判断でループを抜け出す形になっていることです。すなわち、ループ内に複数の exit ポイントがあることによるものです。行番号 214 は、配列間のデータコピーですので、依存性は存在しないため並列化可能です。さて、最も時間が掛かる部分である 220 行からの二重ループは、外側ループで並列化が可能で、221 行の内側ループではベクトル SSE コードの生成が可能であるとコンパイラは認識しています。従って、全く同じ構図で、GPU においても内側ループを SIMD ベクトル化対象として扱い、また外側ループを並列化対象として扱うことが可能です。従って、この計算ブロックは、GPU並列化対象となります。
● サブルーチン jacobi のリスト抜粋 ( 210) do while (k.le.maxit .and. error.gt. tol) ==> 収束判定ループ ( 211) error = 0.d0 ( 212) ( 213) !$acc region ( 214) do j=1,m ( 215) do i=1,n ( 216) uold(i,j) = u(i,j) ( 217) enddo ( 218) enddo ( 219) ( 220) do j = 2,m-1 ( 221) do i = 2,n-1 ( 222) resid = (ax*(uold(i-1,j) + uold(i+1,j)) ( 223) & + ay*(uold(i,j-1) + uold(i,j+1)) ( 224) & + b * uold(i,j) - f(i,j))*1.d0/b ( 225) u(i,j) = uold(i,j) - omega * resid ( 226) error = error + resid*resid ( 227) end do ( 228) enddo ( 229) !$acc end region ( 230) ( 231) * Error check ( 232) K = k + 1 ( 233) error = sqrt(error)/dble(n*m) ( 234) * ( 235) enddo ! End iteration loop jacobi: 210, Loop not vectorized/parallelized: multiple exits 214, Parallel code generated with block distribution for inner loop if trip count is greater ... 220, Parallel code generated with block distribution if trip count is greater than or equal to 50 221, Generated 4 alternate loops for the loop Generated vector sse code for the loop
下記の error_chack ルーチンは、上記で述べた initialize と似た形となっています。262 行ループがベクトル化できない理由は、整数・浮動小数点の混合演算が含まれていると言う上記と全く同じ理由です。さて、この二重ループは、そもそも並列化できるのでしょうか？メッセージを見ると、コンパイラは 261 行ループは並列化できると言っています。外形的に見て、このループ内には複数のスカラ変数とただ一つの配列データ u(i,j) が存在します。参照のみの配列 u(i,j) は存在していますが、左辺部はスカラ変数 temp であり配列ではないため、265行目までの演算式の状態では並列化は困難です。繰り返しになりますが、GPU上での並列化とは、添字要素毎に並列計算が可能となる「データ並列」が基本形ですので、外形的には左辺部に配列が置かれ、各要素並列に値がストアできることが基本です。266 行目を見て下さい。これは、 error と言うスカラ変数の総和を求めるもので、こう言った演算形式を「リダクション」と言います。「リダクション」形式のものとして、その他に、最大値、最小値を求める演算なども含みます。一般的なシミュレーション・アルゴリズムでは、こうしたリダクション形式の演算が必ず出現するために、実は、これを並列化するためのロジックを PGI コンパイラが内部的に用意しています。PGIアクセラレータでは、コンパイラがリダクション演算部を自動検出して、GPU上の並列演算をできるようにしています。また、OpenMP のようなスレッド並列化の場合もコンパイラ機能で並列化します。MPI のような並列ライブラリでは、こうしたリダクションを行うルーチンを明示的に提供しています。しかしながら、CUDA Fortranや CUDA C のようなネイティブな言語体系では、ユーザが明示的に並列リダクション・コードを記述しなければなりません。NVIDIA GPUのリダクション・アルゴリズムに関しては、CUDA SDK の並列リダクション用のリファレンス・ルーチンを参考にすれば良いのですが、実は一般のユーザにとって、これを自分のコードの中に組み込むことは非常に難しいと言うことも事実です。ともかく、この計算ブロックもGPU並列化対象となります。
● サブルーチン error_check のリスト抜粋 ( 260) !$acc region ( 261) do j = 1,m ( 262) do i = 1,n ( 263) xx = -1.0d0 + dx * dble(i-1) ( 264) yy = -1.0d0 + dy * dble(j-1) ( 265) temp = u(i,j) - (1.0-xx*xx)*(1.0-yy*yy) ( 266) error = error + temp*temp ( 267) enddo ( 268) enddo ( 269) !$acc end region error_check: 261, Parallel code generated with block distribution if trip count is greater than or equal to 50 262, Loop not vectorized: mixed data types
PGI アクセラレータでコンパイルする（-ta=nvidia オプション）
以上のように、三つのルーチンの中の（多重）ループはベクトル化ならびに並列化もできそうだと言うことが分かりました。それでは、アクセラレータ用の -ta=nvidia オプションを指定してコンパイルしてみましょう。この場合は、少なくとも、ベクトル・並列化可能となるループ（構造化ブロック）を「アクセラレータ計算領域」として定義してコンパイルする必要があります。すなわち、当該ループの前に !$acc region を、その後に !$acc end region ディレクティブを挿入します。上記で示したコンパイルリストには、予め、これらのディレクティブを挿入してある状態を示しています。
以下のリストは、提供したソースファイル jacobi1.F を -ta=nvidia -Minfo=accel オプションを付けてコンパイルした際のメッセージリストです。Make ユーティリティを使用して、コマンド make jacobi1.exe を実行すると、コンパイル・メッセージの出力と実行モジュールの作成ができます。コンパイラは、自動でベクトル化、並列化対象ループを特定します。また、ホスト側とGPU 間で必要となるデータ転送に関しても自動解析します。以下のように、三つのルーチンは、GPU並列コードが生成されたことが分かります。
[kato@photon29 DOUBLE]$ make jacobi1.exe
pgfortran -o jacobi1.exe jacobi1.F -fastsse -Minfo=accel -ta=nvidia,cuda3.0
initialize:
147, Generating copyout(f(1:n,1:m))
Generating copyout(u(1:n,1:m)) ==>GPUからホストへの配列データ転送
Generating compute capability 1.3 binary
Generating compute capability 2.0 binary
148, Loop is parallelizable
149, Loop is parallelizable （146,129行目のループは並列化可能）
Accelerator kernel generated ==> GPU実行用の Kernel を生成した。
148, !$acc do parallel, vector(16)
149, !$acc do parallel, vector(16) ==> [16x16] のスレッド・ブロックを並列に実行
CC 1.3 : 13 registers; 24 shared, 64 constant, 0 local memory bytes; 100 occupancy
CC 2.0 : 19 registers; 8 shared, 88 constant, 0 local memory bytes; 100 occupancy
jacobi:
213, Generating copyin(f(2:n-1,2:m-1)) ==>ホストからGPUへの配列データ転送
Generating copyin(u(1:n,1:m))
Generating copyout(u(2:n-1,2:m-1)) ==>GPUからホストへの配列データ転送
Generating copyout(uold(1:n,1:m))
Generating compute capability 1.3 binary
Generating compute capability 2.0 binary
214, Loop is parallelizable
215, Loop is parallelizable
Accelerator kernel generated ==> GPU実行用の Kernel を生成した。
214, !$acc do parallel, vector(16)
215, !$acc do parallel, vector(16) ==> [16x16] のスレッド・ブロックを並列に実行
CC 1.3 : 6 registers; 24 shared, 132 constant, 0 local memory bytes; 100 occupancy
CC 2.0 : 10 registers; 8 shared, 156 constant, 0 local memory bytes; 100 occupancy
220, Loop is parallelizable
221, Loop is parallelizable
Accelerator kernel generated ==> GPU実行用の Kernel を生成した。
220, !$acc do parallel, vector(16)
221, !$acc do parallel, vector(16) ==> [16x16] のスレッド・ブロックを並列に実行
Cached references to size [18x18] block of 'uold'
CC 1.3 : 31 registers; 2616 shared, 152 constant, 0 local memory bytes; 50 occupancy
CC 2.0 : 26 registers; 2600 shared, 156 constant, 0 local memory bytes; 66 occupancy
226, Sum reduction generated for error
error_check:
260, Generating copyin(u(1:n,1:m)) ==>ホストからGPUへの配列データ転送
Generating compute capability 1.3 binary
Generating compute capability 2.0 binary
261, Loop is parallelizable
262, Loop is parallelizable
Accelerator kernel generated ==> GPU実行用の Kernel を生成した。
261, !$acc do parallel, vector(16)
262, !$acc do parallel, vector(16) ==> [16x16] のスレッド・ブロックを並列に実行
CC 1.3 : 14 registers; 2072 shared, 72 constant, 0 local memory bytes; 100 occupancy
CC 2.0 : 17 registers; 2056 shared, 84 constant, 0 local memory bytes; 100 occupancy
266, Sum reduction generated for error
ステップ 1 （単純に、ループをアクセラレータ計算領域ディレクティブで囲む）
jacobi1.F のソースは、単純に並列対象ループを「アクセラレータ計算領域ディレクティブ」で囲んだものです。jacobi1.F を実行した結果を以下に示します。見ての通り、ホスト側 1 コアを使用して実行した結果 11.6 秒よりも時間が掛かっています。この理由は、一番時間を消費する jacobi ルーチンの中で、「収束判定を行う do while ループ」の中で、「アクセラレータ計算領域」が指定されているからです。すなわち、アクセラレータ計算領域ディレクティブ（!$acc region) の指定場所では ホスト～GPU間のデータ転送を行うため、その収束判定ループの繰り返し数分、GPUへのデータ転送が毎回行われています。以下の「Accelerator Kernel Timing data」にも、data=14752538 が記録されており、14.7秒の時間がデータ転送だけに費やされています。こうしたことを回避するには、「収束判定を行う do while ループ」の外側で1 回のみ、ホスト～GPU間のデータ転送を行うようにすることです。これは、以前のコラムで説明した「アクセラレータデータ領域ディレクティブ」を活用すれば実現できます。これをステップ 2 で実施します。
●jacobi1.exe 実行結果 [kato@photon29 DOUBLE]$ jacobi1.exe Input n,m - grid real*8 in x,y direction N= 5120 M= 5000 Input alpha - Helmholts constant Alpha= 1.000000000000000 Input relax - Successive over-relaxation parameter relax= 0.5000000000000000 Input tol - error tolerance for iterative solver error tolerance= 1.0000000000000000E-013 Input mits - Maximum iterations for solver Max iterations= 100 Time measurement accuracy : .10000E-05 initialize nvidia Total Number of Iterations 101 Residual 3.8512793897632957E-011 Solution Error : 1.0538681005932478E-004 Elpased Time (Initialize + Jacobi solver + Chack) : 18.303 FORTRAN STOP Accelerator Kernel Timing data /home/kato/GPGPU/OpenMP/DOUBLE/jacobi1.F error_check 260: region entered 1 time time(us): total=47701 kernels=7983 data=35522 262: kernel launched 1 times grid: [313x320] block: [16x16] time(us): total=7716 max=7716 min=7716 avg=7716 266: kernel launched 1 times grid: [1] block: [256] time(us): total=267 max=267 min=267 avg=267 /home/kato/GPGPU/OpenMP/DOUBLE/jacobi1.F jacobi 213: region entered 100 times time(us): total=18066200 init=27 region=18066173 kernels=2163193 data=14752538 ==>ホストからGPUへのデータ転送時間が 14.7 秒も w/o init: total=18066173 max=221807 min=180129 avg=180661 215: kernel launched 100 times grid: [313x320] block: [16x16] time(us): total=312758 max=3133 min=3106 avg=3127 221: kernel launched 100 times ==> 221行ループのカーネル実行時間が最大である grid: [313x320] block: [16x16] time(us): total=1823725 max=18289 min=17585 avg=18237 ==>平均 18.2 msec 226: kernel launched 100 times grid: [1] block: [256] time(us): total=26710 max=270 min=266 avg=267 /home/kato/GPGPU/OpenMP/DOUBLE/jacobi1.F initialize 147: region entered 1 time time(us): total=188541 kernels=5471 data=175107 149: kernel launched 1 times grid: [313x320] block: [16x16] time(us): total=5471 max=5471 min=5471 avg=5471 acc_init.c acc_init 41: region entered 1 time time(us): init=567958
ステップ 2 （ホスト～GPU間のデータ転送の最小化）
ステップ 2（jacobi2.F)は、「アクセラレータデータ領域ディレクティブ」を「収束判定を行う do while ループ」の外で指定したものです。具体的には、以下に示すようなディレクティブの指定を行っています。また、データの転送に関しては、copy, copyin 等のクローズを使用して、明示的に当該配列・変数を指定しています。さらに、二重ループ内の並列分割とベクトル長を指示するため、「ループマッピング・ディレクティブ」も明示的に指定しました。
( 210) !$acc data region copy(u(1:n,1:m)) copyin(f(1:n,1:m)) ( 211) !$acc+ local(uold(1:n,1:m)) local(resid) ( 212) ==> 収束判定ループの前にデータを GPU へ明示的に転送する ( 213) error = 10.0d0 * tol ( 214) k = 1 ( 215) do while (k.le.maxit .and. error.gt. tol) ==> 収束判定ループ ( 216) error = 0.d0 ( 217) 以下のアクセラレータ計算領域内は、GPUのデバイスメモリ内のデータを使用する ( 218) !$acc region copyin( dx, dy,omega,b ), copy(error) ( 219) !$acc do parallel ( 220) do j=1,m ( 221) !$acc do vector(256) ( 222) do i=1,n ( 223) uold(i,j) = u(i,j) ( 224) enddo ( 225) enddo ( 226) ( 227) !$acc do parallel ( 228) do j = 2,m-1 ( 229) !$acc do vector(256) この内側ループはベクトル・モードだけで実行する ( 230) do i = 2,n-1 ( 231) resid = (ax*(uold(i-1,j) + uold(i+1,j)) ( 232) & + ay*(uold(i,j-1) + uold(i,j+1)) ( 233) & + b * uold(i,j) - f(i,j))*1.d0/b ( 234) u(i,j) = uold(i,j) - omega * resid ( 235) error = error + resid*resid ( 236) end do ( 237) enddo ( 238) !$acc end region ( 239) ( 240) * Error check ( 241) K = k + 1 ( 242) error = sqrt(error)/dble(n*m) ( 243) * ( 244) enddo ! End iteration loop ( 245) !$acc end data region ==> !$acc end data region の指定ポイントで GPU からホストへデータを戻す ●コンパイル・メッセージ jacobi: 210, Generating local(resid) ==> 210行目でデータ転送が行われることが分かる Generating local(uold(1:n,1:m)) Generating copyin(f(:n,:m)) Generating copy(u(:n,:m)) 220, Loop is parallelizable 222, Loop is parallelizable Accelerator kernel generated 220, !$acc do parallel 222, !$acc do vector(256) 228, Loop is parallelizable 230, Loop is parallelizable Accelerator kernel generated 228, !$acc do parallel 230, !$acc do vector(256) Cached references to size [258x3] block of 'uold' ==> uold を shared memory 上に 235, Sum reduction generated for error キャッシュ化した
実際に実行してみると、実行時間は、一気に 1.77 秒まで短縮されました。実行プロファイルデータを見ると、jacobi ルーチンでの 210 行ループのホスト～GPU間のデータ転送は 1 回のみで、0.1秒に過ぎません。一方、jacobi ルーチンの処理時間は、230 行ループが大勢を占めています。GPU内で kernel として 100 回起動されていますが、1 回当たりの平均時間は 11.3 msec であることが分かります。前に述べた「ステップ 1」の時のこのループでの時間は、18.2 msec でしたので、ステップ 2 の方のループスケジューリングの方が良いと言うことも分かります。今後、kernel scheduling（ループマッピング）のチューニングを行う場合は、この時間を指標として試行錯誤を行えば良いことになります。
●jacobi2.exe 実行結果 [kato@photon29 DOUBLE]$ jacobi2.exe Input n,m - grid real*8 in x,y direction N= 5120 M= 5000 Input alpha - Helmholts constant Alpha= 1.000000000000000 Input relax - Successive over-relaxation parameter relax= 0.5000000000000000 Input tol - error tolerance for iterative solver error tolerance= 1.0000000000000000E-013 Input mits - Maximum iterations for solver Max iterations= 100 Time measurement accuracy : .10000E-05 initialize nvidia Total Number of Iterations 101 Residual 3.8512793897632873E-011 Solution Error : 1.0538681005932478E-004 Elpased Time (Initialize + Jacobi solver + Chack) : 1.777 FORTRAN STOP Accelerator Kernel Timing data (GPUの実行プロファイル情報) /home/kato/GPGPU/OpenMP/DOUBLE/jacobi2.F error_check 270: region entered 1 time time(us): total=42066 kernels=2554 data=35484 274: kernel launched 1 times grid: [5000] block: [256] time(us): total=2534 max=2534 min=2534 avg=2534 278: kernel launched 1 times grid: [1] block: [256] time(us): total=20 max=20 min=20 avg=20 /home/kato/GPGPU/OpenMP/DOUBLE/jacobi2.F jacobi 218: region entered 100 times time(us): total=1427184 init=2 region=1427182 kernels=1419469 data=0 w/o init: total=1427182 max=14492 min=14216 avg=14271 222: kernel launched 100 times grid: [5000] block: [256] time(us): total=285850 max=2863 min=2857 avg=2858 230: kernel launched 100 times ==> 230行ループのカーネル実行時間が最大である grid: [4998] block: [256] ==>ブロックサイズとグリッドサイズ time(us): total=1131619 max=11372 min=11271 avg=11316 ==>平均 11.3 msec 235: kernel launched 100 times grid: [1] block: [256] time(us): total=2000 max=20 min=20 avg=20 /home/kato/GPGPU/OpenMP/DOUBLE/jacobi2.F jacobi 210: region entered 1 time ==> データ転送は 1 回のみで 0.1 秒 time(us): total=1548195 data=109480 /home/kato/GPGPU/OpenMP/DOUBLE/jacobi2.F initialize 147: region entered 1 time time(us): total=186724 kernels=4326 data=174435 151: kernel launched 1 times grid: [5000] block: [256] time(us): total=4326 max=4326 min=4326 avg=4326 acc_init.c acc_init 41: region entered 1 time time(us): init=574953
ステップ 3 （ループスケジューリングを少し調整してみる）
PGIアクセラレータを使用して、性能をさらにチューニングしようとした場合、最後は、PGIアクセラレータの「ループスケジューリング・クローズ」の調整と言うことになります。以下にステップ 2 と 3 の「ループスケジューリング・クローズ」の違いを示しました。この部分は、最も計算時間の掛かる jacobi サブルーチンの 227 行目の2重ループです。この二重ループの「ループスケジューリング」を変更してみました。ステップ 3 の方は、i と j ループのどちらも、ベクトル長（SIMD長）16 で分割をし、16 x 16 の 2 次元のスレッド・ブロックの生成を指示しています。さらに j ループは、全体を 16 分割する、大きな単位のグリッド構成を行うように指示しています。実行時間は、ステップ 2 に較べ若干、速くなっています。
【ステップ 2 jacobiルーチン227行目以降】 【ステップ 3 jacobiルーチン227行目以降】 !$acc do parallel !$acc do parallel(16) vector(16) do j = 2,m-1 do j = 2,m-1 !$acc do vector(256) !$acc do parallel vector(16) do i = 2,n-1 do i = 2,n-1 resid = (ax*(uold(i-1,j) + uold(i+1,j)) +-----------------+ & + ay*(uold(i,j-1) + uold(i,j+1)) | | & + b * uold(i,j) - f(i,j))*1.d0/b | （左と同じ） | u(i,j) = uold(i,j) - omega * resid | | error = error + resid*resid +-----------------+ end do end do enddo enddo !$acc end region !$acc end region
[kato@photon29 DOUBLE]$ jacobi3.exe Input n,m - grid real*8 in x,y direction N= 5120 M= 5000 Input alpha - Helmholts constant Alpha= 1.000000000000000 Input relax - Successive over-relaxation parameter relax= 0.5000000000000000 Input tol - error tolerance for iterative solver error tolerance= 1.0000000000000000E-013 Input mits - Maximum iterations for solver Max iterations= 100 Time measurement accuracy : .10000E-05 initialize nvidia Total Number of Iterations 101 Residual 3.8512793897632873E-011 Solution Error : 1.0538681005932478E-004 Elpased Time (Initialize + Jacobi solver + Chack) : 1.753 FORTRAN STOP Accelerator Kernel Timing data /home/kato/GPGPU/OpenMP/DOUBLE/jacobi3.F error_check 270: region entered 1 time time(us): total=42111 kernels=2529 data=35549 274: kernel launched 1 times grid: [5000] block: [256] time(us): total=2509 max=2509 min=2509 avg=2509 278: kernel launched 1 times grid: [1] block: [256] time(us): total=20 max=20 min=20 avg=20 /home/kato/GPGPU/OpenMP/DOUBLE/jacobi3.F jacobi 218: region entered 100 times time(us): total=1408566 init=12 region=1408554 kernels=1400869 data=0 w/o init: total=1408554 max=14229 min=14008 avg=14085 222: kernel launched 100 times grid: [5000] block: [256] time(us): total=286170 max=2866 min=2858 avg=2861 230: kernel launched 100 times grid: [16x320] block: [16x16] time(us): total=1112699 max=11262 min=11051 avg=11126 ==>平均 11.1 msec 235: kernel launched 100 times grid: [1] block: [256] time(us): total=2000 max=20 min=20 avg=20 (以下、省略）
ステップ 4 （その他のチューニング：ループ内の除算の排除）
その他、ソースレベルでチューニングすべきところがあるか調べます。このプログラムで最も時間が掛かるループ（230行目）中に、1.d0/b と言う除算演算が含まれます。除算を行うためのコストは、従来の x64 系のプロセッサでも GPU のハードウェアにおいても非常に重いため、できるだけループ内では使用しないことが肝要です。今回の場合は、ループ内で定数である b の逆数を得るために除算を行っていますので、これは、ループの外で逆数の値を計算させて、これをループ内で乗算で使うように修正することができます。
【jacobiルーチン230行目以降】 ( 229) !$acc do vector(16) ( 230) do i = 2,n-1 ( 231) resid = (ax*(uold(i-1,j) + uold(i+1,j)) ( 232) & + ay*(uold(i,j-1) + uold(i,j+1)) ( 233) & + b * uold(i,j) - f(i,j))*1.d0/b ( 234) u(i,j) = uold(i,j) - omega * resid ( 235) error = error + resid*resid ( 236) end do 【jacobiルーチン変更部】 ( 207) b = -2.0/(dx*dx)-2.0/(dy*dy) - alpha ! Central coeff ( 208) b1b= 1.d0/b (中略） ( 229) !$acc do parallel vector(16) ( 230) do i = 2,n-1 ( 231) resid = (ax*(uold(i-1,j) + uold(i+1,j)) ( 232) & + ay*(uold(i,j-1) + uold(i,j+1)) ( 233) & + b * uold(i,j) - f(i,j))*b1b ( 234) u(i,j) = uold(i,j) - omega * resid ( 235) error = error + resid*resid ( 236) end do ( 237) enddo ( 238) !$acc end region
実行結果は、以下のようになります。1.374 秒と大幅に性能が向上しました。jacobi ルーチンの 230 行ループの 1 回当たりの平均時間は 7.33 msec となり、GPU 内での除算コストはかなり大きいことが分かります。
[kato@photon29 DOUBLE]$ jacobi4.exe Input n,m - grid real*8 in x,y direction N= 5120 M= 5000 Input alpha - Helmholts constant Alpha= 1.000000000000000 Input relax - Successive over-relaxation parameter relax= 0.5000000000000000 Input tol - error tolerance for iterative solver error tolerance= 1.0000000000000000E-013 Input mits - Maximum iterations for solver Max iterations= 100 Time measurement accuracy : .10000E-05 initialize nvidia Total Number of Iterations 101 Residual 3.8512793897632866E-011 Solution Error : 1.0538681005932478E-004 Elpased Time (Initialize + Jacobi solver + Chack) : 1.374 FORTRAN STOP Accelerator Kernel Timing data /home/kato/GPGPU/OpenMP/DOUBLE/jacobi4.F error_check 270: region entered 1 time time(us): total=42074 kernels=2529 data=35521 274: kernel launched 1 times grid: [5000] block: [256] time(us): total=2509 max=2509 min=2509 avg=2509 278: kernel launched 1 times grid: [1] block: [256] time(us): total=20 max=20 min=20 avg=20 /home/kato/GPGPU/OpenMP/DOUBLE/jacobi4.F jacobi 218: region entered 100 times time(us): total=1029078 init=10 region=1029068 kernels=1021441 data=0 w/o init: total=1029068 max=10912 min=10171 avg=10290 222: kernel launched 100 times grid: [5000] block: [256] time(us): total=286177 max=2864 min=2858 avg=2861 230: kernel launched 100 times grid: [16x320] block: [16x16] time(us): total=733264 max=7799 min=7217 avg=7332 ==> 平均 7.33 msec 235: kernel launched 100 times grid: [1] block: [256] time(us): total=2000 max=20 min=20 avg=20 (以下、省略）
