第9回：速くならない？とりあえずライブラリに頼ろう！
gdep前回に引き続き、OpenACCでうまく行かないケースとその解決策について学びます。
特に計算が複雑になればなるほど、思ったほどの性能が出ないケースが増えてきます。
そんな時はどうしたら良いのでしょうか？
速くならない？とりあえずライブラリに頼ろう！
OpenACCを使ってみたはいいものの、思ったほど計算が速くならない…そういうケースは多々あります。
代表的な例は行列積演算です。行列と行列の積を求める計算（は行列）は、シミュレーションなどでも頻出する計算です。
OpenACCを使えばとても簡単にGPUで実行できます（図１）。
行列サイズを1024×1024として、東大のReedbush-HスパコンのP100 GPUで計測してみましょう。
P100の理論性能は5300 GFlops（1秒間に5.3兆回の倍精度浮動小数点演算が可能）ですが、この行列積の結果は…
305 GFlops！
スペックの十分の一も出てませんね！
この結果については、別にOpenACCが悪いわけではありません。
実は速い行列積プログラムを実装するというのは大変で、CPUでもGPUでもかなりの最適化が必要です。仮にCUDAを使ったとしても、図１と同程度に単純な実装であれば、OpenACCと似たような性能になるでしょう。
ではOpenACCだと何が大変なのか。
それは最適化です。
この行列積を最適化して速くしようと思っても、OpenACCの制約の中ではここからほとんど何もやりようがありません！
これはGPUのshared memoryと呼ばれる高速なメモリをOpenACCから自由に使うことができないためです（一応、cache指示文というものを使えばshared memoryを利用すること自体はできますが、スレッド間の同期がとれないので十全には使えません）。
であればどうすればよいのでしょうか。ライブラリを呼びましょう！
OpenACCからも、CUDA向けに作られたcublasなどのライブラリを利用できるのです！
倍精度の行列-行列積なら、cublasDgemmという関数（図2）が使えます！
cublasDgemmは、(Cはm*n, Aはm*k, Bはk*nの行列)を計算する関数です。
alpha = 1, beta = 0としたら、図1と同じ行列積ができます。
lda, ldb, ldcは、行列の一次元目の長さなので、この場合ならlda = m, ldb = k, ldc = mとなります。
transaとtransbは、A, Bを転置するかどうかを表す変数で、CUBLAS_OP_Nが転置しない場合、CUBLAS_OP_Tが転置する場合を表します。
Handleはcublas側で利用するもので、使い方は図3を見てください。
図3は、OpenACCでcublasDgemmを使った利用例です。
ソースコードはこのリンクからダウンロードできます。
Cublasの利用方法に従い、handleなどの変数を用意していきますが、ここで一つ困ったことが起きます。
CublasはCUDA向けに作られたBLASライブラリですから、配列A,B,CはGPU上に確保されてなければなりません。
加えて、cublasDgemmの引数として渡すのは、A,B,CのGPU上のアドレスでなければなりません！
OpenACCのdata指示文を使った場合、data指示文で囲まれた範囲で、CPU用の配列とGPU用の配列のペアが生成されるのでした。
細かい話ですがこの際、配列へのアクセスがdata指示文の内側かつkernels指示文の内側で行われたならGPU用の配列が使われ、data指示文の内側であってもkernels指示文の外側で行われたならCPU用の配列が使われます。
つまり普通にcublasDgemmに配列A,B,Cを渡しても、それはCPU上のアドレスになってしまうのです！
ここで登場するのが、#pragma acc host_data という指示文です。z
この指示文を使うと、kernelsの内側でなくてもGPU上のアドレスを得られます。Host_dataに続いてuse_device（変数名）とすることで、host_dataで囲まれた範囲の中では、指定された変数のGPU側のアドレスを使うことができるのです。
この指示文はライブラリの呼び出しに使える他、CUDAと組み合わせる際にも利用できます！
Host_data指示文を使うことで、無事cublasを使うことができました。
その性能は…3,950 GFlops！
10倍以上速くなりましたね。
今回の行列積のように、計算量オーダーの大きい計算パターンでは、OpenACCで十分な性能が得られないケースが多いです。
計算量とは、データ量nに対して必要な計算回数のことです。
例えば行列サイズが10*10=100の行列の積であれば、浮動小数点演算は、足し算が100*100*(100-1) = 990,000回、掛け算が100*100*100 = 1,000,000回、合わせて大体200万回。これはおおよそ、データ量nに対してn3に比例した回数になるので、O(n3)の計算量オーダーと言います。
計算量オーダーの大きい計算パターンを高速化する上で、必須の最適化が高速なメモリの利用であり、GPUの場合はshared memoryの利用なのです。
つまりOpenACCは、設計段階から計算量オーダーの大きい計算パターンを捨てている、と言っても過言ではありません。
その代わりに、今回のhost_data指示文のように、cublasなどのライブラリやCUDAと連携するための機能を、これまた設計段階から用意していたのです。
このOpenACCの設計思想をくみ取れば、
- ひとまず全てOpenACCで並列化
- OpenACCでは十分に性能が得られなかった部分を、ライブラリやCUDAに置き換えて高速化
という手順が王道であることがわかります。
無理にOpenACCで完結するのではなく、積極的にライブラリなどを使ってみてください。
1ヵ月間有効のスパコンお試しアカウント
東京大学情報基盤センターでは、教育の一環として、制限はあるものの一ヵ月の間有効なスパコンアカウントを提供しています。
現在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
