第208回 お試しアカウント付き
並列プログラミング講習会「OpenMPで並列化されたC++プログラムのGPU移植手法」
東京大学 情報基盤センター
三木 洋平

2

講習会概略
• 開催日： 2023年6月29日（木） 13:00--17:00• 形態： ZoomおよびSlackを用いたオンライン講習会• 使用システム： Wisteria/BDEC-01 (Aquarius)• 講習会プログラム：• 13:00--13:30• 13:40--14:20• 14:30--15:10• 15:20--15:50• 16:00--17:00

2023/6/29

GPU向けプログラミング手法の紹介（座学）
OpenACCを用いたGPU化（座学＋演習）
OpenMPを用いたGPU化（座学＋演習）
C++17を用いたGPU化（座学＋演習）
CUDA C++を用いたGPU化（座学＋演習）

講習会：GPU/C++

3

Contents• GPUに関する基礎知識• （主にNVIDIA製の）GPU向け各種プログラミング手法の紹介• OpenACC• OpenMP の target 指示文（OpenMP GPU offloading）• 言語の標準規格を利用してのGPU化（C++17以降）• CUDA C++• 発展的な内容（情報提供のみ）• Independent Thread Scheduling• A100/CUDA 11から導入された機能• CUDA-aware MPI• GPUDirect Storage2023/6/29

講習会：GPU/C++

4

講習会概略
• 開催日： 2023年6月29日（木） 13:00--17:00• 形態： ZoomおよびSlackを用いたオンライン講習会• 使用システム： Wisteria/BDEC-01 (Aquarius)• 講習会プログラム：• 13:00--13:30• 13:40--14:20• 14:30--15:10• 15:20--15:50• 16:00--17:00

2023/6/29

GPU向けプログラミング手法の紹介（座学）
OpenACCを用いたGPU化（座学＋演習）
OpenMPを用いたGPU化（座学＋演習）
C++17を用いたGPU化（座学＋演習）
CUDA C++を用いたGPU化（座学＋演習）

講習会：GPU/C++5

GPU（Graphics Processing Units）• 元々は画像処理を行うために特化していた専用プロセッサを汎用計算にも使えるように拡張• 最近のGPUには（主に深層学習用に）行列積向けユニットも• 高い演算性能，太いメモリバンド幅，消費電力あたり性能も高い• 昔は安かったが，最近はとても高い• 多数のコア（最新GPUは１万越え）を搭載した並列計算機• NVIDIA A100 (SXM, PCIe): 64 × 108 SMs = 6912 cores• NVIDIA H100 (SXM):
128 × 132 SMs = 16896 cores• NVIDIA H100 (PCIe):
128 × 114 SMs = 14592 cores• AMD MI250X:
64 × 110 CUs × 2 GCDs = 14080 cores• AMD MI250:
64 × 104 CUs × 2 GCDs = 13312 cores• スレッド並列（e.g., OpenMP）的考えが分かっていると良い
2023/6/29

講習会：GPU/C++Green500 Ranking (June 2023)
TOP
500
Rank

System

Accelerator

Cores

6

https://www.top500.org/

HPL Rmax[PFlop/s]

Power
[kW]

GFLOPS/W

1

255

Henri, Flatiron Institute, USA

NVIDIA H100 80GB PCIe

8,288

2.88

44

65.3962

34

Frontier TDS, ORNL, USA

AMD Instinct MI250X

120,832

19.20

309

62.684312Adastra, GENCI-CINES, France

AMD Instinct MI250X

319,072

46.10

921

58.021

4

17

Setonix, Pawsey, Australia

AMD Instinct MI250X

181,248

27.16

477

56.983577

Dardel GPU, KTH, Sweden

AMD Instinct MI250X

52,8648.26

146

56.491

6

1

Frontier, ORNL, USA

AMD Instinct MI250X

8,699,904

1,194.00

22,70352.592

7

3

LUMI, EuroHPC/CSC, Finland

AMD Instinct MI250X

2,220,288

309.106,016

51.382

8

483

6,768

1.95

38

51.343970

96,7689.09

19546.543

10

187

25,0563.50

8641.411

7,630,848

442.01

29,89915.41843,008

2.59

409

6.327

49
77

amplitUDE (GPU Partition), University
NVIDIA H100 80GB PCIeof Duisburg-Essen, Germany
Goethe-NHR, Universitaet Frankfurt,
AMD Instinct MI210 64GB
Germany
ATOS THX.A.B, Atos, France

NVIDIA A100 SXM 64GBSupercomputer Fugaku, R-CCS,Japan
Conesus, University of Rochester,311
USA
2023/6/29(CPU-only:
Fujitsu A64FX)(x86 CPU-only:Intel Xeon
Platinum 8480+)
講習会：GPU/C++

2

7

次期スーパーコンピュータ OFP-II• 最先端共同HPC基盤施設（JCAHPC：東京大学情報基盤セン
ターと筑波大学計算科学研究センターが共同で運営）の次期シス
テム（OFP-II）の調達が進行中• 2025年1月運用開始予定• 理論ピーク性能150 PFlop/s 以上を目標• GPU搭載ノードを中心としたシステムとなる予定• GPU移行を支援するためのポータルページを公開• https://jcahpc.github.io/gpu_porting/• 各種情報を追加していく予定• （今日使うN体コードも，このページに置くために作ったもの）• GPU移行相談会も月一ぐらいの頻度で開催中
2023/6/29

講習会：GPU/C++

8（NVIDIAの）GPUの構成（1/2）

NVIDIA A100 Tensor Core GPU Architecture2023/6/29

講習会：GPU/C++

9（NVIDIAの）GPUの構成（2/2）• SM：Streaming Multiprocessor• SMの中での構造は気にしなくても性能が出せる• 注：32スレッドのグループ（ワープ）内で同じ動作をするように意識しておく• SM内ではL1キャッシュ/シェアードメモリを共有（A100では192KB）• 配分は（ある程度）調節可能• シェアードメモリはCUDAでは使えるがOpenACCでは（明示的には）使えない機能• グローバルメモリよりも圧倒的に速い
2023/6/29

講習会：GPU/C++

NVIDIA A100 Tensor Core GPU Architecture

10

GPUプログラミングでの基本思想• このスライドではCUDA用語で説明• 演算器の数よりも多数のスレッドを立てて計算• 各種レイテンシを隠蔽するため• スレッドブロックが基本単位• ブロックあたりのスレッド数は32の倍数，
できれば128以上が推奨• SMに複数ブロックを割り当てるのが普通• ブロック内では32スレッド単位で動作（この組をワープと呼ぶ）• スレッドブロックの集合をグリッドと呼ぶが，意識しなくてOK
2023/6/29

https://developer.nvidia.com/blog/cudarefresher-cuda-programming-model/

講習会：GPU/C++11

他手法での用語紹介（覚える必要なし）• Vector: CUDA でのスレッドに対応• ここの並列度は性能へのインパクト大（32の倍数，できれば128以上）• OpenMPのGPUオフローディングでは SIMD に対応• Worker: CUDAでのスレッドブロックに対応• Vector の集合• Worker 単位でSMに投入される（が，こんなことまで気にしなくてOK）• OpenMPのGPUオフローディングでは thread に対応• Gang: CUDAでのグリッドに対応• Worker の集合• ジョブの全体ということになるが，全く気にしなくてOK• OpenMPのGPUオフローディングでは teams に対応
2023/6/29

講習会：GPU/C++（独断と偏見に基づく不完全な）
GPUプログラミング手法の比較表

12• GPU化は自分のやりたいことに適した手法の選択から始まる• 耳障りの良い言葉ばかり聞いていると後で苦しむことになるので注意
ベース言語など
強み
弱み
Fortran対応
対象GPU

手法

CUDA C++

C++の拡張詳細な最適化可能 記述量が多い
最新機能が使える GPU専用コード

CUDA Fortran NVIDIA限定

OpenACC

指示文

移植コスト低CUDAより遅い
CPUコードと共通 （メモリ律速なら化可能それなり？）

OK

ほぼNVIDIA限定（GCCならAMDも動くが低性能）OpenMP（target指示文）

指示文

移植コスト低CUDAより遅い
CPUコードと共通 （メモリ律速なら化可能それなり？）

OK（flang版はまだ性能いまいち？）NVIDIA, AMD,Intel

言語の標準規格

C++17以降
Fortran 2008

CPUでも同じコー 多くの制約あり
ドが動く
CUDAより遅い

Fortran 2008

NVIDIA（対応GPUが今後増えるかは不明）

HIP

C++の拡張詳細な最適化可能 記述量が多い
ほぼCUDA C++ GPU専用コード

GPUFORTの発展に期待？

AMD, NVIDIA（将来はIntelも？）DPC++

C++，SYCL？

？

Intel, (NVIDIA)（hipSYCLもある）

2023/6/29？
講習会：GPU/C++13

大まかな考え方の違い• 簡易GPU化（OpenACC, OpenMP, C++17, …）• 並列化可能なループであることをコンパイラに伝える（自明な並列度があるループであるので，簡単にGPU化できるループ）• 「どうGPU化するか」はコンパイラ任せ• なるべくコンパイラの自由度を増やしてあげる方が速くなる傾向• 自力でのGPU化（CUDA, HIP, …）• 「どうGPU化するか」を自分で考え，自分で実装• どういう演算命令を使うか，など自分で何でもできる• OpenACC, OpenMP, C++17ではGPU化が難しいループであっても，GPU化できる• 例：ツリー構築のGPU化は指示文ベースでできる気はしないが，CUDAなら可能

2023/6/29

講習会：GPU/C++

14

GPUプログラミングに関する資料• 「UTokyo N-Ways to GPU Programming Bootcamp」講習会資料• https://www.cc.u-tokyo.ac.jp/events/lectures/207/• 講習会の録画も公開されています• ISO標準言語，OpenACC，CUDA• 「GPUプログラミング入門」講習会資料• https://www.cc.u-tokyo.ac.jp/events/lectures/206/• OpenACC• 「OpenACCとMPIによるマルチGPUプログラミング入門」講習会資料• https://www.cc.u-tokyo.ac.jp/events/lectures/209/• OpenACC+MPI• GPU移行に関するポータルサイト• https://jcahpc.github.io/gpu_porting/2023/6/29

講習会：GPU/C++

15

Zoom関連• 「手をあげる」機能• 質問がある際，全体の状況を確認するため使用• ブレークアウトセッション• 画面を共有しながらエラー対応する際に使用• （なるべく口頭でのやりとりやSlackで対応する予定）• https://utelecon.adm.utokyo.ac.jp/zoom/how_to_use

2023/6/29

講習会：GPU/C++

16「手を挙げる」方法
1. Zoomメニュー中の「リアクション」をクリック

2. ポップアップで表示された「挙手」をクリック

2023/6/29

講習会：GPU/C++

17

手が挙がっていることの確認方法
1. Zoomメニュー中の「参加者」をクリックして，参加者一覧を表示

2. 表示された参加者一覧の，自分のところを見ると手が挙がっている

2023/6/29

講習会：GPU/C++

18「手をおろす」方法
1. ポップアップで表示されている「手をおろす」をクリック• もし「手をおろす」が表示されていなければ，「リアクション」の中から探す

2023/6/29

講習会：GPU/C++

19

ブレイクアウトルーム（1/6）• 演習時に使用するかもしれません• 演習中に「ヘルプを求める」ことができます• ホストを招待した後に「画面を共有」することで，皆さんの記述したプログ
ラムを一緒に見ながら問題解決にあたります• Zoomメニュー中の「ブレイクアウトルーム」をクリック

2023/6/29

講習会：GPU/C++

20

ブレイクアウトルーム（2/6）• 進行中のブレイクアウトルームのリストが表示されるので，空いている部屋に「参加」してください• 左の例では5部屋がすべて空室，右の例ではルーム1のみ参加者がいる

2023/6/29

講習会：GPU/C++

21

ブレイクアウトルーム（3/6）• 「参加」をクリックすると確認画面が出てくるので，「はい」を選択すると入室できます

2023/6/29

講習会：GPU/C++

22

ブレイクアウトルーム（4/6）• 再度メニュー中の「ブレイクアウトルーム」をクリックすると，「ヘルプを求める」が増えているのでクリックしてください

2023/6/29

講習会：GPU/C++

23

ブレイクアウトルーム（5/6）• ポップアップで出てきた「ホストを招待」をクリック• ホスト（講師）の参加待ちに移行します（画面はそのまま）• 他の受講者のヘルプ中など，直ちに対応できない場合もあります

2023/6/29

講習会：GPU/C++

24

ブレイクアウトルーム（6/6）• 問題解決後は，Zoomメニュー中の「ルームを退出する」• 「ブレイクアウトルームを退出」が表示されるのでクリックして，元の講習会会場にお戻りください• 間違えて「ミーティングを退出」すると講習会から退出します

2023/6/29

講習会：GPU/C++

25

Slack関連• ブラウザ上で使う場合には：• https://w1590055008-bgo338004.slack.com/• 注：ログインには，事前にお配りしたリンクからの登録が必要です• 質問対応に使用• コードの貼り付け方• スレッドの確認方法• 以下，ブラウザ版で説明しますがアプリ版でも操作は同じです

2023/6/29

講習会：GPU/C++

26

質疑応答チャンネルへの移動• 左側のメニューバーのチャンネル一覧内に「第208回-cpp2gpu」があるので，クリック• 見つからない場合
1. 「チャンネルを追加する」をクリック
2. 「チャンネル一覧を確認する」をクリック
3. 「第208回-cpp2gpu」があるので「参加する」をクリック

2023/6/29

講習会：GPU/C++

27

メッセージの入力方法• 最下部に入力欄があるので，質問内容を記載して Ctrl+Enter• 入力後に右下の「メッセージを送信する」をクリックしても同じ（メッセージ入力前には，「メッセージを送信する」は押せない）
メッセージの入力欄
コードブロックの生成

メッセージを送信する• コードを入力する際には，「コードブロック」がおすすめ• 枠が生成されるので，この中にコピペするのが簡単かつ見やすい• ```（JIS配列ならばShift+@を3連打）しても枠が生成される

2023/6/29

講習会：GPU/C++

28

自分が参加したスレッドのみの表示方法• 左上の「スレッド」をクリックすると，自分が参加しているスレッドの一覧が表示されます• 質問内容には，「スレッドで返信する」形式で回答するので，自分が聞いた内容のみが表示できます

2023/6/29

講習会：GPU/C++

29

Slackメッセージの保存• 講習会slackは無償版なので，90日以上過去のメッセージは表示されなくなります• メッセージの保存方法は複数ありますが，比較的簡単な方法を紹介• SlatickHTML:https://sites.google.com/view/slatickhtml/• Chrome拡張機能なので，Slack Appの作成などが不要• Vivaldi，Edgeなどのブラウザからも同様に使用できる（はずですが，確認していません）• （添付するファイル形式をチェックしておけば）ファイル取得も可能• HTMLファイルとして保存するため，保存後はブラウザを用いて表示可能• デフォルトでは今年のメッセージしか表示されないので注意（上部のメニューで「全て表示」か該当年を選択すれば，過去のメッセージも表示されます）
2023/6/29

講習会：GPU/C++

30

ユーザアカウント• 使用システム： Wisteria/BDEC-01 (Aquarius)• $ ssh USERNAME@wisteria.cc.u-tokyo.ac.jp• 本講習会でのユーザ名（上記コマンドでの USERNAME のこと）
• 利用者番号： tABCDE（ABCDEは，適宜書き換えてください）• 利用グループ： gt00

• 利用期限• 7/29 9:00まで有効• 注:本講習会関連の質問はymiki[at]cc.u-tokyo.ac.jpまで• Slack で質問していただいても結構です• （講習会アカウントでは）公式の相談対応システムは使わないでください

2023/6/29

講習会：GPU/C++

31

講習会概略
• 開催日： 2023年6月29日（木） 13:00--17:00• 形態： ZoomおよびSlackを用いたオンライン講習会• 使用システム： Wisteria/BDEC-01 (Aquarius)• 講習会プログラム：• 13:00--13:30• 13:40--14:20• 14:30--15:10• 15:20--15:50• 16:00--17:00

2023/6/29

GPU向けプログラミング手法の紹介（座学）
OpenACCを用いたGPU化（座学＋演習）
OpenMPを用いたGPU化（座学＋演習）
C++17を用いたGPU化（座学＋演習）
CUDA C++を用いたGPU化（座学＋演習）

講習会：GPU/C++

32

NVIDIA GPU向けに各手法でGPU化• 対象とする計算：N体計算（直接法）• 無衝突系向け（2次Leapfrog, shared time step）• 衝突系向け（4次Hermite，block time step）• 今日試してみるプログラミング手法• OpenACC• OpenMP• C++17• CUDA C++• 今日使うサンプルコード（上記のGPU実装）• https://github.com/ymiki-repo/nbody• CMake のサンプルとしても参照可（HDF5部分は実装が特殊）
2023/6/29

講習会：GPU/C++

33

N体計算（重力多体計算）• 粒子どうしに働く自己重力による系の時間進化を，運動方程式に基づいて計算• データ量: 𝒪(𝑁)• 重力計算: 𝒪(𝑁 2 )• 時間積分: 𝒪(𝑁)ij

• 業界的によく使う用語j• i-粒子: 重力を受ける粒子• j-粒子: 重力を及ぼす粒子jj• （N体計算の詳細は国立天文台CfCAのN体学校を参照）
2023/6/29

講習会：GPU/C++

34

OpenACC/OpenMP での実装方針
1. Unified Memory を使って実装• GPU上のメモリ確保，CPU-GPU 間のデータ転送は全てお任せ• まずは演算部分のGPU実装に注力• （マルチコアCPU向けの）OpenMP実装されていれば，#pragma omp parallel for をGPU向けの指示文に置き換えていく

2. （Unified Memory 実装での性能に満足できない場合）
データ指示文を使ってコードをアップデート• Unified Memory では，一旦メモリを読んで，ページフォルトがあればハードウェアレベルでページ単位で転送という仕組み• 必要なデータ転送は自分で指示した方が必然的に速くなる• Unified Memory ではCPU上のアドレスとGPU上のアドレスを同一視してしまうので，GPUDirect系の機能を活用する際に不利
2023/6/29

講習会：GPU/C++

35

OpenACCでの実装（演算部分）• #pragma omp parallel for を置き換えていく• 性能的に特に重要なものについては，スレッド数の調整も• vector（スレッド数） としてコンパイラに示唆（強制はできない）#pragma omp parallel forfor (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}#pragma acc kernels#pragma acc loop independent vector(NTHREADS)
for (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}#pragma omp parallel forfor (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}#pragma acc kernels#pragma acc loop independentfor (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}

2023/6/29

講習会：GPU/C++

36

OpenACCコードのコンパイル方法など• NVIDIA HPC SDK向けの情報• $ nvc++ -acc=gpu -gpu=cc80 -Minfo=accel,opt• リンク時にも -acc=gpu を指定する• $ nvc++ -acc=gpu -gpu=cc80,managed -Minfo=accel,opt• Unified Memory使用時のコンパイル方法• リンク時にも -acc=gpu -gpu=managed を指定する• デバッグ時などに便利な環境変数• NVCOMPILER_ACC_NOTIFY=1• GPU上でカーネルが実行されるたびに情報を出力• NVCOMPILER_ACC_NOTIFY=3• CPU-GPU間のデータ転送に関する情報も出力• NVCOMPILER_ACC_TIME=1• CPU-GPU間のデータ転送およびGPU上での実行時間を出力
2023/6/29

講習会：GPU/C++

37

計算機へのログインとサンプルの取得• $ ssh USERNAME@wisteria.cc.u-tokyo.ac.jp• $ cd /work/gt00/$USER• $ git clone https://github.com/ymiki-repo/nbody.git• $ cd nbody

2023/6/29

講習会：GPU/C++

38

コンパイル方法• $ module load cmake nvidia/22.7 hdf5
• nvidia/23.3(default) 用のHDF5環境が未構築のため• $ cmake –S . –B build_nvhpc• $ mkdir build_nvhpc; cd build_nvhpc; cmake .. # と同じ• $ cd build_nvhpc• $ ccmake –S .. # デフォルトの設定を調整したい場合• 倍精度で計算したければ，FP_L, FP_Mを64に変更• Hermite法で試したい方は，HERMITE_SCHEMEをONに変更• $ ninja acc_managed• 単に $ ninja だけでもOKですが，他にもたくさんのビルドが走ります2023/6/29

講習会：GPU/C++

39

moduleの指定• コンパイラ・ライブラリ等の環境をセットアップ• Odyssey向け富士通コンパイラ、MPIを使用$ module load fj （fjmpiは自動的にロードされる） または$ module load odyssey• Aquarius向けgcc、cuda, Open MPI(CUDA対応)を使用$ module load gcc cuda ompi-cuda または$ module load aquarius cuda ompi-cuda• 現在指定済みのmoduleを確認するには$ module list

• 困った時（module環境を整理したくなった時）は$ module purge
2023/6/29

講習会：GPU/C++

40

moduleの一覧• 現在利用中の環境で追加できるものを確認$ module avail• 使い方の確認（例はModuleNameというmoduleのヘルプを表示）$ module help ModuleName• 設定されるPATHなどの確認（例はModuleNameというmoduleの場合）$ module show ModuleName• 全ての環境を確認（Wisteria/BDEC-01向けに提供されるコマンド）$ show_module

ApplicationName
ModuleName
Node
BaseCompiler/MPI-------------------------------------------------------------------------Archiconda
archiconda3/0.2.3 odyssey
Arm Forge
forge/21.0.2
aquarius…

2023/6/29

講習会：GPU/C++

41

実習• $ pjsub -x EXEC=bin/acc_managed
sh/wisteria/run_nvidia.sh• sh/wisteria/run_nvidia.sh 内に直接 EXEC=bin/acc_managed と書き込んだ上で $ pjsub sh/wisteria/run_nvidia.sh としても同じ• 正常に終了すると，下記ファイルが出力されるはずです• log/collapse_run.csv• エネルギー保存していることを確認（右端2列の数字が 0.005 より十分に小さい）• dat/collapse_snp[000 – 080].[h5 xdmf]• 本当はここで可視化して計算結果をチェックするところだが，（環境構
築の都合上）今回は省略• Julia スクリプトを使ってエネルギー保存，ビリアル比の時間進化を表示• Julia スクリプトを使って粒子分布の時間進化を表示• VisIt や ParaView を使って粒子分布の時間進化を表示
2023/6/29

講習会：GPU/C++

42

ジョブスクリプトの説明#!/bin/bash#PJM -L rscgrp=lecture-a#PJM -L gpu=1#PJM -L elapse=0:15:00#PJM -g gt00

リソースグループ名: lecture-a
利用GPU数: 1 GPU使用
実行時間制限: 15分

if [ -z "${EXEC}" ]; then
EXEC=bin/basefi

利用グループ名: gt00
実行ファイルの指定
環境設定

module purge
module load nvidia/22.7
module load hdf5

GPUジョブを実行

cd $PJM_O_WORKDIR
numactl --localalloc ${EXEC} ${OPTION}
2023/6/29

講習会：GPU/C++

43

Wisteria/BDEC-01でのジョブ実行• 以下の2通りの実行形態があります
1. バッチジョブ実行• バッチジョブシステムに処理を依頼して実行• 実行したい処理をファイル（ジョブスクリプト）で指示• スパコン環境で一般的• 大規模実行用※講習会アカウントでは
バッチジョブ実行のみ，最大4 GPUs 15分まで• Wisteria/BDEC-01 (Odyssey)では，最大2304ノード（110592コア），24時間
まで• Wisteria/BDEC-01 (Aquarius)では，最大8ノード（64 GPUs），24時間まで

2. インタラクティブジョブ実行• PCでの実行のように，コマンドを入力して実行• スパコン環境では一般的ではない• デバッグ用，大規模実行はできない• 1ノード（48コア）: 30分まで• 12ノード（576コア）: 10分まで

2023/6/29

講習会：GPU/C++

44

バッチ処理とは• スパコン環境では，通常は，インタラクティブ実行（コマンドラインで実行すること）はできません• ジョブはバッチ処理で実行します
バッチキュー
バッチ処理
システムがジョブを取り出す

ジョブの依頼実行
ユーザ
2023/6/29

講習会：GPU/C++

スパコン

45

バッチキューの設定方法• Wisteriaでのバッチ処理は，富士通のバッチシステムで管理• 主要コマンド：• ジョブの投入: pjsub <ジョブスクリプト名>• 自分が投入したジョブの状況確認: pjstat• 投入ジョブの削除: pjdel <ジョブID>• 計算ノードの込み具合を見る: pjstat --rscuse• バッチキューの状態を見る: pjstat --rsc• バッチキューの詳細構成を見る: pjstat --rsc -x• 投げられているジョブ数を見る: pjstat --rsc -b• 過去の投入履歴を見る: pjstat -H• 同時に投入できる数・実行できる数を見る: pjstat --limit
2023/6/29

講習会：GPU/C++参考：インタラクティブ実行のやり方（本講習会アカウントでは使えません）• 1ノード実行$ pjsub --interact -g グループ名 -L rg=interactivea,elapse=01:00• 1 GPU実行$ pjsub --interact -g グループ名 -L rg=shareinteractive,elapse=01:00• インタラクティブ用のノードが全て使われている場合，資源が空くまでログインできません• 本講習会用のアカウントでは使えません
2023/6/29

講習会：GPU/C++

46

47

本お試し講習会でのキュー・グループ名• 本講習会中のキュー名• tutorial-a• 最大15分まで• 最大GPU数は4 GPUs（1ノード内）まで• 本講習会終了後のキュー名• lecture-a• 利用条件tutorial-aと同様• グループ名: gt00
2023/6/29

講習会：GPU/C++

48

計算結果の可視化例（1/2）• https://github.com/ymikirepo/nbody/tree/main/gallery/validation/fig

2023/6/29

講習会：GPU/C++

49

計算結果の可視化例（2/2）

2023/6/29

講習会：GPU/C++

50

OpenACCでの実装（データ転送部分）• Unified Memoryを使わない場合のみ必要• GPU上に置くべきデータ，必要なデータ転送を指定// メモリ確保#pragma acc enter data create(pos_ptr [0:num], vel_ptr [0:num], acc_ptr [0:num])// CPU → GPU のデータ転送#pragma acc update device(pos_ptr [0:num], vel_ptr [0:num])// GPU → CPU のデータ転送#pragma acc update host(acc_ptr [0:num])// メモリ解放#pragma acc exit data delete (pos_ptr [0:num], vel_ptr [0:num], acc_ptr [0:num])2023/6/29

講習会：GPU/C++

51

コンパイルと実行• $ ninja acc_data• $ pjsub -x EXEC=bin/acc_data
sh/wisteria/run_nvidia.sh• 出力されるファイル名を変更したければ，• $ pjsub -x EXEC=bin/acc_data,OPTION=“--file=NAME”
sh/wisteria/run_nvidia.sh• 正常に終了すると，下記ファイルが出力されるはずです• log/collapse_run.csv• エネルギー保存していることを確認（右端2列の数字が 0.005 より十分に小さい）• dat/collapse_snp[000 – 080].[h5 xdmf]

2023/6/29

講習会：GPU/C++

52

講習会概略
• 開催日： 2023年6月29日（木） 13:00--17:00• 形態： ZoomおよびSlackを用いたオンライン講習会• 使用システム： Wisteria/BDEC-01 (Aquarius)• 講習会プログラム：• 13:00--13:30• 13:40--14:20• 14:30--15:10• 15:20--15:50• 16:00--17:00

2023/6/29

GPU向けプログラミング手法の紹介（座学）
OpenACCを用いたGPU化（座学＋演習）
OpenMPを用いたGPU化（座学＋演習）
C++17を用いたGPU化（座学＋演習）
CUDA C++を用いたGPU化（座学＋演習）

講習会：GPU/C++

53

OpenMPを用いてのGPU化• OpenMP 4.0 以降ではアクセラレータへのオフロードをサポート• OpenMP 5.0 では loop 指示節が追加• OpenACC 的な実装ができる• スレッドブロック構造などは全てコンパイラにお任せ• 政治的にはOpenACCよりも優勢• NVIDIA, AMD, Intel が OpenMP でのGPU化をサポート• AMD, Intel は（きっと）OpenACCをサポートしない• Cray コンパイラであればAMD GPU向けのOpenACCもサポート• 自分にとってベンダーロックインはどこまで気にすべきか？
2023/6/29

講習会：GPU/C++

54

OpenMPでの実装（演算部分：distribute）• #pragma omp parallel for を置き換えていく• 性能的に特に重要なものについては，スレッド数の調整も• thread_limit（スレッド数） としてコンパイラに示唆（強制はできない）• NVIDIA HPC SDKではこの示唆が無視されて，128になっている印象• num_teams(チーム数) という方法もある（が，あまり使いたくはない）#pragma omp parallel forfor (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}#pragma omp target teams distribute parallel for
simd thread_limit(NTHREADS)
for (type::int_idx ii = 0U; ii < Ni; ii++) {//ループ内の実装は省略}#pragma omp parallel forfor (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}#pragma omp target teams distribute parallel forsimd
for (type::int_idx ii = 0U; ii < Ni; ii++) {//ループ内の実装は省略}
講習会：GPU/C++

2023/6/29

55

OpenMPでの実装（演算部分：loop）• #pragma omp parallel for を置き換えていく• （先述の）distributeよりも，コンパイラに多くを委ねる実装法• OpenMP 5.0 で導入された記法• スレッド数をコンパイラに示唆する方法はない• NVIDIA HPC SDKの場合は，こちらの方が速い印象（コンパイラが選択したスレッド数が最適であることが前提条件）#pragma omp parallel forfor (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}

2023/6/29#pragma omp target teams loopfor (type::int_idx ii = 0U; ii < Ni; ii++) {//ループ内の実装は省略}

講習会：GPU/C++

56

OpenMPでの実装（データ転送部分）• Unified Memoryを使わない場合のみ必要• GPU上に置くべきデータ，必要なデータ転送を指定// メモリ確保#pragma omp target enter data map(alloc : pos_ptr [0:num], vel_ptr [0:num], acc_ptr [0:num])// CPU → GPU のデータ転送#pragma omp target update to(pos_ptr [0:num], vel_ptr [0:num])// GPU → CPU のデータ転送#pragma omp target update from(acc_ptr [0:num])// メモリ解放#pragma omp target exit data map(delete : pos_ptr [0:num], vel_ptr [0:num], acc_ptr [0:num])2023/6/29

講習会：GPU/C++

57

OpenMPコードのコンパイル方法など• NVIDIA HPC SDK向けの情報• $ nvc++ -mp=gpu -gpu=cc80 -Minfo=accel,opt,mp• リンク時にも -mp=gpu を指定する• $ nvc++ -mp=gpu -gpu=cc80,managed -Minfo=accel,opt,mp• Unified Memory使用時のコンパイル方法• リンク時にも -mp=gpu -gpu=managed を指定する• 参考までにAMD製GPU（MI100）向けにコンパイルする場合• $ amdclang++ -target x86_64-pc-linux-gnu -fopenmp -fopenmptargets=amdgcn-amd-amdhsa -Xopenmp-target=amdgcn-amd-amdhsa-march=gfx908

2023/6/29

講習会：GPU/C++

58

コンパイルと実行• $ ninja omp_dist omp_dist_data omp_loop
omp_loop_data• $ pjsub -x EXEC=bin/omp_dist
sh/wisteria/run_nvidia.sh # 他の実行ファイルも同様• 出力されるファイル名を変更したければ，• $ pjsub -x EXEC=bin/omp_dist,OPTION=“--file=NAME”
sh/wisteria/run_nvidia.sh• 正常に終了すると，下記ファイルが出力されるはずです• log/collapse_run.csv• エネルギー保存していることを確認（右端2列の数字が 0.005 より十分に小さい）• dat/collapse_snp[000 – 080].[h5 xdmf]
2023/6/29

講習会：GPU/C++

59

講習会概略
• 開催日： 2023年6月29日（木） 13:00--17:00• 形態： ZoomおよびSlackを用いたオンライン講習会• 使用システム： Wisteria/BDEC-01 (Aquarius)• 講習会プログラム：• 13:00--13:30• 13:40--14:20• 14:30--15:10• 15:20--15:50• 16:00--17:00

2023/6/29

GPU向けプログラミング手法の紹介（座学）
OpenACCを用いたGPU化（座学＋演習）
OpenMPを用いたGPU化（座学＋演習）
C++17を用いたGPU化（座学＋演習）
CUDA C++を用いたGPU化（座学＋演習）

講習会：GPU/C++

60

C++17 Parallel Algorithms版の実装• #pragma omp parallel for を置き換えていく• GPU向けに特化した実装方法ではない• まったく同じコードをマルチコアCPU向けにコンパイルできる• CUDA実装を混ぜることもできるが，上記メリットが消えるので非推奨• 現時点では先進的（実験的）な部分もある取り組みなので，人柱精神の強い人でないならば様子見で良い（というのが個人的なコメント）#pragma omp parallel forfor (type::int_idx ii = 0U; ii < Ni; ii++) {// ループ内の実装は省略}

2023/6/29#include <algorithm>#include <execution>#include <boost/iterator/counting_iterator.hpp>

std::for_each_n(std::execution::par,boost::iterators::counting_iterator<type::int_idx>(
0U), Ni, [=](const type::int_idx ii) {// for文の中身はそのままでOK});

講習会：GPU/C++

61

標準言語での並列化に関する補足• 前ページの例では，counting_iterator をBoost C++ Librariesから呼んだが，もちろんこれ以外の実装も可能• 自分で実装してもOK（面倒だが..）• NVIDIA HPC SDKの場合は，thrustから呼んでもOK• コンパイラがNVIDIA縛りになるので，標準言語の機能を使う利点は減るが..

• 並列版実装が提供されているもの（std::sort, std::reduceなど）については，実行ポリシーを指定することで並列版アルゴリズムを実行• std::execution::seq （逐次処理）• std::execution::par （マルチスレッド化を許可）• std::execution::par_unseq （マルチスレッド化 and/or ベクトル化を許可）• スレッド数を示唆する方法，データ指示文的な記法は存在しない• この手法は，GPU実装に特化したプログラミング手法ではない• Unified Memory を使うことが前提となっている
2023/6/29

講習会：GPU/C++

62

標準言語でのコンパイル方法• NVIDIA HPC SDK向けの情報• $ nvc++ -stdpar=gpu -gpu=cc80 -Minfo=accel,opt,stdpar• リンク時にも -stdpar=gpu を指定する• $ nvc++ -stdpar=multicore -Minfo=accel,opt,stdpar• マルチコアCPU向けにコンパイルする場合

2023/6/29

講習会：GPU/C++

63

コンパイルと実行• $ ninja stdpar• $ pjsub -x EXEC=bin/stdpar
sh/wisteria/run_nvidia.sh• 出力されるファイル名を変更したければ，• $ pjsub -x EXEC=bin/stdpar,OPTION=“--file=NAME”
sh/wisteria/run_nvidia.sh• 正常に終了すると，下記ファイルが出力されるはずです• log/collapse_run.csv• エネルギー保存していることを確認（右端2列の数字が 0.005 より十分に小さい）• dat/collapse_snp[000 – 080].[h5 xdmf]

2023/6/29

講習会：GPU/C++

64

Fortranユーザ向けに• Fortran 2008 から，標準言語としての並列化がサポート• C++17 の場合と同様に，nvfortran ではGPU化の対象• do文をdo concurrentに直して並列化可能であると明記• C++17 の場合よりも書き換えコストはかなり低い• 各種制限などはC++とFortranで大きく違う• Unified Memoryが前提なのはC++17の場合と同じ• 参考資料• https://www.cc.u-tokyo.ac.jp/events/lectures/202/• OpenACC/OpenMP などについても説明されています

2023/6/29

講習会：GPU/C++

65

各実装手法間での性能比較をしてみる• $ ccmake -S ..
• BENCHMARK_MODE=ON に変更（当該行で Enter を入力）
• [c]onfigure ➔ [g]enerate• $ ninja• $ pjsub -x EXEC=bin/acc_data,OPTION=“--num_min=2097152
num_max=2097152 --num_bin=1 --file=benchmark”
sh/wisteria/run_nvidia.sh• 同様に omp_dist_data, omp_loop_data, stdpar も実行• 性能を確認・比較してみましょう• log/benchmark_run.csv の右端から3列目が演算性能• $ ccmake -S .. から RELAXED_RSQRT_ACCURACY=ON すると？
• $ ccmake -S .. から OVERWRITE_DEFAULT=ON にしたうえで NTHREADS の値を切り替えるとどうなる？
2023/6/29

講習会：GPU/C++

66

実習環境の理論ピーク演算性能• NVIDIA A100 (SXM, 40GB)• 108 SMs (32 FP64 cores per SM, 64 FP32 cores per SM)• GPU Boost Clock: 1410 MHz• 注：テンソルコア（FP64には対応）はN体では使えないので考慮しない• FP64: 9.7 TFlop/s• FP32: 19.5 TFlop/s• 理論ピーク性能はFMA（fused multiply-add）命令が前提なので，N体計算で実際にこの性能が出ることはない• 実習コードでの仮定（逆数平方根が4 Flops相当と仮定）• Leapfrog: 24 Flops per interaction• Hermite: 46 Flops per interaction
2023/6/29

講習会：GPU/C++

67

自分でも実装してみる• 対象のコード： cpp/base/0_base をコピーしただけのもの• OpenACC: cpp/openacc/exercise/nbody_*.cpp• OpenMP: cpp/openmp/exercise/nbody_*.cpp• 実装方針：1. まずは Unified Memory を前提にして計算部分だけGPU化
2. 次に，Unified Memory を外して実装• cpp/[openacc openmp]/exercise/CMakeLists.txt の中で，managed とある部分（2ヶ所）を編集• -gpu=cc${NVIDIA_GPU_HARD},managed ➔ -gpu=cc${NVIDIA_GPU_HARD}• -gpu=managed とある行を削除• $ ccmake -S .. から EXERCISE_MODE=ON して ninja• OpenACC: $ ninja acc_exercise• OpenMP: $ ninja omp_exercise

2023/6/29

講習会：GPU/C++

68

フローチャート• 右図はLeapfrog法の場合• Hermite法も大まかな流れは同様だが，時間積分が予測子・修正子法になっており，また粒子時刻によるソート処理も追加される• 初期条件生成，ファイル出力部分はGPU化しなくてOK• 注：最新の粒子データはGPU上にあるため，
スナップショット出力の直前にCPU上のデータを更新する作業は必要（データ指示文を使った実装の場合）• 性能測定モードにおいては，重力計算だ
けを繰り返す実装になっている
2023/6/29

講習会：GPU/C++

69

講習会概略
• 開催日： 2023年6月29日（木） 13:00--17:00• 形態： ZoomおよびSlackを用いたオンライン講習会• 使用システム： Wisteria/BDEC-01 (Aquarius)• 講習会プログラム：• 13:00--13:30• 13:40--14:20• 14:30--15:10• 15:20--15:50• 16:00--17:00

2023/6/29

GPU向けプログラミング手法の紹介（座学）
OpenACCを用いたGPU化（座学＋演習）
OpenMPを用いたGPU化（座学＋演習）
C++17を用いたGPU化（座学＋演習）
CUDA C++を用いたGPU化（座学＋演習）

講習会：GPU/C++

70

Compute Unified Device Architecture• NVIDIA社のGPUを対象としたプログラミング
環境• C++の拡張（昔はC/C++の拡張だったが…）• GPUの内部構造を（あまり）意識しないで良い• 多くのスレッドを生成• 多数のコアを有効に使うため• グローバルメモリへのアクセスレイテンシを隠蔽• ブロック：スレッドの集合• 典型的には128-512スレッド程度• シェアードメモリ，L1キャッシュを共有• 同期はいくつかの粒度で取れる• 最も遅いスレッドが全体を律速
2023/6/29

講習会：GPU/C++

CUDA C++ Programming Guide (Version 11.2)

71

CUDA C++で実装する際の記述項目• GPUの起動（cudaSetDevice()）• これが必要なのは複数GPUを触りに行く場合• 1 GPU / MPI プロセスであれば，CUDA_VISIBLE_DEVICESが便利• デバイスメモリの確保（cudaMalloc()）• データ転送用メモリの確保（cudaMallocHost()）• CPU⇔GPU間のデータ転送（cudaMemcpy()）• __global__関数の実装（CPUから起動するGPU関数）• __device__関数の実装（GPU関数から呼び出すGPU関数）• カーネル立ち上げ命令の追加（kernel<<<blck, thrd>>>()）• 確保したメモリの開放（cudaFree(), cudaFreeHost() ）2023/6/29

講習会：GPU/C++

72

簡易GPU化との性能比較• OpenACCでCUDAの8-9割を目指す or 達成とは言いますが，• CUDAで全力の最適化を施したコードと比較されていないことが多い• CUDAで実装しなおしました，という程度のコードとの比較が多い• たまにCUDAより速くなったという結果を見せられることもありますが，
CUDA側の最適化が足りていないんでしょ？ と思って聞いている• OpenACC/OpenMP/C++17 で実装したコードと，
CUDA C++で実装したコードとの性能差を確認• 今回の講習会ではCUDA C++全力版までは紹介しませんが，その少し手前ぐらいのコードとの比較まで行います

2023/6/29

講習会：GPU/C++

73

簡易GPU化との思想の違い• GPU上で実行される関数は非同期実行される• GPU上での実行完了を待たずにCPU上で次の命令が実行される• OpenACCなどでは，同期的に実行されていた• cudaDeviceSynchronize() によって完了が保証できる• 特に Unified Memory を使う場合には気をつけないと危険• GPU上で実行される関数などはCUDAストリームに紐付く• デフォルトストリームのみを使っている場合には，意識する必要なし• 同一ストリーム内では実行順序などの一貫性が保証される• cudaDeviceSynchronize() を入れて同期する必要はない• 複数のCUDAストリームを使って処理をオーバーラップすることも可能• cudaStreamSynchronize() などを使って適切に管理する

2023/6/29

講習会：GPU/C++

74

CUDA C++での実装（お手軽版：1/5）• まずはGPU上で動作させたい関数をGPU化• 関数定義の先頭に__global__をつける• GPU上の関数から呼ぶ関数については，__device__をつける• 一番外側のfor文を削除し，代わりに自動設定されるスレッドIDと紐づけ• if(ii < Ni){…} をつけないですむように，スレッド数の定数倍のメモリを確保• 余分に確保したメモリ領域には，質量0の粒子を置いておけば計算結果は正しい• （ツリー法の場合にはもう一工夫必要だが，細かい話なので省略）
void calc_acc(const type::int_idx Ni, const
type::position *const ipos, …) {#pragma omp parallel forfor (type::int_idx ii = 0U; ii < Ni; ii++) {// 関数の中身は省略}}
2023/6/29__global__ void calc_acc_device(const
type::position *const ipos, …) {const type::int_idx ii = blockIdx.x * blockDim.x+ threadIdx.x;// 関数の中身は省略}
講習会：GPU/C++

75

CUDA C++での実装（お手軽版：2/5）• GPU化した関数をCPUから起動する• スレッド数，問題サイズから必要なブロック数を設定（マクロ関数が便利）• スレッド数： NTHREADS• ブロック数： マクロ関数 BLOCKSIZE を使用• <<<ブロック数，スレッド数，動的確保するシェアードメモリ容量，ストリーム>>>• 後ろ2つは省略されることが多い（デフォルト設定をそのまま使用）
constexpr auto BLOCKSIZE(const type::int_idx num, const type::int_idx thread){ return (1U + ((num - 1U) / thread)); }static inline void calc_acc(const type::int_idx Ni, const type::position *const ipos,
type::acceleration *__restrict iacc, const type::int_idx Nj, const type::position*const jpos, const type::flt_pos eps2) {calc_acc_device<<<BLOCKSIZE(Ni, NTHREADS), NTHREADS>>>(ipos, iacc, Nj, jpos, eps2);}
2023/6/29

講習会：GPU/C++

76

CUDA C++での実装（お手軽版：3/5）• Unified Memory を使用する場合• メモリの確保・解放だけ記述すればOK• （CPU-GPU間のデータ転送は自分では何もしない）// 配列サイズを NTHREADS の整数倍にするための細工
auto size = static_cast<size_t>(num);if ((num % NTHREADS) != 0U) {size += static_cast<size_t>(NTHREADS - (num % NTHREADS));}// マネージドメモリの確保
cudaMallocManaged((void **)pos, size * sizeof(type::position));// マネージドメモリの解放
cudaFree(pos);

2023/6/29

講習会：GPU/C++

77

CUDA C++での実装（お手軽版：4/5）• Unified Memoryを使わない場合は，データ転送も記述// GPU上のメモリ確保
cudaMalloc((void **)pos_dev, size * sizeof(type::position));// CPU上の（pinned）メモリ確保（CPU・GPU間のデータ転送高速化のため）
cudaMallocHost((void **)pos_hst, size * sizeof(type::position));// 上記で確保したメモリの解放
cudaFree(pos_dev);
cudaFreeHost(pos_hst);// CPU → GPUのデータ転送
cudaMemcpy(pos_dev, pos_hst, num * sizeof(type::position), cudaMemcpyHostToDevice);// GPU → CPUのデータ転送
cudaMemcpy(acc_hst, acc_dev, num * sizeof(type::acceleration), cudaMemcpyDeviceToHost);
2023/6/29

講習会：GPU/C++

78

CUDA C++での実装（お手軽版：5/5）• もし必要があれば，cudaDeviceSynchronize()を追加• GPU上で関数を起動すると，関数の終了を待たずにCPUに処理が帰る• 複数のCUDAストリームを使った場合には，どこかで同期が必要• 性能測定時など，GPU上の関数の状態を把握すべき場合（下の例）• Unified Memoryを使用した際に，CPUから読み出したデータが不正だった場合
auto timer = util::timer();
cudaDeviceSynchronize();timer.start();calc_acc(num, pos_dev, acc_dev, num, pos_dev, eps2);
cudaDeviceSynchronize();timer.stop();2023/6/29

講習会：GPU/C++

79

CUDA C++でのコンパイル方法• $ nvcc -gencode arch=compute_80,code=sm_80 -Xptxas-v,-warn-spills,-warn-lmem-usage -lineinfo2023/6/29

講習会：GPU/C++

80

コンパイル方法• $ cd ../ # 大元の CMakelists.txt のあるフォルダに移動• $ module purge # load済みのmoduleを一旦全て外す• $ module load cmake cuda gcc hdf5
• $ cmake -DUSE_CUDA=ON –S . –B build_cuda• $ cd build_cuda• $ ccmake –S .. # デフォルトの設定を調整したい場合• 倍精度で計算したければ，FP_L, FP_Mを64に変更• Hermite法で試したい方は，HERMITE_SCHEMEをONに変更• $ ninja cuda_unified_base cuda_memcpy_base• 単に $ ninja だけでもOKですが，他にもたくさんのビルドが走ります2023/6/29

講習会：GPU/C++

81

実行• $ pjsub -x EXEC=bin/cuda_unified_base
sh/wisteria/run_cuda.sh• $ pjsub -x EXEC=bin/cuda_memcpy_base
sh/wisteria/run_cuda.sh• 出力されるファイル名を変更したければ，• $ pjsub -x EXEC=bin/cuda_unified_base,OPTION=“-file=NAME” sh/wisteria/run_cuda.sh• 正常に終了すると，下記ファイルが出力されるはずです• log/collapse_run.csv• エネルギー保存していることを確認（右端2列の数字が 0.005 より十分に小さい）• dat/collapse_snp[000 – 080].[h5 xdmf]
2023/6/29

講習会：GPU/C++

82

CUDA C++コードの最適化• 最適なスレッド数の選択• 後でスレッド数を変えるとどの程度性能が変わるか実験してみてください• できるだけ条件分岐を削除• これはCPUコードの場合も同じ• 高速な近似逆数平方根命令の使用• N体計算の場合は，rsqrtf()による高速化が最重要（NVIDIA GPU）
• NVHPCで RELAXED_RSQRT_ACCURACY=ON すると速くなる理由（nvc++のコンパイルオプションに -Mfprelaxed=rsqrt を追加）
• AMD GPUの場合には，__frsqrtf_rn()がベスト• シェアードメモリの活用• A100では，L2キャッシュの容量が増えた（6 MB➔40 MB）こともあり効果減（=最適化なしの状態でもかなり速くなった）• 命令レベルの並列性の導入（今回はここまではやりません）
2023/6/29

講習会：GPU/C++

83

コンパイルと実行• $ ninja cuda_unified_rsqrt cuda_memcpy_rsqrt• $ pjsub -x EXEC=bin/cuda_unified_rsqrt
sh/wisteria/run_cuda.sh # 他の実行ファイルも同様• 出力されるファイル名を変更したければ，• $ pjsub -x EXEC=bin/cuda_unified_rsqrt,OPTION=“-file=NAME” sh/wisteria/run_cuda.sh• 正常に終了すると，下記ファイルが出力されるはずです• log/collapse_run.csv• エネルギー保存していることを確認（右端2列の数字が 0.005 より十分に小さい）• dat/collapse_snp[000 – 080].[h5 xdmf]

2023/6/29

講習会：GPU/C++

84

シェアードメモリを使った実装の例• 静的に確保するには，GPU関数内で __shared__をつけて宣言しておけばOK• ブロック内の全スレッドから読み書き可能なので，シェアードメモリ上のデー
タを変更する際には__syncthreads() などをつけて制御する• 実行性能はループアンロー
リングの段数にも依存__global__ void calc_acc_device(…) {// 初期化部分省略__shared__ type::position jpos_shmem[NTHREADS];
for (type::int_idx jh = 0U; jh < Nj; jh += NTHREADS) {// load j-particle
const auto pj_tmp = jpos[jh + threadIdx.x];__syncthreads();
jpos_shmem[threadIdx.x] = pj_tmp;__syncthreads();

PRAGMA_UNROLL
for (type::int_idx jj = 0U; jj < NTHREADS; jj++) {const auto pj = jpos_shmem[jj];// 以下省略}}
iacc[ii] = ai;}2023/6/29

講習会：GPU/C++

85

コンパイルと実行• $ ninja cuda_unified_shmem cuda_memcpy_shmem• $ pjsub -x EXEC=bin/cuda_unified_shmem
sh/wisteria/run_cuda.sh # 他の実行ファイルも同様• 出力されるファイル名を変更したければ，• $ pjsub -x EXEC=bin/cuda_unified_shmem,OPTION=“-file=NAME” sh/wisteria/run_cuda.sh• 正常に終了すると，下記ファイルが出力されるはずです• log/collapse_run.csv• エネルギー保存していることを確認（右端2列の数字が 0.005 より十分に小さい）• dat/collapse_snp[000 – 080].[h5 xdmf]

2023/6/29

講習会：GPU/C++

86

各実装手法間での性能比較をしてみる• $ ccmake -S ..
• BENCHMARK_MODE=ON に変更（当該行で Enter を入力）
• [c]onfigure ➔ [g]enerate• $ ninja• $ pjsub -x EXEC=bin/cuda_memcpy_base,OPTION=“-num_min=2097152 num_max=2097152 --num_bin=1 -file=benchmark” sh/wisteria/run_cuda.sh• 同様に cuda_memcpy_rsqrt, cuda_memcpy_shmem も実行• 性能を確認・比較してみましょう• log/benchmark_run.csv の右端から3列目が演算性能• 簡易GPU化したコードとの性能比は？
• $ ccmake -S .. から OVERWRITE_DEFAULT=ON にしたうえで NTHREADS の値を切り替えるとどうなる？
2023/6/29

講習会：GPU/C++

87

CUDA版コードの性能の粒子数依存性• Wisteria/BDEC-01(Aquarius) 上の測定結果• NVIDIA A100 (SXM, 40GB)• CUDA 11.4
• Leapfrog法向け実装• 全てのSMが埋まるまでは，
性能は粒子数Nに比例• コア数の10倍程度の粒子数が必要• 13.6 TFlop/s (FP32) が最高• 全力版は 14.7 TFlop/s まで行く

2023/6/29

講習会：GPU/C++

88

自分でも実装してみる• 対象のコード： cpp/base/0_base をコピーしただけのもの• cpp/cuda/exercise/nbody_*.cu• 実装方針：1. まずは Unified Memory を前提にして計算部分だけGPU化
2. 次に，Unified Memory を外して実装• $ ccmake -S .. から EXERCISE_MODE=ON して make• $ ninja cuda_exercise

2023/6/29

講習会：GPU/C++

89

プロファイラなどの使い方• $ nsys profile --stats=true ./a.out• --stats=true をつけておくと，結果の概要も出力されるので便利• report?.nsys-rep というファイルが生成されるので，（手元で or Xを飛ばして） nsys-ui で開く• $ cuobjdump -sass ./a.out• ディスアセンブルして，発行されている命令を確認したくなった人向け• 時折 PTX の確認方法も紹介されるが，おすすめできません（PTX から実行ファイル生成までにもう一段コンパイルが入るため，
実際に発行されている命令とは違っている場合があります）

2023/6/29

講習会：GPU/C++90

より高度な内容• 以降はGPUを使い始める（た）段階の人向けにはハードルが高めの情報なので，演習などには含めず情報提供のみとしておきます• Volta世代以降のGPUでのワープの挙動• 最近のGPU上で（CUDAから）使える命令・機能• 複数GPU実装に向けての情報• GPUからの直接ファイルIO

2023/6/29

講習会：GPU/C++

91

Independent thread scheduling

2023/6/29

講習会：GPU/C++

92

Independent thread scheduling

2023/6/29

講習会：GPU/C++

93

暗黙の同期の挙動変更に対する処方箋• （暗黙の同期を使わないように）コードを書き換える• NVIDIA的にはこちらが推奨方針
1. ワープシャッフル命令などに _sync を追加（マスクは慎重に！）
2. 暗黙の同期点，if()文直後などに __syncwarp() を追加
3. ブロック同期用の __syncthreads() の挿入位置を再検討• 暗黙の同期を無理矢理使う（勝手にPascalモードと呼ぶ）• 実はこちらの方が速かったりする（今後どうなるかは不明）• Volta: コンパイル時に arch=compute_60,code=sm_70 を指定• Ampere: arch=compute_60,code=sm_80 を指定• Ampereから導入された warp-wide reduction は，Pascalモード
では使えない模様（コンパイルが通らなかった）
2023/6/29

講習会：GPU/C++

94

Warp shuffle命令• 同一ワープ（32スレッドの塊）内にある他スレッドのレジスタの値を（シェアードメモリなどを介さずに）取得できる• maskの値は基本的には 0xffffffff（=32スレッド全員）でOK• 複雑な実装をしている場合は，これではNGの場合もある• 同時に入ってくるスレッド数にデータ依存性があって予測できない場合には，__activemask()命令を使って動的に制御をかける

2023/6/29

講習会：GPU/C++

95

A100/CUDA 11から導入された機能• テンソルコア関連の機能• Sparsity: 小行列の中に0が入っていたらその分をスキップして高速化• TF32: 仮数部10ビット，指数部8ビットのデータ型（実態はFP19）• Asynchronous copy/barrier• グローバルメモリからシェアードメモリに直接（レジスタを介さずに）デー
タを置ける• Warp-wide reduction• ハードウェア的な高速化が効く• L2 cache residency control• （1つの連続な）データ領域をL2キャッシュに固定しておける

2023/6/29

講習会：GPU/C++

96

Asynchronous copy (memcpy_async)• グローバルメモリからシェアードメモリへのdirect memcpy• A100では，hardware accelerated• Svedin et al. (2021):

Svedin et al. 2021

• Pipeline API は Barrier API よりも高速• メモリ律速な問題では性能向上• 演算律速な問題では性能低下• N体問題は残念ながらこちらの場合• H100ではハード的な改良（TMA）が入る
ので，演算律速な場合でも性能低下しない
FLOP / Bytes
2023/6/29

講習会：GPU/C++

97

Warp-wide reduction• Throughputは16（INT32演算が64なので，4演算相当）• Warp shuffleが32（2演算相当）なので，5段（32スレッド）必要な旧実装に比べて（演算コストを無視しても）圧倒的に速い• __reduce_*_sync(unsigned mask, T value)• T: unsigned/intに対してはadd, min, max• 小細工すればfloatに対してmin/maxを返させることは可能• T: unsignedに対してはand, or, xor も可能• Supported by devices of compute capability 8.x or higher• Pascalモード（arch=compute_60,code=sm_80）の場合には，コンパイルが通らなかった（compute_80の指定が必要）• Pascalモード使用による高速化か，Ampereモード＋warp-wide reductionによる高速化か，の競争（tree構築はPascalモードの勝ち）

2023/6/29

講習会：GPU/C++

98

L2 cache residency control• Best Practice Guideからの抜粋"A portion of the L2 cache can be set aside for
persistent accesses to a data region in global
memory"• 1/16（= 2.5 MB）刻みで調整可能• White paper中の記述• CUDAストリームごとに1つの配列中の連続領域を指定• 1/16刻みで調整可能であれば，最大16個の領域を指定できても良さそうだが，残念ながらそういうAPIにはなっていない

2023/6/29

講習会：GPU/C++

99

複数GPU実装に関する情報提供• CUDA-aware MPI が使える環境であれば，実装は超簡単• MPI関数にGPU上のアドレスをそのまま指定すれば実装完了• CPU-GPU 間のデータ転送を自分で書く必要なし• GPUDirect のことを考えると，むしろ書かない方が良い• システムのMPIがGPUDirect RDMA（GDR）などをサポートしていれば，GPU間の（CPUを介しない）直接通信もできる• Open MPI (w/ UCX) や，MVAPICH2 GDR• Unified Memory を使っていると，GDR などは使ってくれない• 一度ホスト側にデータをコピーした上で通信されてしまう• 回避策もあるとのことだが，（現時点での） Unified Memory の弱点の一つ

2023/6/29

講習会：GPU/C++

100

環境変数 CUDA_VISIBLE_DEVICES• この環境変数を用いると，cudaSetDevices() などを使わずにどのGPUを用いるかを外から制御できる• AMD製GPUでも，ROCR_VISIBLE_DEVICESで同様の処理ができる• 例えばMPIプロセスあたりのGPU数を1とする（これがおすすめ）
場合には，• $ mpiexec –n 4 ./wrapper.sh ./a.out• wrapper.sh の中身(chmod +x をお忘れなく)：#!/bin/shexport LOCAL_ID=$OMPI_COMM_WORLD_LOCAL_RANK
export CUDA_VISIBLE_DEVICES=$LOCAL_ID$*• これはOpen MPI の場合の例

2023/6/29

講習会：GPU/C++101

GPUDirect Storage（GDS）• 通常，計算結果はGPU上• 通常のスナップショット生
成手順：
1. CPUにcudaMemcpy()
2. CPUがファイル書き込み• GDSが使えると，CPUを経由せずにGPUから直接
ファイルアクセスできる• PCIeをまたぐ回数が減る• HDF5 1.13+はGDSaware HDF5的に使えるので，実装コストも削減（VFDを使う）
2023/6/29

Memory

CPU1

CPU2

NVMe
SSD

NVMe
SSD

PCIe SW

PCIe SW

InfiniBand

GPUDirect
RDMA

講習会：GPU/C++

MemoryInfiniBand

GPU

GPU

GPU

GPU

GPUDirect
Storage102

最後に（1/2）• 今後GPUを搭載した計算機は増えていきそう• 通常のCPUに比べて消費電力あたりの演算性能が高いため• OFP-II（JCAHPC）はGPUを主としたスパコンになる予定• GPU向けプログラミング手法は多数ある• 性能，移植コスト，可搬性などを考慮して自分に適したものを選択• 今日紹介していないものではKokkos, RAJA などのフレームワークも• 実は最適化も（世間一般で言われているほどは）難しくない• GPU向けにアルゴリズムを設計とか言い出すと話は別だが，その場合はCPU向けであっても大変なはず• GPU移行に関するポータルサイトに情報を記載中• https://jcahpc.github.io/gpu_porting/2023/6/29

講習会：GPU/C++103

最後に（2/2）• アンケートの回答をお願いします• 講習会の改善のために有用な情報なので，ご協力お願いします• この講習会は今回が初めてなので，特にフィードバックが欲しいです• 本講習会アカウントは，7/29（土） 9:00まで使えます• キュー名はlecture-aです（tutorial-aは6/29の17:00以降使えなくなります）• 最大15分まで• 最大GPU数は4 GPUsまで

2023/6/29

講習会：GPU/C++