2021年6月7日（月）10:00-17:30
第15回スーパーコンピュータ「不老」利用型講習会 OpenACC（初級）
Zoomによるオンライン開催

名古屋大学 情報基盤センター

大島聡史（問い合わせ先 ohshima@cc.nagoya-u.ac.jp）

第15回スーパーコンピュータ「不老」利用型講習会
OpenACC（初級）12

プログラムと時間の目安• 9:30 Zoom接続開始• 10:00 – 12:00

イントロダクション、端末設定など– 名古屋大学情報基盤センターの計算機および利用形態– スーパーコンピュータ「不老」の使い方、ジョブの投入方法、実行確認– スーパーコンピュータ「不老」へのログイン• 13:00 - 17:00

並列プログラミングの基本とOpenACCの学習– OpenACCの基礎：仕様と使い方– 並列計算の考え方とOpenACCプログラムの最適化– 並列化演習• 17:00 - 17:30

自由演習、スパコン利用相談会

3

はじめに
• 講師について– 名前：大島 聡史（おおしま さとし）
• 情報基盤センター 准教授– 出身：栃木県塩谷郡（那須と日光の間）– 主な経歴• 出身大学 電気通信大学• → 東京大学 情報基盤センター 助教（2009.09-2017.03）• → 九州大学 情報基盤研究開発センター 助教（2017.04-2019.06）• → 名古屋大学 情報基盤センター 准教授（2019.07-）• スパコンの運用・調達に関わりながら最新の計算機環境の活用について研究– 「GPUコンピューティング（およびアプリケーションのGPU化）」– 「並列数値計算（行列計算、疎行列ソルバー、ライブラリ）」– 「プログラミング環境（言語やライブラリ）」4

サンプルプログラム（ソースコード）について• 資料中で用いているソースコードは /home/center/a49979a/share/20210607.tgz にて公開しています• cpでコピーして tar zxvf 20210607.tgz などで展開して使ってください• コピーコマンド例：cp /home/center/a49979a/share/20210607.tgz .
• 注意：ディレクトリにcdで移動できるようにはしていません。ファイルコピーして使ってください。
（または直接tarファイルを展開してください。）• 見た目の都合等から公開コードと資料中のコードが完全一致しているわけではない点に注意
してください• 資料の1ページに収まるように省略していることなどがあります• 講習会用のアカウントは本日中のみ有効です• 必要なファイルはscpやsftpで手元にコピーしておいてください5（日本語で書かれた）参考書の例• 「計算科学のためのHPC技術1 」– 下司雅章 (編集), 片桐孝洋 , 中田真秀, 渡辺宙志, 山本有作, 吉井範行,
Jaewoon Jung, 杉田 有治, 石村和也, 大石進一, 関根晃太, 森倉悠介, 黒田久泰，著– 大阪大学出版会、ISBN-10: 4872595866、ISBN-13: 978-4872595864、
発売日：2017年4月3日• 「並列プログラミング入門：サンプルプログラムで学ぶOpenMPとOpenACC」– 片桐 孝洋 著– 東大出版会、ISBN-10: 4130624563、ISBN-13: 978-4130624565、発
売日：2015年5月25日• OpenACCの仕様書や参考Webサイトについては講習会の中で紹介

6

ノードの使い分け（ログインノードと計算ノードの違い）• ログインノード– 待たずに利用できる– 多ユーザがリソースを共有している（メモリをたくさん必要とするような処理はしないこと）– 外部ネットワークに直接（ではないが、近く）つながっている– PCI-Express版のTesla V100が搭載されているためGPUの動作確認程度はできる– プログラムのコンパイル・インストールや動作確認、データの整理といった用途を想定– ログインノードからバッチジョブを投入すると計算ノード上で計算が実行される• 計算ノード– ノードが空いていないと利用できない– リソースを占有できる（例外：*-shareの名前のリソースグループ）– 外部ネットワークへのアクセスはログインノード経由– バッチジョブ実行が基本（例外：*-interactiveの名前のリソースグループ）• ホットストレージはどちらからも同じように見える– ホームディレクトリも/data/group1/${USER}も利用可能（後者の利用を推奨）

7

改行コードの確認と変換
• 確認

• 変換

– odコマンドを使うと簡単に確認できる– 使い方：od -c 確認したいファイル

Windows標準の改行コードCR+LFに対応

Linux標準の改行コードLFに対応※ | head は先頭の数行（デフォルトで10行）だけ表示させたいときによく使う方法。確認したいファイルの先頭10、行という意味ではなく、実行結果の先頭10行。

– nkfコマンドを使うと簡単に変換できる– 使い方• nkf -Lu 変換元ファイル > 変換先ファイル– 左側の例は実際に以下のコマンドで変換した結果である• nkf -Lu testjob.sh > testjob2.sh
改行コード
のせいで
→
エラーした例

8

内容• OpenACCを学ぶ前に– 並列計算の基礎– GPUとOpenACC• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• その他、OpenACCに関する話題– 最適化のための一般的なヒント– デバッグとプロファイリング
– 複数GPUの活用 など• まとめ• 参考（演習課題）：CG法のOpenACC化

9

並列計算とは？
• 並列計算とは何か？並列プログラミング（並列化プログラミング）とは何か？
– 並列計算：何らかの計算処理を同時並行的に行うこと– 並列プログラミング（並列化プログラミング）：並列化・並列計算を行うためのプログラミング
逐次計算

2並列計算
CPU 1

CPU

• 並列？ 並行・平行？

CPU 2
time

– 並列
• 同じ処理を同時に行う、ある処理を幾つかのサブ処理に分けて同時並行的に行う

time– 1つのアプリケーションを高速に実行するために分割して同時実行するイメージ– 並行・平行• 何らかの処理を同時に行う、時分割多重のような疑似並列処理も含む– 直接関係のない複数の処理を同時に行うイメージ– 多数のプログラムの動作を調停するシステムソフトウェア（OSなど）の話をする場合はこちら» 直接の関係がない処理同士だが、同じ資源を利用するために調整が必要

10

なぜ並列計算を行うのか？
• 短時間で終わらせたい計算があるから、計算機（CPUやパソコンなど、計算を行うハードウェア）の持つ能力をフル活用したいから– 現代の計算機は並列計算により高い性能を達成、並列計算を行わなければ高い性能を得られない• PC用CPUもスマートフォン用CPUもマルチコアCPUが主流– マルチコアCPU：コア（＝計算をするユニットのかたまり）を複数搭載したCPU• GPUも「超」並列計算向けのプロセッサ– HWの性能を十分に引き出すには並列計算が必須• 逐次計算では搭載された全ての計算コアを使い切れず、性能を活かしきれない

Surface Pro7の公式サイトの例。
プロセッサの選択肢が3つあり、コア数も明記されている。
「コア数＝性能」ではないが、重要な性能の目安ではある。
画像引用元：https://www.microsoft.com/ja-jp/surface/devices/surface-pro-7/tech-specs（2021.04.17の掲載内容）

11

最近のマルチコアCPUの性能の例
CPU名

コア数（ス
レッド数）

クロック周波数×同時実行命令数

おおよその理論演
算性能（FP64）

理論メモリ
バンド幅

備考

SPARC64VIIIfx

8 (8)2.0 GHz × 8

128 GFLOPS

64 GB/s「京」コンピュータに搭載

SPARC64XIfx

32 (32)

2.2 GHz × 16

1.1 TFLOPS

480 GB/s

名大旧システムFX100に搭載

A64FX

48 (48)2.2 GHz × 32

3.3 TFLOPS

1,024 GB/s

「富岳」、「不老」Type Iに搭載Intel Xeon Gold 6230

20 (40)3.00 - 3.70 GHz × 32

1.3 TFLOPS

140.75 GB/s

「不老」Type IIに搭載

IBM POWER922 (88)3.07 GHz × 8

540 GFLOPS

170.7 GB/s

Summit（2018年11月時点の世界最速スパコン）に搭載

Intel Core i9 10900K

10 (20)3.70 - 5.30 GHz × 16

592 GFLOPS

45.8 GB/s

最近のPC向けハイエンド

AMD Threadripper
3990X

64 (128)

2.90 - 4.30 GHz × 16

3.0 TFLOPS

102.4 GB/s

最近のPC向けハイエンド

•
•
•
•

1FLOPS=1秒間に（倍精度）浮動小数点演算（加算・乗算）を1回行える性能
1K=1000, 1M=1000K, 1G=1000M, 1T=1000G, 1P=1000T, 1Exa=1000P
複数のコアを搭載するのが当たり前、コア数は徐々に増加してきている
実際のプログラムの性能はコア数や（理論）演算性能だけでは決まらない点には注意が必要

12

GPUの性能の例• HPC分野でよく利用される（された、されるであろう）GPUの例
コア数

動作周波数

理論演算性能(FP64)

理論メモリ
バンド幅

製品発表時期など

Tesla P100

3,584 (64*56)

1,328 ~ 1,480 MHz

5.3 TFLOPS

720 GB/s

2016.04 製品発表

Tesla V100

5,120 (64*80)~ 1,530 MHz

7.8 TFLOPS

900 GB/s

2017.06 製品発表
「不老」Type IIに搭載

NVIDIA A100

6,912 (64*108)~ 1,410 MHz

9.7 TFLOPS

1,555 GB/s

2020.05 製品発表
普及はこれから• GPUはCPUと比べて……
– 大量の計算コア、低い動作周波数、高い理論演算性能、高いメモリバンド幅– 数字以外の違いについては後述➢ 製品発表時期はイベント等で正式に発表した年月。提供開始時期については、
一部ユーザへの先行提供などがあるためよくわからないことが多い。

13

スーパーコンピュータ（スパコン）とGPU• そもそもスーパーコンピュータとは何か？
– 一般的な計算機システムよりも大幅に高い性能をもつ計算機システム、ただし明確な定義はない– ある程度の基準になりそうなものはある：TOP500リストや外為法の規制対象など– 必要条件ではないが、現実的に多数の計算機を連結したシステム• 近年ではGPUを大量に搭載したGPUスパコンが普及– GPUスパコンに搭載されているGPUは高価なため、個人用途・練習環境ではより安価なGPUを使うことも多い（機能や性能に違いはあるが、ある程度は有用に使える）
NVIDIA社

AMD社

Intel社

オンボードグラフィックス

GeForce

RadeonIntel HD Graphics

ゲーミング

GeForce

RadeonプロフェッショナルCG

QuadroFire GL
Radeon PRO

高性能計算

Tesla（ただし最新製品 Fire Stream
はNVIDIA A100）
Radeon Instinct

組み込み、車載など

Tegra, Xavier, Jetson

IntelもXeという新しいGPUをリリースしてより広い範囲をカバーしようとしているようだが、詳細はまだ不明。

14

TOP500 2020年6月版（前回）• 連立一次方程式の数値解法– 巨大なLU分解• 理研RCCSの「富岳」が2位以下を大きく引き離して1位で初登場– 日本の1位は2011年11月の「京」以来– 「富岳」はTOP500、HPCG、Graph500、HPL-AIの4つのベンチマークで世界一を達成• 上位10システム中4システムが新システム（ ）• 上位10システム中6システムがNVIDIA GPU搭載– TOP100中34、全500中141• アメリカと中国のシステムが多い• 上位システムは総コア数が100万を超える• 消費電力も数十MW級– 10MWで一般家庭300世帯程度
画像引用元：https://www.top500.org/

15

TOP500 2020年11月版（最新版）• 「富岳」が2連覇– 前回に引き続き4項目で世界一• 上位10システム中6システムがNVIDIA GPU搭載– 上位50システム中20システムがNVIDIA GPU搭載• 「不老」はType Iが6.62 PFLOPSで41位
Type IIが4.88 PFLOPSで58位– ※サブシステムを繋げて測定しても性能は上がらない– 最近は機械学習などのおかげでGPUの人気が高いが、同程度の理論性能を持つType IとType IIでこれだけ差が付くことからも「富岳」（A64FX、FX1000）の効率の高さがわかる

画像引用元：https://www.top500.org/

16

ちなみに
• 電力効率（性能を消費電力で割った値）を競う
Green500はGPUの独壇場に近い• 2位のMN-Coreは独自のアクセラレータ• 残り2つは「富岳」とそのプロトタイプ

画像引用元：https://www.top500.org/

17

並列計算機の分類• 並列計算機は共有メモリ型と分散メモリ型に大別される• 共有メモリ型：メモリを共有している• 同じメモリ（データ）に各プロセッサが自由に直接アクセスできる、メモリアクセスへの競合（衝突）が起こる• 分散メモリ型として扱うこともできる
計算
コア

計算
コア

計算
コア

共有メモリ

計算
コア

• 分散メモリ型：メモリを共有していない• 分散メモリモデル：各プロセッサが個別のメモ
リ（データ）を持つ、互いに直接アクセスできない、他のコアの持つメモリにアクセスするには通信が必要• 共有メモリ型として扱うためには、分散メモリ
を共有メモリとして扱う特別な仕組みが必要
メモリ

メモリ

メモリ

メモリ

計算
コア

計算
コア

計算
コア

計算
コア

➢ 大規模環境では混在した（階層的な）構成となる➢ 物理的な構成とプログラミング手法は1対1の対応関係ではない➢ ノードの概念とも1対1ではない（基本的にはノード内が共有、ノード間が分散ではある）

18

並列計算環境を利用するための手段• 並列計算環境が普及した今日では様々な「並列化のための道具」が使われている– バラバラだと提供側・利用側ともに不便なため共通化されることが多い• ある環境向けに書いたものが別の環境でも利用できる（互換性）• 性能まで互換性があるとは限らない点には注意が必要（性能可搬性）– ある環境では有効であった最適化が別の環境では性能低下要因になることも– 環境が変わってもその環境ごとに常に最大の性能が得られるようにするための研究も行われている→自動チューニング• 自動化はできないのか？
– 全く不可能なわけではない、ある程度はコンパイラ等が行ってくれる– 「コンパイラ様のご機嫌を伺う」コードを書く必要がある、コンパイラによって差が大きい、簡単なコードでないとうまくいかないことが多い• 単純な行列積などコードの形状と最適化のパターンが決まっているものは人が書くよりコンパイ
ラやライブラリの方が得意– OpenMPやOpenACCは簡単・少量の記述で効果的な並列高速化を可能とする（良いところ取り、規格化により汎用性・可搬性も担保）

19

並列計算を行う方法（言語など）の例• コンパイラによる自動並列化：SIMD並列化など一部の処理で有用• pthread：スレッド並列処理のための関数群– pthread_createなどのスレッド操作関数を使う– 現在ではアプリケーションコードで直接利用することはあまりない• OpenMP：スレッド並列処理、主にループ並列化向けの指示文規格– #pragma omp parallel for 、!$omp parallel do– 最近の規格ではタスク処理やGPU対応などを拡充• OpenACC：GPU向けのOpenMPのようなもの• CUDA：NVIDIA社のGPU向け、GPUを普及させた最大要因の一つ、
NVIDIA GPUの能力をフル活用できる• OpenCL：「汎用版CUDA」のようなもの、FPGAなどでも利用可能• MPI：プロセス間の通信規格、特に複数ノード利用時に必須– MPI_Send, MPI_Recv, MPI_Gatherなどの通信関数を使う• 必要に応じてこれらを組み合わせて用いる（OpenMP+MPIなど）• ノード内CPU内
スレッド並列化• 共有メモリモデル• GPU内並列化• デバイス内では共有メモリモデル、外部とのやりとりは分散メ
モリモデル• ノード間並列化• 分散メモリモデル

20

並列計算の基本的なイメージ：ループ並列化（C版）• 元となる逐次計算– 単純な繰り返しループ計算• ループそのものを分割し、同時に計算
PE1

for(i=0; i<N/2; i++){A[i] = B[i] + C[i];D[i] = E[i] + F[i];}PE2for(i=N/2; i<N; i++){A[i] = B[i] + C[i];D[i] = E[i] + F[i];}for(i=0; i<N; i++){A[i] = B[i] + C[i];D[i] = E[i] + F[i];}• ループ内の処理を分割し、同時に計算
PE1

for(i=0; i<N; i++){A[i] = B[i] + C[i];}PE2for(i=0; i<N; i++){D[i] = E[i] + F[i];}PE: Processor Element、CPUやGPUなど計算をする「何か」• OpenMPやOpenACCはこちらのイメージ

21

並列計算の基本的なイメージ：ループ並列化（Fortran版）• 元となる逐次計算– 単純な繰り返しループ計算• ループそのものを分割し、同時に計算
PE1

do i=1, N/2
A(i) = B(i) + C(i)D(i) = E(i) + F(i)end doPE2

do i=N/2+1, N
A(i) = B(i) + C(i)D(i) = E(i) + F(i)end dodo i=1, N
A(i) = B(i) + C(i)D(i) = E(i) + F(i)end do• ループ内の処理を分割し、同時に計算
PE1

do i=1, N
A(i) = B(i) + C(i)end doPE2

do i=1, N
D(i) = E(i) + F(i)end do

PE: Processor Element、CPUやGPUなど計算をする「何か」• OpenMPやOpenACCはこちらのイメージ

22

並列計算の基本的なイメージ：タスク並列化• ループによる並列化よりも大きい粒度の並列計算に適する• 大規模なプログラム・複雑なプログラムの並列化に有用• OpenMPにおいても近年サポートが活発– task構文• GPU（OpenACC）には向いていないため
今回の講習会では扱わない

23

並列計算の限界• 並列化できないプログラムも少なくない– 計算順序に制約がある場合は困難• 工夫により可能となることもある

for(i=1; i<N; i++){A[i] = A[i-1] + B[i];C[i] = A[i] + D[i];}例：ループイタレーション間にも、同じループイタレーション内にも、依存関係がある ⇒ 並列化が困難• （プログラム高速化全般に言えることだが）速くした部分しか速くならない
逐次処理A

一部並列 処理A

処理B
B

time
頑張って処理Bを並列高速化したが、全体の実行時間を考えると……

– 探索問題などでは分割数（並列度）よりずっと速くなることもある★解• 探索対象• 逐次探索★
• 並列探索★

プログラムのどこを並列化させるか、GPUに担当させるかを
しっかり考える必要がある

24

性能向上率と並列化の限界• アムダールの法則– 並列化できる範囲の割合と、並列化により得られる性能向上の関係• そもそも対象問題（プログラム）に十分な並列度がなければならない• 並列化できない部分が大きいといくら並列化しても時間が短くならない• 並列化に伴い通信などのオーバーヘッドが増える可能性がある（もし並列化可能な部分がプログラム全体の＊＊％を占めていたら、性能は最大で＊＊倍まで向上しうる。）– 逐次実行では不要であった通信が増える– 共有メモリモデルであるOpenMPに通信はない
が、スレッドの生成・破棄やスレッド間の同期は必要でありオーバーヘッドとなる※グラフはWikipedia「アムダールの法則」から引用

25

この講習会の目的：GPUを活用する方法の基礎を学ぶ• 何故GPUを活用する必要があるのか？
– GPUは非常に高い性能をもつハードウェアであり、うまく活用できれば大変強力な研究の道具となる：成果を得るための加速装置– GPUは国内外の様々な計算環境に導入されているため、利用スキルを得ておくことはきっとプラスになる• 名古屋大学 情報基盤センターでもType IIサブシステムにGPUを搭載• どうやって使う？
1. 並列化されたソフトウェア（アプリケーション）を使う
2. （並列化されていないプログラムから）並列化されたライブラリやフレームワークを使う
3. 自分で並列計算を実装する• 自分が扱いたい任意の問題（アプリ）をGPU化できる• 扱いたい問題が既にGPU化されているならそれを使っても良いが、GPUのことを理解しているかどうかでさらに高い性能が得られるかどうかが決まる、こともある– 例：特定の形状・大きさの行列しか扱わない26

GPU (Graphics Processing Unit)• 画像処理用のハードウェア– CPUやマザーボードに組み込まれたチップ、または拡張スロットに搭載するビデオカードとして広く普及– 本来の役割：高速・高解像度描画、3D描画処理（透視変換、陰影・照明）、画面出力
3次元画像描画の手順① (2, 2)② (8, 3)③ (5, 7)
オブジェクト単位、頂点単位、ピクセル単位で同時（並列）処理が可能– 同時にできることが多いため、それにあわせたハードウェアへと進化– 数値シミュレーション（科学技術計算）を高速に行えるハードウェアと両立することがわかり、GPUの重要な市場の一つとなってきた– 現在では特にビッグデータ、機械学習、AI処理などで性能を発揮

27

CPUとGPUの違い• GPUはCPUと比べて単純な処理を（順序を問わず）たくさん行うことが求められるため、それに適したハードウェアへと進化してきた– HW構成バランスのイメージ
マルチコアCPU（低性能多コアな）
メニーコアプロセッサ

GPU

計算ユニット
メモリ
キャッシュなど
制御部など– GPUの特徴：たくさんの計算ユニット、高速なメモリ、OSレス• OSを動かし様々な仕事をせねばならないCPUとは大きく異なる• 単体で利用できない、CPUによるサポートが必要• 現代の計算加速装置（アクセラレータ）の代表格• CPUよりスゴイ、というよりも、得意とする仕事が違うことが重要

28

GPUの中身（もう少し詳しく）• GPUとはどのようなハードウェアなのだろうか？
– CPUとはバランスが違う、たくさんの計算コアと高速なメモリを搭載、ということは既に紹介したが……

• 具体的な例（Tesla V100の例）– 以下、NVIDIA TESLA V100 GPU ARCHITECTUREより引用– https://images.nvidia.com/content/volta-architecture/pdf/volta-architecture-whitepaper.pdf• 日本語版の資料はこちら• https://images.nvidia.com/content/pdf/tesla/Volta-Architecture-Whitepaper-v1.1-jp.pdf

29

全体構成

30

全体構成
PCIe
GigaThread Engine

GPC

メ
モ
リ
コ
ン
ト
ロ
ー
ラ
群

GPC＝計算するところ
L2キャッシュ（合計6144KB）

GPC

NVLinkGPC

NVLinkNVLinkGPC

NVLinkNVLinkNVLink

HBM2

HBM2

メ
モ
リ
コ
ン
ト
ロ
ー
ラ
群

GPC

HBM2

HBM2

GPC

GPCの中身• GPC: GPU Processing Cluster 31• TPC: Texture Processing Cluster• SM: Streaming Multiprocessor• 1GPUに6GPC搭載• 1GPCに7TPCs搭載、各TPCに2SMs搭載• 1SMに……
• 64 FP32 cores• 64 INT32 cores• 32 FP64 cores• 8 Tensor cores• 4 Texture units• GPU全体では6*14=84SMs搭載、合計で……
• 5,120 FP32 cores• 5,120 INT32 cores• 2,560 FP64 cores• 640 Tensor cores• 320 Texture units• これら大量のコアが最大1,530MHzで動作• GPU全体で32GBのHBM2を搭載• SMあたり96KBまでの高速共有メモリも搭載• PCIeやNVLinkでホストCPUや他のGPUなど
と接続

32

SMの中身
• 多数のコアを単一のスケジューラが調整している点も大きな特徴の一つ• 同じタイミングで32スレッドが同じ計算を行う（異なるデータに対して同一の計算を実行）ことに注意が必要• Voltaで変更が入り細かい挙動が変わったため、Pascal以前とVolta以降
では最適化の仕方が異なることがある• Voltaで新しく搭載された計算コア• 4x4行列[FP16or32] ＝ 4x4行列[FP16] × 4x4行列[FP16] ＋ 4x4行列[FP16or32]という計算だけを高速に行うことができる、低精度小密行列積演算加速装置

33

SIMT(Single Instruction Multiple Thread)• CPUで用いられているのはSIMD• NVIDIAの示している例– Single Instruction Multiple Data
– 1命令で複数のデータを扱う– PascalまでのSIMTの例• SIMTは1命令で複数のスレッドが動作する– 言い換えれば、複数のスレッドが同時に同じ命令を実行する• 命令分岐はマスク処理される– 実行はしないがプログラムカウンタは遷移す
る、全スレッドが全分岐を実行するのと同程度の時間がかかる• データ転送分が省かれたりするだけマシではあるが– Volta以降のSIMTの例

34

もっと単純に言うと？
• 『多数の計算コアを搭載した「計算コア群」』が多数搭載されたハードウェア• 合計で数千コアが搭載されている、ただし「計算コア群」の中の計算コアの動作には制限があり、完全にバラバラに計算できるわけではない• スレッド並列化 + SIMD並列化をイメージするとわかりやすい、かもしれない– CPU：16コア・AVX512 → 16コアが個別に動作、各コア内で512bit（64bit*8）がまとめて動作– GPU：84SMs・32FP64cores → 84SMsが個別に動作、各SM内で32のFP64coresがまとめて動作– （ただし細かい動作モデルなどは異なる）• GPU全体の構成やSM内の構成はGPUの世代によって変わるため、具体的な数字を一生懸命覚える必要はない• 実際のプログラミングにおいてはある程度抽象化して扱う– もちろん、最適化の際には具体的な数字を考えることになることもある

35

GPUを用いた並列計算の注意点（性能を得るコツ）• 多数のコア全てを満たすに十分な仕事を与える– 数千のコア全てをフル稼働させる並列度が必要– コア数の数倍以上の仕事があることが望ましい• メモリアクセスを待つ間に別の仕事ができる• WARPを意識する– 内部的には32演算器単位で動作しており、分岐の際などに注意が必要• だったのだが、Voltaから変わってしまった• 変わってしまったが、まずはWARPの動作をイメージしてプログラムを作成すると良い性能が得られやすい• 連続メモリアクセスについて考える– コアレスなメモリアクセスが高速• 詳細は終盤（OpenACCプログラムの最適化）にて解説する– 基本的には、並列度が高くて分岐の少ない単純なプログラムが良い36

参考：CUDAプログラム• 単純な配列コピーの例__global__ void gpukernel(int N, float* C, float* A){
int tid = blockIdx.x*blockDim.x + threadIdx.x;
int nt = gridDimd.x * blockDim.x;
for(int i=tid; i<N; i+=nt){
C[i] = A[i];}}GPU上で行われる処理（GPUカーネル）••

同時にたくさん実行される
自分のIDを元に計算するべき範囲を特定する

int main(int argc, char **argv){
int N = 100000;
float *A, *C;
CPU上で行われる処理
float *d_A, *d_C;
A = (float*)malloc(sizeof(float)*N);C = (float*)malloc(sizeof(float)*N);• 専用の関数などを用いて様々な処理を行う
cudaMalloc((void**)&d_A, sizeof(float)*N);
cudaMalloc((void**)&d_C, sizeof(float)*N);
cudaMemcpy(d_A, A, sizeof(float)*N, cudaMemcpyHostToDevice);
cudaMemcpy(d_C, C, sizeof(float)*N, cudaMemcpyHostToDevice);gpukernel<<<4,4>>>(N, d_C, d_A);
cudaMemcpy(C, d_C, sizeof(float)*N, cudaMemcpyDeviceToHost);
cudaFree(d_A);
return 0;
慣れてしまえばそれほど難しいものではないが、記述量がそ}れなりにあり「カジュアルなGPU利用」にはあまり向かない37

GPU（CUDA）向けライブラリ・フレームワーク• GPUの性能を活用できる様々なライブラリが公開されている– https://developer.nvidia.com/gpu-accelerated-libraries– 例• Deep Learningライブラリ– cuDNN, TensorRT, DeepStream SDK• 数値計算・数学ライブラリ– cuBLAS, cuSPARSE, cuRAND, cuFFT, ...

• 通信、C++クラス、etc.
– NCCL, Thrust, OpenCV, MAGMA, ...

– NVIDIA社自ら手がけるものは「CUDA-X」と呼ばれるようになった– 各自が行いたい処理とマッチしていれば容易にGPUの性能を活用することができる– OpenACCとの連携についてはとくに書かれていないものが大半を占めると思うが、ペナルティなしに連携して利用できるものもある（ユーザがcudaMallocで確保したメモリのポインタを引き渡すタイプだと簡単）

38

内容• OpenACCを学ぶ前に– 並列計算の基礎– GPUとOpenACC• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• その他、OpenACCに関する話題– 最適化のための一般的なヒント– デバッグとプロファイリング
– 複数GPUの活用 など• まとめ• 参考（演習課題）：CG法のOpenACC化

39

OpenACCとは？
• GPU（などのアクセラレータ）向けのプログラムを簡単に記述することができる並列化プ
ログラミング言語（指示文規格）– GPU向けのOpenMPのようなもの、C/C++やFortranで書かれたプログラムに簡単な指示文を追加するだけでGPU対応が可能• コンパイラ向けのコメントを記述する• 最適化を行うにはある程度GPUの知識があった方が良い• マルチGPU・マルチノードについてはMPIなどと組み合わせて利用• 最近はOpenMPのGPU対応も進んでいる、今後どちらが主流になるかはわからない– OpenMPに集約されるという話はあるが、なかなか進んでいない– いますぐに試しやすいのはOpenACC– 幾つかの会社が独自に開発していたものが共通規格として集約されたもの• 初登場が2011年、ようやく10周年を迎える新しい言語– OpenMPと似た部分も多いため、どちらかを学んであるともう一方の学習も容易• 少なくとも指示文を用いた並列化という基本的な部分は共通しているため馴染みやすいはず– AMDのGPUでも使えるという話はあるが、HWを有していないため未調査

40

OpenACCとCUDA• OpenACCは「どんなプログラムでもGPU化できる」「どんなプログラムでも高速化できる」わけではないが、多くのプログラムに対して有用– 自作コードのGPU化を行うことができる上で最も簡単でコストパフォーマンスの高い（労力に対して十分な性能が得られる）方法• CUDAでしか書けない処理も多い– 「とにかく高い並列度で一気に計算すれば良い」という典型的なGPU向けプログラムではOpenACCでも十分高い性能が得られる– CUDAを使うべきプログラム• GPU上の高速共有メモリやシャッフル命令を意識したアルゴリズム• インスタンスIDを意識したアルゴリズム– CUDAはThreadIDなど「ID」を意識して並列処理を記述する– OpenACCは「ID」の概念そのものがない• その他、最新のハードウェア機能をフル活用したい場合– Tensor core、RT core、半精度演算

41

OpenACCに関する情報、ツール、etc.
• 規格などの情報– https://www.openacc.org/• 主な対応コンパイラ– 商用：PGI（→NVIDIA HPC SDKへ移行）、Cray、国家超级计算无锡中心– 研究：Omni、OpenARC、OpenUH、ROSEACC– OSS：GCC• プロファイラやデバッガなど– allinea MAP/DDT、PGI pgprof、NVIDIA nvprof/nvvp• 日本語で書かれた資料– 旧PGIコンパイラ技術情報 | HPC WORLD https://hpcworld.jp/pgi-compiler-archive/• 旧PGIコンパイラの代理店であるソフテック社のWebサイトに掲載されていた資料のアーカイブ• 言語仕様が大きく変わってしまわない限りは今後も有用だろう

42

PGIコンパイラとHPC SDK• PGIは独立したコンパイラ開発会社であったが、2013年にNVIDIAに買収された。PGIコン
パイラは継続して提供（販売、一部無償公開）されていたが、2020年6月末にNVIDIA HPC
SDKに統合された。
– なお、AMDのGPUへの対応はHPC SDKへ統合されるずっと前に打ち切られている• これからOpenACCを始める場合にはHPC SDKを使うのが基本となると思われるため、本資料も全てHPC SDKの利用を前提として作成してある• NVIDIA HPC SDK– Webページ https://developer.nvidia.com/hpc-sdk– 「不老」には対応するmodulefileが用意されているため簡単に利用できる– （SDK自体はユーザ権限でもインストールできるため、最新版をはやく使いたい場合は自分で導入してください）43

OpenACCの想定する典型的なGPU計算環境• 1つの計算ノードに1つのCPUと1つのGPUが搭載された環境• OpenACCを用いてある程度まともな性能を得るために必要なGPUに対する最低限の理解– GPUは多数の計算コア（複数の計算コア群）が搭載されたプロセッサ– CPUとGPUは個別に計算用のメモリを持っており、相手側のメモリを読み書きするには
CPU-GPU間の通信が必要• （CPU-GPU間の通信はあまり速くないため、あまり頻繁に・大量に通信したくない）

計算ノード
CPU

メインメモリ（DDR*）

通信が必要（あまり速くない）
GPU

デバイスメモリ（GDDR*、HBM2）

多数の計算コア（複数の計算コア群）

44

OpenACCの基本的な使い方
1. 専用の指示文を含むソースコードを用意する
特別な拡張子などは定められていない–

2. 専用のコンパイラでコンパイルする
OpenACC向けのオプションを指定する–

3. 実行する
CPU向けの実行可能ファイルと同様にそのまま実行可能• ./a.out– （少なくともHPC SDKであれば）特別なプログラムを介したりする必要はない• LD_LIBRARY_PATHなどの対応は必要–––「不老」ではmodule loadの時点で設定されている

想定されているGPUが利用できないと実行時エラー• nvaccelinfoプログラムでGPU情報が見えていることが前提条件

45$ nvaccelinfo

GPU環境の確認• nvaccelinfoコマンドでGPUの存在が確認できる状態である必要がある
Type IIサブシステム用ログイン
ノードで実行した場合の例 →CUDA Driver Version:
NVRM version:

11020
NVIDIA UNIX x86_64 Kernel Module

Device Number:
0
Device Name:
Tesla V100-PCIE-32GB
Device Revision Number:
7.0
Global Memory Size:
34089730048
Number of Multiprocessors:
80
Concurrent Copy and Execution: Yes
Total Constant Memory:65536
Total Shared Memory per Block: 49152
Registers per Block:65536
Warp Size:
32
Maximum Threads per Block:1024
Maximum Block Dimensions:
1024, 1024, 64
Maximum Grid Dimensions:
2147483647 x 65535 x 65535
Maximum Memory Pitch:
2147483647B
Texture Alignment:
512B
Clock Rate:
1380 MHz
Execution Timeout:
No
Integrated Device:No
Can Map Host Memory:
Yes
Compute Mode:default
Concurrent Kernels:
Yes
ECC Enabled:
Yes
Memory Clock Rate:
877 MHz
Memory Bus Width:4096 bits
L2 Cache Size:
6291456 bytes
Max Threads Per SMP:
2048
Async Engines:
7
Unified Addressing:YesManaged Memory:
Yes
Concurrent Managed Memory:
Yes
Preemption Supported:
Yes
Cooperative Launch:
Yes
Multi-Device:
Yes
Default Target:
cc70

460.32.03

Sun Dec 27 19:00:34 UTC 2020

1GPU分の情報（下にもう1つのGPUの分の情報がつづく）

OpenACCコンパイル時に指定するべきGPU情報46

単純なOpenACCプログラムの例（配列の定数倍計算）• C (vector1.c)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

• Fortran (vector1.f90)#include <stdio.h>
int main(int argc, char **argv){
int i;
int n=10;
double v1[10], v2[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);v2[i] = 0.0;}

OpenACC並列化対象＝GPU上で実行される#pragma acc kernels
for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}

for(i=0; i<n; i++)printf(" %.2f", v1[i]);printf("¥n");for(i=0; i<n; i++)printf(" %.2f", v2[i]);printf("¥n");return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

program main
implicit none
integer :: i, n=10
double precision :: v1(10), v2(10)do i=1, nv1(i) = dble(i)v2(i) = 0.0d0
enddo!$acc kernels
do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels➢ 単純なプログラムであればkernels指示文で対象を指定するだけでGPU化が可能

do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""
do i=1,nwrite(*,'(1H F8.2)',advance="NO")v2(i)enddowrite(*,*)""
end program main

47

OpenACCの実行モデル• 指定した部分だけがGPU上で実行される
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

#include <stdio.h>
int main(int argc, char **argv){
int i;
int n=10;
double v1[10], v2[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);v2[i] = 0.0;}

• 全体としてCPUが「主」、GPUが「従」の関係• 指定されていない部分はCPU上で実行される
C
P
U
G
P
U

#pragma acc kernels
for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}

for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");Cfor(i=0; i<n; i++){printf(" %.2f", v2[i]);}P
printf("¥n");return 0;}

U• CPUからGPUに対して計算指示が行われる• GPU上のメモリを確保、CPUからGPUへ必要なデータを転送
GPU上の計算コアにより並列計算される• CPUはGPUの計算終了を待つ• GPUからCPUへ結果データを転送、GPU上のメモリを解放

48「不老」Type IIサブシステムでOpenACCを使う方法（準備）• ログインノード上でmodule loadコマンドを実行して準備を行う• moduleコマンドの使い方（よく使うオプションとその意味）
module avail

利用可能なモジュールの一覧を表示

module load <modulename>

対象モジュールのload

module unload <modulename>

対象モジュールのunload

module list

load済みモジュールの一覧を表示

module purge

load済み全モジュールのunload• ログイン時は何もloadされていない状態• hpc_sdk/21.2 をloadするとコンパイラ等が利用可能になる– module load hpc_sdk/21.2– /21.2 は省略可能。省略した場合、今後より新しいものをインストールして提供する場合は新しい
ものが自動的に読み込まれるようになる。

•

コンパイル例• nvcまたはnvfortranでコンパイル•49-gpu（対象GPU）• nvaccelinfoの「Default Target」に対応• 例：Tesla P100ならcc60、V100ならcc70• 搭載されているGPUがV100であるためcc70を指定-tp（対象CPU）• hostにすればコンパイル環境が参照される
互換性のある古い環境を指定しても動くはずだが、コン
パイラによる最適化の度合いが下がる。

– -Minfo=accel OpenACC化に関する情報を出力 •– -acc OpenACC指示文を有効化– -gpu OpenACCの対象ハードウェア（対象GPUの種類）を指定– -tp 対象CPUを指定$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o v1_c_acc vector1.c• コンパイル例 main:
データ転送やGPU上の計算コアの使い方は
16, Loop is parallelizable
コンパイラがある程度判断してくれるが、
Generating Tesla code
16, #pragma acc loop gang, vector(32) /*最適でない（問題が起きる）こともある
blockIdx.x threadIdx.x */※出力情報の具体的な
16, Generating implicit copyout(v2[:]) [if not already present]
読み方は後述する
Generating implicit copyin(v1[:]) [if not already present]$ nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o v1_f_acc vector1.f90
main:
12, Generating implicit copyin(v1(:)) [if not already present]
Generating implicit copyout(v2(:)) [if not already present]
13, Loop is parallelizable
Generating Tesla code
13, !$acc loop gang, vector(32) ! blockidx%x threadidx%x

50

実行例
• 正しく実行されている限り、GPU上で計算されていることは意識できない– あえていうなら性能
C版の場合$ ./v1_c_acc
1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 20.00

Fortran版の場合$ ./v1_f_acc
1.00
2.00
2.00
4.00

3.00
6.00

4.00
8.00

5.00
10.00

6.00
12.00

7.00
14.00

8.00
16.00

←元データ
←計算結果

9.00
18.00

10.00
20.00

←元データ
←計算結果

51

GPUが仕事をしていることを確認する 1/2• 幾つかの環境変数を用いることでGPUの利用状況を確認可能• 環境変数 NVCOMPILER_ACC_TIME– export NVCOMPILER_ACC_TIME=1 等と指定してから実行すると、計算や通信の回数や時間が確認できる– ログインノード上での実行例 →GPU上での計算に関する情報
計算回数、計算時間
CPU-GPU間の通信に関する情報
通信回数、通信時間$ export NVCOMPILER_ACC_TIME=1$ ./v1_c_acc
1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 20.00

Accelerator Kernel Timing data/home/center/a49979a/share/20210607/vector1.cmain NVIDIA devicenum=0time(us): 30
16: compute region reached 1 time
16: kernel launched 1 time
grid: [1] block: [32]elapsed time(us): total=26 max=26 min=26 avg=26
16: data region reached 2 times
16: data copyin transfers: 1
device time(us): total=14 max=14 min=14 avg=14
18: data copyout transfers: 1
device time(us): total=16 max=16 min=16 avg=16

52

GPUが仕事をしていることを確認する 2/2• NVCOMPILER_ACC_NOTIFYも有効– 1, 2, 4, 8, 16のビット組み合わせ（論理和）で指定、以下の情報が出力される• 1: GPUカーネル起動• 2: データ転送• 4: regionのentry/exit• 8: wait/sync• 16: デバイスメモリの確保と破棄– 例：export NVCOMPILER_ACC_NOTIFY=3• 3=1と2の論理和、GPUカーネル起動情報とデータ転送情報が出力される

53

NVCOMPILER_ACC_NOTIFY情報の例
1: GPUカーネル起動
launch CUDA kernel file=/path/to/source.c function=main line=16 device=0 threadid=1 num_gangs=1
num_workers=1 vector_length=32 grid=1 block=32

並列実行形状（後述）も確認できる

2: データ転送
upload CUDA data file=/path/to/source.c function=main line=16 device=0 threadid=1 variable=v1 bytes=80
download CUDA data file=/path/to/source.c function=main line=18 device=0 threadid=1 variable=v2 bytes=80
4: regionのentry/exit
Enter enter data construct file=/path/to/source.c function=main line=16 device=0 threadid=1
Leave enter data construct file=/path/to/source.c function=main line=16 device=0 threadid=1
Enter compute region file=/path/to/source.c function=main line=16 device=0 threadid=1
Leave compute region file=/path/to/source.c function=main line=16 device=0 threadid=1
Enter exit data construct file=/path/to/source.c function=main line=16 device=0 threadid=1
Leave exit data construct file=/path/to/source.c function=main line=18 device=0 threadid=1
8: wait/sync
Implicit wait file=/path/to/source.c function=main line=16 device=0 threadid=1
Implicit wait file=/path/to/source.c function=main line=16 device=0 threadid=1
Implicit wait file=/path/to/source.c function=main line=18 device=0 threadid=1

16のデバイスメモリ確保/破棄は省略

54

演習• サンプルプログラム（vector1.c または vector1.f90）をコンパイルし、まずはそのままロ
グインノードで実行してみる– ログインノードにもPCI Express版のV100が搭載されているためGPU向けプログラムを実行可能– （実行前にexport NVCOMPILER_ACC_TIME=1しておくと、GPUが本当に動いていることが確認できる）• さらに、バッチジョブとして実行してみる• バッチジョブとして実行できたら、さらに環境変数NVCOMPILER_ACC_TIMEや
NVCOMPILER_ACC_NOTIFYをセットして実行してみる• 注意– ログインノードは数が限られる（誰かがGPUを使っていると実行できない）ため、最初の動作確認など一部の場合を除いてはバッチジョブとして実行すること実際にやってみましょう55

OpenACCプログラムの構成• 指示文：directive• 指示節（指示句）：clause
この資料中ではまとめて指示文と呼ぶことにする• 指示文により全てを記述する– 指示文：コンパイラに対して指示を行う特殊なコメント• C/C++：#pragma acc ~• Fortran：!$acc ~• 基本的には「無視してしまっても問題が起きない文」– コンパイラが対応していない場合もコンパイルと実行自体は可能、もちろんGPUは使えない– 対応環境向けコードとそれ以外を別のコードに分けなくても良いので管理が楽» 対応の有無によって別のアルゴリズムを使い時などはどうにもならないが• 具体的な指示文の例– 並列計算の方法を指示するもの• kernels, parallel• loop, seq, collapse• gang/num_gangs, worker/num_workers, vector/vector_length– データの移動について指示するもの• data, enter/exit data, copy{,in,out}, present, update, create, delete

56

並列化対象範囲の指定：kernelsとparallel• 「この範囲内をGPU上で並列実行したい」ことを示す– kernelsとparallelではコンパイラによる解釈の仕方が異なる• parallel：基本的に利用者が細かく指定する• kernels：ある程度コンパイラが判断、手動で調整（上書き）可能• 最適化をしていくと結局同じようなコードになる、はずである• 範囲の途中で離脱するような構造は不可（forループのbreakなど）– 利用する指示節にも違いが生じる
kernelsと共に利用するもの• async / wait
• device_type• if• default• copy系

parallelと共に利用するもの• OpenACCの仕様の元となった指示文規格• async / waitの違いに起因しているようだ• device_type• kernelsとparallelのどちらを用いても良• if
いが、本講義ではkernelsを用いる• default• copy系• num_gangs / num_workers / vector_length• reduction• private57

再掲：単純なOpenACCプログラムの例• C (vector1.c)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

• Fortran (vector1.f90)#include <stdio.h>

1 program main2
implicit none
int main(int argc, char **argv)
3
integer :: i, n=10{
4
double precision :: v1(10), v2(10)int i;
5
int n=10;
6
do i=1, ndouble v1[10], v2[10];7v1(i) = dble(i)8v2(i) = 0.0d0
for(i=0; i<n; i++){9enddo
v1[i] = (double)(i+1);10v2[i] = 0.0;
11}
12 !$acc kernels
13
do i=1, n
14v2(i) = v1(i) * 2.0d0#pragma acc kernels
15
enddo
for(i=0; i<n; i++){16 !$acc end kernels
v2[i] = v1[i] * 2.0;
17}18
do i=1,n19
write(*,'(1H F8.2)',advance="NO")v1(i)for(i=0; i<n; i++)printf(" %.2f", v1[i]);20
enddo
printf("¥n");21
write(*,*)""
for(i=0; i<n; i++)printf(" %.2f", v2[i]);22do i=1,n（OpenMPと同様に）printf("¥n");23
write(*,'(1H F8.2)',advance="NO")v2(i)24
enddo➢ C/C++では指示文直後のループや{ write(*,*)""} で括った部分（構造化ブロック）がreturn 0;25
指示文の適用対象となる}
26 end program main

➢ Fortranではendで閉じる必要がある58

コンパイル時のメッセージ再確認• C$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o v1_c_acc vector1.c
main:
コンパイラの判断によって
16, Loop is parallelizablecopyin, copyoutという命令が生成されたGenerating Tesla code
16, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */16, Generating implicit copyout(v2[:]) [if not already present]
Generating implicit copyin(v1[:]) [if not already present]• Fortran

Tesla(NVIDIA GPU)向けのコードが生成された、ループの並列化が行われた$ nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o v1_f_acc
vector1.f90
main:
12, Generating implicit copyin(v1(:)) [if not already present]
Generating implicit copyout(v2(:)) [if not already present]
13, Loop is parallelizable
Generating Tesla code
13, !$acc loop gang, vector(32) ! blockidx%x threadidx%x※implicit：暗黙的なプログラムには書かれていなかったが、コンパイラが判断した• どのように判断・処理されたのかを確認することは非常に重要※ if not already present とblockidx や threadidxについては後述59

手動による適用：C版• C (vector1.c)
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

• C (vector2.c)

int main(int argc, char **argv){
int i;
int n=10;
double v1[10], v2[10];3 int main(int argc, char **argv)
4 {5
int i;
6
int n=10;
7
double v1[10], v2[10];8for(i=0; i<n; i++){9for(i=0; i<n; i++){v1[i] = (double)(i+1);10v1[i] = (double)(i+1);v2[i] = 0.0;11v2[i] = 0.0;
}
12
}
13
14 #pragma acc kernels copyin(v1[:]) copyout(v2[:])#pragma acc kernels
15 #pragma acc loop gang, vector(32)for(i=0; i<n; i++){16
for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;17v2[i] = v1[i] * 2.0;
}
18
}
19
for(i=0; i<n; i++)printf(" %.2f", v1[i]);20for(i=0; i<n; i++)printf(" %.2f", v1[i]);printf("¥n");
21
printf("¥n");for(i=0; i<n; i++)printf(" %.2f", v2[i]);22for(i=0; i<n; i++)printf(" %.2f", v2[i]);printf("¥n");
23
printf("¥n");24
return 0;
25
return 0;
コンパイラの判断により、「copyin」「copyout」}
26 }「loop gang, vector(32)」が自動的に挿入されていたと思えば良い60

手動による適用：Fortran版• Fortran (vector1.f90)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

program main
implicit none
integer :: i, n=10
double precision :: v1(10), v2(10)• Fortran (vector2.f90)

1 program main2
implicit none
3
integer :: i, n=10
4
double precision :: v1(10), v2(10)5do i=1, n6
do i=1, nv1(i) = dble(i)7v1(i) = dble(i)v2(i) = 0.0d08v2(i) = 0.0d0
enddo
9
enddo
10
11 !$acc kernels copyin(v1(:)) copyout(v2(:))!$acc loop kernels
12 !$acc loop gang, vector(32)do i=1, n13
do i=1, nv2(i) = v1(i) * 2.0d0
14v2(i) = v1(i) * 2.0d0※end loopは不要
enddo
15
enddo
（書いてもエラーには!$acc end kernels
16 !$acc end kernels
ならないようだ）
17
do i=1,n18
do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)19write(*,'(1H F8.2)',advance="NO")v1(i)
enddo
20
enddo
write(*,*)""
21
write(*,*)""
do i=1,n22
do i=1,nwrite(*,'(1H F8.2)',advance="NO")v2(i)23write(*,'(1H F8.2)',advance="NO")v2(i)コンパイラの判断により、「copyin」「copyout」
enddo
24
enddo
「loop gang, vector(32)」が自動的に挿入されていたと思えば良いwrite(*,*)""
25
write(*,*)""
end program main
26 end program main61

参考：指示文の継続行• 指示文行を次の行に継続させることも可能– 長くなってしまったときなどに• C/C++とFortranで少し違うので注意※これはOpenACCの都合というよりCとFortranの仕様の問題（OpenMPでも同様に違いがある）#pragma acc kernels copyin(v1[:]) copyout(v2[:])#pragma acc loop gang, vector(32)for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}!$acc kernels copyin(v1(:)) copyout(v2(:))!$acc loop gang, vector(32)do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels#pragma acc kernels ¥copyin(v1[:]) copyout(v2[:])#pragma acc loop gang, vector(32)for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}!$acc kernels &!$acc copyin(v1(:)) copyout(v2(:))!$acc loop gang, vector(32)do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels• 末尾の¥（バックスラッシュ）で継続• 継続行の先頭にはpragmaは不要• 末尾の&で継続• 継続行の先頭にも!$accが必要

62

CPU-GPU間のデータ転送• CPUとGPUは個別のメモリを持っており、直接相手側のメモリにアクセスできない– 例外あり、詳しくは後述• 適切なデータ送受信を行わねば正しい計算が行えない– GPUカーネル起動時：メインメモリからデバイスメモリへのデータ転送– GPUカーネル終了後：デバイスメモリからメインメモリへのデータ転送• 単純なプログラムでは自動的にデータ転送を行ってくれるが、ある程度複雑な場合には明示する必要がある– 特にC/C++では、コンパイラが配列の長さを把握できないために生じる実行時エラーをくらいやすいため注意が必要– GPUカーネルが生成されなかったり、実行時にエラーしたりする原因となるCPU

メインメモリ（DDR*）

GPU

デバイスメモリ（GDDR*, HBM*）

63

データ転送に関する指示節• kernels指示文に追加して配列のデータ転送を明示する– GPUカーネル実行前に、デバイスメモリを確保し、ホストからデバイスへコピーする• copyin– GPUカーネル終了後に、デバイスからホストへ書き戻し、デバイスメモリを破棄する• copyout※OpenACC2.5からは常に「present_or_*」の挙動となり、存在していれば使い回してくれる。実際にどう扱われるかは– copyin + copyout
コンパイル時のメッセージやNVCOMPILER_ACC_NOTIFYを• copy用いれば確認できる。コンパイル時に[if not already– デバイスメモリを確保するのみ
present]と出ていたのはこのため。
• create– 既にデバイスメモリに存在していることをコンパイラに伝える• present– 存在していない場合のみcopy{,in,out}する• present_or_copy{,in,out}※データを使い回す、該当するものが無ければ実行時エラー

64

データ転送範囲の指定• 配列全体ではなく一部のみを送受信することも可能• 注意：C/C++とFortranでは部分配列の指定方法が異なる– C/C++：先頭と長さを指定する#pragma acc kernels copy(A[head:length])– Fortran：開始点と終了点を指定する!$acc kernels copy(A(begin:end))– A[:N], A(:N) のような省略表記も可能（先頭からN要素が送られる）

65

繰り返しデータ転送を行う場合の問題• GPUカーネルを何度も実行する場合はどうなるだろうか？
int main(int argc, char **argv){
int i, j, n=10;
double v1[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);}

for(j=0; j<10; j++){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}

(vector3.c)

program main
implicit none
integer :: i, j, n=10
double precision :: v1(10)do i=1, nv1(i) = dble(i)enddo

do j=1, 10!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels
enddo(vector3.f90)

66

繰り返しデータ転送を行う場合の問題• GPUカーネルを何度も実行する場合はどうなるだろうか？
int main(int argc, char **argv){
int i, j, n=10;
double v1[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);}

for(j=0; j<10; j++){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}

(vector3.c)

program main
implicit none
integer :: i, j, n=10
double precision :: v1(10)do i=1, nv1(i) = dble(i)enddo

do j=1, 10!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels
enddo(vector3.f90)

デバイスメモリの生成と破棄を繰り返してしまうため、本来は不要なはずの余計な時間がかかる

繰り返し
デバイスメモリの生成
計算
デバイスメモリの破棄

67

データ転送速度• メインメモリやデバイスメモリの転送速度に対してCPU-GPU間のデータ転送速度はずっと
低速、頻繁な通信は避けたい– CPU - メインメモリ• DDR4、100GB/s/socket程度（STREAM Triad実測）– GPU - デバイスメモリ• HBM2、P100で550GB/s程度、V100では800GB/s弱（STREAM Triad実測）– CPU - GPU~200 GB/s• PCI Express Gen.3 x16– 理論性能でも最大16GB/s（×双方向）CPU– GPU - GPU• NVLink、 20(Pascal) or 25(Volta)GB/s×双方向– 1GPUあたり4(Pascal) or 6(Volta)本のNVLink、GPU数によってつなぎ方はかわる– Power系CPUではCPU-GPU間でもNVLinkが使える

メインメモリ（DDR*）

PCIe、~32GB/s
デバイスメモリ（GDDR*, HBM*）GPU~1 TB/s

68

データ転送のみを行う指示文• data指示文– ループ並列化のタイミング以外で、データのみを操作できる#pragma acc data copyin(A) copyout(B)
構造化ブロック!$acc data copyin(A) copyout(B)
構造化ブロック!$acc end data

• enter/exit data指示文– 構造化ブロックを囲まずに自由な位置で送受信を行うことも可能#pragma acc enter data copyin(A)!$acc enter data copyin(A)#pragma acc exit data copyout(B)!$acc exit data copyout(B)– どちらを使っても良い• enter/exit data指示文の方が便利だが、プログラムの見通しが悪くならないように注意が必要– GPU化範囲の前でとにかく全部送信したいとき？
– 複数ソースコードにプログラムが分割されているとき？

69

data指示文によるデバイスメモリの生存期間の最適化• GPUカーネルを何度も実行する場合などにdata指示文が有効
int main(int argc, char **argv){
int i, j, n=10;
double v1[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);}#pragma acc data copy(v1[:])for(j=0; j<10; j++){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}

program main
implicit none
integer :: i, j, n=10
double precision :: v1(10)do i=1, nv1(i) = dble(i)enddo!$acc data copy(v1(:))
do j=1, 10!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels
enddo!$acc end data

デバイスメモリの生成繰り返し
計算
デバイスメモリの破棄
生成と破棄は最初と最後にのみ行われ、無駄がない70

data指示文によるデバイスメモリの生存期間の最適化• GPUカーネルを何度も実行する場合などにdata指示文が有効
int main(int argc, char **argv){
int i, j, n=10;
double v1[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);}#pragma acc data copy(v1[:])for(j=0; j<10; j++){#pragma acc kernels present(v1)for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}

• 配列のアクセスに間接参照がある場合など

program main
はコンパイラが判断に失敗して繰り返し
implicit none
integer :: i, j, n=10
データコピーをしてしまうこともある
double precision :: v1(10)• kernelsにpresent節を加えるとコンパイラの
do i=1, nv1(i) = dble(i)enddo

判断を上書きできる• acc kernels present(v1)!$acc data copy(v1(:))
do j=1, 10!$acc kernels present(v1)do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels
enddo!$acc end data

デバイスメモリの生成繰り返し
計算
デバイスメモリの破棄
生成と破棄は最初と最後にのみ行われ、無駄がない71

実習• vector3.c, vector3.f90を元に、data指示文を挿入したプログラムvector4.c, vector4.f90を作成し、実行時間を比較する– 環境変数NVCOMPILER_ACC_TIMEをセットして実行すると容易に比較が可能
よく見ると、実行回数も実行時間も大きく異なることが分かるはずである

実際にやってみましょう72

実習：実行例
Accelerator Kernel Timing data/home/center/a49979a/share/20210607/vector3.cmain NVIDIA devicenum=0time(us): 117
15: compute region reached 10 times
15: kernel launched 10 times
grid: [1] block: [32]elapsed time(us): total=131 max=28 min=11 avg=13
15: data region reached 20 times
15: data copyin transfers: 10
device time(us): total=55 max=14 min=4 avg=5
17: data copyout transfers: 10
device time(us): total=62 max=16 min=5 avg=6

Accelerator Kernel Timing data/home/center/a49979a/share/20210607/vector3.f90
main NVIDIA devicenum=0time(us): 110
12: compute region reached 10 times
13: kernel launched 10 times
grid: [1] block: [32]elapsed time(us): total=127 max=26 min=10 avg=12
12: data region reached 20 times
12: data copyin transfers: 10
device time(us): total=53 max=13 min=4 avg=5
16: data copyout transfers: 10
device time(us): total=57 max=13 min=4 avg=5

Accelerator Kernel Timing data/home/center/a49979a/share/20210607/vector4.cmain NVIDIA devicenum=0time(us): 30
13: data region reached 2 times
13: data copyin transfers: 1
device time(us): total=13 max=13 min=13 avg=13
18: data copyout transfers: 1
device time(us): total=17 max=17 min=17 avg=17
15: compute region reached 10 times
15: kernel launched 10 times
grid: [1] block: [32]elapsed time(us): total=116 max=26 min=9 avg=11

Accelerator Kernel Timing data/home/center/a49979a/share/20210607/vector4.f90
main NVIDIA devicenum=0time(us): 29
10: data region reached 2 times
10: data copyin transfers: 1
device time(us): total=14 max=14 min=14 avg=14
18: data copyout transfers: 1
device time(us): total=15 max=15 min=15 avg=15
12: compute region reached 10 times
13: kernel launched 10 times
grid: [1] block: [32]elapsed time(us): total=148 max=34 min=10 avg=14
12: data region reached 20 times

73

多次元配列の送受信• 多次元配列の送受信も可能– ただし、連続したメモリの範囲しか扱えない• 低次元分は全て転送する必要がある➢ C言語の二次元（多次元）配列は
行方向優先
0x**00[0,0]0x**01[0,1]0x**02[0,2]#pragma acc data copyin(A[head:length][0:N])– C/C++• C/C++は右側の次元が低次元

0x**03[1,0]

0x**04[1,1]

0x**05[1,2]

0x**06[2,0]

0x**07[2,1]

0x**08[2,2]– Fortran !$acc data copyin(A(1:N, begin:end))• Fortranは左側の次元が低次元➢ Fortranの二次元（多次元）配列は
列方向優先➢ 物理的なメモリアドレスは一次元
0x**000x**010x**020x**030x**040x**050x**060x**070x**080x**00[0,0]0x**03[1,0]

0x**06[2,0]

0x**01[0,1]0x**04[1,1]

0x**07[2,1]

0x**02[0,2]0x**05[1,2]

0x**08[2,2]

74

動的な配列のデータ送受信• 部分転送時の範囲指定には変数を用いても良い• 動的に確保した配列などコンパイル時に配列の長さが分からないものは長さを明示する必要がある– 間接参照をしているときなどに明示的な指定が必要となりやすい(allocate1.c)(allocate1.f90)

double *v1, *v2;
int *index;
v1 = (double*)malloc(sizeof(double)*n);v2 = (double*)malloc(sizeof(double)*n);index = (int*)malloc(sizeof(int)*n);#pragma acc kernels
for(i=0; i<n; i++){v2[i] = v1[index[i]] * 2.0;}

double precision, allocatable :: v1(:), v2(:)
integer,allocatable :: index(:)
allocate(v1(n), v2(n), index(n))!$acc kernels
do i=1, nv2(i) = v1(index(i)) * 2.0d0
enddo!$acc end kernels➢ v1の範囲（長さ）がうまく認識できず、コンパイルはできるが実行時エラーや正しくない結果になることがある（特にC言語版）➢ copyin(v1[n]) および copyin(v1(n)) を加えると正しく動作する75

動的な配列のデータ送受信• 部分転送時の範囲指定には変数を用いても良い• 動的に確保した配列などコンパイル時に配列の長さが分からないものは長さを明示する必要がある– 間接参照をしているときなどに明示的な指定が必要となりやすい(allocate2.c)(allocate2.f90)

double *v1, *v2;
int *index;
v1 = (double*)malloc(sizeof(double)*n);v2 = (double*)malloc(sizeof(double)*n);index = (int*)malloc(sizeof(int)*n);#pragma acc kernels copyin(v1[n])for(i=0; i<n; i++){v2[i] = v1[index[i]] * 2.0;}

double precision, allocatable :: v1(:), v2(:)
integer,allocatable :: index(:)
allocate(v1(n), v2(n), index(n))!$acc kernels copyin(v1(n))
do i=1, nv2(i) = v1(index(i)) * 2.0d0
enddo!$acc end kernels➢ v1の範囲（長さ）がうまく認識できず、コンパイルはできるが実行時エラーや正しくない結果になることがある（特にC言語版）➢ copyin(v1[n]) および copyin(v1(n)) を加えると正しく動作する76

コンパイル時の出力の比較
allocate1.c

allocate1.f90$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host allocate1.c
main:
21, Accelerator restriction: size of the GPU copy of v1 is unknown
Complex loop carried dependence of v1-> prevents parallelization
Loop carried dependence of v2-> prevents parallelization
Loop carried backward dependence of v2-> prevents vectorization
Accelerator serial kernel generated
Generating Tesla code
21, #pragma acc loop seq
21, Generating implicit copyin(index[:10]) [if not already present]
Generating implicit copyout(v2[:10]) [if not already present]
Generating implicit copyin(v1[:]) [if not already present]Complex loop carried dependence of v1-> prevents parallelization
Loop carried dependence of v2-> prevents parallelization
Loop carried backward dependence of v2-> prevents vectorization$ nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host allocate1.f90
main:
14, Generating implicit copyin(v1(:)) [if not already present]
Generating implicit copyout(v2(1:10)) [if not already present]
Generating implicit copyin(index(1:10)) [if not already present]
15, Loop is parallelizable
Generating Tesla code
15, !$acc loop gang, vector(32) ! blockidx%x threadidx%x

allocate2.c

allocate2.f90$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host allocate2.c
main:
21, Generating copyin(v1[:n]) [if not already present]Complex loop carried dependence of v1-> prevents parallelization
Loop carried dependence of v2-> prevents parallelization
Loop carried backward dependence of v2-> prevents vectorization
Accelerator serial kernel generated
Generating Tesla code
21, #pragma acc loop seq
21, Generating implicit copyout(v2[:10]) [if not already present]
Generating implicit copyin(index[:10]) [if not already present]
Loop carried dependence of v2-> prevents parallelization
Loop carried backward dependence of v2-> prevents vectorization$ nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host allocate2.f90
main:
14, Generating implicit copyin(index(1:10)) [if not already present]
Generating implicit copyout(v2(1:10)) [if not already present]
Generating copyin(v1(:n)) [if not already present]
15, Loop is parallelizable
Generating Tesla code
15, !$acc loop gang, vector(32) ! blockidx%x threadidx%x

配列v1の長さがわからない模様※ただしFortranについてはこのままでも（結果的に）問題ない77

実行結果の比較
配列v1初期値
配列v2初期値
index（C言語版）
index（F90版）

：
：
：
：

1.00
0.00
0.00
1.00

2.00
0.00
1.00
2.00

3.00
0.00
2.00
3.00

4.00
0.00
3.00
4.00

正しいv2の結果

：

2.00

4.00

6.00

8.00 10.00 12.00 14.00 16.00 18.00 20.00

allocate1.cのv2結果 ：allocate1.f90のv2結果 ：

2.00
2.00

0.00
4.00

0.00
6.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00
8.00 10.00 12.00 14.00 16.00 18.00 20.00

allocate2.cのv2結果 ：allocate2.f90のv2結果 ：

2.00
2.00

4.00
4.00

6.00
6.00

8.00 10.00 12.00 14.00 16.00 18.00 20.00
8.00 10.00 12.00 14.00 16.00 18.00 20.00

5.00
0.00
4.00
5.00

6.00
0.00
5.00
6.00

7.00
0.00
6.00
7.00

8.00
0.00
7.00
8.00

9.00 10.00
0.00 0.00
8.00 9.00
9.00 10.00• allocate1.cでは配列長をうまく認識できていないためか計算結果がおかしい。
• どの時点で問題が起きているのか詳細を確認することはできるだろうか？

78

NVCOMPILER_ACC_NOTIFYをセットして動作を確認$ export NVCOMPILER_ACC_NOTIFY=2←2に設定するとCPU-GPU間の通信の情報が出力される

allocate1.cではv1が8byte、つまり1要素
しか送信されていないことがわかる$ ./a1_c_acc
upload CUDA data file=/home/center/a49979a/share/20210607/allocate1.c function=main line=21 device=0 threadid=1 variable=index bytes=40
upload CUDA data file=/home/center/a49979a/share/20210607/allocate1.c function=main line=21 device=0 threadid=1 variable=v1 bytes=8
download CUDA data file=/home/center/a49979a/share/20210607/allocate1.c function=main line=23 device=0 threadid=1 variable=v2 bytes=80
1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

allocate1.f90は正しく10要素(80byte)転送できたようだ$ ./a1f_f_acc
upload CUDA data file=/home/center/a49979a/share/20210607/allocate1.f90 function=main line=14 device=0 threadid=1 variable=v1 bytes=80
upload CUDA data file=/home/center/a49979a/share/20210607/allocate1.f90 function=main line=14 device=0 threadid=1 variable=index bytes=40
download CUDA data file=/home/center/a49979a/share/20210607/allocate1.f90 function=main line=18 device=0 threadid=1 variable=v2 bytes=80
1.00
2.00
3.00
4.00
5.00
6.00
7.00
8.00
9.00
10.00
allocate2.cも正しく10要素
2.00
4.00
6.00
8.00
10.00
12.00
14.00
16.00
18.00
20.00

(80byte)転送できたようだ$ ./a2c_c_acc
upload CUDA data file=/home/center/a49979a/share/20210607/allocate2.c function=main line=21 device=0 threadid=1 variable=v1 bytes=80
upload CUDA data file=/home/center/a49979a/share/20210607/allocate2.c function=main line=21 device=0 threadid=1 variable=index bytes=40
download CUDA data file=/home/center/a49979a/share/20210607/allocate2.c function=main line=23 device=0 threadid=1 variable=v2 bytes=80
1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 20.00$ ./a2f_f_acc
upload CUDA data file=/home/center/a49979a/share/20210607/allocate2.f90 function=main line=14 device=0 threadid=1 variable=v1 bytes=80
upload CUDA data file=/home/center/a49979a/share/20210607/allocate2.f90 function=main line=14 device=0 threadid=1 variable=index bytes=40
download CUDA data file=/home/center/a49979a/share/20210607/allocate2.f90 function=main line=18 device=0 threadid=1 variable=v2 bytes=80
1.00
2.00
3.00
4.00
5.00
6.00
7.00
8.00
9.00
10.00
2.00
4.00
6.00
8.00
10.00
12.00
14.00
16.00
18.00
20.00

79

動的な要素を持つ集合的な要素の転送には注意が必要• 完全なDeep copyができない（動的な要素を持つ集合的な要素をまとめて転送できない）– C/C++：動的に確保された配列をメンバとして持つ構造体やクラスを一括コピーできない– Fortran：allocatable属性やpointer属性を持つメンバを含む派生型を一括コピーできない– 複雑なデータ構造を扱う場合に注意が必要• 解決方法– 従来からの解決方法（手間はかかるが確実）• 必要な分だけ手動でcopyする– 最近の解決方法（簡単だが新しめの機能のため注意して使う必要あり）• コンパイル時の -gpu オプションに deepcopy を追加– Fortranのみ対応、完全なDeep copyができる（はず）• Unified memoryを使う– メインメモリとデバイスメモリを同一に扱う技術– コンパイル時の -gpu オプションに managed を追加、制限もあるため注意が必要– シンプルなデータ構造にコードを書き換える

80

Deep copy問題（1/2）• 動的な要素を持つ集合的な要素を扱う際には注意が必要– C/C++：動的に確保された配列をメンバとして持つ構造体やクラスを一括コピーする場合– Fortran：allocatable属性やpointer属性を持つメンバを含む派生型を一括コピーする場合– ポインタ・アドレスのみがコピーされてその中身がコピーされず、中身を読めないために実行時エラー（エラーしないとしても計算結果があわない）という問題が生じる– メンバ長が固定であれば問題は起きない
struct Sobj{
double values[10];};struct Sobj obj;

type :: Tobj
double precision :: values(10)end type Tobj
type(Tobj) :: obj

for(i=0; i<10; i++){obj.values[i] = (double)(i+1);}

do i=1, n
obj%values(i) = dble(i)enddo#pragma acc kernels
for(i=0; i<10; i++){obj.values[i] = obj.values[i] * 2.0;}!$acc kernels
do i=1, n
obj%values(i) = obj%values(i) * 2.0d0
enddo!$acc end kernels(deepcopy1.c)(deepcopy1.f90)

81

Deep copy問題（2/2）• メンバに動的要素を持つコードの例
struct Sobj{
double *values;};
struct Sobj obj;

←コンパイル時にサイズが未定

type :: Tobj
double precision,allocatable :: values(:)
end type Tobj
↑コンパイル時にサイズが未定
type(Tobj) :: obj! 実行時引数からnを得る処理は省略

実行時の引数からnを取得してmalloc

if(argc<2)n=10; else n=atoi(argv[1]);obj.values = (double*)malloc(sizeof(double)*n);for(i=0; i<n; i++){obj.values[i] = i+1;}

(deepcopy2.c)#pragma acc kernels
for(i=0; i<n; i++){obj.values[i] = obj.values[i] * 2.0;}

allocate(obj%values(n))
do i=1, n
obj%values(i) = dble(i)enddo!$acc kernels
do i=1, n
obj%values(i) = obj%values(i) * 2.0d0
enddo!$acc end kernels(deepcopy2.f90)• Fortran版についてはdeepcopyオプションをつけると正しく転送してくれるようになる• ……はずだったのだが、現在はオプションなしでも正しく動いてしまうようだ– 複雑なプログラムでもうまくいくかは不明$ nvfortran -Minfo=accel -acc -gpu=cc70-tp=host -o d2f_f_acc
deepcopy2.f90$ nvfortran -Minfo=accel -acc -gpu=cc70,deepcopy -tp=host -o d2f_f_acc_dc deepcopy2.f90

82

Deep copy問題の解決方法：Unified Memoryの活用• managedオプションを付けるとC/C++, Fortranいずれも正しく動作する• Unified MemoryというCPU(ホスト)とGPUのメモリをまとめて1つに見せる仕組みを用いる– 便利だが転送効率は下がるため、Unified Memoryを使わず手動で実装した場合より長い時間がかかる可能性がある– 対象プログラムによってはあまり性能が落ちないこともあるため、プログラム移植の際には、ひと
まずmanagedで実装してみてからデータ転送の最適化を行う、という手順は検討に値する

CPU

GPU

CPU

GPU

MEM

MEM

MEM

MEM

UnifiedMemory83

実行例• managedオプションを付けた場合
C$ nvc -Minfo=accel -acc -gpu=cc70,managed -tp=host -o d2_c_acc_m deepcopy2.c
main:
25, Loop is parallelizable
Generating Tesla code
25, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */25, Generating implicit copy(obj.values[:n]) [if not already present]
Generating implicit copyin(obj) [if not already present]

計算結果は正しいが、転送状況を見てもよくわからなくなった……

$ ./d2_c_acc_m
upload CUDA data file=/home/center/a49979a/share/20210607/deepcopy2.c function=main line=25 device=0 threadid=1 variable=obj bytes=8
before1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
after
2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 20.00

Fortran$ nvfortran -Minfo=accel -acc -gpu=cc70,managed -tp=host -o d2f_acc_m deepcopy2.f90
main:
32, Generating implicit copyin(obj) [if not already present]
Generating implicit copy(obj%values(1:n)) [if not already present]
33, Loop is parallelizable
Generating Tesla code
33, !$acc loop gang, vector(128) ! blockidx%x threadidx%x$ ./d2_f_acc_m
upload CUDA data file=/home/center/a49979a/share/20210607/deepcopy2.f90 function=main line=32 device=0 threadid=1 variable=obj bytes=144
upload CUDA data file=/home/center/a49979a/share/20210607/deepcopy2.f90 function=main line=32 device=0 threadid=1 variable=.attach. bytes=8
upload CUDA data file=/home/center/a49979a/share/20210607/deepcopy2.f90 function=main line=32 device=0 threadid=1 variable=descriptor bytes=128
before
1.00
2.00
3.00
4.00
5.00
6.00
7.00
8.00
9.00
10.00
after
2.00
4.00
6.00
8.00
10.00
12.00
14.00
16.00
18.00
20.00

84

Deep copy問題の解決方法：手動でフルコピーする• 大規模なプログラムでは手間がかかるが、全て手動で対応することも可能– enter dataとupdateの組み合わせ（updateについては次頁以降で説明）#pragma acc enter data copyin(obj) copyin(obj.values[n])#pragma acc kernels
for(i=0; i<n; i++){obj.values[i] = obj.values[i] * 2.0;}
#pragma acc update host(obj.values[n])#pragma acc exit data delete(obj.values[n], obj)!$acc enter data copyin(obj) copyin(obj%values(n))!$acc kernels
do i=1, n
obj%values(i) = obj%values(i) * 2.0d0
enddo!$acc end kernels!$acc update host(obj%values)!$acc exit data delete(obj%values, obj)

(deepcopy3.c)(deepcopy3.f90)– data節で対応する場合はcreateとpresentを使うと良い• （createをcopyにするとNG）#pragma acc data create(obj) copy(obj.values[n]){#pragma acc kernels present(obj.values[n])for(i=0; i<n; i++){obj.values[i] = obj.values[i] * 2.0;}}

(deepcopy4.c)!$acc data create(obj) copy(obj%values(n))!$acc kernels present(obj%values(n))
do i=1, n
obj%values(i) = obj%values(i) * 2.0d0
enddo!$acc end kernels!$acc end data

(deepcopy4.f90)

85

データの更新• data指示文内のGPUカーネル間でデータの確認や更新を行いたい場合はどうすれば良いだ
ろうか？ #pragma acc data copy(v1)!$acc data copy(v1){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}••

並列化ループの途中で値を確認したい・更新したい
たとえばここでv1の値を出力したら何が出力されるだろうか？

#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}•!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels!$acc end data

data範囲のあとであれば計算結果が全てメインメモリに書き戻されているのだが……？

for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");(update1.c)

do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""(update1.f90)

86

データの更新• data指示文内のGPUカーネル間でデータの確認や更新を行いたい場合はどうすれば良いだ
ろうか？ #pragma acc data copy(v1)!$acc data copy(v1){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}••

並列化ループの途中で値を確認したい・更新したい
たとえばここでv1の値を出力したら何が出力されるだろうか？

#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}•!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels

計算前の値が出力されてしまう!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels!$acc end data

data範囲のあとであれば計算結果が全てメインメモリに書き戻されているのだが……？

for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");(update1.c)

do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""(update1.f90)

87

データの更新：update指示文• デバイスメモリの生成・破棄を伴わないデータ転送（更新）にはupdateを用いる– ホストからデバイス：update device– デバイスからホスト：update host または update self
– 一部のみの更新も可能、範囲の指定方法はdata指示文と同様#pragma acc data copy(v1){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}!$acc data copy(v1)!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels#pragma acc update host(v1)for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");!$acc update host(v1)do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels!$acc end data

(update2.c)(update2.f90)

88

内容• OpenACCを学ぶ前に– 並列計算の基礎– GPUとOpenACC• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• その他、OpenACCに関する話題– 最適化のための一般的なヒント– デバッグとプロファイリング
– 複数GPUの活用 など• まとめ• 参考（演習課題）：CG法のOpenACC化

89

ループ並列化方法の指定• 前回はkernels節を用いて並列化対象の指定を行った• 現実のプログラム（複雑なプログラム）ではコンパイラにはループ並列化が可能であるかを判断できないことがある• kernels指示文でループを囲んでも、コンパイラが並列化できないと判断することがある• プログラマが並列化の判断を明示すれば並列化させることができる• 逆に、並列化させないこともできる90

loop指示文による並列化判断の指定• loop指示文– 並列化対象ループを指定する• loop指示文に加えて以下の指示節を使うことで判断・動作を制御・修正することができる➢ independent指示節とseq指示節• 対象ループを並列実行するか逐次実行するかを明示する• 強制力があり、コンパイラによる判断は行われなくなる➢ collapse(n)指示節：collapse(2), collapse(3)など• 多重ループをまとめて並列化する• 並列度の低い階層ループの並列化などに極めて重要➢ reduction指示節：reduction(+:a）など• 計算結果の集約などを行う• 多くの場合はコンパイラが正しく判断してくれるため書く必要はないが、コンパイラがどのように判断したか確認できるように読み方を覚えておくと良い

91

三重ループ構造による単純な行列積プログラムの例• C (matmul1.c)• Fortran (matmul1.f90)

double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保
36 #pragma acc kernels
37
for(i=0; i<n; i++){38
for(j=0; j<n; j++){39
for(k=0; k<n; k++){40c[i][j] += a[i][k] * b[k][j];
41
}
42
}
43
}

n=10で実行して計算時間を比較してみたC
37: compute region reached 1 time

double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n)
36 !$acc kernels
37 do i=1, n38
do j=1, n39
do k=1, n
40c(j,i) = c(j,i) + a(k,i) * b(j,k)
41
enddo
42
enddo
43 enddo
44 !$acc end kernels

37: kernel launched 1 time
grid: [1] block: [1]elapsed time(us): total=80 max=80 min=80 avg=80

Fortran 36: compute region reached 1 time
39: kernel launched 1 time
grid: [1] block: [128]elapsed time(us): total=28 max=28 min=28 avg=28

どうやらCとFortranで実行時間に大きな差があるようだ、何故だろ
うか？

92

三重ループ構造による単純な行列積プログラムの例• Fortran (matmul1.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n)
36 !$acc kernels
37 do i=1, n38
do j=1, n39
do k=1, n
nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o mm1_f_acc
matmul1.f90
40
c(j,i) = c(j,i) + a(k,i) * b(j,k)
main:
41already present]
enddo
36, Generating implicit copyin(a(1:n,1:n)) [if not
Generating implicit copy(c(1:n,1:n)) [if not 42
already enddo
present]
43alreadyenddo present]
Generating implicit copyin(b(1:n,1:n)) [if not
44 !$acc end kernels
37, Loop is parallelizable
38, Loop is parallelizable➢ copy関係はコンパイラの判断で39, Complex loop carried dependence of c prevents parallelization
特に問題はない
Loop carried dependence of c prevents parallelization
Loop carried backward dependence of c prevents vectorization
Inner sequential loop scheduled on accelerator
Generating Tesla code
37, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
38,! blockidx%x threadidx%x auto-collapsed➢ 37, 38, 39行目が3重ループ
39, !$acc loop seq
39, Complex loop carried dependence of c prevents parallelization
seq = sequential
ループが逐次実行されることを意味する➢ 3重ループの外側2つは並列化され（auto-collapsedについては後述）、最内ループは逐次実行。c(j,i)が競合しない妥当な判断。

93

三重ループ構造による単純な行列積プログラムの例• C (matmul1.c)
double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保
36 #pragma acc kernels
37
for(i=0; i<n; i++){38
for(j=0; j<n; j++){39
for(k=0; k<n; k++){40c[i][j] += a[i][k] * b[k][j];
41
}
42
}
nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o mm1_c_acc matmul1.c
43
}
main:
37, Complex loop carried dependence of a->->,b->->,c->-> prevents parallelization
Accelerator serial kernel generated
Generating Tesla code
➢ 依存関係があり並列化できず、逐次実行コードが37, #pragma acc loop seq
生成された旨が出力されている➢ 37, 38, 39行目が38, #pragma acc loop seq➢ 正しく実行はできるが、逐次実行のため低速
3重ループ
39, #pragma acc loop seq
37, Generating implicit copy(c[:n][:n]) [if not already present]➢ copy関係はコンパイラの判断
Generating implicit copyin(a[:n][:n],b[:n][:n]) [if not already present]で特に問題はない
Complex loop carried dependence of a->->,c->-> prevents parallelization
38, Complex loop carried dependence of b->->,a->->,c->-> prevents parallelization
39, Complex loop carried dependence of b->->,a->->,c->-> prevents parallelization
Loop carried dependence due to exposed use of c[i1][i2] prevents parallelization

94

loop指示文の追加• C (matmul2.c)
double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保
36 #pragma acc kernels
37 #pragma acc loop independent38
for(i=0; i<n; i++){
並列実行する39 #pragma acc loop independent40
for(j=0; j<n; j++){
41 #pragma acc loop seq
並列実行しない42
for(k=0; k<n; k++){43c[i][j] += a[i][k] * b[k][j];
44
}
45
}
46
}

n=10で実行してみた
38: compute region reached 1 time
42: kernel launched 1 time
C版が劇的に高速化
grid: [1] block: [128]elapsed time(us): total=16 max=16 min=16 avg=16• Fortran (matmul2.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n)
36 !$acc kernels
37 !$acc loop independent
38 do i=1, n
並列実行する39 !$acc loop independent40
do j=1, n
41 !$acc loop seq
並列実行しない42
do k=1, n43c(j,i) = c(j,i) + a(k,i) * b(j,k)
44
enddo
45
enddo
46 enddo
47 !$acc end kernels
36: compute region reached 1 time
42: kernel launched 1 time
grid: [1] block: [128]elapsed time(us): total=26 max=26 min=26 avg=26• 一般的に、Fortranプログラムの方がコンパイラによる並列化判断が適切に働く• C/C++はポインタ参照の都合で不具合が起きないよう保守的な判断がされやすい• loop指示文で指定すればコンパイラの判断を上書きできる95

コンパイラによる判断の比較• C言語、loop independent指定なし(matmul1.c)
37, Complex loop carried dependence of a->->,b->->,c->-> prevents parallelization
Accelerator serial kernel generated
Generating Tesla code
37, #pragma acc loop seq
38, #pragma acc loop seq
各ループに依存性があり並列化できないことが示されている
39, #pragma acc loop seq
37, Generating implicit copy(c[:n][:n]) [if not already present]
Generating implicit copyin(a[:n][:n],b[:n][:n]) [if not already present]Complex loop carried dependence of a->->,c->-> prevents parallelization
38, Complex loop carried dependence of b->->,a->->,c->-> prevents parallelization
39, Complex loop carried dependence of b->->,a->->,c->-> prevents parallelization
Loop carried dependence due to exposed use of c[i1][i2] prevents parallelization• C言語、loop independent指定あり(matmul2.c)
38, Loop is parallelizable
Generating implicit copyin(b[:n][:n]) [if not already present]
Generating implicit copy(c[:n][:n]) [if not already present]
Generating implicit copyin(a[:n][:n]) [if not already present]
40, Loop is parallelizable
指示に従って各ループを並列化したことが示されている
42, Generating Tesla code
38, #pragma acc loop gang, vector(128) collapse(2) /* blockIdx.x threadIdx.x */40,/* blockIdx.x threadIdx.x auto-collapsed */42, #pragma acc loop seq

96collapse版• C (matmul3.c)
double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保
36 #pragma acc kernels
37 #pragma acc loop independent collapse(2)38
for(i=0; i<n; i++){39
for(j=0; j<n; j++){40 #pragma acc loop seq
41
for(k=0; k<n; k++){42c[i][j] += a[i][k] * b[k][j];
43
}
44
}
45
}

• Fortran (matmul3.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n)
36 !$acc kernels
37 !$acc loop independent collapse(2)38 do i=1, n39
do j=1, n
40 !$acc loop seq
41
do k=1, n
42c(j,i) = c(j,i) + a(k,i) * b(j,k)
43
enddo
44
enddo
45 enddo
46 !$acc end kernels• ループを融合してから並列化する• 短いループがネストしている際に有用• 実行時間が短くなるかはループの構造と長さ次第
C

38, #pragma acc loop gang, vector(128) collapse(2) /* blockIdx.x threadIdx.x */39,/* blockIdx.x threadIdx.x collapsed */41, #pragma acc loop seq

Fortran 38, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
39,! blockidx%x threadidx%x collapsed
41, !$acc loop seq

97

collapseのイメージ#pragma acc kernels#pragma acc loop independent collapse(2)for(i=0; i<n; i++){for(j=0; j<n; j++){for(k=0; k<n; k++){mat[i][j] = ...;}}}

#pragma acc kernels#pragma acc loop independentfor(ij=0; ij<n*n; ij++){i = ij/n; j = ij%n;for(k=0; k<n; k++){mat[i][j] = ...;}}}

• 短いループでも多数の計算コアを使い切れるようになる• OpenMPのcollapseと同様の効果

98kernels loop
• kernelsの中にloopがひとつあるだけであれば、kernelsとloopをまとめて書くことも可能– （OpenMPでいうところのparallel for/parallel doと同様）#pragma acc kernels#pragma acc loopfor(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}!$acc kernels!$acc loop
do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels(loop1.c)#pragma acc kernels loopfor(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}

(loop1.f90)!$acc kernels loop
do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels loop

end側は書かなくてもOK（対象が1つのloopであることが自明なため省略できる）

99

OpenACCにおける変数や配列の共有・非共有の扱い• これまでのコードは多重ループを含むが、そのループカウンタなどが競合する可能性について特に気にしていなかった• OpenACCにおける変数や配列の扱い– スカラ変数はfirstprivateとなる• 並列化範囲外の値を引き継ぎ、互いに干渉しあわない– 配列はデバイスメモリにて共有される• 互いに干渉する、private指示節により挙動を変更することも可能• 参考：OpenMPの場合– private/shared指示節で指定– 何も指定しないとsharedとなりスレッド間で干渉する– Cでは並列化範囲内で宣言される局所変数はprivate扱い– Fortranのみ、並列化範囲内の逐次ループのループカウンタはprivate扱い

100

並列実行形状の指定• ここまで、ループをどのようにGPU上の計算コアに割り当てるかは「おまかせ」だった• 簡単なプログラムでは特に問題ないことが多いが、明示的に調整したい場合もある– ネストしたループ（多重ループ）はどのように計算コアに割り当たっている？
– ハードウェアの特徴にあわせたループ構造にしたつもりだが、コンパイラはそれに合わせた実行をしてくれているのか？

• 実はコンパイル時のメッセージにどのように割り当てるかが出力されていた– gang, vector と blockIdx, threadIdx という概念が存在することが確認できる
nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o v1_c_acc vector1.c
16, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o v1_f_acc vector1.f90
13, !$acc loop gang, vector(32) ! blockidx%x threadidx%x– 数字を指定することもできるようだ• vector(32)101gang, worker, vectorとblockIdx, threadIdx• OpenACCではハードウェアに階層的な並列性があることを想定しており、上位階層から順にgang, worker, vectorとなっている• これらをどのように組み合わせて実行するかを指定できる• CUDAの並列実行モデルが階層的になっているため、それにあわせて言語設計された、と考えるのが適切– CUDAはgrid, threadblock, threadの三階層構造– OpenACCのおおまかな並列実行モデルの対応付け
gangvectorblockId 0
threadId 0,1,2…

blockId 1
threadId 0,1,2…

blockId 2
threadId 0,1,2…

102

ハードウェアとのおおまかな対応付け
gang : block
vector : thread• Streaming Multiprocessor(SM)内の並列性はvector、SM単位の並列性はgang• おおまかには• 連続メモリアクセスする最内側のループはvector• より外側のループはgangくらいのイメージ• HWの制約上実際には32コアなどの単位で動作していることを覚えておくと良い性能が出ることもある• ()で数字を与えた場合はその数単位で割り当てられる• ループ長やプロセッサ数にあわせて調整する余地がある、程度に考えておこう
もう少し詳細に言えば……
• gangはHWレベルで同期できない単位の粗粒度並列性（SM単位）• workerはHWレベルで同期できる単位の細粒度並列性（SM内のWARP群）• vectorはworker内部でのSIMDやベクトル並列処理（WARPの中）※外側のループほど上位（gang側）でなければならない：OpenACC 2.0以降で
厳密化

103

演習：単純な行列積プログラム• C (matmul1.c)
double **a=NULL, **b=NULL, **c=NULL;
a = (double**)malloc(sizeof(double*)*n);b = (double**)malloc(sizeof(double*)*n);c = (double**)malloc(sizeof(double*)*n);#pragma acc kernels
for(i=0; i<n; i++){for(j=0; j<n; j++){for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}}• Fortran (matmul1.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n))
allocate(b(n,n))
allocate(c(n,n))!$acc kernels
do i=1, n
do j=1, n
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
!$acc end kernels

matmul1, matmul2, matmul3を実際にコンパイルして実行してみよう➢ 並列化されたか？
➢ 実行時間は短くなったか？
➢ independentやseqを変更してみるとどうか？
➢ collapse指定を変更するとどうか？
➢ 問題サイズを大きくしてみるとどうか？

104

演習：行列積を2回実行してみる• 単純に2回記述すると、その間で全データの転送が発生するはずである。data指示文で通信をなくしてみよう。（もちろん、independentなども指定したうえで。）
double **a=NULL, **b=NULL, **c=NULL;
a = (double**)malloc(sizeof(double*)*n);b = (double**)malloc(sizeof(double*)*n);c = (double**)malloc(sizeof(double*)*n);#pragma acc kernels
for(i=0; i<n; i++){for(j=0; j<n; j++){for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}}#pragma acc kernels
for(i=0; i<n; i++){for(j=0; j<n; j++){for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}}double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n))
allocate(b(n,n))
allocate(c(n,n))!$acc kernels
do i=1, n
do j=1, n
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
NVCOMPILER_ACC_TIME=1!$acc end kernels
や!$acc kernels
NVCOMPILER_ACC_NOTIFY=2
do i=1, nを設定して通信の状況を確認してみよう
do j=1, n
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
実際にやってみましょう!$acc end kernels

105

内容
• OpenACCを学ぶ前に– 並列計算の基礎– GPUとOpenACC• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• その他、OpenACCに関する話題– 最適化のための一般的なヒント– デバッグとプロファイリング
– 複数GPUの活用 など ※やや難しめ、初心者向けとは言い難い話題• まとめ• 参考（演習課題）：CG法のOpenACC化

106

最適化のためのヒント（1／2）• GPUは大量の計算コアをフル活用することで高い合計性能を得るハードウェアであることを忘れてはならない– GPUは多数の計算コアによる並列計算によって高性能を得ているため、並列化対象ループに十分な長さがあるようにする• vector(threadIdx)は32以上、gang(blockIdx)はSM数以上• 短いループはcollapseで結合させるなどする• 問題サイズを大きくしてみる– とにかく、単純な計算を大量に行うことを心がける• CPU-GPU間のデータ通信は最小限にする– 転送回数も転送量も減らす• 実はGPUカーネルの起動と終了にも少し時間がかかっている– 大規模なプログラムでは実行時間に影響する、まとめられる計算カーネルはまとめた方が良い107

最適化のためのヒント（2／2）• 連続メモリアクセスできるようにする– 最内ループでvectorによる連続メモリアクセスを行う– コアレスなメモリアクセスを狙う– （CPUプログラミングでも連続メモリアクセス自体は常識だが）#pragma acc kernels#pragma acc loop gangfor(i=0;i<N;i++){#pragma acc loop vector
for(j=0; j<N; j++){○ a[i][j] = b[i][j] + c[i][j];× a[j][i] = b[j][i] + c[j][i];}}!$acc kernels!$acc loop gang
do i=1, N!$acc loop vector
do j=1, N× a(i,j) = b(i,j) + c(i,j);○ a(j,i) = b(j,i) + c(j,i);end doend do!$acc end kernels• 無理なものは無理：演算性能・メモリ転送性能を超える性能は出ない、CPUに勝てるとは
限らない

108

コアレスなメモリアクセス• 同一SM内の複数CUDAコアによる同時メモリアクセスのアクセス先が近い場合にまとめてアクセスしてくれる機能– 詳細な条件はGPUの世代によって異なる、最新世代ほど条件が緩い• アクセスがバラバラな（遠い）場合
コア0 コア1 コア2 コア3

4回のメモリアクセスが行われる

GlobalMemory• アクセスが揃っている（近い）場合
コア0 コア1 コア2 コア3

1回のメモリアクセスに纏めて処理される

GlobalMemory• 基本的にはvectorループが連続メモリアクセスになるようにすれば良い

109

OpenACC並列化部における関数の呼び出し• kernels/parallel内部（OpenACC並列化対象内部＝GPU上）で関数を実行したい場合– 呼び出される関数にroutine指示文（acc routine）を書いておく必要がある– 例（スペースの都合上、変数の宣言などは省略）#pragma acc routine
void calc(int n, double *v);
int main(int argc, char **argv){#pragma acc data copy(v1){#pragma acc kernels{
どちらか片方でも良い
calc(n, v1);}
（両方書いても問題はない）}}#pragma acc routine
void calc(int n, double *v){
int i;
for(i=0; i<n; i++){v[i] = v[i] * 2.0;}
(routine0.c)}

module modcontains
subroutine calc(n, v)!$acc routine
do i=1, nv(i) = v(i) * 2.0d0
enddo
end subroutine calc
end module modprogram mainuse mod!$acc data copy(v1)!$acc kernels
call calc(n,v1)!$acc end kernels!$acc end data
end program main(routine0.f90)➢ Cでは関数名の直前に、Fortranでは関数名の直後に指示文を挿入➢ しかし、これらの書き方では呼び出された関数内のループは逐次実行になってしまう（calc関数内のループにloop指示文を加えても逐次実行になってしまう）

110

関数内のループを並列実行する方法• routineにgangなどの並列実行形状情報を追加する
void calc(int n, double *v);

#pragma acc routine gang
void calc(int n, double *v);

int main(int argc, char **argv){

int main(int argc, char **argv){#pragma acc data copy(v1){#pragma acc kernels{
calc(n, v1);}}

#pragma acc data copy(v1){#pragma acc kernels{
calc(n, v1);
}
}

どちらに
書いても
良い

}

}

#pragma acc routine gang
void calc(int n, double *v){
int i;#pragma acc loop gang vector
for(i=0; i<n; i++){v[i] = v[i] * 2.0;}}
(routine1.c)void calc(int n, double *v){
int i;#pragma acc loop gang vector
for(i=0; i<n; i++){v[i] = v[i] * 2.0;}}
(routine2.c)関数内部の繰り返し部にloop指定がない場
合はコンパイラ任せになるため、確実に並列化したい場合は指定すること
module modcontains
subroutine calc(n,v)!$acc routine gang
integer :: ndouble precision :: v(*)!$acc loop gang vector
do i=1, nv(i) = v(i) * 2.0d0
enddo
end subroutine calc
end module modprogram mainuse mod
implicit none
integer :: i, j, n=10
double precision :: v1(10)!$acc data copy(v1)!$acc kernels
call calc(n,v1)!$acc end kernels!$acc end data
end program main(routine1.f90)

111

関数内を並列実行する方法：複数ソースファイルの場合（C）• routine指示文を適切に使う必要がある
routine3a.c#pragma acc routine gang
void calc(int n, double *v){
int i;#pragma acc loopfor(i=0; i<n; i++){v[i] = v[i] * 2.0;}}

routine3b.c#pragma acc routine gang
void calc(int n, double *v);
int main(int argc, char **argv){#pragma acc data copy(v1){#pragma acc kernels{
calc(n, v1);}}

• C言語版ではプロトタイプと実体の両方にroutine指示文が必要• コンパイル・リンクは普通に行えばよい（それぞれのコードに-gpu=cc70等を付ける）••••

nvc -Minfo=accel -acc -gpu=cc70 -tp=host -c routine3a.c
nvc -Minfo=accel -acc -gpu=cc70 -tp=host -c routine3b.c
nvc -Minfo=accel -acc -gpu=cc70 -tp=host-o r3_c_acc routine3a.o routine3b.o（最後はリンクするだけなので-Minfo=accelをつけても得られる情報はないはず、つけたままでも特に問題はない）

112

関数内を並列実行する方法：複数ソースファイルの場合（Fortran）• routine指示文を適切に使う必要がある
routine3a.f90
module modcontains
subroutine calc(n,v)!$acc routine gang
integer :: ndouble precision :: v(*)!$acc loop
do i=1, nv(i) = v(i) * 2.0d0
enddo
end subroutine calc
end module modroutine3b.f90
program mainuse mod!$acc data copy(v1)!$acc kernels
call calc(n,v1)!$acc end kernels!$acc end data
end program main

• Fortran版ではただ分割するだけで良い• コンパイルも普通に行える（それぞれのコードに-gpu=cc70等を付ける）•••

nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -c routine3a.f90
nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -c routine3b.f90
nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host-o r3_f_acc routine3a.o routine3b.o113

デバッグと性能解析• デバッグ– 単純なプログラムのOpenACC化でデバッグに困ることはあまりないと思われる– 意図した通りにGPU化やデータ転送ができているか、コンパイラの出力に注意が必要• 性能解析
– NVIDIA社の提供するnvprofやnvvp、nsys(Nsight Systems)などが利用可能• 従来はnvprofとnvvpが使われてきたが、今後はnsysに移行するとNVIDIAが宣言している– GUIが表示されるものはX転送が必要• MovaXtermは標準対応、Putty単体では対応不能、CygwinはX Window関係のパッケージが必要• sshコマンド実行に-Yオプションを付けておく必要がある– デバッグ・プロファイリングに関係するオプションとしてはコンパイラオプション -g や lineinfoがあるが、（少なくともプロファイリングについては）付けても付けなくても変わらないようである（必要なケースがあるのかもしれない）$ nvc -Minfo=accel -g -acc -gpu=cc70,lineinfo -tp=host -o m2_c_acc_gl matmul2.c※実際には-gとlineinfoは同時に指定できない、指定するとlineinfoが無視される

114

nvprof の利用例
nvprof ./a1_c_acc の出力結果例==13056== NVPROF is profiling process 13056, command: ./a1_c_acc
1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00==13056== Profiling application: ./a1_c_acc==13056== Profiling result:Type Time(%)
Time
Calls
Avg
Min
GPU activities:
51.15% 4.9920us
1 4.9920us 4.9920us25.25% 2.4640us
2 1.2320us 1.0560us23.61% 2.3040us
1 2.3040us 2.3040us
API calls:
94.21% 391.79ms
1 391.79ms 391.79ms
5.46% 22.695ms
1 22.695ms 22.695ms
0.20% 842.41us
1 842.41us 842.41us
0.05% 205.17us
4 51.293us 3.0360us
0.01% 60.081us
1 60.081us 60.081us
0.01% 57.935us
1 57.935us 57.935us
0.01% 34.941us
1 34.941us 34.941us
0.01% 34.115us
2 17.057us 4.7160us（中略）
OpenACC (excl):
93.18% 146.14ms
2 73.068ms 10.280us
6.58% 10.324ms
2 5.1619ms 22.703us
0.08% 119.67us
1 119.67us 119.67us
0.05% 85.814us
1 85.814us 85.814us
0.03% 46.566us
2 23.283us 5.6830us
0.02% 39.116us
1 39.116us 39.116us
0.02% 38.304us
3 12.768us 1.7200us

実は内部ではCUDAに関係する関数が呼ばれており
OpenACCがCUDAの上で動いているような関係にある、ということを（再）確認することができる
Max
4.9920us1.4080us2.3040us391.79ms
22.695ms
842.41us181.27us60.081us57.935us34.941us29.399us

Name
OpenACCコードの何行目
main_21_gpuにどれだけ時間がかかっ[CUDA memcpy HtoD]
たかがわかる[CUDA memcpy DtoH]
cuDevicePrimaryCtxRetain
cuMemHostAlloc
cuMemAllocHost
cuMemAlloc
cuMemcpyDtoHAsync
cuModuleLoadDataEx
cuStreamCreate
cuMemcpyHtoDAsync

146.13ms
10.301ms
119.67us85.814us40.883us39.116us34.149us

acc_enter_data@allocate1.c:21
acc_exit_data@allocate1.c:21
acc_enqueue_download@allocate1.c:23
acc_device_init@allocate1.c:21
acc_enqueue_upload@allocate1.c:21
acc_enqueue_launch@allocate1.c:21 (main_21_gpu)
acc_wait@allocate1.c:23

115

nvvp（NVIDIA Visual Profiler）• 測定方法（バッチジョブ内で）：nvprof -o filename ./a.out– -oで出力先ファイル名を指定、これをnvvpで読む• 閲覧方法– nvvpコマンドで起動– メニューのfileからimportを選び、Select an import source:でNvprofを選びNext、Single processでNext、
Timeline data file:で出力された
結果ファイルを選びFinishを選択すれば結果を見ることができる– ログインノード上で実行可能（レスポンスがあまり良くないため
ローカルPCにインストールして実行した方が快適かもしれないが、現在HPC SDKはLinux向けしか存在しない）• 画面例 →

116

Nsight• 新しいプロファリングツール– Nsight Systems, Nsight Compute,
Nsight Graphicsから構成される– CUIとGUIそれぞれで利用できる• CUIによるプロファイル情報の取得と表示– nsys profile ./a.out– --stats=trueオプションをつけておくとその場で統計情報も得られる• 統計情報の一部 →– デフォルトの出力ファイル名はreport1.qdstrm– Nsight Systemsで閲覧するには-oでqdrepファイルを指定しておく必要があるようだ

Generating CUDA API Statistics...
CUDA API Statistics (nanoseconds)
Time(%)
Total Time
Calls
Average
Minimum
Maximum Name------- -------------- ---------- -------------- -------------- -------------- --------------------77.1
23989788
1
23989788.0
23989788
23989788 cuMemHostAlloc
13.2
412383514123835.0
4123835
4123835 cuModuleLoadDataEx
3.8118156310118156.3
2924
1110527 cuMemAlloc_v2
2.9898283
1
898283.0898283898283 cuMemAllocHost_v2
1.44235942021179.75039
52341 cuLaunchKernel
0.6181128286468.92476
20178 cuStreamSynchronize
0.3
94849
7
13549.94444
25066 cuMemcpyDtoHAsync_v2
0.3
82336
6
13722.7
3396
27629 cuMemsetD32Async
0.2
53309
7
7615.63497
29734 cuMemcpyHtoDAsync_v2
0.1
37604
7
5372.01915
13222 cuEventRecord
0.1
28981
1
28981.0
28981
28981 cuStreamCreate
0.1
22506
7
3215.1
1366
5176 cuEventSynchronize
0.0
6956
2
3478.01764
5192 cuEventCreate

Generating CUDA Kernel Statistics...
Generating CUDA Memory Operation Statistics...
CUDA Kernel Statistics (nanoseconds)
Time(%)
Total Time Instances
Average
Minimum
Maximum Name------- -------------- ---------- -------------- -------------- -------------- --------------------12.2
4448
22224.01856
2592 main_159_gpu
11.2
4064
2
2032.01696
2368 main_132_gpu
10.8
3936
2
1968.01760
2176 main_132_gpu__red
10.3
3744
2
1872.01760
1984 main_171_gpu__red
10.2
3712
2
1856.01760
1952 main_191_gpu__red
9.8
3584
2
1792.01664
1920 main_171_gpu
9.5
3456
2
1728.01312
2144 main_123_gpu
9.4
3424
2
1712.01664
1760 main_191_gpu
9.1
3296
2
1648.01504
1792 main_181_gpu
3.8
1376
1
1376.01376
1376 main_142_gpu
3.8
1376
1
1376.01376
1376 main_150_gpu

117

Nsight SystemsのGUIによる結果の表示• nsysで出力したqdrepファイルを読み込んでGUI表示する• WindowsやMacのクライアントがあるため、qdrepファイルをダウンロードして閲覧するのがおすすめ• 例：nsys profile -o report1.qdrep --stats=true ./a.out• Nsight Systemsを起動しFile - Openでqdrepファイルを開いて閲覧• カーネル内の性能調査にはさらにNsight Computeとの連携も必要– クライアントからログインノードに接続して使えると良いのだが、いまのところうまく動かせていない

118

OpenACC API関数の利用• OpenACCは指示文を用いて使うものであるが、その多くの機能は関数によっても提供されている– OpenMPのomp_で始まる関数のようなものと思えば良い– 例• 利用可能なGPUの数を得る– acc_get_num_devices• CPU-GPU間のデータコピー– acc_copyin, acc_copyout– C言語では#include "openacc.h"、Fortranではuse openaccして使う– 特に使わなくても良いが、複数GPUを制御する際などに有用なこともある（後述）

119

外部連携
• OpenACCは配列（メモリ）がホスト上にあるのかデバイス上にあるのかをあまり意識しな
くて良い作りになっている– 例外：update指示文による更新– プログラム中で用いる配列はホストとデバイス両方を指すような概念になっている• 一方で、明示的に一方のメモリを指したくなるときもあるため、そのための記述方法が提供されている– host_data, use_device, deviceptr– うまく使うことでCUDAやMPIとの連携が便利に行える• ホストメモリかデバイスメモリかを具体化することで、そのままcudaMemcpyやMPI通信関数の引数に与えられるようになる（少し後でOpenACC + MPIのより具体的な利用例を示す）#pragma acc data copy(buf){……
#pragma acc host_data use_device(buf)
MPI_Send(buf, n, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD)

120

複数GPUの活用• 1ノードに複数のGPUを搭載した環境も増えている– 「不老」Type IIサブシステムも1ノードに4GPU• 複数のGPUを使うことで、より高速・より大規模なプログラムの実行が可能になる可能性がある• OpenACC内部で利用するGPUを切り替えるにはacc_set_device_num関数を使う• OpenACCプログラムが認識できるGPUを調整・制御するには環境変数
CUDA_VISIBLE_DEVICESを使う– その名の通りCUDAでも有効（NVIDIA GPUを使ったGPUプログラム全般で有効な可能性が高い）• GPUスパコンの構成（たとえば各ノードに搭載されたCPUとGPUの数）は色々考えられるが、1MPIプロセスが1GPUを担当するような実行形態を考えておくと汎用性があって使い
やすい

121

ノード内複数GPUの活用：OpenACCだけを利用する例• 4つの配列を用意し、それぞれ異なるGPUに処理させることを考える#include <stdio.h>#include <openacc.h>
int main(int argc, char **argv){
計算対象の配列を二次元化
int i, j;
int n=10;
double v1[4][10], v2[4][10];for(j=0; j<4; j++){for(i=0; i<n; i++){v1[j][i] = (double)(j+1)+(double)(i+1)/100.0;
v2[j][i] = 0.0;}
GPUカーネルの直前で対象GPUを切り替え}
for(j=0; j<4; j++){acc_set_device_num(j, acc_device_nvidia);#pragma acc kernels
for(i=0; i<n; i++){v2[j][i] = v1[j][i] * 2.0;}}
for(j=0; j<4; j++){printf("GPU %d¥n", j);for(i=0; i<n; i++)printf(" %8.3f", v1[j][i]);printf("¥n");for(i=0; i<n; i++)printf(" %8.3f", v2[j][i]);printf("¥n");}return 0;}

multigpu1.c

program mainuse openacc
integer :: i, j, n=10
double precision :: v1(10,4), v2(10,4)
do j=1, 4
do i=1, nv1(i,j) = dble(j) + dble(i)/100.0v2(i,j) = 0.0d0
enddo
enddo
do j=1, 4
call acc_set_device_num(j, acc_device_nvidia)!$acc kernels
do i=1, nv2(i,j) = v1(i,j) * 2.0d0
enddo!$acc end kernels➢ NVCOMPILER_ACC_TIME=1を設定して実行
enddo
するとGPU別に情報が表示されるため、複数
do j=1, 4
GPUが使えていることが簡単に確認できる
write(*,*)"GPU",jdo i=1,nwrite(*,'(1H F8.3)',advance="NO")v1(i,j)
enddo
write(*,*)""
do i=1,nwrite(*,'(1H F8.3)',advance="NO")v2(i,j)
enddo
write(*,*)""
enddo
end program main

multigpu1.f90

122

ノード内複数GPUの活用：OpenACC + OpenMPの例• 前頁の例では4つのGPUに対する処理が逐次的になってしまう• 各GPUを並列に動作させたい場合はOpenMPやMPIと組み合わせるのが良い

j=0;{
acc_set_device_num(j, acc_device_nvidia);#pragma acc kernels
for(i=0; i<n; i++){v2[j][i] = v1[j][i] * 2.0;
}
}
j=1;
{
acc_set_device_num(j, acc_device_nvidia);#pragma acc kernels
for(i=0; i<n; i++){v2[j][i] = v1[j][i] * 2.0;
}
}
j=2;
{
……

multigpu2.c

multigpu2.f90#pragma omp parallel for private(i)for(j=0; j<4; j++){acc_set_device_num(j, acc_device_nvidia);#pragma acc kernels
for(i=0; i<n; i++){v2[j][i] = v1[j][i] * 2.0;}}

#pragma omp parallel for private(i)for(j=0; j<4; j++){acc_set_device_num(j, acc_device_nvidia);#pragma acc kernels
for(i=0; i<n; i++){v2[j][i] = v1[j][i] * 2.0;
}
}

time

time

• 動作イメージ

j=0;
j=1;
j=2;
j=3;
{
{
{
{
acc_set_device_num(j,acc_set_device_num(j,acc_set_device_num(j,acc_device_nvidia);acc_set_device_num(j,acc_device_nvidia);
acc_device_nvidia);
acc_device_nvidia);#pragma acc#pragmakernelsacc#pragmakernelsacc#pragmakernels
acc kernels
for(i=0; i<n;for(i=0;i++){i<n;for(i=0;i++){i<n;for(i=0;i++){i<n; i++){v2[j][i] =v2[j][i]v1[j][i]=v2[j][i]*v1[j][i]2.0; =v2[j][i]*v1[j][i]2.0; = *v1[j][i]2.0;* 2.0;
}
}
}
}
}
}
}
}

123

ノード内複数GPUの活用：OpenACC + MPIの例（実行方法の例）• MPIプログラム内で利用するGPUを指定するには環境変数CUDA_VISIBLE_DEVICESを使う– CUDA_VISIBLE_DEVICESはプロセスから見えるGPUとその順序を定める環境変数– acc_set_device_num関数でも指定が可能、組み合わせて利用できるが混乱しないように注意– 使用例1job1.shexport CUDA_VISIBLE_DEVICES=0,1./a.out– 使用例2job2.shmpirun -n 4 ./run.shrun.shこの実行方法ではa.outの中では2GPUのみ確認できるようになる。ちなみにCUDA_VISIBLE_DEVICES=3,2と指定するとGPU3とGPU2が順番に見えることになるが、NVCOMPILER_ACC_TIMEの出力上のdevicenumは0と1になるので注意が必要。
1ノード内の4GPUを利用する例。OpenMPIを使う際は環境
変数OMPI_COMM_WORLD_RANK にランク番号がセットされるので、これをCUDA_VISIBLE_DEVICESに与えて各プロセスが異なるGPUを叩くように仕向けている。各プロセスには個別の1GPUのみが見える状態になるため、プログラム内
でacc_set_device_numを指定する必要はない。

#!/bin/bash
export CUDA_VISIBLE_DEVICES=${OMPI_COMM_WORLD_RANK}./a.out

124

ノード内複数GPUの活用：OpenACC + MPIの例（プログラムの例1）
元となるMPIプログラムの概形

OpenACC+MPIプログラムの概形1// 配列rbufとsbufは0.0で初期化されているものとする// 配列rbufとsbufは0.0で初期化されているものとする// プロセスごとに送信用データを生成
for(i=0; i<10; i++)sbuf[i] = (double)rank + (double)(i+1)/100.0;// GPUへ配列をコピー#pragma acc data copy(sbuf[10], rbuf[10]){#pragma acc kernels{// 各プロセスはGPU上でsbufを生成#pragma acc loopfor(i=0; i<10; i++)sbuf[i] = (double)rank + (double)(i+1)/100.0;}// MPI_Reduceで値を集約
MPI_Reduce(sbuf, rbuf, 10, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);}// ここでsbufの値をprintf出力（コードは省略）// 各プロセスのsbufをMIP_Reduceでプロセス0のrbufに集約（加算）MPI_Reduce(sbuf, rbuf, 10, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);// ここでrbufの値をprintf出力（コードは省略）// ここでrbufの値をprintf出力（コードは省略）

プロセス0のrbufに全プロセスが生成したsbufの総和が格納されているため、最後のprintでは正しく総和が表示される。
参考：multigpu_mpi0.c

MPI_Reduceを実行する時点でメインメモリ上のsbufには生成した値が入っていない（GPU上のsbufにのみ生成した値が入っている）ため、
MPI_Reduceで総和を求めても0.0であり、最後のprintfでは0.0が表示される。
参考：multigpu_mpi1.c

125

ノード内複数GPUの活用：OpenACC + MPIの例（プログラムの例2）
OpenACC+MPIプログラムの概形2OpenACC+MPIプログラムの概形3// 配列rbufとsbufは0.0で初期化されているものとする// 配列rbufとsbufは0.0で初期化されているものとする// GPUへ配列をコピー#pragma acc data copy(sbuf[10], rbuf[10]){#pragma acc kernels{// 各プロセスはGPU上でsbufを生成#pragma acc loopfor(i=0; i<10; i++)sbuf[i] = (double)rank + (double)(i+1)/100.0;}#pragma acc update host(sbuf[10])// MPI_Reduceで値を集約
MPI_Reduce(sbuf, rbuf, 10, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);#pragma acc update device(rbuf[10])}// GPUへ配列をコピー#pragma acc data copy(sbuf[10], rbuf[10]){#pragma acc kernels{// 各プロセスはGPU上でsbufを生成#pragma acc loopfor(i=0; i<10; i++)sbuf[i] = (double)rank + (double)(i+1)/100.0;}#pragma acc host_data use_device(sbuf,rbuf)// MPI_Reduceで値を集約
MPI_Reduce(sbuf, rbuf, 10, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);}// ここでrbufの値をprintf出力（コードは省略）// ここでrbufの値をprintf出力（コードは省略）

MPI_Reduceを実行する前にsbufをupdateでホストに持ってくれば、
最後のprintでは正しく総和が表示される。（rbufをupdate deviceしておかないとdata copyによりGPU上のrbufが送られてきて破壊されてしまう点に注意。）

host_data, use_deviceを使えば、MPI_Reduceを実行する際にメインメ
モリ上の配列ではなくGPU上の配列を使ってMPI通信が行われる（CUDA-Aware MPI）ため、最後のprintfでは正しく総和が表示される。
性能も向上することが期待される（が、向上するとは限らず、調査中……）

参考：multigpu_mpi2.c
参考：multigpu_mpi3.c
MPI通信の度にホストメモリに書き戻すのは手間、性能的にも悪影響。 同様にNCCLにも応用可能→参考：multigpu_mpi4.c

126

複数ノード複数GPUの活用：OpenACC + MPI• プログラムとしては「OpenACC+MPIプログラムの概形3」の例のままでよい• ジョブスクリプトに多少の工夫が必要• 2ノード×4GPU（合計8MPIプロセス）の例
job.shrun.sh（chmodで実行権を与えておくこと）

省略#!/bin/bash
exportCUDA_VISIBLE_DEVICES=${OMPI_COMM_WORLD_LOCAL_RANK}
numactl -l $1

module load hpc_sdk/21.2
mpirun -n 8 -machinefile $PJM_O_NODEINF -map-by ppr:2:socket ./run.sh ./a.out

複数ノードにプロセスを適切に配置するためのオプション（詳しくはWebの「スーパーコンピュータ「不老」基本マニュアルおよび関連資料」に掲載している「Type IIサブシステム向けのプロセス・スレッド配置方法」を参照）

127

複数ノード実行に関するソースコードとジョブスクリプト例一覧• job_multigpu_mpi1.sh, job_multigpu_mpi1_a.sh– 1ノード4GPUの実行例（multigpu_mpi0.cからmultigpu_mpi4.cまで比較）• job_multigpu_mpi2.sh, job_multigpu_mpi2_a.sh– multigpu_mpi2の2ノード×4GPU版• job_multigpu_mpi3.sh, job_multigpu_mpi3_a.sh– multigpu_mpi3の2ノード×4GPU版• job_multigpu_mpi4.sh, job_multigpu_mpi4_a.sh– multigpu_mpi4の2ノード×4GPU版• 各ソースコードの冒頭にコンパイルの例も記載– 動作確認はしましたが、最低限度の簡単なプログラムのみで、性能の良し悪しは見ていません。十
分に高い性能を得るためにはシステム側の設定や環境変数の追加が必要かも知れません。新たな情報はセンターのWebにて紹介していきます。

128

内容
• OpenACCを学ぶ前に– 並列計算の基礎– GPUとOpenACC• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• その他、OpenACCに関する話題– 最適化のための一般的なヒント– デバッグとプロファイリング
– 複数GPUの活用 など• まとめ• 参考（演習課題）：CG法のOpenACC化

129

まとめ
• OpenACCの仕様や使い方について紹介した– わずか数行の指示文でGPUが利用できる– OpenMPと似ているため、両方使えると色々と活用できる– GPUが利用できる環境がある場合は是非使ってみて欲しい• 残りの時間は演習時間とします– 演習問題としてCG法の逐次コードを提供しているので、指示文を挿入して並列化してみましょう– その他、任意のプログラムのGPU化やプロファイラのテストなどをしてもらって構いません

130

内容
• OpenACCを学ぶ前に– 並列計算の基礎– GPUとOpenACC• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• その他、OpenACCに関する話題– 最適化のための一般的なヒント– デバッグとプロファイリング
– 複数GPUの活用 など• まとめ• 参考（演習課題）：CG法のOpenACC化

131

CG法プログラムのOpenACC化• 単純なCG法の計算カーネル（反復計算部）を題材に、プログラムのOpenACC化を考える– 簡単にするため、行列は密行列、前処理は対角スケーリング• CG法（共役勾配法、Conjugate Gradient Method）– 対称正定値行列を係数とする連立一次方程式Ax=bを解く手法• 行列A、既知のベクトルb、未知のベクトルx– 基本アルゴリズム• Wikipediaから引用、前処理なし• 利用するコードは計算順序が変更されている版• 単純な行列Aとランダム行列xxからbを求めておき
Ax=bを解いてxとxxが（ほぼ）一致することを確認する、という構造にしてある（すぐに収束してしまう
ため、固定回数ループで時間を測定すると良い）• 時間測定を簡単に書くためOpenMP関数を利用– コンパイル時に-mpオプションを加える必要あり

132

今回用いるCG法コードの計算手順
Compute r (0) = b-[A]x (0)for i= 1, 2, …
solve [M]z (i-1) = r (i-1)ρi-1 = r (i-1) z (i-1)
if i=1p (1) = z (0)else
βi-1 = ρi-1 /ρi-2
p (i) = z (i-1) + βi-1 z (i)endif
q (i) = [A]p (i)αi = ρi-1 /p (i) q (i)x (i) = x (i-1) + αi p (i)r (i) = r (i-1) - αi q (i)check convergence |r|
end

•
•
•
•

初期値、x(0)は適当な値（今回は0ベクトル）
収束するまで繰り返す
前処理、今回はrを対角要素で割るだけ
リダクション• コピー• リダクション• ベクトル積和•••••

行列ベクトル積＝一番時間がかかる処理
リダクション
ベクトル積和ベクトル積和収束判定（中身はリダクションと平方根）※上付き文字は反復回数に対応• Ax=b：Aは行列、xとbはベクトル• z,r,p,qはベクトル• α・β・ρはスカラー（ベクトルのリダクション結果）

for(iter=1; iter<=maxiter; iter++){printf("iter %d ", iter);// {z} = [Minv]{r}for(i=0;i<N;i++){z[i] = dd[i]*r[i]; }// {rho} = {r}{z}※対角要素分の1だけのベクトルddを用意済
rho = 0.0;
for(i=0;i<N;i++){rho += r[i]*z[i]; }// {p} = {z} if iter=1// beta = rho/rho1 otherwise
if(iter==1){
for(i=0;i<N;i++){p[i] = z[i];}}else{
beta = rho/rho1;
for(i=0;i<N;i++){p[i] = z[i] + beta*p[i];}}// {q} = [A]{p}for(i=0;i<N;i++){★行列ベクトル積
q[i] = 0.0;
行ごとの計算を並列に行える
for(j=0;j<N;j++){q[i] += A[i*N+j]*p[j];}}// alpha = rho / {p}{q}pq = 0.0;
for(i=0;i<N;i++){pq += p[i]*q[i]; }alpha = rho / pq;// {x} = {x} + alpha*{p}// {r} = {r} - alpha*{q}for(i=0;i<N;i++){x[i] += + alpha*p[i];r[i] += - alpha*q[i];}// check converged
dnrm = 0.0;
for(i=0;i<N;i++){dnrm += r[i]*r[i]; }resid = sqrt(dnrm/bnrm);
if(resid <= cond){break;}if(iter == maxiter){break;}rho1 = rho;
}

★

主要コード

★

• 行列やベクトルに対する単純な計算ばかりで構成されているため、並列化・OpenACC化は容易– ★：行列要素同士のコピーや四則演算– ★：集約演算（dot product, reduction）– これらを全てGPUカーネルにしてしまえば良さそうである– ※複雑な前処理を適用する場合は難易度が上がる– 具体的にはどのような手順でOpenACC化すれば良いだろうか？
※ すぐに収束してしまうため、ここに妨害用のコードを追加してある– 各説明に対応するコードを提供しているが、まずはcg1.cまたはcg1.f90だけをみて並列高速化を目指してみよう

★

★

★

★

★

133

134

CG法のOpenACC化： 1.各計算部の並列化：行列ベクトル積（1／2）• まずは一番実行時間が長い処理である行列ベクトル積を並列化してみる
151 // {q} = [A]{p}152 #pragma acc kernels loop
153
for(i=0;i<N;i++){154
q[i] = 0.0;155
for(j=0;j<N;j++){156
q[i] += A[i*N+j]*p[j];
157
}
158
}
cg1.c

143 ! {q} = [A]{p}※cg0.c, cg0.f90は指示文が144 !$acc kernels loop
145 do i=1, N
まったく入っていないコード
146
q(i) = 0.0147
do j=1, N
148
q(i) = q(i) + A(j,i) * p(j)149
end do150 end do
cg1.f90• 実はkernelsとloopはまとめて1行で指定することができる• Fortran版の閉じる指示文も不要• コンパイルしたところ、Fortran版は外側ループが並列化されたが、C言語版は両ループともseq扱いされ、配列Aの長さも不明と判断されてしまった
153, Accelerator restriction: size of
the GPU copy of A is unknown
153, #pragma acc loop seq
155, #pragma acc loop seq

145, Loop is parallelizable
Generating Tesla code
145, !$acc loop gang, vector(128) ! blockidx%x threadidx%x147, !$acc loop seq

135

CG法のOpenACC化： 1.各計算部の並列化：行列ベクトル積（2／2）• independent/seqとcopyinを指定
151 // {q} = [A]{p}152 #pragma acc kernels copyin(A[N*N])
153 #pragma acc loop independent154
for(i=0;i<N;i++){155
q[i] = 0.0;
156 #pragma acc loop seq
157
for(j=0;j<N;j++){158
q[i] += A[i*N+j]*p[j];
159
}
160
}
cg2.c

143 ! {q} = [A]{p}144 !$acc kernels
145 !$acc loop independent
146 do i=1, N
147
q(i) = 0.0
148 !$acc loop seq
149
do j=1, N
150
q(i) = q(i) + A(j,i) * p(j)151
end do152 end do
cg2.f90
153 !$acc end kernels• C言語の場合はFortranと比べてcopyやindependentを指定せねば適切に並列化されないことが多いことに気を付けること• コンパイラの出力をしっかり確認することが重要
154, Loop is parallelizable
Generating Tesla code
154, #pragma acc loop gang, vector(128) /* 省略 */157, #pragma acc loop seq

146, Loop is parallelizable
Generating Tesla code
146, !$acc loop gang, vector(128) ! 省略
149, !$acc loop seq

136

CG法のOpenACC化： 1.各計算部の並列化：各計算の並列化• 同様にして、反復計算部の全ての行列・ベクトル計算を並列化していく• コンパイラのメッセージを確認し、狙い通りに並列化されているかを確認する– 特にCの場合はindependent/copy*/present節をしっかり挿入する– reductionが適切に生成されているかを確認する• （このような簡単なプログラムで失敗することはないと思って良いが）
cg3.c

134 // {rho} = {r}{z}135rho = 0.0;
136 #pragma acc kernels
137 #pragma acc loop independent138
for(i=0;i<N;i++){139
rho += r[i]*z[i];}140137, Loop is parallelizable
Generating Tesla code
137, #pragma acc loop gang, vector(128) /* 省略 */Generating implicit reduction(+:rho)

cg3.f90

127
128
129
130
131
132
133
134

! {rho} = {r}{z}rho = 0.0d0!$acc kernels!$acc loop
do i=1, N
rho = rho + r(i) * z(i)end do!$acc end kernels

131, Loop is parallelizable
Generating Tesla code
131, !$acc loop gang, vector(128) ! 省略
Generating implicit reduction(+:rho)

137

状況の確認
cg3.cの例• 反復計算部の全ての行列・ベクトル計算を130: data region reached 20 times
並列化してみた
130: data copyin transfers: 20
device time(us): total=91 max=15 min=3 avg=4• 簡単な問題ではあるが132: data copyout transfers: 10
device time(us): total=62 max=15 min=5 avg=6（問題サイズが小さければ時間もとても短い）
138: compute region reached 10 times
環境変数 NVCOMPILER_ACC_TIME=1
138: kernel launched 10 times
grid: [1] block: [128]を設定して実行してみると
device time(us): total=30 max=3 min=3 avg=3
elapsed time(us): total=178 max=20 min=17 avg=17
頻繁にデータ転送していることがわかる– NVCOMPILER_ACC_NOTIFY=2 を設定すると
kernelsの度に配列全体をコピーしている
こともわかる• 本当にこんなに頻繁に転送する必要はあるのだろうか？

138: reduction kernel launched 10 times
grid: [1] block: [256]device time(us): total=30 max=3 min=3 avg=3
elapsed time(us): total=166 max=19 min=16 avg=16138: data region reached 20 times
138: data copyin transfers: 30
device time(us): total=126 max=12 min=3 avg=4140: data copyout transfers: 10
device time(us): total=68 max=20 min=5 avg=6147: compute region reached 1 time
147: kernel launched 1 time
grid: [1] block: [128]device time(us): total=2 max=2 min=2 avg=2
elapsed time(us): total=18 max=18 min=18 avg=18

138

cg3コンパイル時のメッセージの再確認• 確かに頻繁にコピーが行われている可能性がある
nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o cg3_c_acc cg3.c
main:130, Loop is parallelizable
Generating Tesla code
130, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */130, Generating implicit copyout(z[:N]) [if not already present]
Generating implicit copyin(r[:N],dd[:N]) [if not already present]
138, Loop is parallelizable
Generating Tesla code
138, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */Generating implicit reduction(+:rho)nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o cg3_f_acc cg3.f90
main:120, Generating implicit copyin(dd(1:n),r(1:n)) [if not already present]
Generating implicit copyout(z(1:n)) [if not already present]
122, Loop is parallelizable
Generating Tesla code
122, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
129, Generating implicit copyin(z(1:n)) [if not already present]
Generating implicit copy(rho) [if not already present]
Generating implicit copyin(r(1:n)) [if not already present]131, Loop is parallelizable
Generating Tesla code
131, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
Generating implicit reduction(+:rho)

139

データ転送の最適化を考える• 現在のデータ転送状況のイメージ
Compute r (0) = b-[A]x (0)for i= 1, 2, …
solve [M]z (i-1) = r (i-1)ρi-1 = r (i-1) z (i-1)
if i=1p (1) = z (0)else
βi-1 = ρi-1 /ρi-2
p (i) = z (i-1) + βi-1 z (i)endif
q (i) = [A]p (i)αi = ρi-1 /p (i) q (i)x (i) = x (i-1) + αi p (i)r (i) = r (i-1) - αi q (i)check convergence |r|
end

CPU

GPU

計算をする度に、関係する配列全てを送受信している↓
更新があったデータだけを送受信すれば良いはずでは？

140

CG法のOpenACC化手順例： 2.データ転送の最適化• data節を用いてデータ転送を削減してみる– 送受信が必要なデータはどれだろうか？
#pragma acc data copyin(？) copyout(？)!$acc data copyin(？) copyout(？)

Compute r (0) = b-[A]x (0)for i= 1, 2, …
solve [M]z (i-1) = r (i-1)ρi-1 = r (i-1) z (i-1)
if i=1p (1) = z (0)else
βi-1 = ρi-1 /ρi-2
p (i) = z (i-1) + βi-1 z (i)endif
q (i) = [A]p (i)αi = ρi-1 /p (i) q (i)x (i) = x (i-1) + αi p (i)r (i) = r (i-1) - αi q (i)check convergence |r|end※リダクション結果のスカラ変数は自動的にCPU側に送られるため気にしなくてよい– 途中の計算結果がおかしくないかを確認したい場合には、update節を用いて途中の配列データを覗いてみると良い

141

実は……
• 最初に全てのデータをGPUへ転送し、最後に結果ベクトルを回収するだけで良かった– わかっていればいきなりデータ通信の最適化を行っても良かった– 少しずつの最適化が行いやすい点はプログラム最適化においてとても助かる#pragma acc data copyin(z[N], dd[N], r[N],p[N], q[N], A[N*N]) copy(x[N]){for(iter=1; iter<=maxiter; iter++){!$acc data copyin(z(N), dd(N), r(N), p(N),q(N), A(N,N)) copy(x(N))
do iter=1, maxiter……

for自体を囲む必要あり……

A, p, q全てGPU上にある状
態で行列ベクトル積が開始される// {q} = [A]{p}#pragma acc kernels#pragma acc loop independentfor(i=0;i<N;i++){q[i] = 0.0;
for(j=0;j<N;j++){q[i] += A[i*N+j]*p[j];
}
cg4.c
}

! {q} = [A]{p}!$acc kernels!$acc loop
do i=1, N
q(i) = 0.0
do j=1, N
q(i) = q(i) + A(j,i) * p(j)end doend docg4.f90!$acc end kernels

cg5.c と cg5.f90 はさらに念のためpresent節を加えたもの)

142

初期データ
行列A（サイズN×N）

ベクトルx（サイズN）
全て0.0

対角はN
対角以外は1

ベクトルxx（サイズN）
ランダム
ベクトルb（サイズN）➢ A, x, xxを初期化し、b = A×xxを計算しておいたうえで、cg法でAx = bを計算しxを求めるというプログラムにしてある。
➢ これなら前処理のない単純なCG法でも簡単に収束する。
➢ むしろ簡単過ぎてすぐに収束して終了してしまうため、意図的に値を破壊して収束しないようにしてある。
（行列サイズNと同じ回数分だけ計算して収束しないまま終わる。）

143

実行例
A:
9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000
x:
8.400000
8.700000
7.800000
1.600000
9.400000
3.600000
8.700000
9.300000
5.000000
2.200000b:131.900000
134.300000
127.100000
77.500000
139.900000
93.500000
134.300000
139.100000
104.700000
82.300000
iter 1 1.089293e-01 1.000000e-16
iter 2 2.331098e-16 1.000000e-16
iter 3 5.450800e-18 1.000000e-16
time: 0.000010 sec , 0.000003 sec/iter
result:8.400000
8.700000
7.800000
1.600000
9.400000
3.600000
8.700000
9.300000
5.000000
2.200000
diff:
1 8.400000 8.400000: 3.552714e-15
2 8.700000 8.700000: 5.329071e-15
3 7.800000 7.800000: -8.881784e-16
4 1.600000 1.600000: -8.881784e-16
5 9.400000 9.400000: 0.000000e+00
6 3.600000 3.600000: -1.776357e-15
7 8.700000 8.700000: 0.000000e+00
8 9.300000 9.300000: 0.000000e+00
9 5.000000 5.000000: -8.881784e-16
10 2.200000 2.200000: -1.776357e-15←既知の行列A•••←ベクトルX（求める答え）

第一引数は問題サイズN
第二引数を0にすると初期データや計算結果を出力しなくなる
第三引数に0以外を指定すると収束しても計算を続行する（すぐに終わってしまうのを防ぐもう1つの策）←既知のベクトルb←反復計算の履歴、反復計算終了時に所要時間も表示←計算結果ベクトルX←計算結果ベクトルXと最初に設定したベクトルXの比較