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

スケジュール• 09：30 - 10：00 受付• 10：00 - 12：00 Reedbush-Hログイン、スパコンの使い方• 13：30 - 14：20 GPUのアーキテクチャ• 14：30 - 15：20 OpenACC入門• 15：30 - 16：20 OpenACC演習• 16：30 - 18：00 FDTD法による電磁波伝搬計算
2

スケジュール（改）• 09：30 - 10：00 受付• 10：00 - 12：00 Reedbush-Hログイン、スパコンの使い方• 13：30 - 14：30 GPU・OpenACC入門 （座学）
• 14：45 - 18：00 OpenACC演習 （適宜休憩）– 簡単なプログラムのOpenACC並列化– 3次元拡散方程式– FDTD法による電磁波伝搬計算– 質問など3

本講習会について• 本講習会の趣旨
ü 主に自分のアプリケーションのGPU上での実行を考えている人を対象とした、GPUプログラミングの入門
ü 入門に相応しいOpenACCを用いたプログラミング技能の習得• その他の講習会
https://www.cc.u-tokyo.ac.jp/support/kosyu/
ü 5/31 にはディープラーニングに特化した講習会が開催されます。
ü 後期には複数GPUの利用に関する講習会も開催予定です。（時期未定）• スパコンイベント情報メール配信サービス
https://regist.cc.u-tokyo.ac.jp/announce/
ü 講習会や研究会の案内、トライアルユースの実施のお知らせなどを配信しています。
4

事前準備
1.
2.

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

Reedbush 利用上の注意（1）• ディレクトリの扱いについて– ログイン時のディレクトリ（/home/gt00/txxxxx）にはログイン時に最低限必要なファイルのみを置く– 作業には/lustre以下のディレクトリ（/lustre/gt00/txxxxx）を使う– /homeに置いたファイルは計算ノードから参照できない– cdw コマンドでLustreファイルシステムに移動できます。

• 実行中のジョブの確認– ジョブ投入はqsubで行うが、ジョブ確認は qstat ではなく rbstat

6

Reedbush 利用上の注意（２）• コンパイルおよび実行のための環境準備
ü コンパイルおよび実行のための環境を準備するために module コマンドを使用する。これによって様々な環境を簡単に切り替えて使用できる。
$ module load <module_name>

モジュール名 <module_name> のモジュールをロードして環境を準備。環境変数PATHなどが設定される。
$ module avail

使用可能なモジュール一覧を表示する。
$ module list

使用中のモジュールを表示する。

7

モジュールの切り替え• PGIコンパイラ（OpenACCやCUDA Fortran）を使う場合$ module load pgi• CUDA開発環境を使う場合$ module load cuda

• Intelコンパイラを使う場合$ module load intel• MPIを使う場合$ module load openmpi/gdr/2.1.1/{gnu,intel,pgi}$ module load mvapich2/gdr/2.2/{gnu,intel,pgi}

ü PGIなどのコンパイラに追加して load

• モジュールはジョブ実行時にもコンパイル時と同じものを load する。
• 組み合わせて利用できる
8

サンプルコードのコンパイル•

Reedbush へのログイン$ ssh -Y txxxxx@reedbush.cc.u-tokyo.ac.jp

ü txxxxx 各自の利用者番号（アカウント）に置き換えてください。
ü -Y をつけてください。

•

cdw コマンドで Lustreファイルシステムへ移動する。
$ cdw

•

自分のディレクトリにサンプルコードをコピーする。
$ cp /lustre/gt00/share/openacc_samples.tar.gz .

•サンプルコードを展開する。
$ tar zxvf openacc_samples.tar.gz•

サンプルコードへ移動する。
$ cd openacc_samples•

モジュールをロードする。
$ module load pgi•

コンパイルする。
$ cd openacc_hello/01_hello_acc$ make•

実行ファイルがてきていることを確認する。
$ ls9プログラムの実行• ジョブとして投入し、実行する。
$ qsub ./run.sh• 投入されたジョブを確認する。
$ rbstat• 実行が終了すると、以下のファイルが生成される。
run.sh.o??????
run.sh.e??????

（?????? は数字）

• 上記の標準出力ファイルの中身を確認する。
$ cat run.sh.o??????

• 必要に応じて、上記のエラー出力ファイルの中身を確認する。
$ cat run.sh.e??????

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

11

東大センターのスパコン
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 ?

12

2システム運用中• Reedbush (SGI, Intel BDW + NVIDIA P100 (Pascal))– データ解析・シミュレーション融合スーパーコンピュータ– 3.361 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で9位 (2017年11月) (日本で2位)– Omni-Path アーキテクチャ, DDN IME (Burst Buffer)13東京大学情報基盤センター
スパコン（1/2）
Reedbush （SGI Rackable クラスタシステム ）Reedbush-U (2016/7/1 ~ )• 理論性能：508TFlops• ノード数：420• ノード構成：Intel Xeon Broadwell x2
Reedbush-H (2017/3/1 ~ )• 理論性能：1418TFlops• ノード数：120• ノード構成：Intel Xeon Broadwell x2 + NVIDIA P100 GPU x2
Reedbush-L (2017/10/1 ~ )• 理論性能：1435TFlops• ノード数：64• ノード構成：Intel Xeon Broadwell x2 + NVIDIA P100 GPU x4
14

15

東京大学情報基盤センター
スパコン（1/2）
筑波大学計算科学研究センターと共同運用
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日正式運用開始

15

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


128GB

DDR4
DDR4
DDR4
DDR4

G3 x16

76.8GB/sIntel Xeon
E5-2695 v4(BroadwellEP)

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

DDR4
DDR4
DDR4
DDR4


128GB

76.8GB/s

15.7 GB/s

IB EDR
HCA

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

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
18

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

DDR4
DDR4
DDR4
DDR4

76.8GB/s
QPIIntel 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/sIntel Xeon
E5-2695 v4(BroadwellEP)

20 GB/s
NVLinK
NVLinKNVIDIA
PascalIB FDR
HCA

20 GB/s
EDR switch
EDR
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 shown20

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





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

東大情報基盤センターOakforest-PACSスーパーコン
ピュータシステムの料金表（2018年4月1日）• パーソナルコース（年間） 個人研究者（大学・公共機関）のみ– 1口 100,000円 = 17,280トークン(2ノード年相当)– 基準ノード数※ = 8ノード
– 最大ノード数 = 2,048ノード• グループコース– １口 400,000円 (企業 480,000円) = 69,120トークン(8ノード年相当)– 基準ノード数※ = 口数×8ノード– 最大ノード数 = 2,048ノード※トークン消費について• 基準ノードまでは、トークン消費係数が1.0• 基準ノードを超えると、超えた分は、消費係数が2.0になる• 大学等のユーザはReedbushとの相互トークン移行も可能

22

東大情報基盤センターReedbushスーパーコンピュータ
システムの料金表（2018年4月1日 ）• パーソナルコース（年間） 個人研究者（大学・公共機関）のみ– 1口 150,000円 = 17,280トークン• グループコース– １口 300,000円 (企業 480,000円) = 69,120トークン
最大ノード数

基準ノード数

トークン消費係数(基準ノード超分は2.0倍)

RB-U

128

口数×4

1.0

RB-H

32

口数×1

2.5

RB-L

16

口数×1

4.0

• 大学等のユーザはReedbushとの相互トークン移行も可能• ノード固定もあり
23

トライアルユース制度について•安価に当センターの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回実
施）• 双方のコースともに、簡易な利用報告書の提出が必要
24

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

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

What’s GPU ?
ü Graphics Processing Unit
ü もともと PC の3D描画専用の装置
ü パソコンの部品として量産されてる。

= 非常に安価

3D GameComputer Graphics

GPU
http://www.nvidia.co.jp

27

GPUコンピューティング•

GPUはグラフィックスやゲームの画像計算のために進化を続けている。

•

CPUがコア数が2-12個程度に対し、GPUは1000以上のコアがある。

•

GPUを一般のアプリケーションの高速化に利用することを「GPUコンピューティング」「GPGPU (General Purpose computation on GPU)」などという。

•

2007年にNVIDIA社のCUDA言語がリリースされて大きく発展•

ここ数年、ディープラーニング（深層学習）、機械学習、AI（人工知能）などでも注目を浴びている。

28

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

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

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

Tesla：最初のHPC向けGPU、TSUBAME1.2など
Fermi：2世代目、TSUBAME2.0など

3.

Kepler：現在HPCにて多く利用、TSUBAME2.5など

4.
5.

Maxwell：コンシューマ向けのみ
Pascal：最新GPU、Reedbush-Hに搭載

6.

Volta：次世代GPU、アメリカの国立研究所で大規模運用予定• ECCメモリ、FMA演算、L1 L2 キャッシュ• シャッフル命令、Dynamic Parallelism、Hyper-Q• HBM2、半精度演算、 NVLink、倍精度atomicAdd など• Tensor Coreなど

29GPUの特徴• なぜGPUが使われるのか？
ü 性能高いü NVIDIA P100 (Reedbush-H) 5,304 GFlops
ü Intel Broadwell (Reedbush-U) 604.8 GFlops
ü Intel Xeon Phi (Oakforest-PACS) 3,046.4 GFLopsü 消費電力低い
ü スパコンに搭載されている• コンピュータに取り付ける増設ボード
ü 単体では動作できず、CPUからの指令が必要。
ü CPUとGPUはメモリが異なるため、メモリ間のデータ移動のプログラミングも必要。

• 多数の小さなコア（1000以上）を搭載。これを有効に活用するため「並列計算」が必須。
ü CPUは大きなコア：分岐予測、パイプライン処理などできる。逐次処理が得意。
ü GPUは小さなコア：CPUの持つ機能がほとんどないか、全くない。
30

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

T
T

/ 

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

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

88.8%

●



9/3=3

●

9/2=4.5 ≠6

32

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

GPUは苦手

仕事１
仕事２
仕事３

仕事４

仕事５

時間
33

データ並列• データを分割することで並列化する。
ü データは異なるが計算の手続きは同じ。

• データ並列の例：手分けをして算数ドリルを解く
ü 数字だけ異なるが計算の手続きは同じ。

2+1=
3 + 19 =4-6=
5+3=-2 + 10 =1+2=

12 – 88 =- 20 + 29 =4 + 10 =1 + 10 =1 + 13 =0+0=

34 + 3 =1+2=-3 + 2 =-10 – 10 =
2+1=
1 -10

2+4=
99 – 72 =2 + 10 =3 – 11 =12 – 34 =45 + 19 =GPUではこれが原則！
34

抑えておくべきGPUの特徴• 超並列計算が必須• CPUと独立のGPUメモリ• 階層的スレッド管理とコミュニケーション• Warp 単位の実行• コアレスドアクセス

35

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

P100 whitepaper より

36

参考：NVIDIA Tesla P100 の SM

P100 whitepaper より

37

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

メモリread終了

1core=1スレッドのとき
1core=Nスレッドのとき
38

CPUと独立のGPUメモリ

相対的に遅い１．必要なデータを送る
ノードの外へ
~20GB/s

バス (PCIe など)
CPU
OSが動いているGPU~32GB/s

OSは存在しない

２．計算を行う

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

メインメモリデバイス
メモリ

物理的に独立39

階層的スレッド管理とコミュニ
ケーション

スレッドブロック• 階層的なコア/スレッド管理– P100は56 SMを持ち、1 SMは64 CUDA
coreを持つ。トータル3584 CUDA core
– 1 SMが複数のスレッドブロックを担当
し、1 CUDA core が複数スレッドを担当• スレッド間のコミュニケーション– 同一スレッドブロック内のスレッドは高速コミュニケーション可能– 異なるスレッドブロックに属するスレッ
ド間はコミュニケーションが低速• いったんメモリに書き出したり、CPUに処理を戻さなくてはならない

cited from : http://cudaprogramming.blogspot.jp/2012/12/thread-hierarchy-in-cudaprogramming.html

40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

…

OK！

9

配列 B

3

1

…

NG！
41

Warp内分岐

CUDA 8 以前のバージョン（本講習会はCUDA 8 準拠）• Divergent Branch– Warp 内で分岐すること。Warp単位の分岐ならOK。
Divergent branchなし::
if ( TRUE ) {::} else {::}::else 部分は実行せずジャンプ

Divergent branchあり::
if ( 奇数スレッド ) {::} else {::}::一部スレッドを眠らせて全分岐を実行
42
最悪ケースでは32倍のコストコアレスドアクセス• 同じWarp内のスレッド（連続するスレッド）は近いメモリアドレスへアクセスすると効率的ü コアレスドアクセス（coalesced access）と呼ぶ
ü メモリアクセスは128 Byte 単位で行われる。128 Byte に収まれば1回のア
クセス、超えれば128 Byte アクセスをその分繰り返す。

128 byte x 1回のメモリアクセス
Address 128

160

192

224

256

...

15 16

...

31

160

192

224

256

...

15 16

...

31

Thread0

128 byte x 2回のメモリアクセス
Address 128

Thread043

ストライドアクセスがあるとどうなるか• GPUはストライドアクセスに弱い！

Memory Bandwidth [GB/sec]

600

void AoS_STREAM_Triad(STREAM_TYPE scalar){ssize_t i,j;#pragma omp parallel for private(i,j)#pragma acc kernels present(a_aos[0:STREAM_ARRAY_SIZE] ¥,b_aos[0:STREAM_ARRAY_SIZE],c_aos[0:STREAM_ARRAY_SIZE])#pragma acc loop gang vector independentfor (j=0; j<STREAM_ARRAY_SIZE/STRIDE; j++)for (i=0; i<STRIDE; i++)a_aos[j*STRIDE+i] = b_aos[j*STRIDE+i]+scalar*c_aos[j*STRIDE+i];}

ストライドアクセス付き stream triad

500
400
300
200
100
0
0

20

40

60

80
Stride

100

120

140
44

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

45

GPUコンピューティングの方法• ライブラリの利用（CUFFT, CUBLAS など）

簡単

ü GPU用ライブラリを呼ぶだけで、すぐに利用できる。
ü ライブラリ以外の部分は高速化されない。

• 指示文ベース（OpenACC）
ü 指示文（ディレクティブ）を挿入するだけである程度高速化。
ü 既存のソースコードを活用できる。

• プログラミング言語（CUDA、OpenCLなど）
ü GPUの性能を最大限に活用。
ü プログラミングにはGPGPU用言語を使用する必要あり。

難しい

46

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

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

47

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

本講習会ではPGIコンパイラを用いる
48

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

OpenMP

OpenACC

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

}

CPU

CPU

CPU

デバイス
49

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

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

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

計
算
コ
ア

• 共有メモリ

計算コ
ア

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

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

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

51

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

OpenACC と OpenMP の比較
デフォルトでの変数の扱い• OpenMP– 全部 shared• OpenACC– スカラ変数： firstprivate or private– 配列： shared• プログラム上のparallel/kernels構文に差し掛かった時、OpenACC
コンパイラは実行に必要なデータを自動で転送する• 正しく転送されないこともある。自分で書くべき• 構文に差し掛かるたびに転送が行われる(非効率)。後述のdata
指示文を用いて自分で書くべき• 配列はデバイスに確保される (shared的振る舞い)• 配列変数をprivateに扱うためには private 指示節使う
53

OpenACC の指示文

54

はじめてのOpenACCコード
openacc_hello/01_hello_acc

CPU

int main(){const int n = 1000;
float *a = malloc(n*sizeof(float));float *b = malloc(n*sizeof(float));float c = 2.0;
for (int i=0; i<n; i++) {a[i] = 10.0;}

a b#pragma acc data copyin(a[0:n]), copyout(b[0:n])#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;}
double sum = 0;
for (int i=0; i<n; i++) {sum += b[i];}fprintf(stdout, "%f¥n", sum/n);free(a); free(b);return 0;

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

GPUから
copyout

解放}55

はじめてのOpenACCコード
openacc_hello/01_hello_acc

CPU

int main(){const int n = 1000;
float *a = malloc(n*sizeof(float));float *b = malloc(n*sizeof(float));float c = 2.0;
for (int i=0; i<n; i++) {a[i] = 10.0;}

a b#pragma acc data copyin(a[0:n]), copyout(b[0:n])#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;}
double sum = 0;
for (int i=0; i<n; i++) {sum += b[i];}fprintf(stdout,"%f¥n", sum/n);
コード上同じa, b であっても、原則として
free(a); free(b);ホストコードはホストメモリで確保された
a, b 、GPUで実行さ
return0;
れる並列領域（カーネル）はデバイスメモリで確保された a, b}

を参照しにいく。

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

GPUから
copyout

解放56

参考：OpenACC 化とCUDA化の比較// OpenACC// CUDA

void calc(int n, const float *a,const float *b, float c, float *d){#pragma acc kernels present(a, b, d)#pragma acc loop independentfor (int i=0; i<n; i++) {d[i] = a[i] + c*b[i];}}__global__
void calc_kernel(int n, const float *a, const float *b,
float c, float *d){const int i= blockIdx.x * blockDim.x + threadIdx.x;

kernelint main(){...
#pragma acc data copyin(a[0:n], b[0:n]) copyout(d[0:n]){calc(n, a, b, c, d);}...
}ü kernels 指示文でGPUでの実行領域を指定。

if (i < n) {d[i] = a[i] + c*b[i];}}void calc(int n, const float *a, const float *b, float c,float *d){dim3 threads(128);
dim3 blocks((n + threads.x - 1) / threads.x);
calc_kernel<<<blocks, threads>>>(n, a, b, c, d);cudaThreadSynchronize();}
int main(){...
float *a_d, *b_d, *d_d;
cudaMalloc(&a_d, n*sizeof(float));cudaMalloc(&b_d, n*sizeof(float));cudaMalloc(&d_d, n*sizeof(float));ü loop 指示文で並列処理の最適化。

cudaMemcpy(a_d, a, n*sizeof(float), cudaMemcpyDefault);
cudaMemcpy(b_d, b, n*sizeof(float), cudaMemcpyDefault);
cudaMemcpy(d_d, d, n*sizeof(float), cudaMemcpyDefault);

ü data 指示文でデータ転送を制御。
kernels 指示文でデータ転送を自動的に行うこともできる。

calc(n, a_d, b_d, c, d_d);cudaMemcpy(d, d_d, n*sizeof(float), cudaMemcpyDefault);
...
}

57

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

ホスト

デバイス

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

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

59

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

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 program61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でやってください」
62

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

async
wait
device_typeifdefault(none)
copy…

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

async
wait
device_typeifdefault(none)
copy…
num_gangs
num_workers
vector_length
reduction
privatefirstprivate
63kernels/parallel 指示文：指示節
kernels

parallel
非同期実行に用いる。

実行デバイス毎にパラメータを調整

データ指示文の機能を使える

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

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

async
wait
device_typeifdefault(none)
copy…
num_gangs
num_workers
vector_length
reduction
privatefirstprivate
64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#pragma acc kernels copy(src[0:N] ¥dis[0:N])for(i = 0;i < N;i++){dis[i] = src[i];}end subroutine copy}

65

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 の領域が解放される66

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

67

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

ホスト~32GB/s

デバイス
２．計算を行う

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

データ移動最適化指示文：data, enter/exit data
• デバイス側で必要なデータと範囲を指定– Allocate, Memcpy, Deallocate を行う• data 指示文 (推奨)– Allocate + Memcpy (Host ⇄ Device) + Deallocate– 構造ブロックに対してのみ適用可• コードの見通しが良い• Enter data 指示文– Allocate + Memcpy (Host → Device)– Exit data とセット。構造ブロック以外にも使える• Exit data 指示文– Memcpy (Host ← Device) + Deallocate– enter data とセット。構造ブロック以外にも使える
69データ移動最適化指示文が必要なとき
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++){#pragma acc kernels copy(src[0:N] ¥dis[0:N])for(i = 0;i < N;i++){dis[i] = dis[i] + src[i];}}}Kernels をループで囲
むとどうなるか…

70

データ移動最適化指示文が必要なとき
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++){#pragma acc kernels copy(src[0:N] ¥dis[0:N])for(i = 0;i < N;i++){dis[i] = dis[i] + src[i];}}}Kernels をループで囲
むとどうなるか…

71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
#pragma 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が構造ブロックになってるのでなくても大丈夫だが)72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++){２．計算を行う
３．計算結果を返す#pragma 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が構造ブロックになってるのでなくても大丈夫だが)73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);}74

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

if
copy
copyincopyout
createpresent
present_or_...
deviceptr– CUDAなどと組み合わせる時に利用。
cudaMallocなどで確保済みのデータを指定し、OpenACCで扱い可とする

enter data
•
•
•
•
•
•

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

exit data
•
•
•
•
•

if
async
wait
copyout
delete75

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 と同一
76

データ移動指示文：データ転送範囲指定• 送受信するデータの範囲の指定– 部分配列の送受信が可能– 注意：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 data77

データ移動指示文：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

update
•
•
•
•
•
•
•

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

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

ホスト

デバイス

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

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

79

ループ指示文• loop 指示文– parallel/kernels中でしか使えない– 並列化したいループに指定– ループマッピングのパラメータ調整#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;}

配列の1要素が1スレッドで処理されるイメージ• 粒度(gang, worker, vector, seq)の指定• ある程度はコンパイラが自動で決定してくれるため、さほど気にしなくていい– データの独立性の保証 (independent)• コンパイラは保守的に解析するため、並列性を見切れないことがある。C言語ではほぼ必須。

– リダクション処理

80

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

NVIDIA GPUの構成
GPU– 最近は2階層目(SIMD)がある

デバイスメモリ• CUDAは block と thread の2階層
– NVIDA GPUも2階層

SM• 1 SMX に複数CUDA coreを搭載• 各コアはSMXのリソースを共有• OpenACCは3階層– gang：workerの集合 一番大きな単位– worker：vectorの集合

CUDA
コア

– vector：スレッドに相当する一番小さ
い処理単位
81

81

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

collapsegang
workervector
seq
auto
tile
device_type
independent
privatereduction82

ループ指示文：指示節
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

並列化するにはループ長の短すぎるループに使う83ループ指示文：指示節
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….
数値の指定は難しいので、最初はコンパイラ任せでいい84

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

collapsegang
workervector
seq
auto
tile
device_type
independent
privatereduction

ループがデータ独立であることを明示する
コンパイラが並列化できないと判断したときに使用#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;}

並列化可能（データ独立）なので、independent を指定（コンパイラは並列化可能とは判断してくれなかった）n

データ独立でない（並列化可能でない）例// これは正しくない#pragma acc kernels#pragma acc loop independentfor (int i=1; i<n; i++) {d[i] = d[i-1];}

85

ループ指示文：指示節
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)
演算子
対象とする変数
制限：スカラー変数のみ

86参考：リダクションで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
87

参考：リダクションで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;}88

参考：リダクションで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に任せた方が速い！
89

関数呼び出し指示文：routine• parallel/kernels領域内から関数を呼び出す場合、routine指
示文を使う– 特にcaller, calleeでファイルが違う場合に必要とされる#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++){max = vecsum(A[i*N]);}#pragma acc routine vector
double vecsum(double *A){
double x = 0;
#pragma acc loop reduction(+:x)
for(int j = 0;j < N;j++){x += A[j];}return x;}90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

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

91

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

領域内ではデバイス側のアドレスが使われる

92

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

93OpenACC と Unified Memory• Unified Memory とは…
– 物理的に別物のCPUとGPUのメモリをあたかも一つのメモリのように扱う機能– Pascal GPUではハードウェアサポート• ページフォルトが起こると勝手にマイグレーションしてくれる• OpenACC と Unified Memory– OpenACCにUnified Memoryを直接使う機能はない• PGIコンパイラではオプションを与えることで使える• pgfortran –acc –ta=tesla,managed– 使うとデータ指示文が無視され、代わりにUnified Memoryを使う94Unified Memoryのメリット・デメリット• メリット– データ移動の管理を任せられる– 複雑なデータ構造を簡単に扱える!$acc kernels!$acc loop independent
do il=1,kt!$acc loop independent
do it=1,ndt; itt=it+nstrtt-1
zbu(il)=zbu(il)+st_leafmtxp%st_lf(ip)%a1(it,il)*zu(itt)
enddo
enddo
!$acc end kernels• 本来はメモリ空間が分かれているため、ディープコピー問題が発生する↑こう書くだけで正しく動くのは、従来のCUDA
ユーザーからすると革命的• デメリット– ページ単位で転送するため、細かい転送が必要な場合には遅くなる– CPU側のメモリ管理を監視しているので、allocate, deallocateを繰り返すアプリではCPU側が極端に遅くなる• 今研究で使っているコードでは20倍近く遅くなった95アプリケーションの移植方法

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

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

2.

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

parallel/kernels指示文適用

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

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

既にOpenMP化されているアプリケーションのOpenACC化手順
1. !$omp parallel を!$acc kernelsに機械的に置き換え
2. (Unified Memory を使い、とりあえずGPU上で実行)– 本講習会では扱いません

3. データ指示文を用いて転送の最適
4. コンパイラのメッセージを見ながら、OpenACCカーネルの最適化
5. カーネルのCUDA化など–

スレッド間の相互作用が多いアプリケーションでは、sharedmemory や shuffle 命令を自由に使えるCUDAの方が圧倒的に有利
98

データ指示文による最適化手順
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化を始める

99

データ指示文による最適化手順
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

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

100

データ指示文による最適化手順
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

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

101

データ指示文による最適化手順
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

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

Q&A

103

実習
※今回の実習の例は、全てPGIコンパイラ17.1を使った際の例です

104

実習概要
• 簡単なOpenACC プログラムの実行– コンパイラのメッセージの読み方を学ぼう– OpenACCのプログラムを実行しよう• ３次元拡散方程式のOpenACC化– 自分でOpenACCの指示文を挿入してみよう• FDTD法による電磁波伝搬計算– より実際のアプリケーションに近いプログラム
をOpenACC化しよう

105

OpenACCコードのコンパイル•

PGIコンパイラによるコンパイル
ü ReedbushではOpenACCはPGIコンパイラで利用できます。
$ module load pgi/17.1$ pgcc -O3 -acc -Minfo=accel -ta=tesla,cc60 -c main.c-acc: OpenACCコードであることを指示-Minfo=accel:

OpenACC指示文からGPUコードが生成できたかどうか等のメッセージを出力する。
このメッセージがOpenACC化では大きなヒントになる。
-ta=tesla,cc60:

ターゲット・アーキテクチャの指定。NVIDIA GPU Teslaをターゲットとし、computecapability 6.0 (cc60) のコードを生成する。

•

Makefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load pgi/17.1$ make106簡単なOpenACCコード• サンプルコード: openacc_basic/
ü OpenACC指示文 kernels, data, loop を利用したコード
ü 計算内容は簡単な四則演算• ソースコードopenacc_basic/01_original

CPUコード。

openacc_basic/02_kernels

OpenACCコード。上にkernels指示文を追加。

openacc_basic/03_loop

OpenACCコード。上にloop指示文を追加

openacc_basic/04_data

OpenACCコード。上にdata指示文を明示的に追加。

openacc_basic/05_present

OpenACCコード。上でpresent指示節を使用。

openacc_basic/06_reduction OpenACCコード。上にreduction指示節を使用。

107

配列のインデックス計算
calc関数（C言語）

calcサブルーチン（Fortran）

void calc(unsigned int nx, unsigned int ny, const float*a, const float *b, float *c){for (unsigned int j=0; j<ny; j++) {for (unsigned int i=0; i<nx; i++) {const int ix = i + j*nx;
c[ix] += a[ix] + b[ix];}}}subroutine calc(nx, ny, a, b, c)
integer,intent(in) :: nx,ny
real,dimension(:),intent(in) :: a,breal,dimension(:),intent(out) :: cinteger :: i,j,ix
do j = 1, ny
do i = 1, nx
ix = i + (j-1) * nx
c(ix) = c(ix) + a(ix) + b(ix)end doend doend subroutine calcnyjix = j*nx + i

i
nx
108

簡単なOpenACC: CPUコード• CPUコードのコンパイルと実行ü 配列の平均値と実行時間が出力されています。
$ cd openacc_basic/01_original$ make$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
19.886 [sec]? の数字はジョブごとに変わります。
答えは常に3000.0
openacc_basic/01_original• 計算内容ü 配列 a、b、cをそれぞれ 1.0, 2.0, 0.0 で初期化
ü calc関数内で c += a *b を nt(=1000)回実行。
ü この実行時間を測定

109

簡単なOpenACC: kernels 指示文（１）• 02_kernelsコード: calc関数

openacc_basic/02_kernels

ü CPUコードにkernels 指示文の追加
C言語

Fortranvoid calc(unsigned int nx, unsigned int ny, const float*a, const float *b, float *c){#pragma acc kernels copy(a[0:n], b[0:n], c[0:n])for (unsigned int j=0; j<ny; j++) {for (unsigned int i=0; i<nx; i++) {const int ix = i + j*nx;
c[ix] += a[ix] + b[ix];}}}subroutine calc(nx, ny, a, b, c)
integer,intent(in) :: nx,ny
real,dimension(:),intent(in) :: a,breal,dimension(:),intent(out) :: cinteger :: i,j,ix
allocate, H -> D!$acc kernels copy(a,b,c)do j = 1, ny
do i = 1, nx
ix = i + (j-1) * nx
c(ix) = c(ix) + a(ix) + b(ix)end doend do!$acc end kernels
D->H, deallocate
end subroutine calc

ü kernels 指示文では data 指示文の機能が使えるü 上の場合は、copy を指定
ü カーネル前後でGPUとCPU間のメモリ転送が行われる。

110

簡単なOpenACC: kernels 指示文（２）• 02_kernelsコード:初期化

openacc_basic/02_kernels

ü CPUコードにkernels 指示文の追加
C言語
int main(int argc, char *argv[]){...
#pragma acc kernels copyout(b[0:n], c[0:n]){for (unsigned int i=0; i<n; i++) {b[i] = b0;}
for (unsigned int i=0; i<n; i++) {c[i] = 0.0;
}
}
...
}

b0 はスカラ変数のため自動的に各スレッドへコピーが渡される。

Fortranprogram main…
!$acc kernels copyout(b,c)
do i = 1, nb(i) = b0
end dodo i = 1, nc(i) = 0.0
end do!$acc end kernels…
end program main

allocate

D->H, deallocate

Fortranは配列のサイズ情報を持つので、大抵の場合サイズ指定はいらない。
111

簡単なOpenACC: kernels 指示文（３）• コンパイル
ü データの独立性がコンパイラにはわからず、並列化されない。
$ makepgcc -O3 -acc -Minfo=accel -ta=tesla,cc60 -c main.c
calc:
13, Generating copy(a[:n],c[:n],b[:n])14, Complex loop carried dependence of a-> prevents parallelization
Loop carried dependence due to exposed use of c[:n] prevents parallelization
Complex loop carried dependence of c->,b-> prevents parallelization
Accelerator scalar kernel generated
Accelerator kernel generated
Generating Tesla code
14, #pragma acc loop seq
15, #pragma acc loop seq
15, Complex loop carried dependence of a->,c->,b-> prevents parallelization
Loop carried dependence due to exposed use of c[:i1+n] prevents parallelization
main:
43, Generating copyout(c[:n],b[:n])45, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
45, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */48, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
48, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */pgcc -O3 -acc -Minfo=accel -ta=tesla,cc60
main.o-o run112簡単なOpenACC: loop 指示文（１）• 03_loopコード
ü 02_kernelsコードにloop independent の追加
void calc(unsigned int nx, unsigned int ny, const float *a,const float *b, float *c){const unsigned int n = nx * ny;#pragma acc kernels copy(a[0:n], b[0:n], c[0:n])#pragma acc loop independentfor (unsigned int j=0; j<ny; j++) {#pragma acc loop independentfor (unsigned int i=0; i<nx; i++) {const int ix = i + j*nx;
c[ix] += a[ix] + b[ix];}}}// main 関数内#pragma acc kernels copyout(b[0:n], c[0:n]){#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {b[i] = b0;}#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {c[i] = 0.0;}}

openacc_basic/03_loop

subroutine calc(nx, ny, a, b, c)
integer,intent(in) :: nx,ny
real,dimension(:),intent(in) :: a,breal,dimension(:),intent(out) :: cinteger :: i,j,ix!$acc kernels copy(a,b,c)!$acc loop independent
do j = 1, ny!$acc loop independent
do i = 1, nx
ix = i + (j-1) * nx
c(ix) = c(ix) + a(ix) + b(ix)end doend do!$acc end kernels
end subroutine calc!$acc kernels copyout(b,c)!$acc loop independent
do i = 1, nb(i) = b0
end do!$acc loop independent
do i = 1, nc(i) = 0.0
end do!$acc end kernels

113簡単なOpenACC: loop 指示文（２）• コンパイル

openacc_basic/03_loop

ü ループが並列化され、カーネルが生成された。
$ makepgcc -O3 -acc -Minfo=accel -ta=tesla,cc60 -c main.c
calc:
13, Generating copy(a[:n],c[:n],b[:n])15, Loop is parallelizable
17, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
15, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */17, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */main:
45, Generating copyout(c[:n],b[:n])48, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
48, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */52, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
52, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */114

簡単なOpenACC: loop 指示文（３）openacc_basic/03_loop• 03_loopコードの実行ü 答えは正しいが、実行時間が大変長い。
$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
70.414 [sec]ü ソースコードをみると、calc関数でカーネル前後にGPUとCPU間のデー
タ転送が発生する。この関数はmain関数から1000回呼ばれる。これが性能低下させている。
// main 関数内
const unsigned int nt = 1000;…
for (unsigned int icnt=0; icnt<nt; icnt++) {calc(nx, ny, a, b, c);}

1000回呼ばれる

void calc(unsigned int nx, unsigned int ny, const float *a, const float *b,
float *c){const unsigned int n = nx * ny;#pragma acc kernels copy(a[0:n], b[0:n], c[0:n]) allocate, H -> D#pragma acc loop independentfor (unsigned int j=0; j<ny; j++) {#pragma acc loop independentfor (unsigned int i=0; i<nx; i++) {const int ix = i + j*nx;
c[ix] += a[ix] + b[ix];D->H, deallocate
}
}
}

115

簡単なOpenACC: data指示文（１）openacc_basic/04_data• 04_dataコード
ü 03_loopにdata指示文追加// main関数内#pragma acc data copyin(a[0:n]) create(b[0:n]) copyout(c[0:n]){#pragma acc kernels copyout(b[0:n], c[0:n]){#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {b[i] = b0;}#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {c[i] = 0.0;}}

present として振舞う。

}

for (unsigned int icnt=0; icnt<nt; icnt++) {calc(nx, ny, a, b, c);}! mainプログラム内!$acc data copyin(a) create(b) copyout(c)!$acc kernels copyout(b,c)!$acc loop independent
do i = 1, nb(i) = b0
end do!$acc loop independent
do i = 1, nc(i) = 0.0
end do!$acc end kernels
do icnt = 1, nt
call calc(nx, ny, a, b, c)
end do!$acc end dataa: allocate, H -> D
b: allocate
c: allocate

a: deallocate
b: deallocate
c: D->H, deallocate

ü copy/copyin/copyout/create は既にデバイス上確保されているデータに対しては何もしない。present として振舞う。（OpenACC2.5以降）
ü 配列 a, b, c は利用用途に合わせた指示節を指定。
116

簡単なOpenACC: data指示文（２）• 04_dataコードの実行

openacc_basic/04_data

ü 答えは正しく、速度が上がった。
$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
1.174 [sec]117簡単なOpenACC: present指示節• 05_presentコード
ü 04_dataコードで present 指示節を使用

openacc_basic/05_present

void calc(unsigned int nx, unsigned int ny, const float *a, const float *b, float *c){const unsigned int n = nx * ny;#pragma acc kernels present(a, b, c)#pragma acc loop independentfor (unsigned int j=0; j<ny; j++) {#pragma acc loop independentfor (unsigned int i=0; i<nx; i++) {const int ix = i + j*nx;
c[ix] += a[ix] + b[ix];}present へ変更// main 関数内#pragma acc kernels{#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {b[i] = b0;}

指示節を削除ü データ転送の振る舞いは変化しないため、性能変化はなし。
ü present ではメモリ確保、データ転送をしないため、配列サイズの指定は不要。
ü コードとしては見通しがよい。

118

リダクション計算（１）• リダクション計算
ü 配列の全要素から一つの値を抽出
ü 総和、総積、最大値、最小値など
ü 出力が一つのため、並列化に工夫が必要（CUDAでの実装は煩雑）
double sum = 0.0;
for (unsigned int i=0; i<n; i++) {sum += array[i];
}

スレッド1

スレッド2 スレッド3

1

4

2

3

6

5
15

6

7

8
24

9

1.
2.
3.

6

15 24
4.

45

各スレッドが担当する
領域をリダクション
一時配列に移動
一時配列をリダクショ
ン
出力を得る
119

リダクション計算(２）• loop 指示文に reduction 指示節を指定
ü reduction 演算子と変数を組み合わせて指定
double sum = 0.0;
#pragma acc kernels#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += array[i];}• Reduction 指示節
ü acc loop reduction(+:sum)
ü 演算子と対象とする変数（スカラー変数）を指定する。

• 利用できる主な演算子と初期値
ü 演算子: +, 初期値: 0
ü 演算子: *, 初期値: 1
ü 演算子: max, 初期値: leastü 演算子: min, 初期値: largest
120簡単なOpenACC: reduction指示節（１）• 06_reductionコード
ü 05_presentコードで reductionを使用// main 関数内
for (unsigned int icnt=0; icnt<nt; icnt++) {calc(nx, ny, a, b, c);}#pragma acc kernels#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += c[i];}ü data 指示文で c を create に変更。
ü data指示文の範囲に注意。
ü リダクションコードが生成された。

openacc_basic/06_reduction

do icnt = 1, nt
call calc(nx, ny, a, b, c)
end dosum = 0!$acc kernels!$acc loop reduction(+:sum)do i = 1, n
sum = sum + c(i)end do!$acc end kernels!$acc end data

$ makepgcc -O3 -acc -Minfo=accel -ta=tesla,cc60 -c main.c
(省略)
main:
(省略)
67, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
67, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */Generating reduction(+:sum)121簡単なOpenACC: reduction指示節（２）• 06_reductionコードの実行

openacc_basic/06_reduction

ü 答えは正しく、速度が上がった。
ü 配列 c の転送が削減されたこと、リダクションがGPU上で行われることによる性能向上。
$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
1.089 [sec]122OpenACC化のステップのまとめ• OpenACC化のための3つの指示文の適用
ü kernels 指示文を用いてGPUで実行する領域を指定
ü data 指示文を用い、ホスト-デバイス 間の通信を最適化ü loop 指示文を用い、並列処理の指定#pragma acc data copyin(a[0:n]) create(b[0:n], c[0:n]){#pragma acc kernels{#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {b[i] = b0;}#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {c[i] = 0.0;}}
for (unsigned int icnt=0; icnt<nt; icnt++) {calc(nx, ny, a, b, c);}#pragma acc kernels#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += c[i];}}!$acc data copyin(a) create(b) copyout(c)!$acc kernels!$acc loop independent
do i = 1, nb(i) = b0
end do!$acc loop independent
do i = 1, nc(i) = 0.0
end do!$acc end kernels
do icnt = 1, nt
call calc(nx, ny, a, b, c)
end do!$acc kernels!$acc loop reduction(+:sum)do i = 1, n
sum = sum + c(i)end do!$acc end kernels!$acc end data

openacc_basic/06_reduction123

3次元拡散方程式のOpenACC化• 3次元拡散方程式のCPUコードにOpenACC の kernels, data,loop 指示文を追加し、GPUで高性能で実行しましょう。
– サンプルコード： openacc_diffusion/01_originalfor(int k = 0; k < nz; k++) {diffusion.c, diffusion3d 関数内
for (int j = 0; j < ny; j++) {for (int i = 0; i < nx; i++) {const int ix = nx*ny*k + nx*j + i;const int ip = i == nx - 1 ? ix : ix + 1;
const int im = i == 0
? ix : ix - 1;
const int jp = j == ny - 1 ? ix : ix + nx;
const int jm = j == 0
? ix : ix - nx;
const int kp = k == nz - 1 ? ix : ix + nx*ny;
const int km = k == 0
? ix : ix - nx*ny;
fn[ix] = cc*f[ix]+ ce*f[ip] + cw*f[im]+ cn*f[jp] + cs*f[jm]+ ct*f[kp] + cb*f[km];}}}openacc_diffusion/01_original124

拡散現象シミュレーション（１）• 拡散現象ü コップの中に赤インクを落とすと水中で拡がる
ü 次第に拡散し赤インクは拡がり、最後は均一な色になる。

• 拡散方程式のシミュレーション
ü 各点のインク濃度の時間変化を計算する
125

拡散現象シミュレーション（２）• データ構造ü 計算したい空間を格子に区切り、一般に配列で表す。
ü 計算は3次元であるが、C言語では1次元配列として確保することが一般的。
ü 2ステップ分の配列を使い、タイムステップを進める（ダブルバッファ）。

t=n

t = n+1

t = n+2
t = n+3• サンプルコードは、

時間ステップを進めるü 計算領域: nx * ny * nz （3次元）
ü 最大タイムステップ: ntとなっている。
126

拡散現象シミュレーション（３）• 2次元拡散方程式の離散化の一例

平均後の
自分自身の値

上下左右の値

1

0

0

0
1

0

0

2

0

0
2

0
8

0
2

0

j 3

0
1

0
8

64
32
20

0
8

0
1

4

0

0
2

0
8

0
2

0

5

0

0

0
1

0

0

1

2

3

4

5

i

自分自身の値の４倍

最初の状態
１回目の平均後２回目の平均後繰り返し平均化を行うと、インクが拡散します。
127127

拡散現象シミュレーション（４）• 2次元拡散方程式の計算例

128

CPUコード• CPUコードのコンパイルと実行$ cd openacc_diffusion/01_original$ make$ qsub ./run.sh# cat run.sh.o??????
time(0) = 0.00000
time( 100) = 0.00610
time( 200) = 0.01221...
time(1000) = 0.06104
time(1100) = 0.06714
time(1200) = 0.07324
time(1300) = 0.07935
time(1400) = 0.08545
time(1500) = 0.09155
time(1600) = 0.09766
Time =
20.564 [sec]Performance=
2.17 [GFlops]
Error[128][128][128] = 4.556413e-06

実行性能
解析解との誤差• OpenACCコードでは、どのくらいの実行性能が達成できるで
しょうか？
129

OpenACC化(0): Makefile の修正• Makefile に OpenACC をコンパイルするよう –acc などを追加
しましょうCC= pgccCXX = pgc++
GCC = gcc
RM = rm -f
MAKEDEPEND = makedepend
CFLAGS
GFLAGSCXXFLAGSLDFLAGS...

= -O3 -acc -Minfo=accel -ta=tesla,cc60= -Wall -O3 -std=c99= $(CFLAGS)=

130

OpenACC化(1): kernels• diffusion3d関数に kernelsを追加しましょう#pragma acc kernels copyin(f[0:nx*ny*nz]) copyout(fn[0:nx*ny*nz])
for(int k = 0; k < nz; k++) {for (int j = 0; j < ny; j++) {for (int i = 0; i < nx; i++) {const int ix = nx*ny*k + nx*j + i;const int ip = i == nx - 1 ? ix : ix + 1;
const int im = i == 0
? ix : ix - 1;
const int jp = j == ny - 1 ? ix : ix + nx;
const int jm = j == 0
? ix : ix - nx;
const int kp = k == nz - 1 ? ix : ix + nx*ny;
const int km = k == 0
? ix : ix - nx*ny;
fn[ix] = cc*f[ix]+ ce*f[ip] + cw*f[im]+ cn*f[jp] + cs*f[jm]+ ct*f[kp] + cb*f[km];}}}return (double)(nx*ny*nz)*13.0;

diffusion.c, diffusion3d 関数内}make して実行してみましょう。
131

OpenACC化(2): loop
• diffusion3d関数に loopを追加しましょう#pragma acc kernels copyin(f[0:nx*ny*nz]) copyout(fn[0:nx*ny*nz])#pragma acc loop independentfor(int k = 0; k < nz; k++) {#pragma acc loop independentfor (int j = 0; j < ny; j++) {#pragma acc loop independentfor (int i = 0; i < nx; i++) {const int ix = nx*ny*k + nx*j + i;const int ip = i == nx - 1 ? ix : ix + 1;
const int im = i == 0
? ix : ix - 1;
const int jp = j == ny - 1 ? ix : ix + nx;
const int jm = j == 0
? ix : ix - nx;
高速化よりも、
const int kp = k == nz - 1 ? ix : ix + nx*ny;
const int km = k == 0
? ix : ix - nx*ny; まずは正しい計
fn[ix] = cc*f[ix]+ ce*f[ip] + cw*f[im]+ cn*f[jp] + cs*f[jm]+ ct*f[kp] + cb*f[km];}}}return (double)(nx*ny*nz)*13.0;}

算を行うコードを保つことが大事です。
末端の関数から修正を進めます。

diffusion.c, diffusion3d 関数内make してジョブ投入 qsub ./run.shしてみましょう。
132
遅いですが実行できます。

OpenACC化(3): データ転送の最適化(1)• diffusion3d関数で present とし、main関数で data を追加#pragma acc kernels present(f, fn)diffusion.c, diffusion3d 関数内#pragma acc loop independentfor(int k = 0; k < nz; k++) {#pragma acc loop independentfor (int j = 0; j < ny; j++) {#pragma acc loop independentfor (int i = 0; i < nx; i++) {const int ix = nx*ny*k + nx*j + i;const int ip = i == nx - 1 ? ix : ix + 1;
const int im = i == 0
? ix : ix - 1;
const int jp = j == ny - 1 ? ix : ix + nx;
const int jm = j == 0
? ix : ix - nx;
const int kp = k == nz - 1 ? ix : ix + nx*ny;
const int km = k == 0
? ix : ix - nx*ny;
fn[ix] = cc*f[ix]+ ce*f[ip] + cw*f[im]+ cn*f[jp] + cs*f[jm]+ ct*f[kp] + cb*f[km];}}}return (double)(nx*ny*nz)*13.0;}

なお、 present にしなくても期待通りに動作します。133

OpenACC化(4): データ転送の最適化(2)• diffusion3d関数で present とし、main関数で data を追加#pragma acc data copy(f[0:n]) create(fn[0:n]){
start_timer();main.c, main 関数内

for (; icnt<nt && time + 0.5*dt < 0.1; icnt++) {if (icnt % 100 == 0)
fprintf(stdout, "time(%4d) = %7.5f¥n", icnt, time);
flop += diffusion3d(nx, ny, nz, dx, dy, dz, dt, kappa, f, fn);swap(&f, &fn);time += dt;}
elapsed_time = get_elapsed_time();}copy/create など適切なものを選びます。
make して実行してみましょう。 どのくらいの実行性能が出ましたか？
OpenACC化の例は、 openacc_diffusion/02_openacc
134

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

PGIコンパイラを利用する場合、OpenACCプログラムがどのように実行されているか、環境変数PGI_ACC_TIMEを設定すると簡単に確認することができる。
Linuxなどでは、環境変数PGI_ACC_TIME を1に設定し、プログラムを実行する。
$ export PGI_ACC_TIME=1$ ./run•Reedbush でジョブに環境変数PGI_ACC_TIME を設定する場合は、ジョブスクリプト中に記載する。
$ cat run.sh...
. /etc/profile.d/modules.shmodule load pgi/17.1
export PGI_ACC_TIME=1./run

サンプルコードは、 openacc_diffusion/03_openacc _pgi_acc_time
135

PGI_ACC_TIME によるOpenACC 実行の確認• ジョブ実行が終わると、標準エラー出力にメッセージが出力される。
$ cat run.sh.e??????
Accelerator Kernel Timing data
/lustre/pz0115/z30115/lecture/lecture_samples/openacc_diffusion/03_open
acc_pgi_acc_time/main.c
main NVIDIA devicenum=0time(us): 6,359
38: data region reached 2 times
データ移動の回数
38: data copyin transfers: 1
device time(us): total=3,327 max=3,327 min=3,327 avg=3,327
55: data copyout transfers: 1
device time(us): total=3,032 max=3,032 min=3,032 avg=3,032/lustre/pz0115/z30115/lecture/lecture_samples/openacc_diffusion/03_open
acc_pgi_acc_time/diffusion.c
diffusion3d NVIDIA devicenum=0time(us): 101,731
19: compute region reached 1638 times
25: kernel launched 1638 times
起動したスレッド
grid: [4x128x32] block: [32x4]
カーネル
device time(us): total=101,731 max=64 min=62 avg=62
elapsed time(us): total=136,255 max=540 min=81 avg=83
実行時間
19: data region reached 3276 times
136

参考：プロファイラ（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コマンドを実行

137

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

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

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

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

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

ただし1GB以上あるので、今やるとおそらくネットワークがパンクします…
138①Fileをクリック⓪ nvvp コマンドにより起動後の画面

139

②Importをクリック③ Nvprof を選択④ Next140

⑤ Multiple processesに変更⑦ Browse⑥ Next141

⑧ ファイル選択画面が開くので、nvvpで作成したプロファ
イリングデータを選択し、OK⑨ 選択できたら Finish
142ここで拡大この辺がメインの計算部分naïve な diffusion 計算のプロファイルを取った結果のタイムライン表示143↓データ転送に時間がかかっている←この細いところが計算計算部分をクリックすると、カーネルの情報が得られる

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

GPUを用いた
FDTD法による電磁波伝搬計算
145電磁波の方程式• 真空での電場Eと磁場Hの時間発展
Maxwell 方程式の一部

1
⇥E
=
⇥t
(

: 誘電率)

H

H
=
t

1
µ

E

( µ : 透磁率)

この方程を、2次元FDTD法（Finite-difference time-domain 法）*を用いて解いて行きます。
* K. S. Yee, "Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic
Media,'' IEEE Trans. on Antennas and Propagat., vol. 14, pp. 302-307, May 1966.

146

FDTD法（１）• EとHの時間発展

En
H

En 1
1
=
t

n+ 12

H

n

1
2

t

H
=

1
2

n

1
µ

En

変形して、

E =E
n

H

n+ 21

n 1

+

=H

n

t
1
2

H
t
µ

n

1
2

En
147

FDTD法（２）• 例えば、148

2次元FDTD法の変数配置
Ey14
mgn ( = 4 )

13
12
11

プログラムでは、配列は1次元的に確保します。

9

whole_global.length[0]= inside_global[0] +2*mgn + 18であることに注意

10
inside̲global.length[1] ( = 7 )

7
6
5
4
3

mgn ( = 4 )

2
1
0

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
whole̲global.length[1]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
= inside̲global.length[1] 0+ 2mgn + 1 ( = 16)
mgn ( = 4 )inside̲global.length[0] ( = 7 )mgn ( = 4 )HzExwhole̲global.length[0]= inside̲global.length[0]+ 2mgn + 1 ( = 16)149

ソースコード（１）• サンプルコード: openacc_fdtd/ü OpenACCを利用したFDTD法（電磁波解析）openacc_fdtd/01_original

MPI並列化されたCPUコード。

openacc_fdtd/02_openacc1

calc_ex_ey, pml_boundary_ex, pml_boundary_ey,
がOpenACC。

openacc_fdtd/03_openacc2

時間更新ループ全体が OpenACC。

openacc_fdtd/04_openacc3

初期化を含め OpenACC。

openacc_fdtd/05_openacc4

データ移動の最適化。

※本プログラムはMPI化されていますが、本講習会では扱いません

150ソースコード（２）• それぞれのファイルの内容
main.cプログラムのメインコード

fdtd2d.{c, h}

2次元 FDTD の 計算コード

fdtd2d_sources.{c, h}

入射光設定のための関数

setup.c

計算条件の設定と変数の初期化

config.{c, h}

物理定数の定義

output.{cc, h}

計算結果出力のための関数

bitmap*

BMPファイル作成のための関数

本講習では、“main.c”、“fdtd2d.c”、” fdtd2d_sources.c”、”setup.c” のソー
スコードを追記・修正していきます。

151

計算条件
• 2次元波動伝搬

電磁波の境界での非物理的な反射を防ぐための吸収境界条件 (PML)ü 成分: Ex、Ey、Hzü y 方向下側から平面波を入射
mgn

dx = lx/nx
dy = ly/ny
デフォルト設定:
nx = 512ny = 512
mgn = 8lnx = 529
lny = 529ny0

真空と異なる
媒質を設定

プログラム中では下記の変数が使われているので注意

物体を設置

入射光

mgn

inside_global.length[0] = nx0
inside_global.length[1] = ny0
whole_global.length[0]= nx0 + 2*mgn + 1
whole_global.length[1]= ny0 + 2*mgn + 1

mgn

nx0

mgn

152

コード全体の流れ（main.c 内）
計算領域の設定
Range 構造体
inside_global, whole_global,
inside, whole
CPUメモリに配列を確保
ex, ey, hz, ...

媒質と物体の設定（ユーザ定義）
set_object_er 関数 (main.cで定義)初期化と計算条件の設定
init_vars, set_initial_condition 関数
init_pml_vars 関数
set_pml_initial_condition 関数
set_pml_rer 関数
計算結果出力MPI_Gather収集, write_bmp関数

電場 E の時間発展
calc_ex_ey 関数
電場 E の境界条件
pml_boundary_ex 関数
pml_boundary_ey 関数
電場の入射（入射光）
plane_wave_incidence 関数
磁場 Ex の袖領域更新
MPI_Send/MPI_Recv

FDTD計算繰り返し時間発展

媒質と物体の初期化
init_relative_permittivity関数
init_object関数

磁場 H の袖領域更新
MPI_Send/MPI_Recv

磁場 H の時間発展
calc_hz 関数
磁場 H の境界条件
pml_boundary_hz 関数
計算結果出力MPI_Gather,write_bmp関数

153

計算領域の設定（１）• Range 構造体
ü 計算領域の始点と大きさを保持// config.hstruct Range {int length[2];int begin [2];};// main.c
全領域の中心領域
const struct Range inside_global = { { atoi(argv[1]), atoi(argv[2]) },{ 0, 0 } };const struct Range whole_global = { { inside_global.length[0] + 2*mgn + 1, 全領域のinside_global.length[1] + 2*mgn + 1},{ inside_global.begin[0] - mgn, 全体領域
inside_global.begin[1] - mgn} };const struct Range inside

const struct Range whole

分割領域の= { { inside_global.length[0],inside_global.length[1]/nsubdomains }, 中心領域
{ 0,
inside_global.length[1]/nsubdomains * rank } };= { { inside.length[0] + 2*mgn + 1,
inside.length[1] + 2*mgn + 1},
分割領域の{ inside.begin[0] - mgn,
全体領域
inside.begin[1] - mgn
} };

154

• Range 構造体

mgn

ny

inside.length[0] = nx
inside.length[1] = ny

mgn

mgn

プログラムでは
下記の通り= { { inside_global.length[0],inside_global.length[1]/nsubdomains },{ 0,
inside_global.length[1]/nsubdomains * rank } };= { { inside.length[0] + 2*mgn + 1,
inside.length[1] + 2*mgn + 1}, 分割領域の{ inside.begin[0] - mgn,
全体領域
inside.begin[1] - mgn
} };
mgn

const struct Range whole

whole.length[0]nx0= nx + 2*mgn + 1
whole.length[1] = ny + 2*mgn + 1
mgn

分割領域の
中心領域

mgn

座標(inside.begin[0],inside.begin[1])座標(whole.begin[0],whole.begin[1])

ny

ny0

const struct Range insidemgn

struct Range {int length[2];int begin [2];
};

mgn

mgn

ü 計算領域の始点と大きさを保持

ny

mgn

計算領域の設定（２）

mgn

nx

mgn

155

配列の確保• 物理変数配列は main.c で確保// main.cconst int
nelems= whole.length[0] * whole.length[1];const int
nelems_x= whole.length[0];const int
nelems_y= whole.length[1];const size_t size
= sizeof(FLOAT)*nelems;
const size_t size_x= sizeof(FLOAT)*nelems_x;
const size_t size_y= sizeof(FLOAT)*nelems_y;
const size_t size_global = sizeof(FLOAT)* whole_global.length[0] * whole_global.length[1];
FLOAT *ex
= (FLOAT *)malloc(size); // 電場 Ex
FLOAT *ey
= (FLOAT *)malloc(size); // 電場 Ey
FLOAT *hz
= (FLOAT *)malloc(size); // 磁場 Hz...
// For output
FLOAT *ex_global = (FLOAT *)malloc(size_global);
FLOAT *ey_global = (FLOAT *)malloc(size_global);
FLOAT *hz_global = (FLOAT *)malloc(size_global);– 多くの配列は whole.length[0] * whole.length[1]– ex_global, ey_global, hz_global はファイル出力に使うため、
whole_global.length[0] * whole_global.length[1]156時間発展（１）• 前半
ü 電場Eの時間発展（calc_ex_ey）、境界条件(pml_boundary_...)
ü 入射光（plane_wave_incidence）
while (icnt < nt) {MPI_Status status;
const int tag = 0;
const int nhalo= whole.length[0];const int inside_end1 = inside.begin[1] + inside.length[1];const int src_hz
const int dst_hz= whole.length[0] * (inside_end1- whole.begin[1] - 1);= whole.length[0] * (inside.begin[1] - whole.begin[1] - 1);

MPI_Send(&hz[src_hz], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD);MPI_Recv(&hz[dst_hz], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD, &status);
calc_ex_ey(&whole, &inside, hz, cexly, ceylx, ex, ey);pml_boundary_ex(&whole, &inside, hz, cexy, cexyl, rer_ex, ex, exy);
pml_boundary_ey(&whole, &inside, hz, ceyx, ceyxl, rer_ey, ey, eyx);
const int j_in = 0;
plane_wave_incidence(&whole, &inside, time, j_in, wavelength, ex, ey);time += 0.5*dt;
(後半へ）
157

時間発展（２）• 後半
ü 磁場Hの時間発展（calc_hz）、境界条件(pml_boundary_hz)（前半から）
const int src_ex
const int dst_ex= whole.length[0] * (inside.begin[1] - whole.begin[1]);= whole.length[0] * (inside_end1- whole.begin[1]);MPI_Send(&ex[src_ex], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD);MPI_Recv(&ex[dst_ex], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD, &status);
calc_hz(&whole, &inside, ey, ex, chzlx, chzly, hz);
pml_boundary_hz(&whole, &inside, ey, ex, chzx, chzxl, chzy, chzyl, hz, hzx, hzy);time += 0.5*dt;
icnt++;(出力など)}

158

GPUを用いた
FDTD法による電磁波伝搬計算の実習
159

プログラムのコンパイルと実行（１）• CPUコードのコンパイルと実行$ cd openacc_mpi_fdtd/01_original$ make$ qsub ./run.sh$ cat run.sh.o??????
Rank 0: hostname = a090
Rank 1: hostname = a090
Rank 2: hostname = a091
Rank 3: hostname = a091
Calculation condition
nx_global=
512

openacc_fdtd/01_original? の数字はジョブごとに変わります。
利用したノード（省略）
icnt = 4900, time = 2.3115e-14 [sec]icnt = 5000, time = 2.3587e-14 [sec]-----------------------------Domain= 512 x 512
nsubdomains = 4
output_file = 1
Time=
4.103535 [sec]------------------------------

計算領域サイズ、領域
分割数、出力の有無、
計算時間

なお、qsub ./run_no_out.sh すると出力なしで実行する。性能測定用。
160

プログラムのコンパイルと実行（２）• プログラムの実行時オプション$ cat run.sh#! /bin/sh#PBS -q h-tutorial#PBS -l select=1:mpiprocs=1:ompthreads=0

openacc_fdtd/01_original（省略）
mkdir -p sim_run
cd sim_run
nprocs=1
mpirun -np $nprocs ../run 512 512 $nprocs 5000 50
mpirun -np <nprocs> ../run <nx> <ny> <nprocs> <nt> <nout>

nprocs: 全ランク数（=分割数） ※今回は１
nx, ny: 計算領域サイズ
nt: 全時間ステップ
nout: 出力を行うタイムステップ数。50 の場合、50ステップに1回出力する。0 を指定すると出力しない。

161

計算結果の表示• 計算結果は sim_run に BMP として出力される$ cd sim_run/
openacc_mpi_fdtd/01_original• 計算結果の表示ü 1枚のBMPを見る$ display e05000.bmp

ü 複数のBMPファイルをアニメーションで表示$ animate *.bmp

なおssh -Y txxxxx@reedbush.cc.u-tokyo.ac.jpと –Y をつけていないと表示されない。うまく表示できない場合は画像を手元にコピーして表示してください。

162

計算結果の例• 出力されたBMPファイルの一例

ny0

mgn

ü Ex （電場の x 成分）の出力

y

mgn

入射光

mgn

z

x

nx0

mgn

163

実習1
– calc_ex_ey, pml_boundary_ex, pml_boundary_ey を OpenACC化しましょう。
– Makefileü コンパイルオプションの修正– main.cü OpenACCヘッダーの追加
ü data 指示文の追加

実行速度が遅くても、動くプログラムである
状態を保ちながらOpenACC化します。
末端の関数から
OpenACC化するのがよいでしょう。

– fdtd2d.cü kernels 指示文、loop 指示文の追加

解答例は、 openacc_fdtd/02_openacc1
164

kernels, loop指示文• fdtd2d.c 内の関数
void calc_ex_ey(const struct Range *whole, const struct Range *inside,
const FLOAT *hz, const FLOAT *cexly, const FLOAT *ceylx, FLOAT *ex, FLOAT *ey){
const int nx= inside->length[0];const int ny= inside->length[1];const int mgn[] = { inside->begin[0] - whole->begin[0],inside->begin[1] - whole->begin[1] };const int lnx= whole->length[0];#pragma acc kernels present(hz, cexly, ex)#pragma acc loop independentfor (int j=0; j<ny+1; j++) {#pragma acc loop independentfor (int i=0; i<nx; i++) {const int ix = (j+mgn[1])*lnx + i+mgn[0];const int jm = ix - lnx;//ex[ix] += cexly[ix]*(hz[ix]-hz[jm]) - cexlz[ix]*(hy[ix]-hy[km]);ex[ix] += cexly[ix]*(hz[ix]-hz[jm]);
}
}
(省略）
}

165

実習2
– main 関数内の while 内をすべて OpenACCにしましょう。
– main.cü data 指示文の移動と copyin などの最適化– fdtd2d.cü 残りの関数にkernels 指示文、loop 指示文の追加– fdtd2d_sources.cü kernels 指示文、loop 指示文の追加

解答例は、 openacc_fdtd/03_openacc2
166data 指示文
• main関数のwhile 外に data を移動#pragma acc data ¥copyin(ex[0:nelems], ey[0:nelems], hz[0:nelems])¥
copyin(cexly[0:nelems], ceylx[0:nelems], chzlx[0:nelems], chzly[0:nelems]) ¥copyin(exy[0:nelems], eyx[0:nelems], hzx[0:nelems], hzy[0:nelems]) ¥copyin(cexy[0:nelems_y], ceyx[0:nelems_x], chzx[0:nelems_x], chzy[0:nelems_y]) ¥copyin(cexyl[0:nelems_y], ceyxl[0:nelems_x], chzxl[0:nelems_x], chzyl[0:nelems_y]) ¥copyin(obj[0:nelems], er[0:nelems]) ¥copyin(rer_ex[0:nelems], rer_ey[0:nelems]){
while (icnt < nt) {MPI_Status status;
const int tag = 0;
const int nhalo= whole.length[0];const int inside_end1 = inside.begin[1] + inside.length[1];const int src_hz
const int dst_hz= whole.length[0] * (inside_end1- whole.begin[1] - 1);= whole.length[0] * (inside.begin[1] - whole.begin[1] - 1);#pragma acc host_data use_device(hz){MPI_Send(&hz[src_hz], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD);MPI_Recv(&hz[dst_hz], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD, &status);
}

167

実習3
– 初期化を含めて全てOpenACCにします。ただし、set_object_er がCPU上のユーザ定義関数のため、これ以降の初期化関数をOpenACCにします。
– main.cü data 指示文の移動と最適化（多くが create になるはずです）– setup.cü kernels 指示文、loop 指示文の追加

解答例は、 openacc_fdtd/04_openacc3
168

実習4
– 計算領域のサイズなどを変更して性能測定してみましょう。
– OpenACCコードをさらに最適化しましょう。
ü PGI_ACC_TIMEも活用しましょう。
ü 実は単純に fdtd2d.c に kernels と loop を入れても、いくつかの関数で暗黙の copyin が発生します。これも修正していきましょう。
$ makecalc_ex_ey:
25, Generating present(ex[:],cexly[:])Generating implicit copyin(mgn[:])Generating present(hz[:])27, Loop is parallelizable
29, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
27, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */29, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */37, Generating present(ey[:],ceylx[:])Generating implicit copyin(mgn[:])

解答例は、 openacc_fdtd/05_openacc4

169

Q&A
• アカウントは1ヶ月間有効です。
• 資料のPDF版はWEBページにございます。
– https://www.cc.u-tokyo.ac.jp/events/lectures/97/lec097.php• アンケートへの協力をお願いします。

170