第12回：OpenACCを使ったICCG法の高速化つづき

前回に引き続き、OpenACCを使ったICCG法の高速化手法について考えてみます。

​オリジナルのプログラムは東大情報基盤センターの講習会、「OpenMPによるマルチコア・メニィコア並列プログラミング入門」のページから入手できます。

ダウンロードしたmulticore-c.tarを解凍した後に出来る、マルチカラー並列化されたICCGソルバーのmulticore-c/L2/src/solver_ICCG_mc.c をベースに解説します。

OpenACCを使ったICCG法の高速化つづき

前回は、図1 の solve [M]z(i-1)=r(i-1) 部分のプログラムのOpenACC化を行いました。

以降のループも引き続き、

1. OpenACCで並列化できるのかどうか
2. CPU-GPU間のデータ転送を如何に最小化するか

​の2点を意識しながらOpenACC化を進めましょう。

図1：ICCGソルバーのアルゴリズム
（OpenMPによるマルチコア・メニィコア並列プログラミング入門より）​

図2：図1の4-10行目のプログラムのOpenACC化

図2 の7行目から始まるループは、典型的なリダクションループです。

5,6行目の様にkernels, loop 指示文を挿入します。

次はIF文で分岐しているものの、IF文の中身は非常に素直なループです。

17, 24行目から始まるいずれのループも、ループのi番目が配列のi番目を更新するデータ独立なループですから、#pragma acc loop independentを挿入して並列化します。
（なおこの時、W[P][0:N]とW[Z][0:N]がエイリアス、つまり重複する範囲を持っていると、並列化できないので注意。例えば、W[P] = &W[Z][1] のようにW[P]の指すポインタを書き換えると、W[P][0] == W[Z][1] となってしまい、データ独立ではなくなってしまいます。このようなコードは滅多にありませんが、頭の片隅に入れておきましょう。）

​無事ループの並列化ができたら、次は如何にデータ転送を最小化するかです。

図3：データ指示文の範囲変更

前回、solve [M]z(i-1)=r(i-1)の部分をOpenACC化した時のように、data指示文の適用範囲を徐々に広げていき(図3)、最終的には図1のfor文の外側に追いやります。

図3 のようにData 指示文で囲んでしまえば、図2の5, 15, 22行目においてW[R], W[Z], W[P]の各配列はGPUのメモリ上に確保済みであるとみなされ、copy, copyin, copyout指示子は無視されます。

なお前回説明し忘れましたが、図3のデータ指示文中に出てくる変数、ALMAX, AUMAX, itemLMAX, itemUMAXは、オリジナルのプログラムには存在しません。

データ指示文を使って初めて配列のサイズ情報が必要となるからです。そのため、関数のプロトタイプ宣言などを、例えば図4, 5 のように変更しなくてはなりません。

配列のサイズ情報を組み込み関数で確認できるFortranと違い、C言語を使う際にはこの手の変更がしばしば必要となります。​その過程で、図5のようにOpenACCと関係のないところでバグが顕在化することもあります。

​必ず結果を確認しながら、少しずつOpenACC化を進めていくことが、結果として早道になるのです。​

図4：solver_ICCG_mc.c, solver_ICCG_mc.hの変更（赤字が変更部分）。

図5：main.cの変更（赤字が変更部分）。配列のサイズ情報を付け加える過程でバグが顕在化する。

図6：図1の11行目、疎行列・ベクトル積のOpenACC化

引き続きOpenACC化を進めます。

​図6 はq(i)= [A]p(i)の疎行列・ベクトル積部分です。
まず11行目のループが並列化可能かどうかを考えます。

ループのi番目がそれぞれ配列W[Q]に書き込んでおり、読み込まれる配列にW[Q]はないため、データ依存のない独立なループと言えます。
よって10行目のように#pragma acc loop independentを挿入します。

13, 16行目のループについても考えてみると、変数VALに足し込みが行われている形なので、典型的なリダクションループであり並列化可能です。よって13, 17行目のように#pragma acc loop reductionを挿入します。

ただ13, 16行目のループを本当に並列化すべきかどうかは考える必要があります。

図7：データ指示文の範囲変更

本記事中ではあまり詳しく説明していませんが、

GPUは32本のスレッドを一塊として動かすため、ループ長が32以下のループを並列化した場合、何もしないスレッドが存在することとなり、むしろ遅くなってしまうことがあるのです。

13, 16行目のループ長は実は3～6（疎行列の非ゼロ要素数に相当）であるため、#pragma acc loop seqとした方が高性能であると考えられます。
​後で比較してみましょう。
​
データ転送の範囲についても図7 のように変更します。

実はW[Z]とW[Q]は同じ配列を指しているので、冗長な書き方をしています。
​（無視されるため2回転送されるなどの悪影響はありません。）

図8：図1の12-14行目のプログラムのOpenACC化

最後に残りの部分を並列化しましょう。

​図8 の7, 24行目のループは典型的なリダクションループですので、#pragma acc loop reductionを挿入します。

17行目のループは典型的なデータ独立なループですので、#pragma acc loop independentを挿入します。
これで晴れて、全てのループを並列化することができました。
よっていよいよ、データ指示文の範囲を図1 の1行目のループの外側に追い出すことができます。

しかしこのICCGソルバーの実装では、図9 のようにループを抜けるためにgoto文を利用しています。

図9：データ転送の最小化

このような実装の場合、data指示文は利用できない（data指示文の終わりの“}”部分を飛び越えてしまう可能性がある）仕様なので、代わりにenter data 指示文とexit data指示文を使います。

Enter data とexit data のペアでdata指示文と同等の効果があり、enter data はmalloc + copyin, exit data はcopyout + freeとして使うことができます。

見やすさの点でdata指示文をお勧めしていますが、実際のアプリケーションで使う際にはenter/exit dataの方が便利なことが多いです。

これでようやく、最低限のOpenACC実装が完成しました。
CPUと比較してみましょう。

PGIコンパイラでコンパイルするために、Makefileは CC = pgcc, OPTFLAGS = -O3 –acc –ta=tesla,cc60 –Minfo=accel と書き換えています。Makeするとmulticore-c/L2/solver/runに実行ファイルができます。

INPUT.datの一行目が問題サイズを表すので128, 128, 128とし、./L2-solを実行します。

​最初に並列化を行うにあたっての色数とカラーリング手法を聞かれるため、ここでは-10（CMRCMによりカラーリングし、色数は10の意味）を入力します。結果は以下でした。
 
・OpenACC版のGPU実行：　　　 0.962792秒 (CPU-GPUデータ転送を除くと0.749099)
・オリジナルのCPU1コア実行：　24.847288秒

図10: gang, vector による、スレッド割り当ての指定

また、不完全コレスキー分解を行う図1の3行目部分と、疎行列ベクトル積を行う図6 の部分の内側ループのリダクションはコンパイラに無視されていました。

並列化可能なループが複数ある場合、どのループにどのようにスレッドを割り当てるのかは、書いてなければコンパイラ任せになりますから、これはこれで正しい挙動ではあります。

リダクションしない方が効率良いと判断されたのでしょう。賢いですね。

なお内側のループが無視されないようにするためには、図6 から図10 のように書き換えます。
​​実行結果は以下でした。

​・OpenACC版のGPU実行（リダクションする場合）：22.823339秒 (CPU-GPUデータ転送を除くと22.593328)
​すさまじく遅くなりましたね！

今回で初級編は終了となりますが、いかがでしたでしょうか。

プログラムをOpenACC化する際に考えることというのは、どんなプログラムでも大体変わらず、

1. OpenACCで並列化できるのかどうか
2. CPU-GPU間のデータ転送を如何に最小化するか

という2点が基本です。

しかし今回の例のように、データサイズを指示文に与えるためにプログラムの変更が必要だったり、goto文があるために基本であるdata指示文が使えなかったりと、対象とするプログラムによって、新たに対処しなくてはならないことが色々出てきます。

次回以降は中級編ということで、さらに実践的なアプリケーションのOpenACC化を行いつつ、様々なケースでのOpenACCでの対処方法などを紹介出来ればと思います。

ここまで読んでいただきありがとうございました。
今回作ったOpenACC版の実装はこちらから入手できます。

---

＜ １ヵ月間有効のスパコンお試しアカウント ＞
東京大学情報基盤センターでは、教育の一環として、制限はあるものの一ヵ月の間有効なスパコンアカウントを提供しています。
現在3つのスパコンが運用されていますが、そのうちReedbushと呼ばれるスパコンには、一世代前のものではありますがGPUが搭載されていて、OpenACCを使える環境も整っています。
自分でどんどん自習したい場合は、ご利用を考えてみてください。

トライアルアカウント申し込みページ
https://www.cc.u-tokyo.ac.jp/guide/trial/free_trial.php

＜ 過去の講習会の資料やプログラム公開中 ＞
東大センターが行った過去のOpenACCに関する講習会の資料やプログラムも公開されていますので、自習する場合にはぜひご利用ください。

講習会ページ　※オンライン講習会 定期開催中！
https://www.cc.u-tokyo.ac.jp/events/lectures/

講習会で用いているプログラム
https://www.dropbox.com/s/z4fmc4ibdggdi0y/openacc_samples.tar.gz?dl=0​