第75回お試しアカウント付き
並列プログラミング講習会「GPUプログラミング入門」
星野哲也 (助教) hoshino@cc.u-tokyo.ac.jp
2016年4月12日 (水)東京大学情報基盤センター

スケジュール• 09：30 - 10：00 受付• 10：00 - 12：00 Reedbush-Hログイン、GPU入門• 13：30 - 15：00 OpenACC入門• 15：15 - 18：00 OpenACC演習

2

事前準備
1.
2.

Reedbushにログインする（別資料を参照）※ログインノードは-U/-H共通
利用上の注意3

2. Reedbush 利用上の注意• ディレクトリの扱いについて– ログイン時のディレクトリ（/home/gt00h/txxxxx）にはログイン時に最低限必要なファイルのみを置く– 作業には/lustre以下のディレクトリ（/lustre/gt00h/txxxxx）を使う– /homeに置いたファイルは計算ノードから参照できない• コンパイルの為の準備– moduleコマンドを使ってコンパイルのための準備をする• モジュールをロード：module load 対象モジュール名– 環境変数PATHなどが書き換わり、コンパイラやライブラリが利用可能に– バージョンを指定する場合：module load 対象モジュール名/バージョン• 使用可能なモジュールの一覧を表示：module avail• 使用中のモジュールを確認：module list

詳細は午後説明します4moduleの切替• CUDA開発環境を使う場合– module load cuda

• PGIコンパイラを使う場合– （OpenACCやCUDA Fortranを使う場合）– module load pgi• Intelコンパイラを使う場合– module load intel• MPIを使う場合（コンパイラに追加してloadする）– module load mvapich2-gdr/2.2/{gnu,intel,pgi}– module load openmpi-gdr/2.0.2/{gnu,intel,pgi}• ジョブ実行時にも同じmoduleをloadすること• 組み合わせて使用しても良い– 順序に注意、環境変数PATHやLD_LIBRARY_PATHなどを確認する

詳細は午後説明します5

当センターの運用システム

6

東大センターのスパコン
2基の大型システム，6年サイクル
FY
08 09 10 11 12 13 14 15 16 17 18 19 20 21 22
Hitachi SR11K/J2
IBM Power-5+

18.8TFLOPS, 16.4TB

Yayoi: Hitachi SR16000/M1
IBM Power-7

メニーコア型大規模
スーパーコンピュータ

54.9 TFLOPS, 11.2 TB

Hitachi HA8000 (T2K)
AMD Opteron

JCAHPC:
筑波大・
東大

Oakforest-PACS
Fujitsu, Intel KNL

25PFLOPS, 919.3TB

140TFLOPS, 31.3TB

Oakleaf-FX: Fujitsu PRIMEHPC FX10,SPARC64 IXfx

BDEC system

50+ PFLOPS (?)

1.13 PFLOPS, 150 TB

Oakbridge-FX

136.2 TFLOPS, 18.4 TB
データ解析・シミュレーション
融合スーパーコンピュータ

Big Data &Extreme Computing

Reedbush, SGI
Broadwell + Pascal

1.93 PFLOPS
Reedbush-L
長時間ジョブ実行用演算加速装置
Pascal 1.4+PF付き並列スーパーコンピュータ

Peta

K

K computer

Post-K ?

7

5システム運用中• Yayoi (日立 SR16000, IBM Power7)– 54.9 TF, 2011年11月 〜 2017年10月• Oakleaf-FX (富士通 PRIMEHPC FX10)– 1.135 PF, 京コンピュータ商用版, 2012年4月 〜 2018年3月• Oakbridge-FX (富士通 PRIMEHPC FX10)– 136.2 TF, 長時間実行用（168時間）, 2014年4月 〜 2018年3月• Reedbush (SGI, Intel BDW + NVIDIA P100 (Pascal))– データ解析・シミュレーション融合スーパーコンピュータ– 1.93 PF, 2016年7月〜 2020年6月– 東大ITC初のGPUシステム (2017年3月より), DDN IME (Burst Buffer)• Oakforest-PACS (OFP) (富士通、Intel Xeon Phi (KNL))– JCAHPC (筑波大CCS＆東大ITC)– 25 PF, TOP 500で6位 (2016年11月) (日本で1位)– Omni-Path アーキテクチャ, DDN IME (Burst Buffer)8東京大学情報基盤センター
スパコン（1/3）
Fujitsu PRIMEHPC FX10 （FX10スーパーコンピュータシステム）
Total Peak performance: 1.13 PFLOPS
Total number of nodes: 4,800
Total memory: 150TB
Peak performance per node : 236.5 GFLOPSMain memory per node
: 32 GBDisk capacity: 2.1 PBSPARC64 IXfx 1.848GHz

2012年7月~2018年3月（予定）

Oakbridge-FX：長時間ジョブ用のFX10
ノード数：24～576
制限時間：最大168時間（１週間）
9東京大学情報基盤センター
スパコン（2/3-1）
Reedbush-U （SGI Rackable クラスタシステム ）Total Peak performance: 508 TFLOPS
Total number of nodes: 420
Total memory: 105 TB
Peak performance per node : 1209.6 GFLOPSMain memory per node
: 256 GBDisk capacity: 5.04 PB
File Cache system (SSD): 230 TBIntel Xeon E5-2695v4 2.1GHz 18 core x2 socket

2016年7月1日試験運転開始
2016年9月1日正式運用開始

10東京大学情報基盤センター
スパコン（2/3-2）
Reedbush-H （SGI Rackable クラスタシステム ）Total Peak performance: 145 TFLOPS+ 1272 TFLOPS
Total number of nodes: 120
Total memory: 30 TB + 8 TB
Peak performance per node : 1209.6 GFLOPS+ 10.6 TFLOPSMain memory per node
: 256 GB+ 32 GBDisk capacity (shared w/ U) : 5.04 PB
File Cache system (SSD, (shared w/ U): 230 TBIntel Xeon E5-2695v4 2.1GHz 18 core x2 socket+ NVIDIA Tesla P100 with NVLink x 2

2017年3月1日試験運転開始
2017年4月1日正式運用開始

11

12

東京大学情報基盤センター
スパコン（3/3）
筑波大学計算科学研究センターと共同運用
Oakforest-PACS （Fujitsu PRIMERGY CX600）
Total Peak performance: 25 PFLOPS
Total number of nodes: 8,208
Total memory: 897.7 TB
Peak performance per node : 3.046 TFLOPSMain memory per node
: 96 GB (DDR4)+ 16 GB(MCDRAM)
Disk capacity: 26.2 PB
File Cache system (SSD): 960 TBIntel Xeon Phi 7250 1.4 GHz 68 core x1 socket

2016年12月1日試験運転開始
2017年4月3日正式運用開始

12

FX10TOFU
Network





Core
#0

Core
#1

Core
#2

Core
#3

L1

L1

L1

L1





…
: L1

 

20GB/
Core
#12

Core
#13

Core
#14

Core
#15

32KB L1

L1

L1

L1

ICC

L2 (16  12MB)

85GB/
=(8Byte×1333MHz
×8 channel)

MemoryMemoryMemoryMemory4GB ×2

4GB ×2

4GB ×2

4GB ×2 

8GB×432GB
GPUプログラミング入門

DDR3 DIMM

13

FX10の通信網
TOFU



1TOFU


GPUプログラミング入門

6
5GB/


14

FX10の通信網（１TOFU単位間の結合）• ユーザから見ると、X軸、Y軸、Z軸について、奥の1TOFUと、手前の１TOFU
は、繋がってみえます（３次元トーラス接続）• ただし物理結線では– X軸はトーラス– Y軸はメッシュ– Z軸はメッシュまたは、
トーラスになっていますGPUプログラミング入門

15

Reedbush-Uノードのブロック図• メモリのうち、「近い」メモリと「遠い」メモリがある=> NUMA (Non-Uniform Memory Access)(FX10はフラット)


128GB

DDR4
DDR4
DDR4
DDR4

G3 x16

76.8GB/s

76.8GB/s
QPIIntel Xeon
E5-2695 v4(BroadwellEP)

QPI
76.8GB/sIntel Xeon
E5-2695 v4(Broadwell-EP)

DDR4
DDR4
DDR4
DDR4


128GB

76.8GB/s

15.7 GB/s

IB EDRHCA

GPUプログラミング入門

16

Broadwell-EP

Core
#0

L L
1 2

L3

Core
#5

L L
1 2

Core
#1

L L
1 2

L3

Core
#6

L L
1 2

Core
#2

L L
1 2

L3

Core
#7

L L
1 2

Core
#3

L L
1 2

L3

Core
#8

L L
1 2

Core
#4

L L
1 2

L3

DDR4
DIMM

 L1: 2KB, L2: 256KB,
L3: 2.5MB() => L3 45MB

PCIe

QPI x2

MemoryMemory16GB ×2

16GB ×2

L3

Core
#9

L L
1 2

L3

Core
#14

L L
1 2

L3

L3

Core
#10

L L
1 2

L3

Core
#15

L L
1 2

L3

L3

Core
#11

L L
1 2

L3

Core
#16

L L
1 2

L3

L3

Core
#12

L L
1 2

L3

Core
#17

L L
1 2

L3

Core
#13

L L
1 2

L3

  !16GB 8"128GB
GPUプログラミング入門

MemoryMemory16GB ×2
16GB ×2
76.8 GB/=(8Byte 2400MHz 4 channel)
17

Reedbush-Uの通信網• フルバイセクションバンド幅を持つFat Tree網– どのように計算ノードを選んでも互いに無衝突で通信が可能• Mellanox InfiniBand EDR 4x CS7500: 648ポート– 内部は36ポートスイッチ (SB7800)を (36+18)台組み合わせたものと等
価
648

Director 1
...

36
Spine
18

...

36
Leaf 
36

Uplink: 18

Downlink: 18
1

...

18

19

...

36

37

...

54

GPUプログラミング入門

18

Reedbush-Hノードのブロック図

128GB

DDR4
DDR4
DDR4
DDR4

Intel Xeon
E5-2695 v4(Broadwell-EP)

QPI
76.8GB/s

15.7 GB/s

15.7 GB/s


128GB

76.8GB/s

G3
x16

G3
x16

NVIDIA
Pascal

DDR4
DDR4
DDR4
DDR4

PCIe sw

PCIe sw

IB FDR
HCA

G3
x16

G3x
16

76.8GB/s

76.8GB/s
QPIIntel Xeon
E5-2695 v4(BroadwellEP)

20 GB/s
NVLinK
NVLinKNVIDIA
PascalIB FDR
HCA

20 GB/s
EDR switch
EDR
GPUプログラミング入門

19

Oakforest-PACS 計算ノード

Knights Landing: Next Intel® Xeon Phi™ Processor
Intel® Many-Core Processor targeted for HPC and SupercomputingFirst self-boot Intel® Xeon Phi™ processor that is binary
compatible with main line IA. Boots standard OS.

Significant improvement in scalar and vector performance• Intel Xeon Phi (Knights Landing)Integration of Memory on package: innovative memoryarchitecture for high bandwidth and high capacity

Integration of Fabric on package– 1ノード1ソケット
Three products• MCDRAM: オンパッケージの高バンド幅メモリ16GB+ DDR4メモリ
KNL Self-Boot(Baseline)

HotChips27
Knights Landing
Overview
KNLスライドより

KNL Self-Boot w/ Fabric

KNL Card(Fabric Integrated)(PCIe-Card)

MCDRAMMCDRAM

2 x16
1 x4

X4
DMI

PCIe
Gen 3

D
M
I

MCDRAMMCDRAM

Potential future options subject to change without notice.
All timeframes, features, products and dates are preliminary forecasts and subject to change without further notification.

EDC

EDC

3

ソケット当たりメモリ量：16GB×6＝96GB

2 VPU
Core

1MB
L2

C
H
A
N
N
E
L
S

3

EDC

Core

MCDRAMDDR
4

36 Tilesconnected by
2D MeshInterconnect

DDR MC

2 VPU

4

EDC

Tile

D
D
R
4

MCDRAM: 490GB/秒以上 （実測）DDR4: 115.2 GB/秒=(8Byte×2400MHz×6 channel)

EDC

EDC

MCDRAM

misc

Package

DDR MCEDC

MCDRAM

C
H
A
N
N
E
L
S

EDC

MCDRAM

Omni-path not shown

GPUプログラミング入門

20

Oakforest-PACS: Intel Omni-Path Architecture による
フルバイセクションバンド幅Fat-tree網
768 port Director
Switch12(Source by Intel)

22Uplink: 24

48 port Edge Switch
362 

Downlink: 24
1

...

24

25

...

48

49

...

72





コストはかかるがフルバイセクションバンド幅を維持• システム全系使用時にも高い並列性能を実現• 柔軟な運用：ジョブに対する計算ノード割り当ての自由度が高い
GPUプログラミング入門

21

東大情報基盤センターOakforest-PACSスーパーコン
ピュータシステムの料金表（2017年4月1日 ）• パーソナルコース（年間）– コース１： 100,000円 ： 8ノード(基準)、最大16ノードまで– コース２： 200,000円 ： 16ノード(基準)、最大64ノードまで• グループコース– 400,000円 (企業 480,000円) ： 1口 8ノード（基準）、最大128ノードまで– 以上は、「トークン制」で運営• 申し込みノード数×360日×24時間の「トークン」が与えられる• 基準ノードまでは、トークン消費係数が1.0• 基準ノードを超えると、超えた分は、消費係数が2.0になる• 大学等のユーザはFX10、Reedbushとの相互トークン移行も可能

GPUプログラミング入門

22

東大情報基盤センターReedbushスーパーコンピュータ
システムの料金表（2017年4月1日 ）• パーソナルコース（年間）– 150,000円 ： RB-U: 4ノード（基準）、最大16ノードまで
RB-H: 1ノード（基準）、最大2ノードまで• グループコース– 300,000円： 1口 4ノード（基準）、最大128ノードまで、
RB-H： 1ノード（基準）、最大32ノードまで（トークン係数はUの2.5倍）– RB-Uのみ 企業 360,000円 ： 1口 4ノード（基準）、最大128ノードまで– RB-Hのみ 企業 216,000円 ： 1口 1ノード（基準）、最大32ノードまで– 以上は、「トークン制」で運営• 申し込みノード数×360日×24時間の「トークン」が与えられる• 基準ノードまでは、トークン消費係数が1.0• 基準ノードを超えると、超えた分は、消費係数が2.0になる• 大学等のユーザはFX10, Oakforest-PACSとの相互トークン移行も可能• ノード固定もあり
GPUプログラミング入門

23

東大情報基盤センターFX10スーパーコンピュータシス
テムの料金表（2017年4月1日）• パーソナルコース（年間）– コース１：90,000円 ： 12ノード(基準)、最大24ノードまで– コース２： 180,000円 ： 24ノード(基準)、最大96ノードまで• グループコース– 360,000円 (企業 432,000円) ：1口、12ノード、最大1440ノードまで– 以上は、「トークン制」で運営• 申し込みノード数×360日×24時間の「トークン」が与えられる• 基準ノードまでは、トークン消費係数が1.0• 基準ノードを超えると、超えた分は、消費係数が2.0になる• 大学等のユーザはReedbush, Oakforest-PACSとの相互トークン移行も可能
GPUプログラミング入門

24

トライアルユース制度について•安価に当センターのOakleaf/Oakbridge-FX, Reedbush-U/H, Oakforest-PACSシ
ステムが使える「無償トライアルユース」および「有償トライアルユース」制度があります。
– アカデミック利用• パーソナルコース、グループコースの双方（1ヶ月～3ヶ月）– 企業利用• パーソナルコース（１ヶ月～３ヶ月）（FX10: 最大24ノード、最大96ノード,
RB-U: 最大16ノード, RB-H: 最大2ノード、OFP: 最大16ノード, 最大64ノード）
講習会いずれかの受講が必須、審査無• グループコース– 無償トライアルユース：（１ヶ月～３ヶ月）：無料（FX10:最大1,440ノード、
RB-U: 最大128ノード、RB-H: 最大32ノード、OFP: 最大2048ノード）– 有償トライアルユース：（１ヶ月～最大通算９ヶ月）、有償（計算資源は無償と同等）– スーパーコンピュータ利用資格者審査委員会の審査が必要（年2回実
施）• 双方のコースともに、簡易な利用報告書の提出が必要
25

スーパーコンピュータシステムの詳細• 以下のページをご参照ください– 利用申請方法– 運営体系– 料金体系– 利用の手引などがご覧になれます。
http://www.cc.u-tokyo.ac.jp/system/ofp/
http://www.cc.u-tokyo.ac.jp/system/reedbush/
http://www.cc.u-tokyo.ac.jp/system/fx10/
GPUプログラミング入門

26

GPU入門「GPUプログラミング入門」講習会の目的
GPUのアーキテクチャ
GPUの使い方27

GPUって何？
• GPU : Graphics Processing Unit
– いわゆるグラフィックボード(グラボ)、ビデオカード• 画像処理に特化したハードウェア– 高速・高解像度描画、3D描画処理（透視変換、陰影・照明）、画
面出力– ノートPC、ゲームハード、スマホなどに搭載• GPUを汎用計算に用いる手法をGPGPU、GPUコンピューティ
ングなどと呼ぶGPUプログラミング入門

なぜGPUコンピューティング？
• 性能が高いから！
P100

BDW

KNL

動作周波数（GHz）

1.480

2.10

1.40

コア数（有効スレッド数）

3,584

18 (18)68 (272)

理論演算性能（GFLOPS）

5,304

604.8

3,046.4

主記憶容量（GB）

16

128

16

メモリバンド幅（GB/sec., Stream Triad）

534

65.5

490

Reedbush-HのGPU

Reedbush-U/HのCPU

Oakforest-PACSのCPU (Intel Xeon Phi)

備考

GPUプログラミング入門

29

GPUの応用例

http://www.nvidia.com/object/gpu-applications.html (英語 情報量多い)
http://www.nvidia.co.jp/object/gpu-applications-domain-jp.html (日本語)• 数値流体力学• 数値分析
• 医療画像• 天気と気候• データサイエンス• 機械学習– AlphaGo の学習にもGPU

線形代数演算が速いので、広い範囲で応用可能！
GPUプログラミング入門

30

GPUプログラミングは何が難しい？
• CPU： 大きなコアをいくつか搭載– Reedbush-H の CPU : 2.10 GHz 18コア
– 大きなコア… 分岐予測、パイプライン処理、Out-of-Order• 要はなんでもできる– 逐次の処理が得意• GPU： 小さなコアをたくさん搭載– Reedbush-H の GPU： 1.48 GHz 3,584 コア– 小さなコア... 上記機能が弱いまたはない！
– 並列処理が必須

GPUの難しさ
1. 多数のコアを効率良く扱う難しさ
2. 並列プログラミング自体の難しさ
GPUプログラミング入門

31

並列プログラミングとは何か？
• 逐次実行のプログラム（実行時間T ）を、ｐ台の計算機を使って、T/ ｐ にすること。

T
T

/ 

• 素人考えでは自明。
• 実際は、できるかどうかは、対象処理の内容（アルゴリズム）で 大きく 難しさが違う– アルゴリズム上、絶対に並列化できない部分の存在– 通信のなどのオーバヘッドの存在• 通信立ち上がり時間• データ転送時間32

アムダールの法則
並列化できない部分（１ブロック） 並列化できる部分（８ブロック）● 

88.8%

●



9/3=3

●

9/2=4.5 ≠6

33

並列処理の実行形態：タスク並列• タスク（ジョブ）を分割することで並列化する。
• データの操作（＝演算）は異なるかもしれない。
• タスク並列の例：カレーを作る– 仕事１：野菜を切る– 仕事２：肉を切る– 仕事３：水を沸騰させる– 仕事４：野菜・肉を入れて煮込む– 仕事５：カレールゥを入れる●並列化

GPUは苦手

仕事１
仕事２
仕事３

仕事４

仕事５

時間
34

並列処理の実行形態：データ並列• データを分割することで並列化する。
• データの操作（＝演算）は同一となる。
• データ並列の例：行列－行列積

GPUでは
基本これ！

æ 1 2 3 ö æ 9 8 7 ö æç 1* 9 + 2 * 6 + 3 * 3 1* 8 + 2 * 5 + 3 * 2 1* 7 + 2 * 4 + 3 *1 ö÷ç÷ ç÷ ＝ ç 4 * 9 + 5 * 6 + 6 * 3 4 * 8 + 5 * 5 + 6 * 2 4 * 7 + 5 * 4 + 6 *1÷
ç 4 5 6 ÷ ç 6 5 4 ÷ ç 7 * 9 + 8 * 6 + 9 * 3 7 * 8 + 8 * 5 + 9 * 2 7 * 7 + 8 * 4 + 9 *1÷
øç7 8 9÷ ç 3 2 1÷ è
è
ø è
ø
●並列化
CPU０ æ 1

2 3ö
ç÷
CPU１ ç 4 5 6 ÷
CPU２ ç 7 8 9 ÷èø

全CPUで共有

æ9 8 7öæ 1* 9 + 2 * 6 + 3 * 3 1* 8 + 2 * 5 + 3 * 2 1* 7 + 2 * 4 + 3 *1 öç
÷ ＝ çç 4 * 9 + 5 * 6 + 6 * 3 4 * 8 + 5 * 5 + 6 * 2 4 * 7 + 5 * 4 + 6 *1÷÷
ç 6 5 4÷
ç 7 * 9 + 8 * 6 + 9 * 3 7 * 8 + 8 * 5 + 9 * 2 7 * 7 + 8 * 4 + 9 *1÷
èø
ç3 2 1÷èø 並列に計算：初期データは異なるが演算は同一
35

本講習会の目的
並列アルゴリズムはひとまず置いといて…

GPUの多数のコアを効率良く扱う手法を学ぼう！
• その他の講習会 (予定) http://www.cc.u-tokyo.ac.jp/support/kosyu/– 並列アルゴリズムを学びたい方• 7/11,12 OpenMP/OpenACCによるマルチコア・メニィコア並列プログラ
ミング入門– KNL について学びたい方• 5/15

KNL 講習会(仮称)– その他• 5/12
OpenFOAM入門• 5/23,24 「並列有限要素法とハイブリッド並列プログラミング」• 6/1,2 MPI 基礎• 6/26,27 MPI 上級編• 7/18,19 GPUクラスタ講習会(仮称)「こんな講習会が欲しい！」GPUプログラミング入門という要望があればぜひア
ンケートへ！

36

参考：NVIDIA Tesla P100
• 56 SMs• 3584 CUDA
Cores
• 16 GB HBM2

http://pc.watch.impress.co.jp/docs/column/kaigai/751980.html より
日本語によるP100の解説では一番詳しい
GPUプログラミング入門

37

参考：NVIDIA Tesla P100 の SM

http://pc.watch.impress.co.jp/docs/column/kaigai/751980.html より
日本語によるP100の解説では一番詳しい
GPUプログラミング入門

38

参考：NVIDIA社のGPU• 製品シリーズ– GeForce• コンシューマ向け。安価。

– Tesla• HPC向け。倍精度演算器、大容量メモリ、ECCを備えるため高価。

• アーキテクチャ（世代）1.
2.

Tesla：最初のHPC向けGPU、TSUBAME1.2など
Fermi：2世代目、TSUBAME2.0など• ECCメモリ、FMA演算、L1 L2 キャッシュ3.

Kepler：現在HPCにて多く利用、TSUBAME2.5など• シャッフル命令、Dynamic Parallelism、Hyper-Q

4.
5.

Maxwell：コンシューマ向けのみ
Pascal：最新GPU、Reedbush-Hに搭載• HBM2、半精度演算、 NVLink、倍精度atomicAdd など

6.

Volta：次世代GPU39

押さえておくべきGPUの特徴• CPUと独立のGPUメモリ• 性能を出すためにはスレッド数>>コア数• 階層的スレッド管理と同期• Warp 単位の実行• やってはいけないWarp内分岐• コアレスドアクセス

GPUプログラミング入門

40

CPUと独立のGPUメモリ
１．必要なデータを送る
ノードの外へ
~20GB/s

バス (PCIe など)
CPU
OSが動いている~32GB/sGPU
OSは存在しない

２．計算を行う

３．計算結果を返す~200GB/s~1,000GB/s
デバイス
メモリ

メインメモリ• 計算はCPUから始まる• 物理的に独立のデバイスメモリとデータのやり取り必須
GPUプログラミング入門

41

性能を出すためにはスレッド数>>コア数• 推奨スレッド数– CPU： スレッド数=コア数 (高々数十スレッド)– GPU： スレッド数>=コア数*4~ (数万~数百万スレッド)• 最適値は他のリソースとの兼ね合いによる• 理由：高速コンテキストスイッチによるメモリレイテンシ隠し– CPU : レジスタ・スタックの退避はOSがソフトウェアで行う(遅い)– GPU : ハードウェアサポートでコストほぼゼロ• メモリアクセスによる暇な時間(ストール)に他のスレッドを実行
メモリread開始

メモリread終了

1core=1スレッドのとき
1core=Nスレッドのとき
GPUプログラミング入門

42

階層的スレッド管理と同期• コアの管理に対応

値はP100の場合– 1 SM の中に 64 CUDA core、56 SM で 3584 CUDA core
– 1 CUDA core が複数スレッドを担当• スレッド間の同期– 同一SM内のスレッド間は同期できる• 正確には同一スレッドブロック内– 異なるSMのスレッド間は同期できない• 同期するためにはGPUの処理を終了する必要あり• atomic 演算は可能• メモリ資源の共有– L1 cache, shared memory, Instruction cache などはSM内で共有– L2 cache, Device memory などは全スレッドで共有
GPUプログラミング入門

43

Warp 単位の実行• 連続した32スレッドを1単位 = Warp と呼ぶ• このWarpは足並み揃えて動く– 実行する命令は32スレッド全て同じ– データは違ってもいい
スレッド

1

2

3

…

31 32

スレッド

1

2

3

…

31 32

配列 A

4

3

5

…

8

配列 A

4

3

5

…

8

0

÷ × ＋ …

−

×

2

1

9

配列 B

0

× × × …

× ×

2

1

3

1

…

9

配列 B

OK！

3

1

…

NG！
GPUプログラミング入門44

やってはいけないWarp内分岐• Divergent Branch– Warp 内で分岐すること。Warp単位の分岐ならOK。

::
if ( 奇数スレッド ) {::} else {
:
:
}
:
:

:
:
if ( TRUE ) {::} else {::}::else 部分は実行せずジャンプ

一部スレッドを眠らせて全分岐を実行
45
最悪ケースでは32倍のコスト

GPUプログラミング入門

コアレスドアクセス• 同じWarp内のスレッドが近いアドレスに同時にアクセスするのがメモリの性質上効率的– これをコアレスドアクセス(coalesced access)と呼ぶ
スレッド

1 2 3 4 … 32…

メモリアクセスが1回で済む。
デバイスメモリスレッド

1 2 3 4 … 32

32回のメモリアクセスが行われる
128バイト単位でメモリアクセス。Warp内のアクセスが128バイトに収まってれば1回。外れればその分だけ繰り返す。最悪ケース
では32倍のコスト
GPUプログラミング入門

46

OpenACC入門
OpenACC とは
OpenACC の指示文
OpenACC の実用例
OpenACC 演習

47

GPUを利用する方法
簡単• GPU対応ライブラリ– ライブラリを呼び出すだけ– ライブラリとして存在するもののみ• OpenACC– 既存のC/C++/Fortranプログラムに指示文を入れるだけ– 細かいチューニングは不可• CUDA

難しい

本講習会の主な対象– 自由度が最も高い– プログラムの大幅な書き換えを必要とする
48

OpenACC• OpenACCとは… アクセラレータ(GPUなど)向けのOpenMPのよ
うなもの– 既存のプログラムのホットスポットに指示文を挿入し、計算の重たい部分をアクセラレータにオフロード– 対応言語： C/C++, Fortran• 指示文ベース– 指示文：コンパイラへのヒント– 記述が簡便, メンテナンスなどをしやすい– コードの可搬性(portability)が高い• 対応していない環境では無視される
GPUプログラミング入門

C/C++#pragma acc kernels
for(i = 0;i < N;i++) {….
}Fortran!$acc kernels
do i = 1, N….
end do!$acc end kernels

49OpenACC• 規格– 各コンパイラベンダ(PGI, Crayなど)が独自に実装していた拡張を統合し共通規格化 (http://www.openacc.org/)– 2011年秋にOpenACC 1.0 最新の仕様はOpenACC 2.5• 対応コンパイラ– 商用：PGI, Cray, PathScale– 研究用：Omni (AICS), OpenARC (ORNL), OpenUH (U.Houston)– フリー：GCC 6.x• 開発中 (開発状況： https://gcc.gnu.org/wiki/Offloading)• 実用にはまだ遠い
本講習会ではPGIコンパイラを用いる
GPUプログラミング入門

50

OpenACC と OpenMP の実行イメージ
1スレッド

OpenMP

OpenACC

int main() {#pragma …
for(i = 0;i < N;i++) {
…
}

}

CPU

CPU
GPUプログラミング入門

CPU

デバイス
51

OpenACC と OpenMP の比較
OpenMPの想定アーキテクチャ
マルチコアCPU環境• 計算コアがN個

MEMORY– N < 100 程度 (Xeon Phi除く)
CPU(s)
計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

• 共有メモリ

計算コ
ア

一番の違いは対象アーキテクチャの複雑さ

52

OpenACC と OpenMP の比較
OpenACC の想定アーキテクチャ
アクセラレータを備えた計算機環境
MEMORY(ホスト)

CPU(s)MEMORY(デバイス)• 計算コアN個をM階層で管理– N > 1000 を想定– 階層数Mはアクセラレータによる• ホスト-デバイスで独立したメモリ– ホスト-デバイス間データ転送は低速一番の違いは対象アーキテクチャの複雑さ

53

OpenACC と OpenMP の比較• OpenMPと同じもの– Fork-Joinという概念に基づく並列化• OpenMPになくてOpenACCにあるもの– ホストとデバイスという概念• ホスト-デバイス間のデータ転送– 多階層の並列処理• OpenMPにあってOpenACCにないもの– スレッドIDを用いた処理など• OpenMPのomp_get_thread_num()に相当するものが無い• その他、気をつけるべき違い– OpenMPと比べてOpenACCは勝手に行うことが多い• 転送データ、並列度などを未指定の場合は勝手に決定
54

OpenACC と CUDA の違い
OpenACC

CUDA• 指示文ベース• 対象：アクセラレータ全般• 記述の自由度• 言語拡張• 対象：NVIDIA GPU のみ• 記述の自由度– 高レベルな抽象化– ある程度勝手にやってくれる– 低レベルな記述– 書いたようにしかやらない– デバイスに特化した機能は使えない– デバイスの持つ性能を十分に引き出せる• shuffle機能を使えない など

プログラムの可搬性◎
可読性◎
GPUプログラミング入門

性能◎

55

最低限動くプログラムを作るにはOpenACC

attribute(global) subroutine mm_cuda(a, b, c, n)
integer, value :: nreal(8), dimension(n, n) :: a, b, c
integer :: i, j, kreal(8) :: cci = (blockidx%x-1) * blockdim%x + threadidx%x
j = (blockidx%y-1) * blockdim%y + threadidx%y
cc = 0.0
do k = 1, n
cc = cc + a(i, k) * b(k, j)
end doc(i, j) = cc
CUDA
end subroutine mm_cuda

subroutine matmul(a, b, c, n)real(8), dimension(n, n) :: a, b, c
integer :: ninteger :: i, j, kreal(8) :: cc

subroutine matmul(a, b, c, n)real(8), dimension(n, n) :: a, b, c
integer :: ninteger :: i, j, kreal(8) :: cc

do j = 1, n
do i = 1, n
cc = 0do k = 1, n
cc = cc + a(i,k) * b(k,j)
end doc(i,j) = cc
end doend dodo j = 1, n
do i = 1, n
cc = 0do k = 1, n
cc = cc + a(i,k) * b(k,j)
end doc(i,j) = cc
end doend do1. オフロードする領域を決める

end subroutine matmul

end subroutine matmul

56

最低限動くプログラムを作るには
2. オフロード領域の並列化、カーネルコードの記述
OpenACC

CUDA
attribute(global) subroutine mm_cuda(a, b, c, n)
integer, value :: nreal(8), dimension(n, n) :: a, b, c
integer :: i, j, kreal(8) :: cci = (blockidx%x-1) * blockdim%x + threadidx%x
j = (blockidx%y-1) * blockdim%y + threadidx%y
cc = 0.0
do k = 1, n
cc = cc + a(i, k) * b(k, j)
end doc(i, j) = cc
end subroutine mm_cuda

subroutine matmul(a, b, c, n)subroutine matmul(a, b, c, n)real(8), dimension(n, n) :: a, b, creal(8), dimension(n, n) :: a, b, c
integer :: ninteger :: ninteger :: i, j, k
integer :: i, j, kreal(8) :: ccreal(8) :: cc!$acc kernels
do j = 1, n
do j = 1, n
do i = 1, n
do i = 1, n
cc = 0cc = 0do k = 1, n
do k = 1, n
cc = cc + a(i,k) * b(k,j)
cc = cc + a(i,k) * b(k,j)
end doend doc(i,j) = ccc(i,j) = cc
OpenACCはこの時点で実行可能！
end doend doend doend do!$acc end kernels
end subroutine matmul
end subroutine matmul

57

最低限動くプログラムを作るには
3. GPU用のメモリを確保し、明示的にデータ転送する
OpenACC
subroutine matmul(a, b, c, n)real(8), dimension(n, n) :: a, b, c
integer :: ninteger :: i, j, kreal(8) :: cc!$acc data copyin(a, b) copyout(c)!$acc kernels
do j = 1, n
do i = 1, n
cc = 0do k = 1, n
cc = cc + a(i,k) * b(k,j)
end doc(i,j) = cc
end doend do!$acc end kernels!$acc end data
end subroutine matmul

CUDA
attribute(global) subroutine mm_cuda(a, b, c, n)
integer, value :: nreal(8), dimension(n, n) :: a, b, c
integer :: i, j, kreal(8) :: cci = (blockidx%x-1) * blockdim%x + threadidx%x
j = (blockidx%y-1) * blockdim%y + threadidx%y
cc = 0.0
do k = 1, n
cc = cc + a(i, k) * b(k, j)
end doc(i, j) = cc
end subroutine mm_cuda
subroutine matmul(a, b, c, n)real(8), dimension(n, n) :: a, b, c
integer :: ninteger :: i, j, kreal(8) :: ccreal(8), device, allocatable, dimension(:, :) :: &a_dev, b_dev, c_dev
type(dim3) :: dimGrid, dimBlcok
allocate(a_dev(n, n), b_dev(n, n), c_dev(n, n))a_dev(:, :) = a(:, :)b_dev(:, :) = b(:, :)dimGrid = dim3( n/16, n/16, 1)
dimBlock = dim3( 16, 16, 1)
call mm_cuda<<<dimGrid, dimBlock>>>(a, b, c, n)c(:, :) = c_dev(:, :)end subroutine matmul

58最低限動くプログラムを作るには
4. スレッドを割り当てる
OpenACC
subroutine matmul(a, b, c, n)real(8), dimension(n, n) :: a, b, c
integer :: ninteger :: i, j, kreal(8) :: cc!$acc data copyin(a, b) copyout(c)!$acc kernels!$acc loop gang
do j = 1, n!$acc loop vector
do i = 1, n
cc = 0!$acc loop seq
do k = 1, n
cc = cc + a(i,k) * b(k,j)
end doc(i,j) = cc
end doend do!$acc end kernels!$acc end data
end subroutine matmul

CUDA
attribute(global) subroutine mm_cuda(a, b, c, n)
integer, value :: nreal(8), dimension(n, n) :: a, b, c
integer :: i, j, kreal(8) :: cci = (blockidx%x-1) * blockdim%x + threadidx%x
j = (blockidx%y-1) * blockdim%y + threadidx%y
cc = 0.0
do k = 1, n
cc = cc + a(i, k) * b(k, j)
end doc(i, j) = cc
end subroutine mm_cuda
subroutine matmul(a, b, c, n)real(8), dimension(n, n) :: a, b, c
integer :: ninteger :: i, j, kreal(8) :: ccreal(8), device, allocatable, dimension(:, :) :: &a_dev, b_dev, c_dev
type(dim3) :: dimGrid, dimBlcok
allocate(a_dev(n, n), b_dev(n, n), c_dev(n, n))a_dev(:, :) = a(:, :)b_dev(:, :) = b(:, :)dimGrid = dim3( n/16, n/16, 1)
dimBlock = dim3( 16, 16, 1)
call mm_cuda<<<dimGrid, dimBlock>>>(a, b, c, n)c(:, :) = c_dev(:, :)end subroutine matmul

59

OpenACC の指示文

GPUプログラミング入門

60

OpenACC の主要な指示文• 並列領域指定指示文– kernels, parallel• データ移動最適化指示文– data, enter data, exit data, update• ループ最適化指示文– loop• その他、比較的よく使う指示文– host_data, atomic, routine, declare

GPUプログラミング入門

61

並列領域指定指示文
１．必要なデータを送る

ホスト

デバイス

~32GB/s２．計算を行う

３．計算結果を返す~200GB/s~1,000GB/s
デバイス
メモリ
メインメモリ• 1, 2, 3 全てを行う指示文

GPUプログラミング入門

62並列領域指定指示文：parallel, kernels• アクセラレータ上で実行すべき部分を指定– OpenMPのparallel指示文に相当• 2種類の指定方法：parallel, kernels– parallel：(どちらかというと) マニュアル• OpenMP に近い• 「ここからここまでは並列実行領域です。並列形状などはユー
ザー側で指定します」– kernels：(どちらかというと) 自動的• 「ここからここまではデバイス側実行領域です。あとはお任せします」– 細かい指示子・節を加えていくと最終的に同じような挙動になるので、どちらを使うかは好み
63kernels/parallel 指示文
kernels

parallelprogram main

program main!$acc kernels
do i = 1, N
! loop bodyend do!$acc end kernels!$acc parallel num_gangs(N)!$acc loop gang
do i = 1, N
! loop bodyend do!$acc end parallelend program

end program

GPUプログラミング入門

64kernels/parallel 指示文• ホスト-デバイスを意識するのがkernels• 並列実行領域であることを意識するのがparallel

kernels

parallel

ホスト側
デバイス側
program main

program main!$acc kernels
do i = 1, N
! loop bodyend do!$acc end kernels!$acc parallel num_gangs(N)!$acc loop gang
do i = 1, N
! loop bodyend do!$acc end parallelend program

end program「並列数はデバイスに合わせてください」「並列数Nでやってください」
GPUプログラミング入門

65kernels/parallel 指示文：指示節
kernels
•
•
•
•
•
•

async
wait
device_typeifdefault(none)
copy…

parallel
•
•
•
•
•
•
•
•
•
•
•
•

async
wait
device_typeifdefault(none)
copy…
num_gangs
num_workers
vector_length
reduction
privatefirstprivate

GPUプログラミング入門

66kernels/parallel 指示文：指示節
kernels

parallel
非同期実行に用いる。
今回は扱わない。

実行デバイス毎にパラメータを調整
if(0)/if(.false.)などとするとホスト側で実行される
データの自動転送を行わないようにする
データ指示文の機能を使える

parallelでは並列実行領域であることを意識するため、並列数や変数の扱いを決める指示節がある。

•
•
•
•
•
•
•
•
•
•
•
•

async
wait
device_typeifdefault(none)
copy…
num_gangs
num_workers
vector_length
reduction
privatefirstprivate

GPUプログラミング入門

67kernels/parallel 実行イメージ
FortranC言語

subroutine copy(dis, src)
real(4), dimension(:) :: dis, src

void copy(float *dis, float *src) {int i;!$acc kernels copy(src,dis)do i = 1, Ndis(i) = src(i)end do!$acc end kernels#pragam acc kernels copy(src[0:N] ¥dis[0:N])for(i = 0;i < N;i++){dis[i] = src[i];}end subroutine copy}

GPUプログラミング入門

68kernels/parallel 実行イメージ
Fortran(ホスト)(デバイス)

subroutine copy(dis, src)
real(4), dimension(:) :: dis, src

dis, src⓪dis, src の領域が確保されるend subroutine copy

①dis, src の値がコピーされる③dis’_dev,
src’_dev の値

dis_dev, src_dev
②デバイス
上の計算

!$acc kernels copy(src,dis)do i = 1, Ndis(i) = src(i)end do!$acc end kernels

がコピーされるdis’, src’

dis’_dev, src’_dev④dis, src の領域が解放される

GPUプログラミング入門

69

デバイス上で扱われるべきデータについて• プログラム上のparallel/kernels構文に差し掛かった時、OpenACCコンパイラは実行に必要なデータを自動で転送する– 正しく転送されないこともある。自分で書くべき– 構文に差し掛かるたびに転送が行われる(非効率)。後述のdata指示文を用いて自分で書くべき– 自動転送はdefault(none)で抑制できる• スカラ変数は firstprivate として扱われる– 指示節により変更可能• 配列はデバイスに確保される (shared的振る舞い)– 配列変数をスレッドローカルに扱うためには private を指定する
GPUプログラミング入門

70

データ移動最適化指示文
１．必要なデータを送る

ホスト~32GB/s

デバイス
２．計算を行う

３．計算結果を返す~200GB/s~1,000GB/s
デバイス
メモリ
メインメモリ• 1, 3 を最適化するための指示文• ただしデータの一貫性を維持するのはユーザの責任GPUプログラミング入門

71

データ移動最適化指示文：data, enter/exit data
• デバイス側で必要なデータと範囲を指定– Allocate, Memcpy, Deallocate を行う• data 指示文 (推奨)– Allocate + Memcpy (H ⇄ D) + Deallocate– 構造ブロックに対してのみ適用可• コードの見通しが良い• enter data 指示文– Allocate + Memcpy (H → D)
– exit data とセット。構造ブロック以外にも使える• exit data 指示文– Memcpy (H ← D) + Deallocate– enter data とセット。構造ブロック以外にも使える
GPUプログラミング入門

72

データ移動最適化指示文が必要なとき
FortranC言語

subroutine copy(dis, src)
real(4), dimension(:) :: dis, src
do j = 1, M!$acc kernels copy(src,dis)do i = 1, Ndis(i) = dis(i) + src(i)end do!$acc end kernels
end doend subroutine copy

void copy(float *dis, float *src) {int i, j;for(j = 0;j < M;j++){#pragam acc kernels copy(src[0:N] ¥dis[0:N])for(i = 0;i < N;i++){dis[i] = dis[i] + src[i];}}}Kernels をループで囲
むとどうなるか…

GPUプログラミング入門

73

データ移動最適化指示文が必要なとき
C言語
for(j = 0;j < M;j++){１．必要なデータを送る

ホスト~32GB/s

デバイス
２．計算を行う

３．計算結果を返す~200GB/s

メインメモリ}~1,000GB/s
デバイス
メモリ

1, 2, 3 全てが繰り返される！

void copy(float *dis, float *src) {int i, j;for(j = 0;j < M;j++){#pragam acc kernels copy(src[0:N] ¥dis[0:N])for(i = 0;i < N;i++){dis[i] = dis[i] + src[i];}}}Kernels をループで囲
むとどうなるか…

GPUプログラミング入門

74data指示文
FortranC言語

void copy(float *dis, float *src) {subroutine copy(dis, src)
int i, j;real(4), dimension(:) :: dis, src#pragma acc data copy(src[0:N] ¥!$acc data copy(src,dis)
dis[0:N])do j = 1, M{!$acc kernels present(src,dis)for(j = 0;j < M;j++){do i = 1, N
#pragam acc kernels present(src,dis)dis(i) = dis(i) + src(i)for(i = 0;i < N;i++){present: 既に転送済
end dodis[i] = dis[i] + src[i];であることを示す!$acc end kernels}(OpenACC2.5の仕様以降、copyはpresent_or_copyとしend doて扱われることになったので、 }実は書き換えなくても大丈夫
Cの場合、data指示文の範囲は!$acc end data}になった。){}で指定}
end subroutine copy
(この場合はforが構造ブロックになってるのでなくても大丈夫だが)GPUプログラミング入門

75data指示文
C言語
void copy(float *dis, float *src) {int i, j;#pragma acc data copy(src[0:N] ¥１．必要なデータを送る
dis[0:N])for(j = 0;j < M;j++){ {ホスト
デバイス~32GB/sfor(j = 0;j < M;j++){２．計算を行う
３．計算結果を返す#pragam acc kernels present(src,dis)for(i = 0;i < N;i++){~200GB/s~1,000GB/s
デバイス
dis[i] = dis[i] + src[i];メモリ
メインメモリ}}}
Cの場合、data指示文の範囲は
}
2のみが
{}で指定
繰り返される！
}(この場合はforが構造ブロックになってるのでなくても大丈夫だが)GPUプログラミング入門

76

enter data, exit data指示文
void main() {double *q;
int step;for(step = 0;step < N;step++){if(step == 0) init(q);solverA(q);solverB(q);….
if(step == N) fin(q);}}void init(double *q) {q = (double *)malloc(sizeof(double)*M);
q = … ; // 初期化#pragma acc enter data copyin(q[0:M])}
void fin(double *q) {#pragma acc exit data copyout(q[0:M])
print(q); //結果出力
free(q);}GPUプログラミング入門

77data, enter/exit data 指示文の指示節
data
•
•
•
•
•
•
•
•

enter data
•
•
•
•
•
•

if
copy
copyincopyout
createpresent
present_or_...
deviceptr

if
async
非同期転送用
wait
copyincreatepresent_or_...

exit data– CUDAなどと組み合わせる時に利用。
cudaMallocなどで確保済みのデータを指定し、OpenACCで扱い可とする

•
•
•
•
•

if
async
wait
copyout
deleteGPUプログラミング入門

78

data, enter/exit data 指示文の指示節•

copy– allocate, memcpy (H→D), memcpy (D→H), deallocate•

copyin– allocate, memcpy (H→D), deallocate 結果の出力を行わない•

copyout– allocate, memcpy (D→H), deallocate データの入力を行わない•

create– allocate, deallocate コピーは行わない•

present– 何もしない。既にデバイス上にあることを教える•

present_or_copy/copyin/copyout/create (省略形：pcopy)– デバイス上になければ copy/copyin/copyout/create する。なければ何もしない
ただしOpenACC2.5以降では、
copy, copyin, copyout の挙動はpcopy, pcopyin, pcopyout と同一

79

データ移動指示文：データ転送範囲指定• 送受信するデータの範囲の指定– 部分配列の送受信が可能– 注意：FortranとCで指定方法が異なる• 二次元配列Aを転送する例
Fortran版!$acc data copy(A(lower1:upper1, lower2:upper2) )…
fortranでは下限と上限を指定!$acc end data

C版

#pragma acc data copy(A[start1:length1][start2:length2])…
Cでは先頭と長さを指定#pragma acc end data80

データ移動指示文：update 指示文• データ指示文などで既にデバイス上に確保済みのデータを対象とする– Memcpy (H ⇄ D) の機能を持っていると思えば良い!$acc data copy( A(:,:) )do step = 1, N…
!$acc update host( A(1:2,:) )call comm_boundary( A )!$acc update device( A(1:2,:) )…
end do!$acc end data
GPUプログラミング入門

update
•
•
•
•
•
•
•

if
async
wait
device_type
self #host と同義
host # H ← D
device # H → D
81

ループ最適化指示文
１．必要なデータを送る

ホスト

デバイス

~32GB/s２．計算を行う

３．計算結果を返す~200GB/s~1,000GB/s
デバイス
メモリ
メインメモリ• 2の最適化を行う指示文

GPUプログラミング入門

82

階層的並列モデルとループ指示文• OpenACC ではスレッドを階層的に管理– gang, worker, vector の3階層
– gang：workerの塊 一番大きな単位– worker：vectorの塊– vector：スレッドに相当する一番小さい
処理単位• loop 指示文– parallel/kernels中のループの扱いについて指示• パラメータの設定はある程度勝手にやってくれる– 粒度(gang, worker, vector)の指定– ループ伝搬依存の有無の指定

GPUでの行列積の例!$acc kernels!$acc loop gang
do j = 1, n!$acc loop vector
do i = 1, n
cc = 0!$acc loop seq
do k = 1, n
cc = cc + a(i,k) * b(k,j)
end doc(i,j) = cc
end doend do!$acc end kernels

83階層的並列モデルとアーキテクチャ• OpenMPは1階層– マルチコアCPUも1階層•

NVIDIA GPUの構成
GPU– 最近は2階層目(SIMD)がある• CUDAは block と thread の2階層
– NVIDA GPUも2階層

デバイスメモリ
SMX• 1 SMX に複数CUDA coreを搭載• 各コアはSMXのリソースを共有• OpenACCは3階層– 様々なアクセラレータに対応す
るため

CUDA
コア

84

84

ループ指示文：指示節
loop
•
•
•
•
•
•
•
•
•
•
•

collapsegang
workervector
seq
auto
tile
device_type
independent
privatereduction

GPUプログラミング入門

85

ループ指示文：指示節
loop
•
•
•
•
•
•
•
•
•
•
•

collapsegang
workervector
seq
auto
tile
device_type
independent
privatereduction

3つのループが一重化される!$acc kernels!$acc loop collapse(3) gang vector
do k = 1, 10
do j = 1, 10
do i = 1, 10….
end doend doend do!$acc end kernels

並列化するにはループ長の短すぎるループに使う

GPUプログラミング入門

86

ループ指示文：指示節
loop
•
•
•
•
•
•
•
•
•
•
•

collapsegang
workervector
seq
auto
tile
device_type
independent
privatereduction!$acc kernels!$acc loop gang(N)do k = 1, N!$acc loop worker(1)do j = 1, N!$acc loop vector(128)do i = 1, N

vectorはworkerより内側
workerはgangより内側….
!$acc kernels!$acc loop gang vector(128)do i = 1, N

ただし1つのループに複数つけるのはOK….
数値の指定は難しいので、最初はコンパイラ任せでいい
GPUプログラミング入門

87

do j = 1, N

ループ指示文：指示節
collapsegang
workervector
seq
auto
tile
device_type
independent
privatereductionidxI(i) = i; idxJ(j) = j
end do

loop
•
•
•
•
•
•
•
•
•
•
•

do i = 1, N

end do!$acc kernels &!$acc& copyin(A, idxI, idxJ) copyout(B)!$acc loop independent gang
do j = 1, N!$acc loop independent vector(128)do i = 1, N
B(idxI(i),idxJ(j)) = alpha * A(i,j)

Bに間接参照→

end doend do!$acc end kernels

OpenACCコンパイラは保守的。依存関係が生じそうなら並列化しない。きちんと並列化されているかどうか、必ずコンパイラのメッセージを確認する。(やり方は後述)
GPUプログラミング入門

88

ループ指示文：指示節
loop
•
•
•
•
•
•
•
•
•
•
•

collapsegang
workervector
seq
auto
tile
device_type
independent
privatereduction!$acc kernels &!$acc loop reduction(+:val)
do i = 1, N

簡単なものであれば、PGIコ
ンパイラは自動でreduction
を入れてくれる

val = val + 1

利用できる演算子

end do!$acc end kernels(OpenACC2.0仕様書より)

acc reduction (+:val)
演算子
対象とする変数
制限：スカラー変数のみ

GPUプログラミング入門

89

参考：リダクションでOpenACCが好きになる
リダクションが必要な例
sum = 0.0
for(i = 0;i < N;i++)sum += array[i]• そもそもリダクションって？
スレッド1

array

1

2

スレッド2

3

4

5

スレッド3

6

7

8

9
1. 各スレッドが担当領域で縮約6

15

24
2. 一時配列に書き
込む

6

15

24
3. 一時配列を縮約

sum =

45
GPUプログラミング入門

90

参考：リダクションでOpenACCが好きになる•

CUDAで実装しようと思うと…
– 全部自分でやる– これはshuffle機能を使ったリダクション

int main( int argc, char* argv[] ){ホスト側プログラム….
double *tmp,*ans_d;
cudaMalloc((void**)&tmp, sizeof(double) * 896);
cudaMalloc((void**)&ans_d, sizeof(double) * 1); 一時配列の確保
int chunk = (N+895)/896;
dim3 dimGrid_L(896, 1, 1);
dim3 dimBlock_L(128, 1, 1);
dim3 dimGrid_G(1, 1, 1);
dim3 dimBlock_G(1024, 1, 1);

使用するスレッド数の宣言

GPUカーネル呼び出し

reduction_L <<<dimGrid_L, dimBlock_L>>> (N,A_d,tmp,chunk);
reduction_G <<<dimGrid_G, dimBlock_G>>> (tmp,ans_d);}

cudaMemcpy(&sum,ans_d,sizeof(double),cudaMemcpyDeviceToHost);….

結果の書き戻し__global__ void reduction_L(int N, double *A, double *tmp, int chunk){__shared__ double sum_tmp[4];int tx = threadIdx.x; int bx = blockIdx.x;
GPUカーネル1
int st = bx*chunk + tx; int en = min(N,(bx+1)*chunk);
double sum = 0.0;
for(i = st;i < en;i+=128)
sum += A[i];Warp shuffle機能を
sum += __shfl_xor(sum,16); sum += __shfl_xor(sum,8); 使ったWarp内の縮約
sum += __shfl_xor(sum,4); sum += __shfl_xor(sum,2);
sum += __shfl_xor(sum,1);
if(tx % 32 == 0) sum_tmp[tx/32] = sum;__syncthreads();if(tx == 0) tmp[bx] = sum + sum_tmp[1] + sum_tmp[2] + sum_tmp[3];}shared memoryを使ったWarp間の縮約

syncthreads()による同期必須__global__ void reduction_G(double *tmp, double *ans){
GPUカーネル2__shared__ double sum_tmp[32];double sum;int tx = threadIdx.x;
if(tx >= 896)
一時配列も同様に縮約
sum = 0.0;
else
sum = tmp[tx];sum += __shfl_xor(sum,16); sum += __shfl_xor(sum,8);
sum += __shfl_xor(sum,4); sum += __shfl_xor(sum,2);
sum += __shfl_xor(sum,1);
if(tx % 32 == 0) sum_tmp[tx/32] = sum;__syncthreads();if(tx < 32){
sum = sum_tmp[tx];sum += __shfl_xor(sum,16); sum += __shfl_xor(sum,8);
sum += __shfl_xor(sum,4); sum += __shfl_xor(sum,2);
sum += __shfl_xor(sum,1);}__syncthreads();if(tx == 0) ans[0] = sum;}

GPUプログラミング入門

91

参考：リダクションでOpenACCが好きになる• OpenACC なら…

• 性能も…

sum = 0.0
#pragma acc kernels copyin(A[0:N])#pragma acc loop reduction(+:sum)for(i = 0;i < N;i++){sum += array[i]
}
これだけ！

1.

OpenACC のリダクション

2.

一旦CPUに書き戻して1スレッドで計算

3.

CPUに書きもどさず、GPUの1スレッドで計算

4.

CUDA の shuffleを使ったリダクション

array(倍精度)のサイズ100000000の時 (Fortran)

1.
2.
3.
4.

reduction acc time:
copyback and CPU time:
cuda 1 thread time:
reduction cuda time :0.00148 (s)0.63831 (s)9.63381 (s)0.00140 (s)

下手なプログラムを書くくらいなら
OpenACCに任せた方が速い！
GPUプログラミング入門

92

関数呼び出し指示文：routine• parallel/kernels領域内から関数を呼び出す場合、routine指示文を使う#pragma acc routine vector

プロトタイプ宣言にもつける

extern double vecsum(double *A);…
#pragma acc parallel num_gangs(N) vector_length(128)for (int i = 0;i < N; i++){#pragma acc routine vector

max = vecsum(A[i*N]);double vecsum(double *A){}

double x = 0;
#pragma acc loop reduction(+:x)
for(int j = 0;j < N;j++){x += A[j];}return x;
GPUプログラミング入門}

93host_data 指示文• OpenACC のデータ指示文
を使うと…
subroutine copy(dis, src)
real(4), dimension(:) :: dis, src!$acc data copy(src,dis)!-------------------------! CPU上の処理!-------------------------call hoge(dis,src)!$acc kernels present(src,dis)do i = 1, Ndis(i) = dis(i) + src(i)end do!$acc end kernels(ホスト)(デバイス)①両方に作られるdis, src

dis_dev, src_dev②並列領域以外ではホスト側が使われる関数にもホスト側のアドレスが渡る③並列領域ではデバイス側が使われる!$acc end data
end subroutine copy
GPUプログラミング入門

デバイス側のアドレスを渡したい時は？？

94host_data 指示文• 並列領域の外でデバイス側のアドレスを使うための指示文– データ指示文で確保済の配列が対象• デバイス側のアドレスを使いたいケースって？
– CUDAで書いた関数の呼び出し– GPU用のライブラリの呼び出し– GPU Direct によるMPI通信• GPU Direct： ホスト側のメモリを介さず、GPU間で直接MPIによるデー
タ通信をするもの
FortranC!$acc data create(tmp) copy(val) copyin(A)...
!$acc host_data use_device(A,tmp,val)call reduction_cuda_shuffle(val,N,A,tmp)!$acc end host_data...
!$acc end data#pragma acc data create(tmp[0:896]) copy(val){...
#pragma acc host_data use_device(A,tmp,val){
reduction_cuda_shuffle(&val,N,A,tmp);}...
}

GPUプログラミング入門
領域内ではデバイス側のアドレスが使われる

95

atomic 指示文• 並列化領域内で、どうしても並列化できない部分が存在する場合に使う– 全体的には並列に計算できるが、書き込み先が衝突するようなケース• Pascal GPU は倍精度のatomicAddをハードウェアサ
ポートしてるので、それなりに速い
一文ずつ囲む!$acc kernels async(0)!$acc loop independent gang vector
do k= inls,inle
i= IAL(k)X1= Xin(3*i-2)
X2= Xin(3*i-1)
X3= Xin(3*i )WVAL1= AL(k,1)*X1 + AL(k,2)*X2 + AL(k,3)*X3
WVAL2= AL(k,4)*X1 + AL(k,5)*X2 + AL(k,6)*X3
WVAL3= AL(k,7)*X1 + AL(k,8)*X2 + AL(k,9)*X3
i = INL_G(k)i が関節参照!$acc atomic
tmpL(i,1) = tmpL(i,1) + WVAL1 なので、書き!$acc end atomic
込み先が衝!$acc atomic
突する可能性tmpL(i,2) = tmpL(i,2) + WVAL2 がある!$acc end atomic!$acc atomic
tmpL(i,3) = tmpL(i,3) + WVAL3!$acc end atomic
enddo!$acc end kernels

GPUプログラミング入門

96

アプリケーションの移植方法

GPUプログラミング入門

97

アプリケーションのOpenACC化手順
1. プロファイリングによるボトルネック部位の導出
2. ボトルネック部位のOpenACC化
1.

並列化可能かどうかの検討

2.

(OpenACCの仕様に合わせたプログラムの書き換え)3.

parallel/kernels指示文適用

3. data指示文によるデータ転送の最適化
4. OpenACCカーネルの最適化
1 ~ 4 を繰り返し適用。それでも遅ければ、
5. カーネルのCUDA化–

スレッド間の相互作用が多いアプリケーションでは、shared memory
や shuffle 命令を自由に使えるCUDAの方が圧倒的に有利
GPUプログラミング入門

98アプリケーションのOpenACC化手順
int main(){double A[N];sub1(A);sub2(A);
sub3(A);}sub2(double A){
subA(A);
subB(A);
}

ホスト

デバイス

main

sub1

sub2

subA

sub3

subB

subA(double A){}

for( i = 0 ~ N ) {…
}葉っぱの部分から
OpenACC化を始める

GPUプログラミング入門

99

アプリケーションのOpenACC化手順
int main(){double A[N];sub1(A);sub2(A);
sub3(A);}sub2(double A){
subA(A);
subB(A);}
subA(double A){#pragma acc …
for( i = 0 ~ N ) {
…
}
}

ホスト

デバイス

main

sub1

subA

sub2

sub3

data指示文で配列Aをコピー
subB

subA

この状態でも必ず正しい結果を得られるように作る！
この時、速度は気にしない！

GPUプログラミング入門

100

アプリケーションのOpenACC化手順
int main(){double A[N];sub1(A);#pragma acc data{sub2(A);}
sub3(A);}sub2(double A){
subA(A);
subB(A);}
subA(double A){#pragma acc …
for( i = 0 ~ N ) {
…
}
}

ホスト

デバイス

main

sub1

subA

sub2

data指示文で配列Aをコピー
sub3

subB

sub2

subA

subB

徐々にデータ移動を上流に移動する

GPUプログラミング入門

101アプリケーションのOpenACC化手順
int main(){double A[N];#pragma acc data{sub1(A);sub2(A);
sub3(A);}}sub2(double A){
subA(A);
subB(A);}
subA(double A){#pragma acc …
for( i = 0 ~ N ) {
…
}
}

ホスト

main

sub1

sub2

subA

デバイス

data指示文で配列Aをコピー

sub3

subB

sub1

subA

main

sub3

sub2

subB

ここまで来たら、ようやく個別のカーネルの最適化を始める。
※データの転送時間が相対的に十分小さくなればいいので、かならずしも最上流までやる必要はない
GPUプログラミング入門

102

上で説明した指示文を用いてOpenACC化したもの

OpenACC の適用例

GPUプログラミング入門

103

OpenACCを用いたICCG法ソルバーの
Pascal GPUにおける性能評価とKNL
星野哲也・大島聡史 ・塙 敏博 ・中島研吾 ・伊田明弘
1)

東京大学情報基盤センター・最先端共同HPC基盤施設（JCAHPC）
2)

第158回 HPC研究会＠熱海メニーコアプロセッサと消費電力• Green500 Top10(2016/11)は全部メニーコア

P100
P100
PEZY
Sunway• メニーコア向けのアルゴリズム
や最適化手法の開発が重要

KNL
KNL
KNL
KNL
KNL
KNL

105

本研究の概要• 研究目的– 最新世代のGPU，Xeon Phiである、P100とKNLの性能評価• 評価手法
– OpenACC, OpenMPにより並列化したICCG法ソルバーを使用– P100とKNLの性能を比較– 前世代 (Kepler, KNC) と比較– それぞれへの最適化の適用• 成果– P100において、指示文の最適化などにより1.21倍の性能向上を達成– KNLにおいて、同期の削減などにより1.20倍の性能向上を達成
106

発表概要
• 研究背景• 対象アプリケーション• ICCG法・行列格納手法• OpenACCによる並列化• P100 vs KNL ( vs BDW, K20, KNC)
• P100・KNL向け最適化• まとめ・今後の課題

107

対象アプリケーション• 有限体積法，一様場ポアソン方程式ソルバー– [大島他 SWoPP 2014], [中島 HPC-139], [中島 HPC-147], [中島 HPC-157]– 差分格子：データ構造は非構造，7点ステンシル, 問題サイズ1283– 対称正定な疎行列を係数とする連立一次方程式– ICCG法，上下三角成分を別々に格納• 色付け・リオ－ダリング– CM-RCM + Coalesced/Sequential• 行列格納手法– CRS, Sliced-ELL, Sell-C-σ

NZ
z

DZ

y

NY
x

DY
DX

NX

108

ICCG法• 共役勾配法(Conjugate Gradient,
CG)
• 前処理
– 不完全コレスキー分解(Incomplete Cholesky
Factorization, IC)– 並列化にはカラーリングが必要

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|
end

109

ICCG法並列化：データ依存性の解決• 「色づけ＋リオ－ダリング」による並列性抽出• 同色内は並列に実行可能• 色数 → 収束性 並列性 同期コスト
29 61
45
22 46
16 62
11 47
7 63
4 48
2 64
1

29 22 16 11

7

4

2

1

53 36 19

13
37 29
49
30 14
23 30
51
17 15
12 31
53
8 16
55
5 32
3

37 30 23 17 12

8

5

3

7

54 37 20

44 57
41
38 42
31 58
24 43
18 59
13 44
9 60
6

44 38 31 24 18 13

9

6

25

8

9 25
33
50
45 10
39 26
35
32 11
25 27
37
19 12
14 28
39
10

50 45 39 32 25 19 14 10

43 26

5

52

55 53
37
51 38
46 54
40 39
33 55
26 40
20 56
15

55 51 46 40 33 26 20 1561 44 27 10 57 40 23

6

5 21
6 22
7 23
8 24
17
59
56 19
52
47 21
41
34 23
27
21

59 56 52 47 41 34 27 21

14 62 45 28 11 58 41 24

60 34
57 50
53 35
48 51
42 36
35 52
28
62 49
33

62 60 57 53 48 42 35 28

31 15 63 46 29 12 59 42

1 17
2
3
4
64
63 61
3 18
58 54
5 19
49 43
7 20
36

64 63 61 58 54 49 43 36

48 32 16 64 47 30 13 60

MC (Color#=4)
Multicoloring

RCM
Reverse Cuthill-Mckee

2

49 33 173

55 38 219

1

50 34 18
4

56 39 22

51 35

CM-RCM (Color#=4)
Cyclic MC + RCM

110

ICCG法並列化：データ依存性の解決
MC● 並列性良❌悪条件問題に弱いRCM● 収束性良❌ 並列性悪、同期コ
スト高

29 61
45
22 46
16 62
11 47
7 63
4 48
2 64
1

29 22 16 11

7

4

2

1

53 36 19

13
37 29
49
30 14
23 30
51
17 15
12 31
53
8 16
55
5 32
3

37 30 23 17 12

8

5

3

7

54 37 20

44 57
41
38 42
31 58
24 43
18 59
13 44
9 60
6

44 38 31 24 18 13

9

6

25

8

9 25
33
50
45 10
39 26
35
32 11
25 27
37
19 12
14 28
39
10

50 45 39 32 25 19 14 10

43 26

5

52

55 53
37
51 38
46 54
40 39
33 55
26 40
20 56
15

55 51 46 40 33 26 20 1561 44 27 10 57 40 23

6

5 21
6 22
7 23
8 24
17
59
56 19
52
47 21
41
34 23
27
21

59 56 52 47 41 34 27 21

14 62 45 28 11 58 41 24

60 34
57 50
53 35
48 51
42 36
35 52
28
62 49
33

62 60 57 53 48 42 35 28

31 15 63 46 29 12 59 42

1 17
2
3
4
64
63 61
3 18
58 54
5 19
49 43
7 20
36

64 63 61 58 54 49 43 36

48 32 16 64 47 30 13 60

MC (Color#=4)
Multicoloring

RCM
Reverse Cuthill-Mckee

CM-RCM● 並列性・収束性
本発表で使用249 33 173

55 38 219

1

50 34 18
4

56 39 22

51 35

CM-RCM (Color#=4)
Cyclic MC + RCM

111

Coalesced,
Sequential オー
ダリング

Coalesced

Initial Vector

Coloring(5 colors)+Orderingcolor=1

color=1

color=2

color=2

color=3

color=4

color=3

color=4

color=5

color=512345678 12345678 12345678 12345678 12345678

色順に並べ替えた後、• Coalesced– データの順はそのまま– 各スレッドは不連続なメ
モリアクセス•Sequential– データ順を再オーダリ
ング– 計算単位(スレッドやコア)内で連続アクセス

Sequential

Initial Vector

Coloring(5 colors)+Ordering
各スレッド上で不連続なメモリ
アクセス（色の順に番号付け）

color=1

color=1

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

11111 22222 33333 44444 55555 66666 77777 88888
スレッド内で連続に番号付け

112

行列格納方法
一番一般的な格
納方法

0でpadding

非零要素数でソート。複数配列で表現することで無駄な計算を削減CをSIMD長等に合わせ、適宜0padding で整形

C

s

CRS

ELL

Sliced ELL

SELL-C-s(SELL-2-8)
113

実施ケース

略称

Numbering

格納形式

c-CRS

Coalesced

CRS

c-Sliced-ELL

Sliced-ELL (ブロッキングあり)
SELL-C-σc-SELL-C-σs-CRS

Sequentials-Sliced-ELL
s-SELL-C-σ

CRS
Sliced-ELL (ブロッキングあり)
SELL-C-σ

カラーリングはいずれもCM-RCM
問題サイズは1283

114

OpenACC/OpenMP実装• 同一ソースコードにそれぞれの指示文を挿入– HPC157から、不要な一時配列の除去などの最適化– HPC157から、一部最適化を評価のため除外(OpenMP)• OpenACC– 全並列化ループに !$acc kernels を適用– データ転送は計測外で行う– 予稿から最適化レベルを1段階ダウン• async 指示節を除外• OpenMP– 全並列化ループに !$omp parallel do を適用!$omp parallel do private(... )!$acc kernels!$acc loop independent gang
do ip= 1, PEsmpTOT
ip1= (ip-1)*NCOLORtot + ic!$omp simd!$acc loop independent vector
do i= index(ip1-1)+1, index(ip1)
…
enddo
enddo
!$omp end parallel do!$acc end kernels

115

Hardware• P100: NVIDIA Tesla P100 (予稿はPCI-E版)– Reedbush-H• K20: NVIDIA Tesla K20Xm– TSUBAME2.5• KNC: Intel Xeon Phi 5110P (Knights Corner)– KNSC• KNL: Intel Xeon Phi 7250 (Knights Landing)– Oakforest-PACS• BDW: Intel Xeon E5-2695 v4 (Intel Broadwell-EP)– Reedbush-U
116

Hardware spec
略 称

P100

K20

KNL

KNC

BDW

動作周波数（GHz）

1.480

0.732

1.40

1.053

2.10

コア数（最大有効ス
レッド数）

1,792

896

68 (272)

60 (240)

18 (18)

理論演算性能（GFLOPS）

5,304

1,311.7

3,046.4

1,010.9604.8

主記憶容量（GB）

16

6

16

8

128

メモリバンド幅（GB/sec., Stream
Triad）

534

179

490

159

65.5

Reedbush-H

TSUBAME2.5

OakforestPACS

KNSC

Reedbush-U

System117

コンパイラ・環境変数
P100• コンパイラ&オプション– pgfortran 17.1 (予稿は16.10)– -O3 –ta=tesla:cc60• 環境変数– 特になし118

コンパイラ・環境変数
KNL• コンパイラ&オプション

その他
・numactl –membind=1 で MCDRAMを使用
・メモリモデル：Flat
・サブNUMAモード：Quadrant– ifort (IFORT) 17.0.0 20160721– -align array64byte -O3 -xMIC-AVX512 -qopenmp -qopt-streamingstores=always -qopt-streaming-cache-evict=0• ただしCRSを使うときは-qopt-streaming-stores=never (予稿ではalways)• 環境変数– export OMP_STACKSIZE=1G– ulimit -s 1000000– 66 スレッドの時• export OMP_NUM_THREADS=66• export KMP_AFFINITY=granularity=fine,proclist=[2-67],explicit– 132/198/264 スレッドの時• export OMP_NUM_THREADS=132/198/264• export KMP_AFFINITY=compact• export KMP_HW_SUBSET=66c@2,2t/3t/4t
119

P100
1.1

•

Coalesced版が速い•

Sell-C-σ(C=64, σ=1)とSliced-ELLが同程度で最速•

CRSが3割程度遅い•

色数6程度が最速– 原因：同期コスト(カーネル呼び出し)– 予稿と違うのはasync
外したため

Elapsed time [sec]

1.0
0.9
0.8
0.7
0.6
0.5
2

4

6
c-CRS
s-CRS

8

10
12
Color#
c-Sliced-ELL
s-Sliced-ELL

14

16

18

20

c-SELL-64-1
s-SELL-64-1

120

KNL
2.5

Sequential版が速い

2.3•

Sell-C-σ(C=8, σ=1)が
最速

2.1

•

CRSが2.5倍程度遅い•

色数14程度が最速•

66/132/198/264 ス
レッド中 66 スレッドが最速

Elapsed time [sec]

•

1.9
1.7
1.5
1.3
1.1
0.9
0.7
0.5
2

4
c-CRS
s-CRS

6

8

10
12
Color#
c-Sliced-ELL
s-Sliced-ELL

14

16

18

20

c-SELL-8-1
s-SELL-8-1

121

各プロセッサの性能比較•

各プロセッサ最速の疎行列形状の結果
K20 vs KNC– KNCが1.57倍も遅い•P100 vs KNL– KNLが1.10倍しか遅くない•

P100 vs K20– P100はK20の2.71
倍程度

4.0
3.5
Elapsed time [sec]

•

4.5

BDW

3.0

KNC

2.5
2.0

K20

1.5
KNL

1.0
0.5

P100

0.0
2

4

6
P100

8
K20

10
12
14
Color#
KNL
KNC

16

18

20

BDW

122

実行効率 (vs. BDW)
• 色数10の時の比較• メモリバンド幅X倍で性能もX倍になると仮定すると、P100・KNLの実行効率はBDWの61%程度略 称
メモリバンド幅（GB/sec., Stream
Triad）

P100

K20

KNL

KNC

BDW

534

179

490

159

65.5

相対メモリバンド
幅

8.152

2.733

7.481

2.427

1

実行時間

0.64481.714

0.7061

2.708

3.224

相対性能

5.00

1.88

4.57

1.19

1

(相対性能)/(相対メモリバンド幅)

0.613

0.688

0.610

0.490

1

123

最適化(OpenACC)
1. Baseline– 全並列ループに!$acc kernels

2. Async (予稿のBaseline)– !$acc kernels を !$acc kernels async(0) と置き換え

3. Thread– gang, vector などのパラメータを調整

4. Fusion↑指示文のみの変更↓コードの変更– SPMV部分のカーネルを融合

124

最適化 3 Thread• リダクション部分のスレッド数を調整– 目的：ローカルリダクション
で生成される一時配列のサ
イズを小さくする!$omp parallel do private(i) reduction(+:RHO)!$acc kernels async(0)!$acc loop independent gang(448) vector reduction(+:RHO)
do i= 1, N
RHO= RHO + W(i,R)*W(i,Z)
enddo!$acc end kernels• デフォルト設定– ローカルリダクション• gang = (N-1)/128+1• vector = 128– グローバルリダクション• gang = 1• vector = 256↑ 448要素の一時配列ができ、その配列を1 gangでリダクション← 16384要素の一時配列ができ、その配列を1 gangでリダクション

125

最適化 4 Fusion• カーネルを融合し、呼び出しコ
ストを削減• SpMV部分では、色1, 色2~N-1,
色N の時で処理内容が異なる– 色2~N-1のカーネルをN-2回呼び出す!C!C SpMV!C
do ic = 1, NCOLORtot
if( ic == 1 ) then

! 色ループ! 色1の時の処理
else if ( 2 <= ic, ic < NCOLORtot) then
! 色2 ~ 色NCOLORtot-1の時の処理
else• 色2~N-1の処理を、複数のカーネルから1カーネルに融合! 色NCOLORtot の時の処理
endif
enddo

126

P100 最適化結果
0.75• asyncの効果大• Baselineから1.21倍の性能向上• BDW比実行効率
61.3%→74.5%

Elapsed time [sec]– カーネル起動
オーバーヘッド
を隠せる

0.70
0.65
0.60
0.55
0.50
2

4

6

1 Baseline

8
2 Async

10
12
Color#
3 Thread

14

16

18

20

4 Fusion127最適化(OpenMP)
1. Baseline– 全並列ループに !$omp parallel do

2. mvparallel1– !$omp parallel を色ループの外に

3. nowait– !$omp end do nowait を使う(SpMV部分)

4. mvparallel2– !$omp parallelを収束判定ループの外に

5. rmompdo↑指示文のみの変更↓コードの変更– !$omp doを使わず自力で(機械的に)ループ分割 (reduction以外)

6. rmreduction– reductionを自分で書く

7. loopschedule– ループスケジューリングの変更

128

Move !$omp paralell
Baseline
do L= 1, ITR ! 収束判定ループ

2 mvparallel
do L= 1, ITR ! 収束判定ループ

4 mvparallel!$omp parallel
do L= 1, ITR ! 収束判定ループ!$omp parallel
do ic = 1, NCOLORtot !色ループ!$omp parallel dodo ...

do ic = 1, NCOLORtot !色ループ!$omp doend dodo ...

!$omp end parallel doend doend do!$omp end do
end doenddo!$omp end parallelenddo

do ic = 1, NCOLORtot !色ループ!$omp dodo ...
end do!$omp end do
end doenddo!$omp end parallelループスケジューリ
ング

Sequential
Coloring(5 colors)+Ordering• 5 rmompdo• 自力で機械的にstatic分割• 7 loopscheduling• リオーダリングの知識を使い、前進後退代入などと同一領域を担当
Baseline!$omp dodo i= 1, N
W(i) = ....
enddo

各スレッド上で不連続なメモリ
アクセス（色の順に番号付け）

color=1

color=1

color=2

color=2

color=3

color=4

color=5

color=3

color=4

color=512345678 12345678 12345678 12345678 12345678

11111 22222 33333 44444 55555 66666 77777 88888
スレッド内で連続に番号付け

5 rmompdo
ip = omp_get_thread_num()+1

7 loopschedulingnth = omp_get_num_threads()
ls = (N+nth-1)/nth
do i= (ip-1)*ls+1, min(ip*ls,N)!$omp end do

W(i) = ....

!この後、W(i)を前進後退代入の入力

enddo!として使うInitial Vector!$omp barrier

do i = 自スレッドの担当部分
W(i) = ....
enddo
! バリアはいらない自力リダクション!$omp parallel private(i, ALPHA, C1, ....)
ip = omp_get_thread_num()+1

C1= 0.d0

nth = omp_get_num_threads()!$omp do reduction(+:C1)
do i= 1, N

ls = (N+nth-1)/nth

C1= C1 + W(i,P)*W(i,Q)

...

enddo

C1S(ip)= 0.0d0
do i= (ip-1)*ls+1, min(ip*ls,N)!$omp end do

C1S(ip)= C1S(ip) + W(i,P)*W(i,Q)

ALPHA= RHO / C1

enddo
C1= 0.d0• オリジナルはバリア2回• ローカルリダクションとグローバルリダクションの間• !$omp end do の時• nowaitは不可!$omp barrier
do i = 1, PEsmpTOT
C1= C1 + C1S(i)end do
ALPHA= RHO / C1

•
•

バリア1回
各自グローバルリダクション
131最適化(OpenMP)

1イテレーションあたりのOMP指示文数 (12色の時)

1

2

3

4

5

6

7

parallel do

40

5

5

0

0

0

0

parallel

0

3

3

0

0

0

0

4. mvparallel2

do

0

35

23

28

3

0

0

5. rmompdo

do(nowait)

0

0

12

12

0

0

0

reduction
clause

3

3

3

3

3

0

0

barrier(explicit)

0

0

0

1

26

29

26

barrier(implicit)

43

46

34

31

6

0

0

barrier
(合計)

43

46

34

32

32

29

26

1. Baseline
2. mvparallel1
3. nowait

6. rmreduction
7. loopschedule

最終的に、1回の!$omp parallel と!$omp
barrierだけのプログラムになる
132

KNL最適化結果
0.85– KNLのバリアは遅い？

• Baselineから1.20倍の性能向上• BDW比実行効率
61.0%→70.0%

0.80
Elapsed time [sec]• !$omp parallel の移動、バリアフリー
どちらも効果大

0.75
0.70
0.65
0.60
0.55
2

4

6

1 Baseline
4 mvparallel 2

8

10
12
Color#
2 mvparallel1
5 rmompdo

14

16

18

20

3 nowait
6 rmreduction133バリアフリーの効果(KNL, KNC)

0.7

Elapsed time

0.6
0.5
0.4
0.3
0.2
0.1
0.0

KNL
# of color = 15
# of thread = 66

3.0
2.5

18.1%

KNC
# of color = 13
# of thread = 120

27.2%

2.0
Elapsed time

0.8

1.5
1.0
0.5
0.0

134

バリアフリーの効果(BDW)
3.175
3.170
3.165

• BDWではほとんど効果が得られない– “7”では”4”から0.3%性
能低下

1.0%

3.155
Elapsed time

– 最外に!$omp parallelを追いやることで1.0%
性能向上

3.160

BDW
# of color = 20
# of thread = 18

3.150
3.145
3.140
3.135
3.130
3.125
3.120

0ではない

135

P100 vs KNL
0.9
0.8– P100: 534 GB/s
– KNL: 490 GB/s

0.7
Elapsed time [sec]• バンド幅比を考え
ると同程度

0.6
0.5
0.4
0.3
0.2
0.1
0.0
2

4

6

8

10
12
Color#
P100
KNL

14

16

18

20

136

まとめ
• ICCGソルバーにより、P100・KNLを評価• P100– OpenACCの async 節を適切につけることが必須• CUDAはdefaultで非同期– 指示文の最適化などにより1.21倍の高速化• KNL– バリアフリー化が効果大– 一連の最適化により1.20倍の高速化• 計算速度の向上、コアの増加による相対的な同期コスト増により、バリアフリー化が重要– ただし、既存のCPUで性能低下を引き起こす可能性あり• 今後の課題– MPIなどの通信を含むアプリケーションでの評価– 演算律速なアプリケーションでの評価– 同期コストについての詳細な評価
137

Q&A

GPUプログラミング入門

138

実習
※今回の実習の例は、全てPGIコンパイラ17.1を使った際の例ですGPUプログラミング入門

139

実習概要
• OpenACC プログラムのコンパイル– PGIコンパイラのメッセージの読み方• OpenACC プログラムの作成– 行列積、diffusion• OpenACC プログラムの最適化– NVIDIA visual profiler の使い方など

GPUプログラミング入門

140

OpenACCサンプル集•

Reedbush へログイン– $ ssh -Y reedbush.cc.u-tokyo.ac.jp –l txxxxx•

module のロード
– $ module load pgi/17.1– $ module load cuda/8.0.44•

ワークディレクトリに移動– $ cdw•

OpenACC_samples のコピー•

OpenACC_samples へ移動– $ cp /lustre/gt00h/z30108/OpenACC_samples.tar.gz .
– $ tar zxvf OpenACC_samples.tar.gz– $ cd OpenACC_samples– $ lsC/ F/#CとFortran好きな方を選択
GPUプログラミング入門

141

PGIコンパイラによるメッセージの確認• コンパイラメッセージの確認はOpenACCでは極めて重要– OpenMP と違い、
• 保守的に並列化するため、本来並列化できるプログラムも並列化されないことがある• 並列化すべきループが複数あるため、どのループにどの粒度(gang, worker, vector)が割り付けられたかしるため• ターゲットデバイスの性質上、立ち上げるべきスレッド数が自明に決まらず、スレッドがいくつ立ち上がったか知るため– メッセージを見て、プログラムを適宜修正する• コンパイラメッセージ出力方法– コンパイラオプションに -Minfo=accel をつける

GPUプログラミング入門

142

PGIコンパイラによるメッセージの確認

8.
9.
10.
11.
12.

• OpenACC_samples を利用

13.

• $ make acc_compute

15.

14.

16.
17.

ソースコード

18.
19.
20.

コンパイラメッセージ(fortran)

21.

subroutine acc_kernels()
double precision :: A(N,N), B(N,N)
double precision :: alpha = 1.0
integer :: i, j
A(:,:) = 1.0
B(:,:) = 0.0!$acc kernels
do j = 1, N
do i = 1, N
B(i,j) = alpha * A(i,j)
end doend do!$acc end kernels
end subroutine acc_kernels

pgfortran -O3 -acc -Minfo=accel -ta=tesla,cc60 -Mpreprocess acc_compute.f90-o acc_computeacc_kernels:
14, Generating implicit copyin(a(:,:))
Generating implicit copyout(b(:,:))
15, Loop is parallelizable
16, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
15, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
16, !$acc loop gang, vector(32) ! blockidx%x threadidx%x….

GPUプログラミング入門

143

PGIコンパイラによるメッセージの確認

8.
9.
10.
11.
12.

• OpenACC_samples を利用

13.

• $ make acc_compute

15.

14.

16.
17.

サブルーチン名

ソースコード

18.
19.
20.

コンパイラメッセージ(fortran)

21.

subroutine acc_kernels()
double precision :: A(N,N), B(N,N)
double precision :: alpha = 1.0
integer :: i, j
A(:,:) = 1.0
B(:,:) = 0.0!$acc kernels
do j = 1, N
do i = 1, N
B(i,j) = alpha * A(i,j)
end doend do!$acc end kernels
end subroutine acc_kernels

pgfortran -O3 -acc -Minfo=accel -ta=tesla,cc60 -Mpreprocess acc_compute.f90-o acc_compute
配列aはcopyin, bはcopyoutとして扱われます
acc_kernels:
14, Generating implicit copyin(a(:,:))
Generating implicit copyout(b(:,:))
15, Loop is parallelizable
16, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
15, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
16, !$acc loop gang, vector(32) ! blockidx%x threadidx%x….

15, 16行目の２重ループは(32x4)のスレッドでブロック分割して扱います。
GPUプログラミング入門

144

PGIコンパイラによるメッセージの確認• OpenACC_samples を利用• $ make acc_compute

40.
41.
42.

50.
51.
52.
53.
54.

ソースコード(C) 55.
コンパイラメッセージ(C)56.

void acc_kernels(double *A, double *B){
double alpha = 1.0;
int i,j;/ * A と B 初期化 */#pragma acc kernels
for(j = 0;j < N;j++){for(i = 0;i < N;i++){B[i+j*N] = alpha * A[i+j*N];}}}

pgcc -O3 -acc -Minfo=accel -ta=tesla,cc60 -Mcuda acc_compute.c -o acc_computeacc_kernels:
配列aはcopyin,
50, Generating implicit copy(B[:1000000])bはcopyとして扱われます
Generating implicit copyin(A[:1000000])
51, Loop carried dependence of B-> prevents parallelization
Loop carried backward dependence of B-> prevents vectorization
Complex loop carried dependence of B->,A-> prevents parallelization
Accelerator scalar kernel generated
Accelerator kernel generated ループ伝搬依存が見つかったので並列化しませGenerating Tesla codeんの意。ポインタAとBが同じ領域を指していること
51, #pragma acc loop seqを警戒して、並列化しない。
52, #pragma acc loop seq
52, Complex loop carried dependence of B->,A-> prevents parallelization
GPUプログラミング入門

145

PGIコンパイラによるメッセージの確認• OpenACC_samples を利用• $ make acc_compute

59.

60.
61.

69.
70.
71.
72.

ソースコード(C) 73.
コンパイラメッセージ(C)

74.
75.

void acc_kernels(double *restrict A,
double *restrict B){
double alpha = 1.0;
int i,j;/ * A と B 初期化 */#pragma acc kernels
for(j = 0;j < N;j++){for(i = 0;i < N;i++){B[i+j*N] = alpha * A[i+j*N];}}}

acc_kernels_restrict:
69, Generating implicit copy(B[:1000000])
Generating implicit copyin(A[:1000000])
70, Loop carried dependence of B-> prevents parallelization
Loop carried backward dependence of B-> prevents vectorization
71, Loop is parallelizable
2次元配列を1次元化して扱っているため、i+j*N
Accelerator kernel generatedが実は同じ場所を指す可能性を考慮し、70行
Generating Tesla code
目を並列化していない。
70, #pragma acc loop seq
71, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */GPUプログラミング入門

146

PGIコンパイラによるメッセージの確認• OpenACC_samples を利用

78.

79.
80.

88.
89.

• $ make acc_compute

90.
91.
92.

ソースコード(C)コンパイラメッセージ(C)

93.
94.
95.
96.

void acc_kernels(double *restrict A,
double *restrict B){
double alpha = 1.0;
int i,j;/ * A と B 初期化 */#pragma acc kernels#pragma acc loop independentfor(j = 0;j < N;j++){#pragma acc loop independentfor(i = 0;i < N;i++){B[i+j*N] = alpha * A[i+j*N];}}}

acc_kernels_independent:
88, Generating implicit copy(B[:1000000])
Generating implicit copyin(A[:1000000])
90, Loop is parallelizable
92, Loop is parallelizable
Accelerator kernel generated
ようやくまともに並列化
Generating Tesla code
90, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */92, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */GPUプログラミング入門

147

実習1
コンパイラのメッセージを確認しよう• acc_compute.f90 or acc_compute.c のソースと、コンパイル
メッセージを見比べてください– コンパイル ： $ make acc_compute• 注目点– parallel / kernels での違い
– 間接参照があった際のコンパイラメッセージ• acc_kernels_BAD_indirect_reference と
acc_kernels_indirect_reference の比較

GPUプログラミング入門

148

PGI_ACC_TIME によるOpenACC 実行の確認• PGI環境の場合、OpenACC プログラムが実行されているかを確認するには、環境変数PGI_ACC_TIME を使うのが簡単• 使い方 (一般的なLinux環境、またはインタラクティブジョブ実行時)– $ export PGI_ACC_TIME=1– $ (プログラムの実行)• 一般的なスパコン環境では、ジョブの中で環境変数を設定する必要が有る– ジョブスクリプト中に書いてある

GPUプログラミング入門

149

PGI_ACC_TIME による
OpenACC 実行の確認••

OpenACC_samples を利用$ qsub acc_compute.sh– 実行が終わると以下ができる• acc_compute.sh.eXXXXX (標準エラー
出力)• acc_compute.sh.oXXXXX (標準出力)

•

40.
41.
42.

50.
51.
52.
53.
54.
55.
56.

void acc_kernels(double *A, double *B){
double alpha = 1.0;
int i,j;/ * A と B 初期化 */#pragma acc kernels
for(j = 0;j < N;j++){for(i = 0;i < N;i++){B[i+j*N] = alpha * A[i+j*N];}}}$ less acc_compute.sh.eXXXXXPGI_ACC_TIME による出力メッセージ
Accelerator Kernel Timing data
/lustre/pz0108/z30108/OpenACC_samples/C/acc_compute.c
acc_kernels NVIDIA devicenum=0time(us): 149,101
50: compute region reached 1 time
51: kernel launched 1 time
↓カーネル実行時間← 起動したスレッド数
grid: [1] block: [1]device time(us): total=140,552 max=140,552 min=140,552 avg=140,552
elapsed time(us): total=140,611 max=140,611 min=140,611 avg=140,611
50: data region reached 2 times
50: data copyin transfers: 2← データ移動のdevice time(us): total=3,742 max=3,052 min=690 avg=1,871
56: data copyout transfers: 1
回数・時間
device time(us): total=4,807 max=4,807 min=4,807 avg=4,807
GPUプログラミング入門
150

実習2
OpenACCプログラムを実行してみよう• acc_data.f90 or acc_data.c のソースと、コンパイルメッセージを見比べてください– コンパイル ： $ make acc_data• プログラムを実行し、PGI_ACC_TIMEの出力を確認してください– 実行• $ qsub acc_data.sh– バッチジョブ実行が終了すると acc_data.sh.oXXXXXX (標準出力),
acc_data.sh.eXXXXXX (標準エラー出力) の2ファイルが出来る– PGI_ACC_TIME の出力確認• $ less acc_data.sh.eXXXXXX #標準エラーの方• 注目点– acc_data_copy, acc_data_copyinout の実行時間の違い
GPUプログラミング入門

151

実習3
OpenACCプログラムを作ろう• 行列積のOpenACC化• matmul.f90 または matmul.c を用いる。
– acc_matmul ルーチンにOpenACC指示文を加えてください。
– 注意：コンパイラの出力をよく見てください• よく見るエラー例(C)：– Accelerator restriction: size of the GPU copy of C,B is unknown
出力例====== OpenACC matmul program ======
1024 * 1024 matrix
check result...OKelapsed time[sec] :0.00708
FLOPS[GFlops]: 302.85417

GPUプログラミング入門

152

実習4
OpenACCプログラムを速くしよう• 拡散方程式のプログラムの高速化• diffusion.f90 または diffusion.c を用いる。
– 既にOpenACC化されていますが、実行してみると…
出力例(C)• 確認手順(nx, ny, nz) = (128, 128, 128)
elapsed time : 574.247 (s)flops: 0.078 (GFlops)
throughput: 0.096 (GB/s)
accuracy: 4.443942e-06
count: 1638

すごく遅い(特にC)

1. きちんと並列化されているか？ → コンパイラメッセージ
2. 無駄にデータ転送してないか？ → PGI_ACC_TIME
3. どっちもOKだけどなお遅い？ → NVIDIA Visual Profiler (nvvp)•

https://reedbush-www.cc.u-tokyo.ac.jp/session_login.cgi に詳細
資料あり
GPUプログラミング入門

153

4. プロファイラ（nvvp）を使う
NVIDIA CUDA Visual Profiler• （Windowsユーザのみ）準備– Cygwinインストール時にxorg-serverとxinitをインストールしておく– 「XWin Server」を起動する• タスクトレイにアイコンが2つ増える– 緑の線の入った方のアイコンから「Cygwin Terminal」を起動する• 使用– Reedbushへssh接続する際に-Yオプションを付けておく• ssh -Y reedbush.cc.u-tokyo.ac.jp -l txxxxx– module設定を行い、nvvpを起動する• module load cudaを設定したあとでnvvpコマンドを実行

154

Reedbush での NVVP の使い方
1. インタラクティブノードを使う––

nvvpを起動し、nvvpの中からアプリケーションを起動
講習会アカウントでは使えません

2. 計算ノードでデータ収集、ログインノードでnvvp起動––

スパコン環境で一般的なやり方。今回はこれネットワーク越しに画面転送するので遅い

3. 計算ノードでデータ収集、自分の端末にnvvpインストールして起動–––速い。よく使う人にオススメ
scp で収集データを自分の端末に持って来て起動
https://developer.nvidia.com/nvidia-visual-profiler•

ただし1GB以上あるので、今やるとおそらくネットワークがパンクしま
す…
GPUプログラミング入門

155

Reedbush での nvvp の使い方
• ログインして module をロード
– $ ssh -Y reedbush.cc.u-tokyo.ac.jp -l txxxxx– $ module load cuda pgi/17.1– $ cd OpenACC_samples/C または F• nvprof コマンドを使ってデータを収集– ジョブスクリプトの中に書く• nvvp.sh 参照– 今回は予め取っておいたデータを使う• diffusion.nvp (diffusionを10ステップだけ実行したもの)• nvvp の起動– $ nvvp
GPUプログラミング入門

156

①Fileをクリック

GPUプログラミング入門

157

②Importをクリック③ Nvprof を選択④ NextGPUプログラミング入門

158

⑤ Multiple processesに変更⑦ Browse⑥ NextGPUプログラミング入門

159

⑧ ファイル選択画面が開くので、nvvpで作成したプロファ
イリングデータを選択し、OK
ファイルシステム→上の「場所」欄に/lustre/gt00h/ユーザー名/OpenACC_samples/…
とするのが速いか…？

⑨ diffusion.nvpを選択できたら Finish
GPUプログラミング入門

160

ここで拡大この辺がメインの計算部分

GPUプログラミング入門

161

計算部分をクリックすると、カーネルの情報が得られる

ここをクリックするとプロファイラがヒントをくれる

GPUプログラミング入門

162

Q&A
• アカウントは1週間有効です。
• 資料のPDF版はWEBページにございます。
– http://www.cc.u-tokyo.ac.jp/support/kosyu/75/• アンケートへの協力をお願いします。

GPUプログラミング入門

163