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

名古屋大学 情報基盤センター 大島聡史 （問い合わせ先 ohshima@cc.nagoya-u.ac.jp）

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

お詫び（補足）• hpc_sdk 21.11からコンパイラの挙動が変更になり、動的な配列を扱う際は「tripletな指定」が推奨されるようになっていました。
• Cの動的に確保する配列を用いたプログラムでデータのコピーを明示する際は、data[N]のような長さだけの記述ではだめで、data[0:N]のように先頭と長さを明示する必要があるようです。
• Fortran版のallocatableな配列を用いたプログラムでデータのコピーを明示する際は、data(N)のような長さだけの記述ではだめで、data(1:N)のように上界・下界を明示する必要があるようです。
• いくつかのサンプルプログラムが影響を受けているため注意してください。
– 動いてしまったものもあるため気がつきませんでした。
– コンパイル時の-gpuオプションにimplicitsectionsオプションを追加すると従来の挙動になるようですが、今後もこれが引き継がれるのかはわかりません。
• 参考情報：https://docs.nvidia.com/hpc-sdk/archive/22.5/compilers/hpc-compilers-userguide/index.html#gpu

3

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

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

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

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

4

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

サンプルプログラム（ソースコード）について• 資料中で用いているソースコードは /data/group1/a49979a/share/20220627.tgz にて公開しています• コピーして使ってください（ディレクトリへのcdは許可していません）$ cp /data/group1/a49979a/share/20220627.tgz ./$ tar zxvf 20220627.tgz→20220627というディレクトリが作られる

ディレクトリへのcdは許可していません。
cpで持っていってください。

• 見た目の都合等から公開コードと資料中のコードが完全一致しているわけではない点に注意
してください• 資料の1ページに収まるように省略していることなどがあります• 講習会用のアカウントは今日から三日間（明後日の23:59まで）有効です• 必要なファイルはscpやsftpで手元にコピーしておいてください

6

参考書の例
• 「計算科学のための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日• 「Parallel and High Performance Computing」– Robert Robey, Yuliana Zamora 著– Manning、ISBN-10: 1617296465、ISBN-13: 978-1617296468、
発売日：2021年6月22日• OpenACCの仕様書や参考Webサイトについては講習会の中で紹介する

7

イントロダクション、端末設定など• 「講習会資料-共通編_20220421.pdf」および録画済み動画にて紹介– スーパーコンピュータ「不老」の紹介– ログイン方法と基本的な操作の紹介– いくつか内容が古くなっている部分がありますがご了承ください• 動画作成時は「富岳」がTOP500ランキング世界一位でしたが、2022年6月に二位に後退しました• 「不老」の各ランキング情報ももちろん資料中の順位より低下しています• （各種ランキングは基本的にシステム稼働開始時に測定・登録し、後発のより高性能なシステムが抜き去っていくもの）

8「不老」へのログイン• 実際にログインしてみましょう• 全体的な流れ
1. 配布されたユーザIDとパスワード情報を用いて「HPCポータル」にログインする• https://portal.cc.nagoya-u.ac.jp/
2. SSH公開鍵（秘密鍵と公開鍵のペア）を準備する
3. HPCポータル上でSSH公開鍵を登録する
4. SSHアクセスする• 今回の接続先は flow-cx.cc.nagoya-u.ac.jp• 具体的な手順は以下のいずれかのファイルを参照（手順を飛ばさないこと！）– 不老へのログイン_Windows_MobaXterm.pdf Windows利用者向け、おすすめ– 不老へのログイン_Windows_Cygwin.pdf Windows利用者向け、Cygwinユーザ向け– 不老へのログイン_Windows_Putty.pdf Windows利用者向け、Puttyユーザ向け– 不老へのログイン_MacLinux.pdf MacやLinuxを利用している人向け

9

講習会用アカウント（w49～）でSSH接続する際の詳細手順（補足）• 接続先は flow-cx.cc.nagoya-u.ac.jp• ユーザ名は w49～
• 秘密鍵はid_rsa（MobaXtermではid_rsa.ppk）• MobaXtermを使う場合– 「Remote host」に接続先を記入する– 「Specify username」にw49～を指定する– 「Use private key」に秘密鍵（id_rsa.ppkファイル）を指定する• MacやLinux、Cygwinを使う場合– ssh w49～@flow-cx.cc.nagoya-u.ac.jpというコマンドを実行する– ~/.ssh/id_rsa 以外の秘密鍵を使う場合は -i で明示する
ssh -i 鍵ファイルのパス w49～@flow-cx.cc.nagoya-u.ac.jp

ここをクリックするとファイル選択
画面が開くので、保存した秘密鍵（id_rsa.ppkファイル）を選択する

10

バッチジョブの実行と結果の確認• プログラムの実行を伴わない（環境変数や時刻を表示して待機するだけ）のジョブを実行
し、バッチジョブの実行と結果の確認方法を確認する
• 手順
1.
2.
3.
4.

バッチジョブスクリプトを作成する
ジョブを投入する
ジョブの実行状況を確認する
ジョブの実行結果を確認する

11

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

スパコン上でのプログラム（テキスト）編集• sshアクセス時にも利用できるエディタとしてはEmacsやviが有名だが利用には慣れが必要• 「いわゆる普通のテキストエディタ」的なものとしてはnanoが利用可能（次ページ参照）• 手元で作成したファイルを転送しても良いが、スパコン上でもちょっとしたテキスト編集をしたいケースは多いと思うので、sshアクセス時に使えるエディタにも慣れておくと良い• MobaXtermでemacsを起動しようとするとX転送によりウィンドウが別途開かれてしまう– 便利だが動作は遅め、emacs -nw でCUI版のemacsを利用することを推奨• ファイルの転送にはscp、sftp、MobaXtermの左側のパネルなどを使う– Visual Studio CodeのRemote SSHプラグインなども利用可能– 外部のgitリポジトリなどを経由するなどの方法もある• 主にWindowsユーザ向けの注意– 利用するコンパイラなどによっては文字コードや改行コードが合わずにエラーすることがあるが、見た目では判別できないため注意が必要：文字コードUTF-8、改行コードLFなら安全（なはず）– 最近の「メモ帳」はUTF-8で保存できるが、LFで保存できないため変換が必要– バッチジョブスクリプトは改行コードLFでないとエラーしてしまうため特に注意が必要

13

nanoの使い方• 起動：nanoコマンドで起動– 実行時引数でファイル名を渡すと開いた状態で起動– 例：nano hello.c• Shift+カーソルキーで文字列の選択ができる• 画面の一番下にショートカット
コマンドが書かれている^X Exit → Ctrl+Xで終了※保存が必要そうな場合は確認が出る

M-U Undo→

Alt+U、またはEsc→Uでアンドゥ

14

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

• 変換

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

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

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

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

15

計算ノードの使い方：バッチジョブとインタラクティブジョブ
• Type IIサブシステムの計算ノードは複数のリソースグループ（キュー）にわかれており、状況に応じて使いわけることでプログラムを効率良く実行できる– 空いているキューを積極的に利用して欲しい– SSDを使わないならcx-*とcxgfs-*に差はない– インタラクティブジョブは専用キュー（cx-interactive, cxgfs-interactive）• バッチジョブ– ジョブスクリプトをエディタで作成する– pjsubでジョブを投入し、pjstat, pjstat2で実行状況を確認する• インタラクティブジョブ（会話型バッチ処理）– バッチジョブ実行の仕組みを使って計算ノードにログインする仕組みであり、計算ノード上でイン
タラクティブにコマンドを実行できる– pjsub --interact -L rscgrp=cx-interactive– ログインしているだけでポイントを消費する– 主に環境構築・動作確認のための利用を想定、コンテナの構築にも有用

16

バッチジョブスクリプトの記述例#!/bin/bash -x#PJM -L rscgrp=cx-workshop#PJM -L node=1#PJM -L elapse=1:00#PJM -j#PJM -S

←
←
←
←
←
←

-xを付けておくと、実行するコマンドとその結果の両方が出力されるリソースグループ名：今回はcx-workshop
利用ノード数：1
必須ではないが重要
利用時間制限：1分
標準出力と標準エラー出力を統合
お好みで
統計情報を出力date
env
sleep 30
date

←
←
←
←

時刻を表示
環境変数を表示（特に意味はない）
30秒待つ
時刻を再度表示• 赤字部分はバッチジョブの設定に関する情報• 黒字部分は一般的なbashスクリプト処理

作業手順
1. cd /data/group1/${USER}で作業場所を変更（この通り入力してください。${USER}は自動的にw49~に置き換わって解釈されます。）
2. job_test.sh として上記ファイルを作成→ pjsub job_test.sh でジョブ投入

17

バッチジョブ実行と確認の例
赤文字が実行したコマンド、緑文字が実行結果として出力されたもの←pjsubでジョブを投入、ジョブIDが通知される[a49979a@flow-cx02 ~]$ pjsub job_test.sh（設定を間違えた場合などはpjdelにこのIDを指定して削除する）[INFO] PJM 0000 pjsub Job 75544 submitted.
[a49979a@flow-cx02 ~]$←pjstat2でジョブの状況を確認（pjstatでも確認可能詳細は利用手引書を参照）[a49979a@flow-cx02 ~]$ pjstat2
JOB_ID
JOB_NAMESTATUS GROUPRSCGROUP
START_DATE
ELAPSE
NODE
75544
job_test
QUEUED centercx-single--/-- --:--:-00:00:00
1[a49979a@flow-cx02 ~]$↑QUEUEDは実行されていない状態、開始時刻が決まるとSTART_DATEに表示される[a49979a@flow-cx02 ~]$ pjstat2
JOB_ID
JOB_NAMESTATUS GROUPRSCGROUP
START_DATE
ELAPSE
NODE
75544
job_test
RUNNING centercx-single(01/05 14:15)
00:00:00
1[a49979a@flow-cx02 ~]$↑RUNNINGは実行中、START_DATEが括弧付きの場合は開始準備中[a49979a@flow-cx02 ~]$ pjstat2
JOB_ID
JOB_NAMESTATUS GROUPRSCGROUP
START_DATE
ELAPSE
NODE
75544
job_test
RUNNING centercx-single
01/05 14:15:29
00:00:021[a49979a@flow-cx02 ~]$[a49979a@flow-cx02 ~]$ pjstat2←実行中/実行開始待ちのジョブがなくなった状態
No unfinished job found.

18

実行結果ファイルの確認[a49979a@flow-cx02 ~]$ cat job_test.sh.76342.out+ dateThu Jan 7 17:46:31 JST 2021
+ env
PJM_DPREFIX=#PJM
PJM_JOBNAME=job_test.sh←

catコマンド（またはlessなど）でファイルを表示←

dateコマンドが実行された時刻←

ここから環境変数←

最後の時刻表示中略}_=/usr/bin/env+ sleep 30+ dateThu Jan 7 17:46:41 JST 2021

補足
• 結果ファイル名は「ジョブスクリプト名.ジョブID.out」、#PJM -jを付けなかった場合は標準エラー出力が別途「ジョブスクリプト名.ジョブID.err」に出力される ※オプション指定によって変更可能• + で始まる行はコマンドの実行結果ではなくコマンドそのもの（bashに-xオプションを付けると表示される）• Type Iサブシステムでbash -xを指定すると冒頭に環境設定のための処理などがたくさん表示される

19

リソースグループについての補足（消費ポイントについての注意点）• 本講習会ではcx-workshopを指定しているが、実際に利用する際にはリソースグループとGPU数の関係に注意が必要• cx-shareだけは常に1/4ノード、1GPUのみが確保される• cx-share以外は常に1ノード単位で確保され、4GPUが利用される• Type IIサブシステムの消費ポイント計算式は利用GPU数×経過時間[s] × 0.007• cx-share以外では常に利用GPU数＝ノード数×4であることに注意すること– cx-singleのsingleは1ノードの意味でのsingle• GPU数をベースにするかノード数をベースにするかの問題であり、他のサブシステムとの一貫性からノード数をベースにした20

インタラクティブジョブの実行確認

通常はcx-interactiveを使うが、講習会用にcx-workshopでもインタラクティブジョブが実行可能にしてあるのでcxworkshopを指定してください。30分を超えるインタラクティ
ブジョブを実行したい場合はcx-interactiveを利用、特に60分を超える場合はelapse指定を追加してください。

[a49979a@flow-cx02 ~]$ hostname
flow-cx04← ログインノードflow-cx04にいる[a49979a@flow-cx02 ~]$ pjsub --interact -L rscgrp=cx-workshop[INFO] PJM 0000 pjsub Job 76345 submitted.
↑ インタラクティブジョブ用のオプションを付けてpjsub[INFO] PJM 0081 .connected.
[INFO] PJM 0082 pjsub Interactive job 76345 started.
[a49979a@cx120 ~]$ hostname
cx120← 計算ノードcx120にいる[a49979a@cx120 ~]$ここで何かコマンドを入力すれば計算ノード上で動作する[a49979a@cx120 ~]$ exit
← Ctrl+Dやexitコマンドでログアウト
logout[INFO] PJM 0083 pjsub Interactive job 76345 completed.
[a49979a@cx120 ~]$← インタラクティブジョブ終了

インタラクティブジョブ実行時はpjsubに--interactが必要。
-L オプションはジョブスクリプトのPJM行と同じ役割。インタラクティブではないバッチジョブの投入時にオプションを与えることも可能。
• elapse指定をしていない状態と同じである点に注意。もっと長い時間実行したい場合は-L elapse
で指定する。cx-interactiveは24時間まで延長可能。cx-workshopは標準でも最長でも30分まで。
例：pjsub --interact -L rscgrp=cx-interactive -L elapse=2:00:00
または pjsub --interact -L rscgrp=cx-interactive,elapse=2:00:00
実際にやってみましょう• 計算ノードに空きがないと実行できない。表示されるメッセージをよく見ておくこと。
•
•

21

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

22

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

2並列計算
CPU 1

CPU

• 並列？

並行・平行？

CPU 2
time

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

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

23

なぜ並列計算を行うのか？
• 短時間で終わらせたい計算があるから、計算機（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の掲載内容）

24

最近のマルチコア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(名大旧システム)に搭載

Fujitsu 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回行える性能※「富岳」は2.0 GHz
1K=1000, 1M=1000K, 1G=1000M, 1T=1000G, 1P=1000T, 1Exa=1000P
複数のコアを搭載するのが当たり前、コア数は徐々に増加してきている
実際のプログラムの性能はコア数や（理論）演算性能だけでは決まらない点には注意が必要

25

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

動作周波数

理論演算性能(FP64)* 理論メモリバンド幅

発表年月**

Tesla P100

3584 (64*56)

1328 ~ 1480 MHz

5.3 TFLOPS

720 GB/s

2016.04 製品発表

Tesla V100

5120 (64*80)~ 1530 MHz

7.8 TFLOPS

900 GB/s

2017.06 製品発表

NVIDIA A100

6912 (64*108)~ 1410 MHz

9.7 TFLOPS

1555 GB/s

2020.05 製品発表

NVIDIA H100

16,896 (128*132)

未定

30 TFLOPS

3,000 GB/s

2022.03 製品発表

• GPUはCPUと比べて……
– 大量の計算コア、低い動作周波数、高い理論演算性能、高いメモリバンド幅– 数字以外の違いについては後述* 倍精度(64bit)浮動小数点演算のFLOPS値。32bitなど異なるデータ型に対する演算性能に注目することもある。
** イベント等で正式に発表した年月。一部のユーザに先行提供していたり、一般向けに販売開始したのがいつなのかわからないことも多い。

26

スーパーコンピュータ（スパコン）と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）
Instinct

組み込み、車載など

Tegra, Xavier, Jetson

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

27

TOP500 2022年6月版• 連立一次方程式の数値解法（巨大なLU分解）• 初登場のFrontierが首位（ついに1EFLOPS到達）、
理研RCCSの「富岳」はついに2位に後退• 近年はアメリカと中国が上位の多くを占める• 理論性能に対する性能比はあまり高くない– 昔は高かった。現実のアプリではもっと下がることが多い。

• 消費電力も大きい– 大きなものは発電所を1つ占有するレベル– 10MW：一般家庭300世帯程度• 下線部はNVIDIA社とAMD社のGPU– TOP10内の7システムを含む多数のスパコンがGPUを主要な計算装置として活用– 特に今回はAMDの最新GPU Instinct MI250Xが活躍• 基本的にシステム稼働直後に測定した値を登録、後発
画像引用元：https://www.top500.org/に抜かれていくシステム

28

ちなみに
• 電力効率（性能を消費電力で割った値）を競う
Green500の上位はGPUの独壇場に近い– 10システム中9システムがGPUスパコン– （上位4システムは同じHPE Crayのマシンだけど……）– 5位のMN-Coreは独自のアクセラレータ• 「富岳」も超大規模システムとしては非常に高効率だがGPUがそれを上回っている

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

29

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

計算
コア

計算
コア

共有メモリ

計算
コア

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

メモリ

メモリ

メモリ

計算
コア

計算
コア

計算
コア

計算
コア

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

30

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

31

並列計算を行う方法（言語など）の例• コンパイラによる自動並列化：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内並列化• デバイス内では共有メモリモデル、外部とのやりとりは分散メ
モリモデル• ノード間並列化• 分散メモリモデル

32

並列計算の基本的なイメージ：ループ並列化（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はこちらのイメージ

33

並列計算の基本的なイメージ：ループ並列化（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はこちらのイメージ

34

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

35

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

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に担当させるかを
しっかり考える必要がある

36

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

37

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

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

39

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

GPU

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

40

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• 名古屋大学情報基盤センターに設置されているスーパーコンピュータ「不老」のType IIサブシス
テムに搭載されているGPU– ノードあたり4枚、サブシステム全体で884枚搭載– 資料に出てくる様々な数字を全て覚える必要はない（最適化を行う際には気にする必要がある）

41

全体構成

42

全体構成
PCIe
GigaThread Engine （GPUの動作を制御する部分）

GPC(GPU Processing Cluster)

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

NVLinkNVLinkNVLinkNVLinkNVLinkNVLink

HBM2

HBM2

メ
モ
リ
コ
ン
ト
ロ
ー
ラ
群

HBM2

HBM2

メ
モ
リ
コ
ン
ト
ロ
ー
ラ
群

GPCの中身
GPC

TPC
SM

SM

• GPC: GPU Processing Cluster 43• 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などと接続※ Tesla V100では84SMs中の80SMsのみが有効。4個までは不良品が含まれてもOKとすることで歩留まりを向上させている。

44

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

45

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

46

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

47

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

参考：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利用」にはあまり向かない49

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で確保したメモリのポインタを引き渡すタイプだと簡単）

50

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

51

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

52

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

53

OpenACCに関する情報、ツール、etc.
• 規格などの情報– https://www.openacc.org/• 仕様書の他、学習用の資料へのリンクなどもあり– 例 https://github.com/OpenACC/openacc-training-materials• 主な対応コンパイラ– 商用：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サイトに掲載されていた資料のアーカイブ• 言語仕様が大きく変わってしまわない限りは今後も有用だろう

54

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自体はユーザ権限でもインストールできるため、最新版をはやく使いたい場合は自分で導入してください）55

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

計算ノード
CPU

メインメモリ（DDR*）

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

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

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

56

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

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

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

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

57$ nvaccelinfo

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

11040
NVIDIA UNIX x86_64 Kernel Module 470.57.02 Tue Jul 13 16:14:05 UTC 2021

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

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

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

単純な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

59

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上のメモリを解放

60「不老」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/22.2 をloadするとコンパイラ等が利用可能になる– module load hpc_sdk/22.2– /22.2 は省略可能。省略した場合、システムのデフォルト設定のものがloadされる。
– バージョンを省略しておけば今後さらに新しいものをインストールして提供する場合は新しいものが自動的に読み込まれるようになる。ただしコンパイル時と実行時でバージョンが異なると問題が起きることがある点には注意が必要。今回の講習会ではデフォルトの22.2を使用する。

•

コンパイル例• nvcまたはnvfortranでコンパイル•

61-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

62

OpenACCプログラム向けのジョブスクリプトの書き方#!/bin/bash#PJM -L rscgrp=cx-workshop#PJM -L node=1#PJM -L elapse=1:00#PJM -j#PJM -S←←←←←

リソースグループ名：cx-workshop
利用ノード数：1
利用時間制限：1分
標準出力と標準エラー出力を統合
統計情報を出力

module load hpc_sdk/22.2←

module loadでhpc_sdkを使えるようにする

export NVCOMPILER_ACC_TIME=1./v1_c_acc←←

GPUの実行時間内訳を確認するための設定（必須ではない）
プログラムを実行

必須ではないが重要
お好みで※module loadはコンパイル時のみではなく実行時にも必要（ライブ
ラリなどを参照するため）
job_vector1.sh として作成→ pjsub job_vector1.sh でジョブ投入• 赤字部分はバッチジョブの設定に関する情報• 黒字部分は一般的なbash
スクリプト処理63

リソースグループについての補足（消費ポイントについての注意点）• 本講習会ではcx-workshopを指定しているが、実際に利用する際にはリソースグループとGPU数の関係に注意が必要• cx-shareだけは常に1/4ノード、1GPUのみが確保される• cx-share以外は常に1ノード単位で確保され、4GPUが利用される• Type IIサブシステムの消費ポイント計算式は利用GPU数×経過時間[s] × 0.007• cx-share以外では常に利用GPU数＝ノード数×4であることに注意すること• cx-singleのsingleは1ノードの意味でのsingle– GPU数をベースにするかノード数をベースにするかの問題、他のサブシステムとの一貫性からノード数をベースにした64

バッチジョブ実行と確認の例
赤文字が実行したコマンド、緑文字が実行結果として出力されたもの$ pjsub job_vector1.sh←pjsubでジョブを投入、ジョブIDが通知される[INFO] PJM 0000 pjsub Job 75544 submitted.
（設定を間違えた場合などはpjdelにこのIDを指定して削除する）$←pjstat2でジョブの状況を確認（pjstatでも確認可能詳細は利用手引書を参照）$ pjstat2
JOB_ID
JOB_NAMESTATUS GROUPRSCGROUP
START_DATE
ELAPSE
NODE
75544
job_vector QUEUED centercx-single--/-- --:--:-00:00:00
1$↑QUEUEDは実行されていない状態、開始時刻が決まるとSTART_DATEに表示される$ pjstat2
JOB_ID
JOB_NAMESTATUS GROUPRSCGROUP
START_DATE
ELAPSE
NODE
75544
job_vector RUNNING centercx-single(01/05 14:15)
00:00:00
1$↑RUNNINGは実行中、START_DATEが括弧付きの場合は開始準備中$ pjstat2
JOB_ID
JOB_NAMESTATUS GROUPRSCGROUP
START_DATE
ELAPSE
NODE
75544
job_vector RUNNING centercx-single
01/05 14:15:29
00:00:02
1$$ pjstat2←実行中/実行開始待ちのジョブがなくなった状態
No unfinished job found.

65

実行例
• 正しく実行されている限り、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

←元データ
←計算結果

66

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/20211004/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

67

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カーネル起動情報とデータ転送情報が出力される

68

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
Implicit wait
Implicit wait

file=/path/to/source.c function=main line=16 device=0 threadid=1
file=/path/to/source.c function=main line=16 device=0 threadid=1
file=/path/to/source.c function=main line=18 device=0 threadid=1

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

69

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

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

71

並列化対象範囲の指定：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72

再掲：単純な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で閉じる必要がある73

コンパイル時のメッセージ再確認• 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については後述

74

手動による適用：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)」が自動的に挿入されていたと思えば良い75

手動による適用：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

76

参考：指示文の継続行• 指示文行を次の行に継続させることも可能– 長くなってしまったときなどに• 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が必要

77

データ転送に関する指示節• 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}※データを使い回す、該当するものが無ければ実行時エラー

78

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

メインメモリ（DDR*）

GPU

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

79

データ転送速度• メインメモリやデバイスメモリの転送速度に対してCPU-GPU間のデータ転送速度はずっと
低速、頻繁な通信は避けたい– CPU - メインメモリ• DDR4、100GB/s/socket程度（STREAM Triad実測）– より高速なCPUでも理論性能で200GB/s/socket程度※ A64FXなどHBM系を搭載したCPUであれば、同様のメモリを搭載した
GPUと同程度の速度が出る– GPU - デバイスメモリ• HBM2、P100で550GB/s程度、V100では800GB/s弱（STREAM Triad実測）– CPU - GPU
他より低速～数百 GB/s• PCI Express Gen.3 x16
CPU
メインメモリ（DDR*）– 理論性能でも最大16GB/s（×双方向）» Gen.4 でさらに2倍、それでも32GB/s（×双方向）– GPU - GPU• NVLink、 20(Pascal) or 25(Volta)GB/s×双方向– 1GPUあたり4(Pascal) or 6(Volta)本のNVLink、GPU数によって最適なつなぎ方が変わる– Power系CPUではCPU-GPU間でもNVLinkが使える

PCIe、～数十 GB/s
GPU
～数 TB/s

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

80

データ転送範囲の指定• 配列全体ではなく一部のみを送受信することも可能• 注意：C/C++とFortranでは部分配列の指定方法が異なる– C/C++：先頭と長さを指定する#pragma acc kernels copy(A[head:length])– Fortran：開始点と終了点を指定する!$acc kernels copy(A(begin:end))– A[:N], A(:N) のような省略表記も可能（先頭からN要素が送られる）• 単純なプログラムであればコンパイラが勝手に正しい長さを認識してくれることが多いが、長さの認識が間違っているためにプログラムがエラーすることも少なくない。
• コンパイラの出力を確認し、認識が正しくない場合は手動で指定するよう気を付けること。

81

繰り返しデータ転送を行う場合の問題• 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)

82

繰り返しデータ転送を行う場合の問題• 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)

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

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

83

データ転送のみを行う指示文• 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化範囲の前でとにかく全部送信したいとき？
– 複数ソースコードにプログラムが分割されているとき？

• 配列長（転送範囲）の指定方法はkernelsの場合と同様

84

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;}}

(vector4.c)

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

(vector4.f90)

デバイスメモリの生成繰り返し

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

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

デバイスメモリの生成繰り返し

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

実習• vector3.cとvector4.c（もしくはvector3.f90とvector4.f90）を実行して動作を比較してみよう– NVCOMPILER_ACC_TIMEやNVCOMPILER_ACC_NOTIFY=2を設定して実行し、出力を比較してみよう（通信回数や時間を比較してみよう）• プログラムの動作を理解できたら、配列v1の値を2倍するのではなく、結果を配列v2に格納したり、配列v1とv2の加算結果を配列v3に格納するなど、少し異なる動作をするプログラ
ムも作成してみよう– 最適と思われる指示文の入れ方は？
– コンパイラは最適な指示文を挿入してくれるだろうか？

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

実習：実行例
Accelerator Kernel Timing data/home/center/a49979a/share/20211004/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/20211004/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/20211004/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/20211004/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

88

動的な配列のデータ送受信• 部分転送時の範囲指定には変数を用いても良い• 動的に確保した配列などコンパイル時に配列の長さが分からないものは長さを明示する必要がある– 間接参照をしているときなどに明示的な指定が必要となりやすい(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)) を加えると正しく動作する89

動的な配列のデータ送受信• 部分転送時の範囲指定には変数を用いても良い• 動的に確保した配列などコンパイル時に配列の長さが分からないものは長さを明示する必要がある– 間接参照をしているときなどに明示的な指定が必要となりやすい(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)) を加えると正しく動作する90

コンパイル時の出力の比較
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についてはこのままでも（結果的に）問題ない91

実行結果の比較
配列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

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のv2結果 ： 2.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
allocate1.f90のv2結果 ： 2.00 4.00 6.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• allocate1.cでは配列長をうまく認識できていないためか計算結果がおかしい。
• どの時点で問題が起きているのか詳細を確認することはできるだろうか？

92

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/20211004/allocate1.c function=main line=21 device=0 threadid=1 variable=index bytes=40
upload CUDA data file=/home/center/a49979a/share/20211004/allocate1.c function=main line=21 device=0 threadid=1 variable=v1 bytes=8
download CUDA data file=/home/center/a49979a/share/20211004/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/20211004/allocate1.f90 function=main line=14 device=0 threadid=1 variable=v1 bytes=80
upload CUDA data file=/home/center/a49979a/share/20211004/allocate1.f90 function=main line=14 device=0 threadid=1 variable=index bytes=40
download CUDA data file=/home/center/a49979a/share/20211004/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/20211004/allocate2.c function=main line=21 device=0 threadid=1 variable=v1 bytes=80
upload CUDA data file=/home/center/a49979a/share/20211004/allocate2.c function=main line=21 device=0 threadid=1 variable=index bytes=40
download CUDA data file=/home/center/a49979a/share/20211004/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/20211004/allocate2.f90 function=main line=14 device=0 threadid=1 variable=v1 bytes=80
upload CUDA data file=/home/center/a49979a/share/20211004/allocate2.f90 function=main line=14 device=0 threadid=1 variable=index bytes=40
download CUDA data file=/home/center/a49979a/share/20211004/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

93

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

94

ただし……
• C言語ではよく二次元配列を2段階のmallocで確保するが、これをやるとメモリの連続性が失われてしまうa = (double**)malloc(sizeof(double*)*n);for(i=0;i<n;i++){a[i] = (double*)malloc(sizeof(double)*n);}• OpenACCでどう扱われるのか（気を利かせて転送してくれるのか）よくわからないため、本講習会ではC言語で二次元配列（行列）を動的確保する際は一次元配列で実装している– 送ってくれるようなのだが、必ず送れるのか不安が……
a = (double*)malloc(sizeof(double)*n*n);for(i=0; i<n; i++){for(j=0; j<n; j++){a[i*n+j] = (double)(i+1) + (double)(j+1)*0.01;}}
a[i][j] に相当※ Fortranは一命令で二次元配列を
allocateできるから大丈夫だろうと判断
し、二次元配列をそのまま用いている

95

データの更新• 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)

96

データの更新• 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)

97

データの更新：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)

98

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

99

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

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

101

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

double a[N*N], b[N*N], c[N*N];
初期化等は省略
27 #pragma acc kernels
28
for(i=0; i<N; i++){29
for(j=0; j<N; j++){30
for(k=0; k<N; k++){31
c[i*N+j] += a[i*N+k] * b[k*N+j];
32
}
33
}
34
}

double precision :: a(N,N), b(N,N), c(N,N)
初期化等は省略
29 !$acc kernels
30 do i=1, N
31
do j=1, N
32
do k=1, N
33
c(j,i) = c(j,i) + a(k,i) * b(j,k)
34
enddo
35
enddo
36 enddo
37 !$acc end kernels

n=10で実行して計算時間を比較してみたCtime(us): 62

28: compute region reached 1 time
30: kernel launched 1 time
grid: [1] block: [32]device time(us): total=13 max=13 min=13 avg=13
elapsed time(us): total=549 max=549 min=549 avg=549

Fortrantime(us): 45
29: compute region reached 1 time
32: kernel launched 1 time
grid: [1] block: [128]device time(us): total=4 max=4 min=4 avg=4
elapsed time(us): total=491 max=491 min=491 avg=491• Cがやや遅いが、もの凄く遅いという程ではない（blockの数を見ると妥当性があるのだが、そのあたりは後ほど）• （そもそも問題サイズが小さす
ぎる、もっと大きなサイズでチェックすべきではある）102

三重ループ構造による単純な行列積プログラムの例：動的確保の場合•• C (matmul1a.c)• Fortran (matmul1a.f90)

double *a=NULL, *b=NULL, *c=NULL;
mallocなどは省略
34 #pragma acc kernels
35
for(i=0; i<N; i++){36
for(j=0; j<N; j++){37
for(k=0; k<N; k++){38
c[i*N+j] += a[i*N+k] * b[k*N+j];
39
}
40
}
41
}

double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocateなどは省略
36 !$acc kernels
37 do i=1, N38
do j=1, N
39
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

実行して計算時間を比較しようとすると、そもそもC言語版は実行時にエラーしてしまう$ ./mm1a_c_acc
N = 10out = 1
Failing in Thread:1
call to cuStreamSynchronize returned error 700:
Illegal address during kernel execution$ ./mm1a_f_accn =10
result:58.96
59.06
114.06114.26169.16 169.46
中略
444.66445.46499.76500.66554.86555.86

これらはエラーではないので無視して良い

Warning: ieee_inexact is signaling
FORTRAN STOP

59.17
114.47
169.7759.27
114.67170.0759.38
114.88170.3859.48
115.08170.6859.59
115.29170.9959.69
115.49171.2959.80
115.70171.6059.91
115.91
171.91446.27501.57556.87447.07502.47557.87447.88503.38558.88448.68504.28559.88449.49505.19560.89450.29506.09561.89451.10507.00562.90451.91507.90563.90103

コンパイラの出力するメッセージの確認• C (matmul1a.c)• Fortran (matmul1a.f90)$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o mm1a_c_acc matmul1a.c$ nvfortran -g -Minfo=accel -acc -gpu=cc70 -tp=host -o mm1a_f_acc matmul1a.f90
main:
main:
35, Accelerator restriction: size of the GPU copy of c,a is unknown
36, Generating implicit copyin(b(1:N,1:N)) [if not already present]Complex loop carried dependence of b->,a-> prevents parallelization
Generating implicit copy(c(1:N,1:N)) [if not already present]
Loop carried dependence of c-> prevents parallelization
Generating implicit copyin(a(1:N,1:N)) [if not already present]
Loop carried backward dependence of c-> prevents vectorization
37, Loop is parallelizable
Complex loop carried dependence of c-> prevents parallelization
38, Loop is parallelizable
Accelerator serial kernel generated
39, Complex loop carried dependence of c prevents parallelization
Generating NVIDIA GPU code
Loop carried dependence of c prevents parallelization
35, #pragma acc loop seq
Loop carried backward dependence of c prevents vectorization
seq = sequential
36, #pragma acc loop seq
Inner sequential loop scheduled on accelerator
ループが逐次実行される
37, #pragma acc loop seq
Generating NVIDIA GPU code
ことを意味する
35, Generating implicit copyin(a[:]) [if not already present]
37, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
Generating implicit copy(c[:]) [if not already present]
38, ! blockidx%x threadidx%x auto-collapsed
Loop carried backward dependence of c-> prevents vectorization
39, !$acc loop seq
36, Complex loop carried dependence of b->,a->,c-> prevents parallelization
39, Complex loop carried dependence of c prevents parallelizationz
37, Accelerator restriction: size of the GPU copy of b is unknown
Complex loop carried dependence of b->,a->,c-> prevents parallelization
Loop carried dependence due to exposed use of c[(N)*i1+i2] prevents parallelization
Generating implicit copyin(b[:]) [if not already present]• Fortran言語版は配列のサイズが正しく認識できており、
Complex loop carried dependence of c-> prevents parallelization••

C言語版は配列のサイズが認識できていない（データの転送以前にそもそもメモリの確保が正しく行えていないと考えられる）カーネルの並列化もできていない（できても困るが）•データ転送も問題なく行われているように見える
カーネルの並列化もできているようだ➢ 3重ループの外側2つは並列化され（auto-collapsedについては後述）、最内ループは逐次実行。cへの書き込みは並列化す
ると競合するため、妥当な判断。

104

動的確保＋copy指示文の追加• C (matmul2a.c)• Fortran (matmul2a.f90)

34 #pragma acc kernels copy(c[N*N]) copyin(a[N*N],b[N*N])
35
for(i=0; i<N; i++){36
for(j=0; j<N; j++){37
for(k=0; k<N; k++){38
c[i*N+j] += a[i*N+k] * b[k*N+j];
39
}
40
}
41
}

36 !$acc kernels copy(c(N,N)) copyin(a(N,N),b(N,N))37
do i=1, N
38
do j=1, N
39
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$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o mm2a_c_acc matmul2a.c
main:
35, Generating copyin(b[:N*N]) [if not already present]
Generating copy(c[:N*N]) [if not already present]
Generating copyin(a[:N*N]) [if not already present]Complex loop carried dependence of b->,a-> prevents parallelization
Loop carried dependence of c-> prevents parallelization
Loop carried backward dependence of c-> prevents vectorization
Complex loop carried dependence of c-> prevents parallelization
Accelerator serial kernel generated
Generating NVIDIA GPU code
35, #pragma acc loop seq
36, #pragma acc loop seq
37, #pragma acc loop seq
35, Loop carried dependence of c-> prevents parallelization
Loop carried backward dependence of c-> prevents vectorization
Complex loop carried dependence of c-> prevents parallelization
36, Complex loop carried dependence of b->,a->,c-> prevents parallelization
37, Complex loop carried dependence of b->,a->,c-> prevents parallelization
Loop carried dependence due to exposed use of c[(N)*i1+i2] prevents parallelization••$ nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o mm2a_f_acc matmul2a.f90
main:
36, Generating copy(c(N,N)) [if not already present]
Generating copyin(b(N,N),a(N,N)) [if not already present]
37, Loop is parallelizable
38, Loop is parallelizable
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 NVIDIA GPU code
37, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
38,! blockidx%x threadidx%x auto-collapsed
39, !$acc loop seq
39, Loop carried dependence of c prevents parallelization
Loop carried backward dependence of c prevents vectorization

C言語版は配列サイズの認識問題が解決
カーネルの並列化はできていないまま•

Fortran版は特に変わらず

105

動的確保＋copy指示文の追加：実行時間の確認• C (matmul2a.c)• Fortran (matmul2a.f90)$ NVCOMPILER_ACC_TIME=1 ./mm2a_c_acc
N = 10out = 1
result:省略$ NVCOMPILER_ACC_TIME=1 ./mm2a_f_accn =10
result:省略

Accelerator Kernel Timing data
/data/group1/a49979a/share/20220627/matmul2a.cmain NVIDIA devicenum=0time(us): 108
35: compute region reached 1 time
35: kernel launched 1 time
grid: [1] block: [1]device time(us): total=65 max=65 min=65 avg=65
elapsed time(us): total=530 max=530 min=530 avg=530
35: data region reached 2 times
35: data copyin transfers: 3
device time(us): total=31 max=16 min=4 avg=10
41: data copyout transfers: 1
device time(us): total=12 max=12 min=12 avg=12

Warning: ieee_inexact is signaling
FORTRAN STOP
Accelerator Kernel Timing data
/data/group1/a49979a/share/20220627/matmul2a.f90
main NVIDIA devicenum=0time(us): 62
36: compute region reached 1 time
39: kernel launched 1 time
grid: [1] block: [128]device time(us): total=6 max=6 min=6 avg=6
elapsed time(us): total=510 max=510 min=510 avg=510
36: data region reached 2 times
36: data copyin transfers: 3
device time(us): total=40 max=19 min=10 avg=13
44: data copyout transfers: 1
device time(us): total=16 max=16 min=16 avg=16• C版はFortran版と比べてかなり遅い• GPUカーネルの並列化ができていないので当然ではある106

loop指示文の追加• C (matmul3a.c)• Fortran (matmul3a.f90)

34 #pragma acc kernels copy(c[N*N]) copyin(a[N*N],b[N*N])
35 #pragma acc loop independent36
for(i=0; i<N; i++){
並列実行する37 #pragma acc loop independent38
for(j=0; j<N; j++){39 #pragma acc loop seq
並列実行しない40
for(k=0; k<N; k++){41c[i*N+j] += a[i*N+k] * b[k*N+j];
42
}
43
}
44
}

36 !$acc kernels copy(c(N,N)) copyin(a(N,N),b(N,N))
37 !$acc loop independent38
do i=1, N
並列実行する39 !$acc loop independent40
do j=1, N
41 !$acc loop seq
並列実行しない42
do k=1, N
43c(j,i) = c(j,i) + a(k,i) * b(j,k)
44
enddo
45
enddo
46
enddo
47 !$acc end kernels•

データ転送もGPUによる計算も問題なさそう$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o mm3a_c_acc matmul3a.c
main:
36, Generating copyin(b[:N*N]) [if not already present]
Generating copy(c[:N*N]) [if not already present]
Generating copyin(a[:N*N]) [if not already present]
Loop is parallelizable
38, Loop is parallelizable
40, Generating NVIDIA GPU code
36, #pragma acc loop gang, vector(128) collapse(2)/* blockIdx.x threadIdx.x */38, /* blockIdx.x threadIdx.x auto-collapsed */40, #pragma acc loop seq$ nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o mm3a_f_acc matmul3a.f90
main:
36, Generating copyin(b(N,N)) [if not already present]
Generating copy(c(N,N)) [if not already present]
Generating copyin(a(N,N)) [if not already present]
38, Loop is parallelizable
40, Loop is parallelizable
42, Generating NVIDIA GPU code
38, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
40,! blockidx%x threadidx%x auto-collapsed
42, !$acc loop seq

107

loop指示文の追加（実行時間の比較）• C (matmul3a.c)• Fortran (matmul3a.f90)

34 #pragma acc kernels copy(c[N*N]) copyin(a[N*N],b[N*N])
35 #pragma acc loop independent36
for(i=0; i<N; i++){
並列実行する37 #pragma acc loop independent38
for(j=0; j<N; j++){39 #pragma acc loop seq
並列実行しない40
for(k=0; k<N; k++){41c[i*N+j] += a[i*N+k] * b[k*N+j];
42
}
43
}
44
}

36 !$acc kernels copy(c(N,N)) copyin(a(N,N),b(N,N))
37 !$acc loop independent38
do i=1, N
並列実行する39 !$acc loop independent40
do j=1, N
41 !$acc loop seq
並列実行しない42
do k=1, N
43c(j,i) = c(j,i) + a(k,i) * b(j,k)
44
enddo
45
enddo
46
enddo
47 !$acc end kernels

n=10で実行してみた→C版が劇的に高速化、Fortran版とほぼ同じ性能

time(us): 57
36: compute region reached 1 time
40: kernel launched 1 time
grid: [1] block: [128]device time(us): total=6 max=6 min=6 avg=6
elapsed time(us): total=544 max=544 min=544 avg=544•••

time(us): 59
36: compute region reached 1 time
42: kernel launched 1 time
grid: [1] block: [128]device time(us): total=6 max=6 min=6 avg=6
elapsed time(us): total=559 max=559 min=559 avg=559

一般的に、Fortranプログラムの方がコンパイラによる並列化判断が適切に働きやすい
C/C++はポインタ参照の都合で不具合が起きないよう保守的な判断がされやすい
loop指示文で指定すればコンパイラの判断を上書きできる（妥当な性能を引き出せる）

108

collapse版• C (matmul4a.c)• Fortran (matmul4a.f90)

34 #pragma acc kernels copy(c[N*N]) copyin(a[N*N],b[N*N])
35 #pragma acc loop independent collapse(2)36
for(i=0; i<N; i++){37
for(j=0; j<N; j++){38 #pragma acc loop seq
39
for(k=0; k<N; k++){40c[i*N+j] += a[i*N+k] * b[k*N+j];
41
}
42
}
43
} • ループを融合してから並列化する

36 !$acc kernels copy(c(N,N)) copyin(a(N,N),b(N,N))
37 !$acc loop independent collapse(2)38do i=1, N
39
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••

短いループがネストしている際に有用
実行時間が短くなるかはループの構造と長さ次第$ nvc -Minfo=accel -acc -gpu=cc70 -tp=host -o mm4a_c_acc matmul4a.c
main:
36, Generating copyin(b[:N*N]) [if not already present]
Generating copy(c[:N*N]) [if not already present]
Generating copyin(a[:N*N]) [if not already present]
Loop is parallelizable
37, Loop is parallelizable
39, Generating NVIDIA GPU code
36, #pragma acc loop gang, vector(128) collapse(2)/* blockIdx.x threadIdx.x */37,/* blockIdx.x threadIdx.x collapsed */39, #pragma acc loop seq$ nvfortran -Minfo=accel -acc -gpu=cc70 -tp=host -o mm4a_f_acc matmul4a.f90
main:
36, Generating copyin(b(N,N)) [if not already present]
Generating copy(c(N,N)) [if not already present]
Generating copyin(a(N,N)) [if not already present]
38, Loop is parallelizable
39, Loop is parallelizable
41, Generating NVIDIA GPU code
38, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
39,! blockidx%x threadidx%x collapsed
41, !$acc loop seq

109

collapseのイメージ#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 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と同様の効果

110kernels 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, N
v2(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, N
v2(i) = v1(i) * 2.0d0
enddo!$acc end kernels loop

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

111

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

112

並列実行形状の指定• ここまで、ループをどのように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)113

gang, worker, vectorとblockIdx, threadIdx• OpenACCではハードウェアに階層的な並列性があることを想定しており、上位階層から順にgang, worker, vectorとなっている• これらをどのように組み合わせて実行するかを指定できる• CUDAの並列実行モデルが階層的になっているため、それにあわせて言語設計された、と考えるのが適切– CUDAはgrid, threadblock, threadの三階層構造• 最近はThread Block Clustersという概念もある– OpenACCのおおまかな並列実行モデルの対応付け
gangvectorblockId 0
threadId 0,1,2…

blockId 1
threadId 0,1,2…

blockId 2
threadId 0,1,2…

114

ハードウェアとのおおまかな対応付け
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以降で
厳密化

115

演習：単純な行列積のコンパイル確認と性能確認• 配布している行列積プログラムをコンパイルして実行する– コンパイルメッセージを見てデータの転送と並列化の状況を確認する– 問題サイズを大きくしてみて実行時間を確認する• ソースコード
matmul1

3重ループ行列積

matmul1a1の配列動的確保版

matmul2

copy指示節あり

matmul2a2の配列動的確保版

matmul3

copy指示節＋loop指示文あり

matmul3a3の配列動的確保版

matmul4

copy指示節＋loop指示文＋collapse指示節あり

matmul4a4の配列動的確保版

116

演習：行列積の繰り返し実行• 行列積を単純に2回連続で行うと、その間で全データの転送が発生するはずである。data指示文で通信をなくしてみよう。（もちろん、independentなども指定したうえで。）#pragma acc kernels copy(c[N*N]) copyin(a[N*N],b[N*N])#pragma acc loop independent collapse(2)for(i=0; i<N; i++){for(j=0; j<N; j++){#pragma acc loop seq
for(k=0; k<N; k++){c[i*N+j] += a[i*N+k] * b[k*N+j];}}}#pragma acc kernels copy(c[N*N]) copyin(a[N*N],b[N*N])#pragma acc loop independent collapse(2)for(i=0; i<N; i++){for(j=0; j<N; j++){#pragma acc loop seq
for(k=0; k<N; k++){c[i*N+j] += a[i*N+k] * b[k*N+j];}}}(matmul5.c, matmul5a.c)!$acc kernels copy(c(N,N)) copyin(a(N,N),b(N,N))!$acc loop independent collapse(2)do i=1, N
do j=1, N!$acc loop seq
do k=1, Nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
!$acc end kernels!$acc kernels copy(c(N,N)) copyin(a(N,N),b(N,N))!$acc loop independent collapse(2)do i=1, N
do j=1, N!$acc loop seq
do k=1, N
NVCOMPILER_ACC_TIME=1c(j,i) = c(j,i) + a(k,i)* b(j,k)
や
enddo
enddo
NVCOMPILER_ACC_NOTIFY=2
enddoを設定して通信の状況を確認してみよう!$acc end kernels(matmul5.f90, matmul5a.f90)

実際にやってみましょう

117

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

118

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

119

最適化のためのヒント（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に勝てるとは
限らない

120

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

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

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

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

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

Array of Structure(AoS)とStructure of Array(SoA)• 構造体の配列（Array of Structure, AoS）と配列の構造体（Structure of Array, SoA）をGPUで扱うことを考えてみる
構造体の配列（AoS）

配列の構造体（SoA）

struct myStruct{
double a;double b;
double c;};struct myStruct data[N];struct myStruct{
double a[N];double b[N];double c[N];};struct myStruct data;// メンバ変数の値を2倍にする#pragma acc kernels loop copy(data[N])for(i=0; i<N; i++){data[i].a *= 2.0;
data[i].b *= 2.0;
data[i].c *= 2.0;}// メンバ変数の値を2倍にする#pragma acc kernels loop copy(data)for(i=0; i<N; i++){data.a[i] *= 2.0;
data.b[i] *= 2.0;
data.c[i] *= 2.0;}

(aos.c および soa.c)
どちらの実装がGPUに適しているのだろうか？

122

Array of Structure(AoS)とStructure of Array(SoA)：Fortran版• 構造体の配列（Array of Structure, AoS）と配列の構造体（Structure of Array, SoA）をGPUで扱うことを考えてみる
構造体の配列（AoS）

配列の構造体（SoA）

type MyType
double precision :: adouble precision :: b
double precision :: cend type MyType

type MyType
double precision :: a(n)double precision :: b(n)double precision :: c(n)end type MyType

type(myType) :: data(n)type(myType) :: data// メンバ変数の値を2倍にする!$acc kernels loop copy(data(n))
do i=1, ndata(i)%a = data(i)%a * 2.0d0
data(i)%b = data(i)%b * 2.0d0
data(i)%c = data(i)%c * 2.0d0
end do// メンバ変数の値を2倍にする!$acc kernels loop copy(data)
do i=1, n
data%a(i) = data%a(i) * 2.0d0
data%b(i) = data%b(i) * 2.0d0
data%c(i) = data%c(i) * 2.0d0
end do(aos.f90 および soa.f90)
どちらの実装がGPUに適しているのだろうか？

123

構造体の配列（AoS）の挙動• vector(4)による並列計算を考えてみる– （より長いvectorの場合も同様に考えれば良い）
struct myStruct{
double a;double b;
double c;};struct myStruct data[N];#pragma acc kernels loopfor(i=0; i<N; i++){data[i].a *= 2.0;
data[i].b *= 2.0;
data[i].c *= 2.0;}

メモリ上のdata[N]の配置イメージ
data[0]a0bdata[1]c

1a2

CUDAコアdata[2]

b

c

a

b

3

timestep 0
timestep 1
timestep 2

data[3]cabc

CUDAコアへの計算割り当てのイメージ➢ あるタイミングにおける連続した計算コアのアクセス先は不連続➢ メモリアクセスがまとめて行われないため低速

...

124

配列の構造体（SoA）の挙動• vector(4)による並列計算を考えてみる– （より長いvectorの場合も同様に考えれば良い）
struct myStruct{
double a[N];double b[N];double c[N];};struct myStruct data;

メモリ上のdataの配置イメージ（N=4）

a

0

a

a

1

a

2

CUDAコア#pragma acc kernels loopfor(i=0; i<N; i++){data.a[i] *= 2.0;
data.b[i] *= 2.0;
data.c[i] *= 2.0;
}

b

b

b

b

3

timestep 0
timestep 1
timestep 2

c

c

c

c

CUDAコアへの計算割り当てのイメージ➢ あるタイミングにおける連続した計算コアのアクセス先は連続➢ メモリアクセスがまとめて行われるため高速➢ マルチコアCPUで並列計算を行う場合はアクセス競合の都合でAoSの方が高速になりやすい（CPUとGPUで適したコードが異なる）

125

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

126

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)

127

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

128

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129

実行例• 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/20211004/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/20211004/deepcopy2.f90 function=main line=32 device=0 threadid=1 variable=obj bytes=144
upload CUDA data file=/home/center/a49979a/share/20211004/deepcopy2.f90 function=main line=32 device=0 threadid=1 variable=.attach. bytes=8
upload CUDA data file=/home/center/a49979a/share/20211004/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

130

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)

131

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指示文を加えても逐次実行になってしまう）

132

関数内のループを並列実行する方法• 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)

133

関数内を並列実行する方法：複数ソースファイルの場合（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をつけても得られる情報はないはず、つけたままでも特に問題はない）

134

関数内を並列実行する方法：複数ソースファイルの場合（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135

OpenACC API関数の利用• OpenACCは指示文を用いて使うものであるが、その多くの機能は関数によっても提供されている– OpenMPのomp_で始まる関数のようなものと思えば良い– 例• 利用可能なGPUの数を得る– acc_get_num_devices• CPU-GPU間のデータコピー– acc_copyin, acc_copyout– C言語では#include “openacc.h”、Fortranではuse openaccして使う– 使いたいコンパイラの組み合わせ等によっては有益なこともあるかもしれない• OpenACCに対応していないコンパイラと組み合わせて使いたい場合など？（未検証）

136

外部連携
• 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)

137

デバッガとプロファイラ• デバッガ、プロファイラ：基本的にCUDAと共通のツールが利用できる– HPC SDKに含まれるcuda-gdbやnvprof、nsysなどが利用できる– GUIが表示されるものはX転送が必要• MovaXtermは標準対応、Putty単体では対応不能、CygwinはX Window関係のパッケージが必要• sshコマンド実行に-Yオプションを付けておく必要がある– 必要に応じて、コンパイラオプション -gを付けたり、-taオプションにdebugやlineinfoを付ける（はずなのだが、付けなくても普通に情報を見られる気もする……）$ nvc -Minfo=accel -g -acc -gpu=cc70,lineinfo -tp=host -o m2_c_acc_gl matmul2.c※実際には-gとlineinfoは同時に指定できない、指定するとlineinfoが無視される• デバッガ– 単純なプログラムのOpenACC化でデバッグに困ることはあまりないと思われる– 意図した通りにGPU化やデータ転送ができているか、コンパイラの出力に注意が必要• プロファイラ– NVIDIA社が提供：nvprof, nvvp, Nsight– NVIDIA社以外が提供するツールもいくつかあるようだ138

nvprof の利用例
$ nvprof ./a1c_c_acc==96336== NVPROF is profiling process 96336, 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==96336== Profiling application: ./a1_c_acc==96336== Profiling result:Type Time(%)
Time
Calls
GPU activities: 42.33% 4.0640us235.67% 3.4240us1
22.00% 2.1120us1
API calls: 80.70% 353.32ms
1
12.26% 53.667ms
1
5.14% 22.489ms
1
1.59% 6.9831ms
1
0.19% 817.59us1
0.06% 245.00us
4～中略～OpenACC (excl): 98.08% 22.609ms
1
0.59% 135.03us1
0.38% 87.861us1
0.28% 64.660us
1
0.12% 28.616us1
0.11% 25.423us1
0.10% 23.474us
1
0.07% 17.253us1

実はCUDAに関係する関数が呼ばれており、OpenACCがCUDAの上で動いているような関係にあることを確認することができる

OpenACCコードの何行目

Avg
Min
2.0320us 1.6000us3.4240us 3.4240us2.1120us 2.1120us353.32ms 353.32ms
53.667ms 53.667ms
22.489ms 22.489ms
6.9831ms 6.9831ms
817.59us 817.59us61.249us 4.2880usMax
2.4640us
3.4240us
2.1120us353.32ms
53.667ms
22.489ms
6.9831ms
817.59us223.16us

Nameにどれだけ時間がかかっ[CUDA memcpy HtoD]
たかがわかる
main_21_gpu[CUDA memcpy DtoH]
cuDevicePrimaryCtxRetain
cuDevicePrimaryCtxRelease
cuMemHostAlloc
cuMemFreeHost
cuMemAllocHost
cuMemAlloc

22.609ms 22.609ms
135.03us 135.03us87.861us 87.861us64.660us 64.660us28.616us 28.616us25.423us 25.423us23.474us 23.474us17.253us 17.253us22.609ms
135.03us87.861us64.660us28.616us25.423us23.474us17.253us

acc_enter_data@allocate1.c:20
acc_device_init@allocate1.c:20
acc_enqueue_launch@allocate1.c:21 (main_21_gpu)
acc_enqueue_upload@allocate1.c:20
acc_enqueue_download@allocate1.c:25
acc_enqueue_upload@allocate1.c:21
acc_exit_data@allocate1.c:21
acc_enter_data@allocate1.c:21

139

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向けしか存在しない）

• 画面例

→

140

Nsight• 新しいプロファリングツール– Nsight Systems, Nsight Compute,
Nsight Graphicsから構成される– CUIとGUIそれぞれで利用できる• CUIによるプロファイル情報の取得と表示

Generating CUDA API Statistics...
CUDA API Statistics (nanoseconds)
Time(%)------77.1
13.2
3.8
2.9
1.4
0.6
0.3
0.3
0.2
0.1
0.1
0.1
0.0

Total Time-------------23989788
4123835
1181563898283423594
181128
94849
82336
53309
37604
28981
22506
6956

Calls
---------1
1
10
1
20
28
7
6
7
7
1
7
2

Average-------------23989788.0
4123835.0
118156.3
898283.021179.76468.913549.9
13722.7
7615.65372.028981.03215.13478.0

Minimum-------------23989788
41238352924
898283
5039
2476
4444
3396
3497
1915
28981
1366
1764

Maximum-------------23989788
4123835
1110527
898283
52341
20178
25066
27629
29734
13222
28981
5176
5192

Name
--------------------cuMemHostAlloc
cuModuleLoadDataEx
cuMemAlloc_v2
cuMemAllocHost_v2
cuLaunchKernel
cuStreamSynchronize
cuMemcpyDtoHAsync_v2
cuMemsetD32Async
cuMemcpyHtoDAsync_v2
cuEventRecord
cuStreamCreate
cuEventSynchronize
cuEventCreate– nsys profile ./a.out
Generating CUDA Kernel Statistics...
CUDA Memory Operation Statistics...
– --stats=trueオプションをつけておく Generating
CUDA Kernel Statistics (nanoseconds)とその場で統計情報も得られる
Time(%)
Total Time
Instances
Average
Minimum
Maximum Name------- -------------- ---------- -------------- -------------- -------------- --------------------• 統計情報の一部 →
12.2
4448
2
2224.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
report1.qdstrm
9.8
3584
2
1792.01664
1920 main_171_gpu
9.5
3456
2
1728.01312
2144 main_123_gpu– Nsight Systemsで閲覧するには-oで
9.4
3424
2
1712.01664
1760 main_191_gpu
9.1
3296
2
1648.01504
1792 main_181_gpu
qdrepファイルを指定しておく
3.8
1376
1
1376.01376
1376 main_142_gpu
3.8
1376
1
1376.01376
1376 main_150_gpu必要があるようだ

141

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との連携も必要– クライアントからログインノードに接続して使えると良いのだが、いまのところうまく動かせていない

142

複数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を担当するような実行形態を考えておくと汎用性があって使い
やすい

143

ノード内複数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

144

ノード内複数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;}}!$omp parallel dodo 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
enddo

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;
}
}
}
}
}
}
}
}

145

ノード内複数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

146

ノード内複数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

147

ノード内複数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

148

複数ノード複数GPUの活用：OpenACC + MPI• プログラムとしては「OpenACC+MPIプログラムの概形3」の例のままでよい• ジョブスクリプトに多少の工夫が必要• 2ノード×4GPU（合計8MPIプロセス）の例
job.shrun.sh（chmodで実行権を与えておくこと）#!/bin/bash -x#PJM -L rscgrp=cx-workshop#PJM -L node=2#PJM -L elapse=00:5:00#PJM -j#PJM -S#!/bin/bash
export CUDA_VISIBLE_DEVICES=${OMPI_COMM_WORLD_LOCAL_RANK}
numactl -l $1

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

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

149

複数ノード実行に関するソースコードとジョブスクリプト例一覧• 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にて紹介していきます。

150

その他の指示文について• ここまでに扱ってきた指示文だけでも多くのプログラムをOpenACCによりGPU化できる• まだ扱っていない指示文についても少しだけ紹介しておく– async, wait
• 通信や計算が終了するのを待たずに次の処理を行うasyncと、その終了を待つためのwait• うまく使えれば様々な待ち時間を削減し全体の実行時間を短縮させることができる• 例： for(i=...){#pragma acc update device(array) async(i)#pragma acc parallel loop async(i)for(...){...}#pragma acc update host(array) async(i)}#pragma acc wait••

ループカウンタ i ごとに更新と計算を組にして処理を発行終了を待たずに次の i の発行を行い、最後に待ち合わせる– tile• 二次元タイル分割を行うための指示節• キャッシュ利用効率が向上することで性能向上することがある• 指定したサイズが適していないと性能低下することがある• 例：#pragma acc kernels loop tile(64,4)

151

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

152

まとめ
• OpenACCの仕様や使い方について紹介した– わずか数行の指示文でGPUが利用できる• 触れなかった機能・指示文もいくつかあります：tileやasyncなど– OpenMPと似ているため、両方使えると色々と活用できる– 並列計算環境（プログラミング言語など）を取り巻く状況は流動的であり、今後何十年もOpenACCが生き残るかは誰にもわからないが、OpenACCは簡単に高い性能が得られることが多い
ため、是非活用してみて欲しい• OpenMP, oneAPI, SYCL, DPC++などによるGPU利用が進んでくるとどうなるかわからない状況• 残りの時間は演習時間とします– 演習問題としてCG法の逐次コードを提供しているので、指示文を挿入して並列化してみましょう– その他、任意のプログラムのGPU化やプロファイラのテストなどをしてもらって構いません

講習会アカウントは明後日の23:59まで有効です

153

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

154

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

155

今回用いる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だけをみて並列高速化を目指してみよう

★

★

★

★

★

156

157

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

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

対角はN-1対角以外は1

ベクトルxx（サイズN）
ランダム
ベクトルb（サイズN）➢ A, x, xxを初期化し、b = A×xxを計算しておいたうえで、cg法でAx = bを計算しxを求めるというプログラムにしてある。
➢ これなら前処理のない単純なCG法でも簡単に収束する。
➢ むしろ簡単過ぎてすぐに収束して終了してしまうため、意図的に値を破壊して収束しないようにしてある。
（行列サイズN*10回の反復を行い、収束しないまま終わる。）➢ 正しく実装できているかを確認したい場合は妨害コードをコメントアウトすると良い。

実行例

•
•
•

第一引数は問題サイズN
第二引数を0にすると初期データや計算結果を出力しなくなる
第三引数に0以外を指定すると収束しても計算を続行する（すぐに終わってしまうのを防ぐもう1つの策）$ ./cg5_c_acc
N=10, check=1, nocond=0A:
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↑既知の行列A←ベクトルX（求める答え）＝前ページのベクトルXX←既知のベクトルb（右上に続く）iter 1 1.075989e-01 1.000000e-16
iter 2 1.328273e-03 1.000000e-16（省略）
iter 99 9.776105e-01 1.000000e-16
iter 100 9.874046e-01 1.000000e-16
time: 1.848180e-01 sec , 1.848180e-03 sec/iter
result:9.185754
9.742296
8.072670-3.429193
11.040893
0.281085
9.742296
10.855379
2.878280-2.316110
diff:
1 8.400000 9.185754: -7.857539e-01
2 8.700000 9.742296: -1.042296e+00
3 7.800000 8.072670: -2.726703e-01
4 1.600000 -3.429193: 5.029193e+00
5 9.400000 11.040893: -1.640893e+00
6 3.600000 0.281085: 3.318915e+00
7 8.700000 9.742296: -1.042296e+00
8 9.300000 10.855379: -1.555379e+00
9 5.000000 2.878280: 2.121720e+00
10 2.200000 -2.316110: 4.516110e+00↑反復計算の履歴、反復計算終了時に所要時間も表示←計算結果ベクトルX↑計算結果ベクトルXと最初に設定したベクトルXXの比較（測定のため収束を妨害してあるので結果が合わなくて良い）

158

159

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

160

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

161

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)

162

状況の確認
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

163

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)

164

データ転送の最適化を考える• 現在のデータ転送状況のイメージ
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

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

165

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節を用いて途中の配列データを覗いてみると良い

166

実は……
• 最初に全てのデータを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節を加えたもの）