第107回 お試しアカウント付き
並列プログラミング講習会

OpenACCとMPIによるマルチGPUプログラミング入門
東京大学 情報基盤センター
担当：星野哲也, 下川辺 隆史
hoshino@cc.u-tokyo.ac.jp
shimokawabe @ cc.u-tokyo.ac.jp（内容に関するご質問はこちらまで）

1

講習会スケジュール
n 開催日時ü 11月6日（火） 10:00 – 18:00

n プログラムü 09:30 - 10:00 受付ü 10:00 – 10:40 Reedbush-Hへのログインü 10:40 – 12:00 GPUとOpenACC基礎（座学）
ü （昼休み）ü 13:00 – 14:00 OpenACC基礎（演習）ü 14:00 – 14:50 MPI復習（座学、簡単な演習）（休憩）ü 15:00 – 16:30 OpenACCとMPIによるマルチGPUプログラミ
ング（座学、演習）ü 16:40 – 18:00 複数GPUを用いたFDTD法による電磁波伝搬
計算（座学、演習）ü2

講習会についてn 本講習会では
ü GPUとOpenACCプログラミングの基礎
ü OpenACCとMPIを使った複数GPUプログラミング入門を中心に扱います。
n その他の講習会

https://www.cc.u-tokyo.ac.jp/events/lectures/
ü OpenACC、MPIの機能全般について学びたい場合は、他の講習会も受講することをお勧めします。
n スパコンイベント情報メール配信サービス

https://regist.cc.u-tokyo.ac.jp/announce/
ü 講習会や研究会の案内、トライアルユースの実施のお知らせなどを配信しています。
3

事前準備

4

Reedbush にログインする
n 別資料「Reedbush利用の手引き」を参照し、ロ

グインします。

5

Reedbush 利用上の注意（１）n ディレクトリについて（home と lustre)ü ログイン時のディレクトリ（/home/gt00/txxxxx）にはログイン時に必要なファイルのみを置く
ü プログラム作成や実行などに必要なファイルは /lustre 以下のディレクトリ（/lustre/gt00/txxxxx）に置く
ü /home は計算ノードからは参照できないü cdw コマンドで Lustreファイルシステムへ移動できる。
$ cdw

6

Reedbush 利用上の注意（２）n コンパイルおよび実行のための環境準備
ü コンパイルおよび実行のための環境を準備するために module コ

マンドを使用する。これによって様々な環境を簡単に切り替えて使用できる。
$ module load <module_name>

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

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

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

7

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

n CUDA開発環境を使う場合$ module load cudan Intelコンパイラを使う場合$ module load intel

n MPIを使う場合

訂正$ module load openmpi/gdr/2.1.2/{gnu,intel,pgi}$ module load mvapich2/gdr/2.3a/{gnu,pgi}
ü PGIなどのコンパイラに追加して loadn モジュールはジョブ実行時にもコンパイル時と同じもの

を load する。
n 組み合わせて利用できる

8

サンプルコードのコンパイルnReedbush へのログイン$ ssh -Y txxxxx@reedbush.cc.u-tokyo.ac.jp
ü txxxxx 各自の利用者番号（アカウント）に置き換えてください。
ü -Y をつけてください。

n

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

n

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

n

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

n

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

n

モジュールをロードする。
$ module load pginコンパイルする。
$ cd openacc_hello/01_hello_acc$ maken

実行ファイルがてきていることを確認する。
$ ls9プログラムの実行
n ジョブとして投入し、実行する。
$ qsub ./run.sh

n 投入されたジョブを確認する。
$ rbstatn 実行が終了すると、以下のファイルが生成される。
run.sh.o??????
run.sh.e??????

（?????? は数字）n 上記の標準出力ファイルの中身を確認する。
$ cat run.sh.o??????

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

10

東大情報基盤センター
スーパーコンピュータの概略
11

東大センターのスパコン
FY
11

2基の大型システム，6年サイクル(だった)12 13 14 15 16 17 18 19 20 21 22 23 24 25
Yayoi: Hitachi SR16000/M1
IBM Power-7

メニーコア型大規模
スーパーコンピュータ（JCAHPC：筑波大・東大）

54.9 TFLOPS, 11.2 TB

Oakforest-PACS
Fujitsu, Intel KNL

T2K Tokyo

140TF, 31.3TB

Big Data &Extreme Computing

25PFLOPS, 919.3TB

Oakleaf-FX: Fujitsu PRIMEHPC
FX10, SPARC64 IXfx

BDEC System

50+ PFLOPS (?)

1.13 PFLOPS, 150 TB

Oakbridge-FX

Oakbridge-II

136.2 TFLOPS, 18.4 TB
データ解析・シミュレーション
融合スーパーコンピュータ
長時間ジョブ実行用演算加速装置付き並列スーパーコンピュータ

Intel/AMD/P9/ARM CPU only

5-10 PFLOPS

Reedbush, HPE
Broadwell + Pascal

大規模超並列スーパーコンピュータ

1.93 PFLOPS

Reedbush-L
HPE
1.43 PFLOPS

12

2システム運用中
n Reedbush (SGI, Intel BDW + NVIDIA P100 (Pascal))ü データ解析・シミュレーション融合スーパーコンピュータü 3.361 PF, 2016年7月〜 2020年6月
ü 東大ITC初のGPUシステム (2017年3月より), DDN IME (Burst Buffer)
ü Reedbush-U: CPU only, 420 nodes, 508 TF (2016年7月)
ü Reedbush-H: 120 nodes, 2 GPUs/node: 1.42 PF (2017年3月)
ü Reedbush-L: 64 nodes, 4 GPUs/node: 1.43 PF (2017年10月)

n Oakforest-PACS (OFP) (富士通、Intel Xeon Phi (KNL))ü JCAHPC (筑波大CCS＆東大ITC)
ü 25 PF, TOP 500で6位 (2016年11月) (日本で1位)
ü 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

東京大学情報基盤センター スパコン（2/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ノードのブロック図
n メモリのうち、「近い」メモリと「遠い」メモ

リがある=> 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の通信網
n フルバイセクションバンド幅を持つFat Tree網ü どのように計算ノードを選んでも互いに無衝突で通信が可能
n 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

Integration of Memory on package: innovative memoryn Intel Xeon
Phi (Knights Landing)architecture for high bandwidth and high capacity
Integration of Fabric on package

n 1ノード1ソケット

Three products

HotChips27
Knights Landing
Overview
KNLスライドより

KNL Self-Boot KNL Self-Boot w/ Fabric
KNL Cardn MCDRAM:
オンパッケージ(Baseline)(Fabric Integrated)(PCIe-Card)の高バンド幅メモリ16GB+ DDR4メモリ

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.

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

Core

1MB
L2

C
H
A
N
N
E
L
S

EDC

EDC

3

Tile

EDC

Core

MCDRAMDDR
4

36 Tilesconnected by
2D MeshInterconnect

DDR MC

2 VPU

4

EDC

3
D
D
R
4

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

EDC

EDC

MCDRAM

Omni-path not shown

misc

Package

DDR MCEDC

MCDRAM

C
H
A
N
N
E
L
S

EDC

MCDRAM

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





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

東大情報基盤センターOakforest-PACSスーパーコン
ピュータシステムの料金表（2017年4月1日 ）n パーソナルコース（年間）n コース１：まで
n コース２：
ドまで

100,000円：

8ノード(基準)、最大16ノード

200,000円：

16ノード(基準)、最大64ノーn グループコースn400,000円 (企業 480,000円)
まで： 1口 8ノード（基準）、最大128ノード

n 以上は、「トークン制」で運営

360日 24時間の「トークン」が与えられる
n 基準ノードまでは、トークン消費係数が1.0n 基準ノードを超えると、超えた分は、消費係数が2.0になる
n 大学等のユーザはFX10、Reedbushとの相互トークン移行も可能
n 申し込みノード数

22

東大情報基盤センターReedbushスーパーコンピュータシ
ステムの料金表（2017年4月1日 ）n パーソナルコース（年間）n

150,000円： RB-U: 4ノード（基準）、最大16ノードまで
RB-H: 1ノード（基準）、最大2ノードまで

n グループコース
n 300,000円： 1口 4ノード（基準）、最大128ノードまで、

RB-H： 1ノード（基準）、最大32ノードまで（トークン係数はUの2.5倍）

n RB-Uのみ 企業 360,000円： 1口 4ノード（基準）、最大128

ノードまで
n RB-Hのみ 企業 216,000円 ： 1口 1ノード（基準）、最大32ノードまで
n 以上は、「トークン制」で運営
360日 24時間の「トークン」が与えられる
n 基準ノードまでは、トークン消費係数が1.0n 基準ノードを超えると、超えた分は、消費係数が2.0になる
n 大学等のユーザはFX10, Oakforest-PACSとの相互トークン移行も可能
n ノード固定もあり
n 申し込みノード数

23

東大情報基盤センターFX10スーパーコンピュータシステ
ムの料金表（2017年4月1日）
n パーソナルコース（年間）n コース１：90,000円：

12ノード(基準)、最大24ノードまで

n コース２：

180,000円：

24ノード(基準)、最大96ノードまで

n グループコース
n 360,000円 (企業 432,000円)：1口、12ノード、最大1440ノードまで

n 以上は、「トークン制」で運営

360日 24時間の「トークン」が与えられる
n 基準ノードまでは、トークン消費係数が1.0n 基準ノードを超えると、超えた分は、消費係数が2.0になる
n 大学等のユーザはReedbush, Oakforest-PACSとの相互トークン移行も可能
n 申し込みノード数

24

25

JPY (=Watt)/GFLOPS Rate
Smaller is better (efficient)
SystemOakleaf/Oakbridge-FX (Fujitsu)(Fujitsu PRIMEHPC FX10)
Reedbush-U (SGI)(Intel BDW)
Reedbush-H (SGI)(Intel BDW+NVIDIA P100)
Oakforest-PACS (Fujitsu)(Intel Xeon Phi/Knights Landing)

JPY/GFLOPS
125
62.0
17.1
16.5

トライアルユース制度についてn

安価に当センターのOakleaf/Oakbridge-FX, Reedbush-U/H, OakforestPACSシステムが使える「無償トライアルユース」および「有償トライ
アルユース」制度があります。
n アカデミック利用
n パーソナルコース、グループコースの双方（1ヶ月〜3ヶ月）

n

企業利用
n パーソナルコース（１ヶ月〜３ヶ月）（FX10: 最大24ノード、最大96

ノード,
RB-U: 最大16ノード, RB-H: 最大2ノード、OFP: 最大16ノード, 最大64ノード）講習会いずれかの受講が必須、審査無n グループコース
n 無償トライアルユース：（１ヶ月〜３ヶ月）：無料（FX10:最大
1,440ノード、 RB-U: 最大128ノード、RB-H: 最大32ノード、OFP:
最大2048ノード）
n 有償トライアルユース：（１ヶ月〜最大通算９ヶ月）、有償（計算資源は無償と同等）
n スーパーコンピュータ利用資格者審査委員会の審査が必要（年2回実施）
n 双方のコースともに、簡易な利用報告書の提出が必要
26

スーパーコンピュータシステムの詳細
n 以下のページをご参照ください
n 利用申請方法
n 運営体系
n 料金体系
n 利用の手引などがご覧になれます。
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/
27

GPU入門28

What’s GPU ?
n Graphics Processing Unitn もともと PC の3D描画専用の装置
n パソコンの部品として量産されてる。= 非常に安価
3D GameComputer Graphics

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

29

GPUコンピューティング
n
n
n

n
n

GPUはグラフィックスやゲームの画像計算のために進化を続けている。
CPUがコア数が2-12個程度に対し、GPUは1000以上のコアがある。
GPUを一般のアプリケーションの高速化に利用することを「GPUコンピューティング」「GPGPU (General Purpose computation on
GPU)」などという。
2007年にNVIDIA社のCUDA言語がリリースされて大きく発展
ここ数年、ディープラーニング（深層学習）、機械学習、AI（人工知能）などでも注目を浴びている。

30

GPUの特徴n

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

ラミングする必要がある。

n

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

31

NVIDIA Tesla P100

n 56 SMs, 3584 CUDAコア, 16 GByte
Tesla P100 whitepaper より
32

並列計算
n 実行時間 T の逐次処理のプログラムを p 台の計算機で並列計算することで、実行時間を T / p にする。
Tp台
T/pn 実際にできるかどうかは、処理内容（アルゴリズム）による。アルゴリズムによって難易度は異なる。
ü 並列化できないアルゴリズム、通信のオーバーヘッド
ü 部分的にでも並列化できないアルゴリズムがあると、どれだけ並列数を上げても、その時間は短縮されない。
n 並列処理（計算）の種類
ü 「タスク並列」と「データ並列」

33

タスク並列n タスク（仕事）を分割することで並列化する。
n タスク並列の例：カレーを作る
ü 仕事１：野菜を切る
ü 仕事２：肉を切る
ü 仕事３：水を沸騰させる
ü 仕事４：野菜と肉を入れて煮込む
ü 仕事５：カレーのルウを入れる

n 並列化

GPUは苦手
仕事1
仕事2
仕事3

仕事4

仕事5

時間
34

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

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

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の並列計算はこれが原則。
35

GPUプログラミングでの注意点（１）これらはプログラミング言語が CUDA かOpenACCに関わらず、GPUプログラミングでは考慮する必要がある。
n CPUメモリとGPUメモリが独立
ü 特にOpenACCでは、これを意識する必要がある

n 高い性能を得るには、スレッド数 >> コア数（P100 では3584）
ü 数万から数百万スレッド程度
ü メモリアクセスによる暇な時間（ストール）を、他のスレッドが埋めることができる
ü Intel CPU では Hyperthread に相当。ただしこちらは物理コア x 2

36

GPUプログラミングでの注意点（２）n スレッド間での異なる処理（分岐）は避ける
ü 連続した32スレッド（Warp と呼ぶ）が同じ命令を実行ü Warp内で分岐先が異なると、全ての分岐先を実行し、真になる条件の処理のみ採用
ü branch divergence, divergent branch などという。
...
...
if (奇数番スレッド） {

スレッド

処理 A;
} else {
処理 B;
}
...
...
37

GPUプログラミングでの注意点（３）n 同じWarp内のスレッド（連続するスレッド）は近いメモ

リアドレスへアクセスすると効率的ü コアレスドアクセス（coalesced access）と呼ぶ
ü メモリアクセスは128 Byte 単位で行われる。128 Byte に収まれ

ば1回のアクセス、超えれば128 Byte アクセスをその分繰り返す。
128 byte x 1回のメモリアクセス
Address 128

160

192

224

256

...

15 16

...

31

192

224

256

Thread0

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

160

38

Thread
0

...

15 16

...

31

OPENACC入門

39

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

簡単

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

n 指示文ベース（OpenACC）
ü 指示文（ディレクティブ）を挿入するだけである程度高速化。
ü 既存のソースコードを活用できる。
n プログラミング言語（CUDA、OpenCLなど）
ü GPUの性能を最大限に活用。
ü プログラミングにはGPGPU用言語を使用する必要あ

り。

難しい

40

OpenACCn OpenACCとはü 新しいアクセラレータ用プログラミングインターフェースü OpenMP のようなディレクティブ（指示文）・ベース
ü C 言語/C++, Fortran に対応
ü 2011年秋に OpenACC1.0、最新は 2.5
ü コンパイラ：［商用］PGI, Cray,［フリー］GCC（PGIは無料版あり）
ü WEBサイト：http://www.openacc.org/
n 指示文ベースの利点ü 指示文：コンパイラへのヒントü アプリケーションの開発や移植が比較的簡単
ü ホスト（CPU）用コード、複数のアクセラレータ用コードを単一コードとして記述。メンテナンスが容易。高生産性。
FortranC言語#pragma acc directive-name [clause, …] !$acc directive-name [clause, …]{! Fortran code// C code!$acc end directive-name41}

逐次、OpenMP、OpenACCの実行イメージ
1スレッド

OpenMP

OpenACC

int main(){...
#pragma acc ...
for (i=0; i<n; i++) {
...
}
...
}

CPU

CPU

CPU

デバイス
(GPU)
42

はじめての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;
}

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

解放
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から
copyout

43

はじめての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;
}

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

解放
double sum = 0;
for (int i=0; i<n; i++) {sum += b[i];}fprintf(stdout,"%f¥n", sum/n);
コード上同じa, b であっても、原則としてホストコードはホストメモリで確保された a, b 、GPUで実行さfree(a); free(b);れる並列領域（カーネル）はデバイスメモリで確保された
a, b
return 0;を参照しにいく。
}

GPUから
copyout

44

OpenACCの主な指示文
n 並列領域指定指示文
ü kernels, paralleln データ管理・移動指示文
ü data, enter data, exit data, updaten 並列処理の指定
ü loopn その他
ü host_data, atomic, routine, declare

赤字：この講習会で扱うもの

45

OpenACC によるアクセラレータでの実行
n kernels 指示文
ü 指定された領域がアクセラレータで実行されるカーネルへ
ü 一般には、それぞれのループが別々のカーネルへコンパイル
int main() {#pragma acc kernels{for (int i=0; i<n; i++) {A;}
for (int i=0; i<n; i++) {B;}

kernel 1
kernel 2}}

ü 同様な指示文として、領域内が一つのカーネルとして生成される parallel 指示文もある
46

CPUコードのOpenACC化
openacc_hello/01_hello_acc
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;}

n ループのOpenACC 化
ü並列化したいループにkernels,
loop 指示文を追加

for (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;}47

CPUコードのOpenACC化
openacc_hello/01_hello_acc
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;}

n ループのOpenACC 化
ü並列化したいループにkernels,
loop 指示文を追加#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;}48

CPUコードのOpenACC化
openacc_hello/01_hello_acc
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;}#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;}

n ループのOpenACC 化
ü並列化したいループにkernels,
loop 指示文を追加

カーネルとしてコンパイルされ、GPU上で実行される配列の1要素が1スレッドで処理されるイメージ

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

49

OpenACCのデフォルトの変数の扱い
n スカラ変数
ü firstprivate または private
ü ホストからデバイスへコピーが渡され初期化。ホストに戻せない。
n 配列ü sharedü デバイスメモリに動的に確保され、スレッド間で共有。
ü デバイスからホストへコピーすることが可能。

n kernels 構文に差し掛かると、
ü OpenACCコンパイラは実行に必要なデータを自動で転送する。
ü 配列はデバイスメモリに確保され、shared になる。
ü 構文に差し掛かるたびに転送を行う。data 指示文で制御できる。

50

データ管理・移動
n data 指示文
ü デバイス(GPU)メモリの確保と解放、ホスト(CPU)とデバイス(GPU)間のデータ転送を制御
kernels指示文では、データ転送は自動的に行われる。data指示文でこれを制御することで、不要な転送を避け、性能向上できるü CUDA で言うところの cudaMalloc, cudaMemcpy に相当
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;}

openacc_hello/01_hello_acc#pragma acc data copyin(a[0:n]), copyout(b[0:n])#pragma acc kernels#pragma acc loop independent
変数 c はスカラ変数のため、
for (int i=0; i<n; i++) {自動的に デバイスへコピーさ
b[i] = a[i] + c;}
れ、プライベート変数となる。
51

data 指示文の指示節
n
n

n

n

n

n

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
ü 何もしない。既にデバイス上で確保済みであることを伝える。
copy/copyin/copyout/create は既にデバイス上確保されているデータに対しては何もしない。present として振る舞う。（OpenACC2.5以降）
52

データの移動範囲の指定
n ホストとデバイス間でコピーする範囲を指定
n 部分配列の転送が可能
n Fortran と C言語で指定方法が異なるので注意

n 二次元配列A転送する例
n Fortran: 下限と上限を指定!$acc data copy(A(lower1:upper1, lower2:upper2) )...
!$acc end datan C言語: 始点とサイズを指定#pragma acc data copy(A[begin1:length1][begin2:length2])...

53

並列処理の指定
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 independent
loop指示文for (int i=0; i<n; i++) {b[i] = a[i] + c;
}

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

解放
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から
copyout54

並列処理の指定：loop 指示文
n loop 指示文
ü ループマッピングのパラメータの調整（CUDAのスレッドブ

ロック数を指定）
gang, worker, vector を用いて指定する。大まかに以下のように考えると良い。
- gang: CUDA の thread block 数の指定
- vector: CUDA の block 内の threads 数の指定
ü ループがデータ独立であることを指定 (independent clause)
データ独立でないと並列化できない。C言語ではコンパイラがしばしばデータ独立であることを判断できないので、その場合はこれを指定。
ü リダクション処理 (reduction clause)ü 逐次処理 (seq clause)55データの独立性
n independent 指示節 により指定
ü ループがデータ独立であることを明示する
ü コンパイラが並列化できないと判断したときに使用する#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;}n

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

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

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での実行領域を指定。
loop 指示文で並列処理の最適化。
data 指示文でデータ転送を制御。
kernels 指示文でデータ転送を自動的に行うこともできる。

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));cudaMemcpy(a_d, a, n*sizeof(float), cudaMemcpyDefault);
cudaMemcpy(b_d, b, n*sizeof(float), cudaMemcpyDefault);
cudaMemcpy(d_d, d, n*sizeof(float), cudaMemcpyDefault);calc(n, a_d, b_d, c, d_d);cudaMemcpy(d, d_d, n*sizeof(float), cudaMemcpyDefault);
...
}

57

OpenACCコードのコンパイル
n

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をターゲットとし、compute capability 6.0 (cc60) のコードを生成する。

nMakefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load pgi/17.1$ make58

簡単なOpenACCコードn サンプルコード: openacc_basic/
ü OpenACC指示文 kernels, data, loop を利用したコード
ü 計算内容は簡単な四則演算
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];}}}ü ソースコードopenacc_basic/01_original

CPUコード。

openacc_basic/02_kernels

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

openacc_basic/03_data

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

openacc_basic/04_present

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

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

配列のインデックス計算
n サンプルコード: openacc_basic/
ü OpenACC指示文 kernels, data, loop を利用したコード
ü 計算内容は簡単な四則演算
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];}}}ix = j*nx + i

ny

j

i
nx

60

簡単なOpenACC: CPUコードn 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

n 計算内容ü 配列 a、b、cをそれぞれ 1.0, 2.0, 0.0 で初期化
ü calc関数内で c += a *b を nt(=1000)回実行。
ü この実行時間を測定

61

簡単なOpenACC: kernels 指示文（１）n 02_kernelsコード: calc関数
ü CPUコードにkernels 指示文の追加

openacc_basic/02_kernels

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 -> Dfor (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];}}D->H, deallocate}

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

62

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

ü CPUコードにkernels 指示文の追加
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;
}
}
...
}

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

D->H, deallocate

ü kernels 指示文では data 指示文が使える
ü 上の場合は、copyout を指定
ü GPU上で値を初期化するため、CPUからGPUへのコピーは不要。

63簡単なOpenACC: kernels 指示文（３）n コンパイル
ü データの独立性がコンパイラにはわからず、並列化されない。
$ 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 run64簡単なOpenACC: loop 指示文（１）n 03_loopコード
ü 02_kernelsコードにloop independent の追加

openacc_basic/03_loop

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

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

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 */66

簡単なOpenACC: loop 指示文（３）n 03_loopコードの実行

openacc_basic/03_loop

ü 答えは正しいが、実行時間が大変長い。
$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
70.414 [sec]ü ソースコードをみると、calc関数でカーネル前後にGPUとCPU間のデータ転送が発生する。これが性能低下させている。
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}67簡単なOpenACC: data指示文（１）n 04_dataコード
openacc_basic/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;}}

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

present として振舞う。

for (unsigned int icnt=0; icnt<nt; icnt++) {calc(nx, ny, a, b, c);}}a: deallocate
b: deallocate
c: D->H, deallocate

ü copy/copyin/copyout/create は既にデバイス上確保されているデー

タに対しては何もしない。present として振舞う。（OpenACC2.5
以降）
68ü 配列 a, b, c は利用用途に合わせた指示節を指定。

簡単なOpenACC: data指示文（２）n 04_dataコードの実行ü 答えは正しく、速度が上がった。

openacc_basic/04_data$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
1.174 [sec]69簡単なOpenACC: present指示節
n 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 ではメモリ確保、データ転送をしないため、配列サイズの指定は不要。
ü コードとしては見通しがよい。

70

リダクション計算（１）n リダクション計算
ü 配列の全要素から一つの値を抽出
ü 総和、総積、最大値、最小値など
ü 出力が一つのため、並列化に工夫が必要（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

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

リダクション計算(２）n 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];}n Reduction 指示節
ü acc loop reduction(+:sum)
ü 演算子と対象とする変数（スカラー変数）を指定する。
n 利用できる主な演算子と初期値
ü 演算子: +, 初期値: 0
ü 演算子: *, 初期値: 1
ü 演算子: max, 初期値: leastü 演算子: min, 初期値: largest72簡単なOpenACC: reduction指示節（１）n 06_reductionコード
ü 05_presentコードで reductionを使用

openacc_basic/06_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 に変更。

n 06_reductionコード
ü リダクションコードが生成された。
$ 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)73簡単なOpenACC: reduction指示節（２）n 06_reductionコードの実行ü 答えは正しく、速度が上がった。

openacc_basic/06_reduction

ü 配列 c の転送が削減されたこと、リダクションがGPU上で行われることによる性能向上。
$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
1.089 [sec]74OpenACC化のステップのまとめ
n 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];}}openacc_basic/06_reduction75

OPENACC入門実習

76

実習
n 3次元拡散方程式のOpenACC化ü サンプルコード： openacc_diffusion/01_originaln 3次元拡散方程式のCPUコードにOpenACC の kernels, data,loop 指示文を追加し、GPUで高性能で実行しましょう。
for(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_original77

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

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

78

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

t=n

t = n+1

t = n+2
t = n+3
時間ステップを進める

n サンプルコードは、
ü 計算領域: nx * ny * nz （3次元）
ü 最大タイムステップ: ntとなっている。

79

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

平均後の
自分自身の値80

上下左右の値

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

自分自身の値の４倍

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

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

81

CPUコードn 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

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

OpenACC化(0): Makefile の修正
n 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)=

83

OpenACC化(1): kernels
n 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 して実行してみましょう。
84OpenACC化(2): loop
n 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してみましょう。遅いですが実行できます。

85

OpenACC化(3): データ転送の最適化(1)n 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 にしなくても期待通りに動作します。

86OpenACC化(4): データ転送の最適化(2)n 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

87

PGI_ACC_TIME によるOpenACC 実行の確認nn

PGIコンパイラを利用する場合、OpenACCプログラムがどのように実行されているか、環境変数PGI_ACC_TIMEを設定すると簡単に確認することができる。
Linuxなどでは、環境変数PGI_ACC_TIME を1に設定し、プログラムを実行する。
$ export PGI_ACC_TIME=1$ ./runn

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

88

PGI_ACC_TIME によるOpenACC 実行の確認
n ジョブ実行が終わると、標準エラー出力にメッセージが出力される。
$ 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

カーネル
実行時間

89

MPI復習90

MPIn MPI (Message Passing Interface)n メッセージパッシング用のライブラリの規格の一つü C, C++, Fortran に対応
ü コンパイラの規格、特定のソフトウェアやライブラリではない。

n 分散メモリ並列計算機での並列実行に向くü 分散メモリモデル
ü プロセス間のデータ移動はメッセージn 大規模計算が可能
ü 1プロセッサにおけるメモリサイズを超え、計算時間を短縮
ü 利用可能な並列度はノード数を超えられる91

MPIとOpenMP
n MPI

プロセス

ü 分散メモリモデル
ü プロセス間のデータ移動はメッセージ
ü Critical section の代わりにメッセージによる同期
ü 利用可能な並列度はノード数を超えられる
ü 逐次プログラムからプログラム構造が変更することが多い

データを独立に持つ

n OpenMP

スレッド

ü 共有メモリモデル
ü スレッド間のデータ移動は共有変数で行う
ü Critical section による排他処理
ü メモリ共有できるノード内並列に限られる
ü 逐次プログラムに必要な部分に指示文（#pragmaomp）を入れ簡単に高速化できる（可能性がある）
共通のデータにアクセス
92

MPIの実装
n MPICHü 米国アルゴンヌ国立研究所が開発
n MVAPICH
ü 米国オハイオ州立大学で開発
ü MPICHをベース
ü Reedbush では、以下などでロードできる$ module mvapich2/mxm/2.2/{gnu,intel,pgi}n OpenMPI
訂正ü オープンソース
ü Reedbush では、以下などでロードできる$ module openmpi/gdr/1.10.7/{gnu,intel,pgi}

訂正n ベンダMPI
ü Intel MPI など
ü 大抵、上のどれかがベース

93

MPIプロセスとメモリ
n SPMD(Single Program Multiple Data)ü 複数のプロセスが同一プログラムを実行

n プロセスごとに別のメモリ空間
ü すべての変数は各プロセスで別々ü 協調動作させるには明示的に通信する

n プロセスには固有の番号（ランク, rank）がつく
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);// 全プロセス数取得// ランク取得

ü 0 <= rank < nprocs
ü MPI_COMM_WORLD は全プロセスを含むコミュニケータ。
ü コミュニケータは操作対象となるプロセッサ群を定める。
ü メッセージの送信先、受信先は rank で指定
94

MPIプログラムの大きな構造
mpi_basic/01_hello_mpi

Rank

0

1

2

3

#include <stdio.h>#include <mpi.h>
int main(int argc, char *argv[]){
MPI_Init(&argc, &argv);
int nprocs = 1;
int rank= 0;
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
fprintf(stdout, "rank = %d, nprocs = %d¥n", rank, nprocs);
MPI_Finalize();
return 0;}

MPIの初期化: MPI_Init
MPIの終了: MPI_Finalize
をプログラムの先頭と最後で行う。

変数は各プロセスで独立

95

MPIコードのコンパイルnMPIコードのコンパイルü ここではPGIコンパイラとOpenMPIを利用します。
$ module load pgi/17.1$ module openmpi/gdr/1.10.7/pgi$ mpicc -O3 -c main.c
mpicc

MPIコードをコンパイラできる。オプションは元のコンパイラ（今回はPGI）と同じ。

n

Makefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load pgi/17.1 openmpi/gdr/1.10.7/pgi$ make96

MPIプログラムの実行n

ジョブスクリプト$ cat run.sh#! /bin/sh#PBS -q h-tutorial#PBS -l select=1:mpiprocs=2:ompthreads=0#PBS -W group_list=gt00#PBS -l walltime=00:05:00
cd $PBS_O_WORKDIR. /etc/profile.d/modules.shmodule load pgi/17.1
module load openmpi/gdr/1.10.7/pgi
mpirun -np 2./runselect=1:

使用するノード数を指定。この場合は1ノード。
mpiprocs=2:

各ノードに割り当てるMPIプロセス数。この場合は2プロセス。
mpirun -np 2 ./run-np に全並列数（全MPIプロセス数）を指定。select * mpiprocs となることが多い。./run が実行したいプログラム。
97

MPI関数n システム制御
ü MPI_Init, MPI_Comm_rank, MPI_Comm_size, MPI_Finalize,

MPI_Barriern 1対1通信ü ブロックキング通信
ü MPI_Send, MPI_Recv
ü ノンブロッキング通信
ü MPI_Isend, MPI_Irecv, MPI_Wait

n 集団通信ü MPI_Bcast, MPI_Scatter, MPI_Gather, MPI_Allgather, MPI_Alltoall,

MPI_Reduce, MPI_Allreducen 時間計測ü MPI_Wtime

98

C言語とFortranインターフェースの違い（１）n エラーコードの扱い
ü C言語版は、整数の戻り値(err)err = MPI_Xxxx(...);
ü Fortran版は、最後に整数引数 err を渡して取得

call MPI_XXXX(..., err)n システム用配列の確保の方法
ü C言語版

MPI_Status status;
ü Fortran版は、最後に整数引数 err を渡して取得

integer status(MPI_STATUS_SIZE)99

C言語とFortranインターフェースの違い（２）n MPIデータ型
n C言語版
ü MPI_CHAR（文字型）、MPI_INT（整数型）、MPI_FLOAT（単精度実数型）、MPI_DOUBLE（倍精度実数型）
n Fortran版
ü MPI_CHARACTER（文字型）、MPI_INTEGER（整数型）、

MPI_REAL （単精度実数型）、MPI_DOUBLE_PRECISION（倍精度実数型）、MPI_COMPLEX（複素数型）

以降はC言語インターフェースで説明します。
100

基礎的なMPI関数：MPI_Recv
err = MPI_Recv(recvbuf, count, datatype, src, tag, comm, status);src ランクを持つプロセスから送られた識別番号 tag で連続した count 個のメッセージを recvbuf に受信する。メッセージが到着するまで待たされる（ブロックキ
ング）。
n
n
n

n
n
n

n

n

任意
受信領域の先頭アドレス。
count
整数
受信するデータの個数。
datatype
MPIデータ型
受信領域のデータ型。
ü MPI_CHAR（文字型）、MPI_INT（整数型）、MPI_FLOAT(単精度実数型）、MPI_DOUBLE（倍精度実数型）など
src
整数
送信元プロセスのランク。
tag
整数
受信したいメッセージの識別番号。
comm
コミュニケータ コミュニケータ。
ü 通常は MPI_COMM_WORLD を指定する。
statusMPI_Status型
受信状況に関する情報が入る。
ü 必ず専用の型で受け取る。
err
整数
エラーコードが入る。
101
recvbuf基礎的なMPI関数：MPI_Send
err = MPI_Send(sendbuf, count, datatype, dst, tag, comm);

sendbuf から連続した count 個のメッセージを、識別番号 tag をつけて、dst ランクをもつプロセスへ送信する。メッセージ送受信が完了するまで待たされる（ブロックキング）。
nsendbuf

n

count

n

datatype

n

dst

n

tag

n

comm

n

err

任意
整数
MPIデータ型
整数
整数
コミュニケータ
整数

送信領域の先頭アドレス。
送信するデータの個数。
送信領域のデータ型。
送信先プロセスのランク。
送信するメッセージの識別番号。
コミュニケータ。
エラーコードが入る。

なお、MPI_Send/MPI_Recvはブロックキングするため、注意深く順番を並べて関数を呼ばないとプログラムがデッドロックする。

102

1対1通信（ブロッキング）
mpi_basic/02_recv_send
const int n = 16;
int a[n];int b[n];for (int i=0; i<n; i++) {a[i] = i;b[i] = 0;}
if (rank == 1) {MPI_Send(a, n, MPI_INT, 2, 100, MPI_COMM_WORLD);} else if (rank == 2) {MPI_Status status;
MPI_Recv(b, n, MPI_INT, 1, 100, MPI_COMM_WORLD, &status);}rank 0rank 1

rank 2

rank 3

MPI_Senda
ブロックキング
再開

MPI_Recvb103基礎的なMPI関数：MPI_Irecv
err = MPI_Irecv(recvbuf, count, datatype, src, tag, comm, request);

MPI_Recvのノンブロッキング版。送受信が完了しなくても関数は直ちに終了する。
受信を完了したい箇所で request を渡して MPI_Wait を呼ぶ。
n

recvbuf

n

count

n

datatype

n

src

n

tag

n

comm

n

request

n

err

任意
整数
MPIデータ型
整数
整数
コミュニケータ
MPI_Request型
整数受信領域の先頭アドレス。
受信するデータの個数。
受信領域のデータ型。
送信元プロセスのランク。
受信したいメッセージの識別番号。
コミュニケータ。
受信を要求したメッセージにつけられる識別子
エラーコードが入る。

MPI_Irecvが完了する段階では受信が終わっていない可能性があるため、受信情報を保持する MPI_Status は引数に取らない。
送受信が完了しなくても、 MPI_Irecv 直後から送受信と関係ない作業を再開できる。
104

基礎的なMPI関数：MPI_Isend
err = MPI_Isend(sendbuf, count, datatype, dst, tag, comm, request);

MPI_Sendのノンブロッキング版。送受信が完了しなくても関数は直ちに終了する。
送信を完了したい箇所で request を渡して MPI_Wait を呼ぶ。
nsendbuf

n

count

n

datatype

n

dst

n

tag

n

comm

n

request

n

err

任意
整数
MPIデータ型
整数
整数
コミュニケータ
MPI_Request型
整数

送信領域の先頭アドレス。
送信するデータの個数。
送信領域のデータ型。
送信先プロセスのランク。
送信するメッセージの識別番号。
コミュニケータ。
送信を要求したメッセージにつけられる識別子
エラーコードが入る。

送受信が完了しなくても、 MPI_Isend 直後から送受信と関係ない作業を再開できる。

105

基礎的なMPI関数：MPI_Wait
err = MPI_Wait(request, status);

MPI_Isend/MPI_Irecvの送信/受信リクエストを受け取り、その送受信を待つ。
nrequestn

statusnerr

MPI_Request型
MPI_Status型
整数送信を要求したメッセージにつけられる識別子
受信状況に関する情報が入る。
エラーコードが入る。

106

1対1通信（ノンブロッキング）
mpi_basic/03_irecv_isendfor (int i=0; i<n; i++) {a[i] = i + rank * 10;b[i] = 0;}
const int dst = (rank + 1) % nprocs;
const int src = (rank - 1 + nprocs) % nprocs;
MPI_Status sstatus, rstatus;
MPI_Request sreq, rreq;
MPI_Isend(a, n, MPI_INT, dst, 100, MPI_COMM_WORLD, &sreq);
MPI_Irecv(b, n, MPI_INT, src, 100, MPI_COMM_WORLD, &rreq);
MPI_Wait(&rreq, &rstatus);
MPI_Wait(&sreq, &sstatus);rank 0rank 1
直ちに再開

MPI_Isendarank 2
MPI_Irecv

ノンブロックキング

rank 3

直ちに再開

MPI_Wait
MPI_Wait

b

待つ
107

1対1通信（ブロックキング、デッドロック）
mpi_basic/04_recv_send_deadlock//const int n = 16;
const int n = 1024 * 1024;
int a[n];int b[n];for (int i=0; i<n; i++) {a[i] = i + rank * 10;b[i] = 0;}
const int dst = (rank + 1) % nprocs;
const int src = (rank - 1 + nprocs) % nprocs;
MPI_Status status;
MPI_Send(a, n, MPI_INT, dst, 100, MPI_COMM_WORLD);MPI_Recv(b, n, MPI_INT, src, 100, MPI_COMM_WORLD, &status);

mpi_basic/03_irecv_isend を MPI_Send/MPI_Recvで実装すると、通信が循環しているため、デッドロックします。
MPI_Sendがバッファリングされる場合には実行されてしまいますが、されるかどうかは実装依存であるため、それを前提にしてはいけません。
108

基礎的なMPI関数：MPI_Gather
err = MPI_Gather(sendbuf, sendcount, sendtype,recvbuf, recvcount, recvtype, root, comm);

comm 内のランク0から順に送られる sendbuf から連続した sendcount 個のメッセー
ジをランク root の recvbuf へ recvcount個ずつ入れる。
n
n
n
n

n

任意
送信領域の先頭アドレス。
sendcount
整数
送信するデータの個数。
sendtype
MPIデータ型
送信領域のデータ型。
recvbuf
任意
受信領域の先頭アドレス。
ü 原則、送受信領域は別の領域としなければならない。
recvcount
整数
受信するデータの個数。
sendbuf

ここでは各ランクから送られてくる個数（同一の数）を指定する。各ラン
クから異なる数のデータは送れない。一般には sendcount と同数。
recvtype
MPIデータ型
受信領域のデータ型。
root
整数受信するランク。
comm
コミュニケータ コミュニケータ。
err
整数
エラーコードが入る。
ü

n
n
n
n

109

集団通信：収集
mpi_basic/05_gather
const int n = 8;
int a[n];int b[n];for (int i=0; i<n; i++) {a[i] = i + rank * 10;b[i] = 0;}
MPI_Gather(a, 2, MPI_INT, b, 2, MPI_INT, 1, MPI_COMM_WORLD);

array aarray brank 0

0

1

2

3

4

5

6

7

rank 1

10

11

12

13

14

15

16

17

rank 2

20

21

22

23

24

25

26

27

rank 3

30

31

32

33

34

35

36

37

0

1

10

11

20

21

30

31

基礎的なMPI関数：MPI_Scatter
err = MPI_Scatter(sendbuf, sendcount, sendtype,recvbuf, recvcount, recvtype, root, comm);

ランクroot の sendbuf から sendcount個ずつのメッセージを comm 内のランク0から順に全ランクの recvbuf へ recvcount 個入れる。
nsendbufn

sendcount

任意整数

送信領域の先頭アドレス。
送信するデータの個数。

ここでは各ランクへ送る個数（同一の数）を指定する。各ランクへ異なる数のデータは送れない。一般には recvcount と同数。
sendtype
MPIデータ型
送信領域のデータ型。
recvbuf
任意
受信領域の先頭アドレス。
ü 原則、送受信領域は別の領域としなければならない。
recvcount
整数
受信するデータの個数。
recvtype
MPIデータ型
受信領域のデータ型。
root
整数受信するランク。
comm
コミュニケータ コミュニケータ。
err
整数
エラーコードが入る。
ü

n
n

n
n
n
n
n

111

集団通信：分配
mpi_basic/06_scatter
const int n = 8;
int a[n];int b[n];for (int i=0; i<n; i++) {a[i] = i + rank * 10;b[i] = 0;}
MPI_Scatter(a, 2, MPI_INT, b, 2, MPI_INT, 1, MPI_COMM_WORLD);

array aarray brank 0

0

1

2

3

4

5

6

7

10

11

rank 1

10

11

12

13

14

15

16

17

12

13

rank 2

20

21

22

23

24

25

26

27

14

15

rank 3

30

31

32

33

34

35

36

37

16

17

基礎的なMPI関数：MPI_Reduce
err = MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm);

comm 内の全ランクのsendbuf から連続した count 個のメッセージに対して、演算op を適用し、ランク root の recvbuf へ送る。
n
n
n
n
n

n
n
n

任意
送信領域の先頭アドレス。
recvbuf
任意
受信領域の先頭アドレス。
count
整数送受信するデータの個数。
datatype
MPIデータ型
送受信領域のデータ型。
op
MPIオペランド 演算の種類
ü MPI_SUM（総和）、MPI_PROD（積）、MPI_MAX（最大）、MPI_MIN（最小）などroot
整数
受信するランク。
comm
コミュニケータ コミュニケータ。
err
整数
エラーコードが入る。
sendbuf113

集団通信：リダクション
mpi_basic/07_reduce
const int n = 8;
int a[n];int b[n];for (int i=0; i<n; i++) {a[i] = i + rank * 10;b[i] = 0;}
MPI_Reduce(a, b, 2, MPI_INT, MPI_SUM, 1, MPI_COMM_WORLD);

array aarray brank 0

0

1

2

3

4

5

6

7

rank 1

10

11

12

13

14

15

16

17

rank 2

20

21

22

23

24

25

26

27

rank 3

30

31

32

33

34

35

36

37

60

64

基礎的なMPI関数：MPI_Bcast
err = MPI_Bcast(buffer, count, datatype, root, comm);

ランク root の buffer から連続した count 個のメッセージを comm 内の全ランクのbuffer へ送る。
n

n
n
n
n
n

任意
送信および受信領域の先頭アドレス。
ü ランク root のプロセスにとっては送信領域、それ以外には受信領域
count
整数送受信するデータの個数。
datatype
MPIデータ型
送受信領域のデータ型。
root
整数送信するランク。これ以外は受信となる。
comm
コミュニケータ コミュニケータ。
err
整数
エラーコードが入る。
buffer115

集団通信：ブロードキャスト
mpi_basic/08_bcast
const int n = 8;
int a[n];for (int i=0; i<n; i++) {a[i] = i + rank * 10;}
MPI_Bcast(a, 4, MPI_INT, 1, MPI_COMM_WORLD);

array aarray arank 0

0

1

2

3

4

5

6

7

10

11

12

13

4

5

6

7

rank 1

10

11

12

13

14

15

16

17

10

11

12

13

14

15

16

17

rank 2

20

21

22

23

24

25

26

27

10

11

12

13

24

25

26

27

rank 3

30

31

32

33

34

35

36

37

10

11

12

13

34

35

36

37

その他関数と定数
n プロセス間同期
MPI_Barrier(comm);ü コミュニケータ comm 内のすべてのプロセスで同期をとる。

comm 内のプロセスがこれを呼ばないと次へ進まない。
n MPI_PROC_NULLü MPI_Send の dst や MPI_Recv の src に MPI_PROC_NULLを指定すると、何も起こらず（バッファを書き変えずに）、関数が終了する。端領域などで通信をしない場合に、利用できる。

117

OPENACCとMPIによるマルチGPUプログラミング
118

マルチGPUコンピューティングn 複数GPU計算する目的
ü 1個のGPUに搭載されたメモリよりも大きい問題を解きたい。
ü 1個のGPUで計算するよりも高速に計算したい。

n 複数GPU計算の方法
ü 複数GPUをMPIで並列化
ü 複数GPUをOpenMPで並列化
ü 複数GPUを acc_set_device_num (OpenACC) や cudaSetDevice(CUDA) で切り替えながら計算
n 本講習会では、複数ノードに搭載された複数のGPUを活用できる MPI による並列化を行う。

119

Reedbush-H の計算ノード

6.8 GB/sec
16 GB/sec

16 GB/sec

20 GB/sec x 2
120

Reedbush-H の計算ノード模式図
DeviceMemoryGPU 1

DeviceMemoryNVLink

PCI-e

GPU 0

Host
Memory

CPU
Chip
set

DeviceMemoryGPU 1

NVLinkDeviceMemoryHostMemory

GPU 0

CPU

PCI-e

Network
Card

Chip
set
NetworkCardn GPU-GPU通信にはいくつかの種類がある
121

ノード内GPU間通信（１）n ホストメモリ経由のノード内GPU間通信ü

デバイスメモリとホストメモリの間の PCI-eを通る。

DeviceMemoryGPU 1

DeviceMemoryNVLink

PCI-e

GPU 0

Host
Memory

CPU
Chip
set
Network
Card

DeviceMemoryGPU 1

NVLinkPCI-e

DeviceMemoryHostMemory

GPU 0

CPU
Chip
set
NetworkCard122ノード内GPU間通信（２）n NVLink経由のGPU間通信
ü CUDA Inter Process Communication (IPC) を利用し、同じノード内にあるプロセスではホストを経由せずGPU間で直接通信できる。
ReedbushではGPU間をつなぐNVLink経由で高速通信する。
DeviceMemoryGPU 1

DeviceMemoryNVLink

PCI-e

GPU 0

Host
Memory

CPU
Chip
set
Network
Card

DeviceMemoryGPU 1

NVLinkPCI-e

DeviceMemoryHostMemory

GPU 0

CPU
Chip
set
NetworkCard
123

ノード間のGPU間通信（１）n ホストメモリ経由のノード間GPU間通信
ü あるノードのGPU0から別ノードのGPU0へ転送する際に、それぞれのホストメモリを経由して転送する。
DeviceMemoryGPU 1

DeviceMemoryNVLink

PCI-e

GPU 0

Host
Memory

CPU
Chip
set
Network
Card

DeviceMemoryGPU 1

NVLinkPCI-e

DeviceMemoryHostMemory

GPU 0

CPU
Chip
set
NetworkCard124

ノード間のGPU間通信（２）n GPUDirectRDMAによるノード間GPU間通信
ü ホストメモリを経由することなく、GPUとInfiniBand (NetworkCard) 間で直接データ転送 (RDMA) をすることにより異なるノードのGPU間で高速通信する。
DeviceMemoryGPU 1

DeviceMemoryNVLink

PCI-e

GPU 0

Host
Memory

CPU
Chip
set
Network
Card

DeviceMemoryGPU 1

NVLinkPCI-e

DeviceMemoryHostMemory

GPU 0

CPU
Chip
set
Network
Card
125

CUDA-aware MPIn

従来のMPIü MPIの受信領域や送信領域にホストメモリ上のアドレスのみ指定可能。
ü デバイスメモリのデータを転送する際には、一度ホストメモリへコピーが必要。
n

CUDA-aware MPI
ü CUDA（NVIDA GPU向けの開発環境）とMPIによる複数GPU計算では、しばしばデバイスメモリの内容を他ノードのGPUへMPIで転送する。CUDA-aware MPIでは受信領域や送信領域にデバイスメモリ
上のアドレスも指定可能。MPIライブラリ内部ではGPUDirectなどが利用され、MPIライブラリを利用するだけでGPU 間の高速通信技
術を利用することができる。
ü OpenACC からも利用できる。
ü Reedbush では、 CUDA-aware で GPUDirectRDMA(GDR)に対応したOpenMPI (モジュール名: openmpi/gdr/1.10.7/pgi)
MVAPICH2-GDR (モジュール名: mvapich2/gdr/2.2/pgi)が利用できる。
126OpenACCとMPIコードのコンパイルn

OpenACCとMPIコードのコンパイルü ここではPGIコンパイラとOpenMPIを利用します。
$ module load pgi/17.1
GDR対応であることに注意$ module openmpi/gdr/1.10.7/pgi$ mpicc -O3 -acc -Minfo=accel -ta=tesla,cc60 -c main.c-acc: OpenACCコードであることを指示-Minfo=accel:

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

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

nMakefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load pgi/17.1 openmpi/gdr/1.10.7/pgi$ make127簡単なOpenACCとMPIコード（１）n サンプルコード: openacc_mpi_basic/
ü OpenACCとMPIを利用したコード
ü 計算内容は簡単な四則演算と2プロセス間の通信// main.c 内
for (unsigned int i=0; i<n; i++) {a[i] = 3.0 * rank * ny;b[i] = 0.0;}

次のスライドにも説明あります
rank = 1 では a = 3.0 * ny
rank1の a から rank0 の b へ
w * nx = 10nx が転送される。

const int dst_rank = (rank + 1) % nprocs;
const int tag= 0;
if (rank == 0) {MPI_Status status;
MPI_Recv(b, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD, &status);} else {MPI_Send(a, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD);}
double sum = 0.0;
for (unsigned int i=0; i<n; i++) {sum += b[i];}sum/n = (3.0*ny) * (w*nx) / (nx*ny)= 3.0*w= 30.0128簡単なOpenACCとMPIコード（２）n 計算内容ü rank = 1 では、a = 3.0 * ny で初期化
ü rank = 1 の a から rank = 0 の b へ w * nx = 10nx 個の要素が転送される。
ü sum/n = (3.0*ny) * (w*nx) / (nx*ny) = 3.0*w = 30.0 となる
全要素の値: 3.0 * ny

白要素:
0.0

ny

ny
w

w

MPI

nx

nx

Rank 1, array a

Rank 0, array b

赤要素:
3.0 * ny129

簡単なOpenACCとMPIコード（３）n サンプルコード: openacc_mpi_basic/
ü OpenACCとMPIを利用したコード
ü 計算内容は簡単な四則演算と2プロセス間の通信

openacc_mpi_basic/01_original

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

openacc_mpi_basic/02_kernels

MPI+OpenACCコード。上にkernels/loop指示文を追加。単一ノード内2GPU使用。

openacc_mpi_basic/03_update

MPI+OpenACCコード。上にupdate指示文を追加。
単一ノード内2GPU使用。

openacc_mpi_basic/04_cuda_aware

MPI+OpenACCコード。MPI通信にデバイスポインタを渡す。単一ノード内2GPU使用。

openacc_mpi_basic/05_no_gdr

MPI+OpenACCコード。上を2ノードにある2GPU使用。

openacc_mpi_basic/06_gdr

MPI+OpenACCコード。上にGPUDirectRDMAを適用。
2ノードにある2GPU使用。
130

簡単なOpenACCとMPIコード: CPUコード(1)n CPUコードのコンパイルと実行ü 配列の平均値と実行時間が出力されています。
$ cd openacc_mpi_basic/01_original$ make$ qsub ./run.sh$ cat run.sh.o??????
mean = 30.00
Time =
0.078 [sec]? の数字はジョブごとに変わります。
答えは常に30.0
openacc_mpi_basic/01_original131簡単なOpenACCとMPIコード: CPUコード(2)
n 1ノード内の2GPUを利用ü ジョブスクリプトを確認します。
$ cat run.sh#! /bin/sh#PBS -q h-tutorial#PBS -l select=1:mpiprocs=2:ompthreads=0#PBS -W group_list=gt00#PBS -l walltime=00:05:00
cd $PBS_O_WORKDIR. /etc/profile.d/modules.shmodule load pgi/17.1
module load openmpi/gdr/1.10.7/pgi
mpirun -np 2 ./runselect=1:

使用するノード数を指定。この場合は1ノード。
mpiprocs=2:

各ノードに割り当てるMPIプロセス数。この場合は2プロセス。
mpirun -np 2-np に全並列数（全MPIプロセス数）を指定。 select * mpiprocs となること132が多い。./run が実行したいプログラム。

MPIランクとGPU割り当て
n OpenACC関数を呼ぶためヘッダーを追加#include <openacc.h>n ノード上のGPU数の取得
const int ngpus = acc_get_num_devices(acc_device_nvidia);

acc_device_nvidiaを指定することで、NVIDIA GPUを数える。

n あるプロセスに特定のGPUを割り当て
int rank= 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
acc_set_device_num(rank % ngpus, acc_device_nvidia);

device_num を rank % ngpus とすることで、ノード内のプロセスは異なるGPUを利用することとなる。

openacc_mpi_basic/02_kernels
133kernels 指示文追加コード（１）n 02_kernelsコード
ü kernels, loop を利用したコード。MPIはホストメモリを参照。
#pragma acc kernels copyout(a[0:n], b[0:n])#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {a[i] = 3.0 * rank * ny;b[i] = 0.0;}

openacc_mpi_basic/02_kernels
allocate

D->H, deallocate

const int dst_rank = (rank + 1) % nprocs;
const int tag= 0;
if (rank == 0) {MPI(H)MPI_Status status;
MPI_Recv(b, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD, &status);} else {MPI_Send(a, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD);}MPI(H)

double sum = 0.0;
#pragma acc kernels copyin(b[0:n])#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += b[i];}allocate, H->D
deallocate
134kernels 指示文追加コード（２）n 通信方法ü MPI前後で配列の全領域をデバイスメモリとホストメモリ間でコピー
ü ホストメモリ経由のノード内GPU間通信
n コンパイルと実行$ make$ qsub ./run.sh$ cat run.sh.o??????
num of GPUs = 2
mean = 30.00
Time =
0.233 [sec]135

update 指示文
n data 指示文などで既にデバイスメモリ上に確保されたデータに対して、それ自体または対応するホストメモリを対となるメモリの値で更新する。
ü memcpy(H->D), memcpy(D->H) の機能を有すると思えば良い。

n update 指示文の 主な指示節
n hostü memcpy(D->H)n deviceü memcpy(H->D)n update 指示文の 使い方例...
#pragma acc data copy(a[0:nx * ny]){...
#pragma acc update host(a[0:nx])
host_func(a, nx)#pragma acc update device(a[0:nx])
...
}
...

allocate, H->D
D->H
H->D
deallocate

136

update指示文追加コード（１）n 03_update コード

openacc_mpi_basic/03_updateü data, update を追加したコード。MPIはホストメモリを参照。
double sum = 0.0;
#pragma acc data create(a[0:n], b[0:n]){#pragma acc kernels copyout(a[0:n], b[0:n])#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {a[i] = 3.0 * rank * ny;b[i] = 0.0;}

allocate

const int dst_rank = (rank + 1) % nprocs;
const int tag= 0;
if (rank == 0) {MPI(H)MPI_Status status;
MPI_Recv(b, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD, &status);#pragma acc update device(b[0:w *nx])
H->D} else {#pragma acc update host(a[0:w *nx])
D->H
MPI_Send(a, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD);}MPI通信前後に挿入

MPI(H)#pragma acc kernels copyin(b[0:n])#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += b[i];}}deallocate

137

update指示文追加コード（２）n 通信方法ü MPI前後で配列の必要な領域のみを update 指示文でデバイスメ

モリとホストメモリ間でコピー
ü 02_kernels と比較し、メモリの確保・解放とデータのコピー量が削減されている。
ü ホストメモリ経由のノード内GPU間通信
n コンパイルと実行$ make$ qsub ./run.sh$ cat run.sh.o??????
num of GPUs = 2
mean = 30.00
Time =
0.127 [sec]138host_data 指示文
n data 指示文などで既にデバイスメモリ上に確保されたデータに対して、並列領域の外でデバイスメモリ側のアドレスを使うための指示文
n デバイス側のアドレスを使いたい例
ü GPU用のライブラリの呼び出し
ü CUDA で書かれた関数を呼ぶ
ü CUDA-aware MPIによる通信（GPUDirectの利用）...
#pragma acc data copy(a[0:n]){...
#pragma acc host_data use_device(a){cuda_func(a, n)
}
...
}
...

allocate, H->D

host_data 内ではホストコードにも関わらず a はデバイス側のアドレスが使われる。
deallocate

139

CUDA-aware MPIコード（１）n 04_cuda_aware コード

openacc_mpi_basic/ 04_cuda_aware

ü host_data を利用したコード。MPIはデバイスメモリを参照。
double sum = 0.0;
#pragma acc data create(a[0:n], b[0:n]){#pragma acc kernels copyout(a[0:n], b[0:n])#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {a[i] = 3.0 * rank * ny;b[i] = 0.0;}

allocate

const int dst_rank = (rank + 1) % nprocs;
const int tag= 0;
if (rank == 0) {MPI_Status status;#pragma acc host_data use_device(b)MPI(D)
MPI_Recv(b, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD, &status);} else {#pragma acc host_data use_device(a)MPI(D)MPI_Send(a, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD);}#pragma acc kernels copyin(b[0:n])#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += b[i];}}deallocate

140

CUDA-aware MPIコード（２）n 通信方法ü MPI_Send, MPI_Recv に対してhost_data 指示文でデバイス側のアドレスを渡している。CUDA-aware MPI であれば正しく動作する。
ü 03_update と比較し、CUDA IPC が利用され、ホストメモリを介さない通信となる。NVLink経由でGPU間通信する。
n コンパイルと実行$ make$ qsub ./run.sh$ cat run.sh.o??????
num of GPUs = 2
Rank 0: hostname = a086, GPU num = 0
Rank 1: hostname = a086, GPU num = 1
mean = 30.00
Time =
0.131 [sec]n CUDA IPC の無効化
ü OpenMPI では、デフォルトで CUDA IPC は有効となる。これを無効とするには、実行時に btl_smcuda_use_cuda_ipcを0とする。
mpirun -np 2 --mca btl_smcuda_use_cuda_ipc 0 ./run141

2ノードの2GPUでの実行
n 各ノード1GPUで2ノードを利用 openacc_mpi_basic/ 05_no_gdr
ü ジョブスクリプトを確認します。
$ cat run.sh#! /bin/sh#PBS -q h-tutorial#PBS -l select=2:mpiprocs=1:ompthreads=0#PBS -W group_list=gt00#PBS -l walltime=00:05:00
cd $PBS_O_WORKDIR. /etc/profile.d/modules.shmodule load pgi/17.1
module load openmpi/gdr/1.10.7/pgi
mpirun -np 2 ./runselect=2:

使用するノード数を指定。この場合は2ノード。
mpiprocs=1:

各ノードに割り当てるMPIプロセス数。この場合は1プロセス。
mpirun -np 2:-np に全並列数（全MPIプロセス数）を指定。この場合は2プロセス。

142

GPUDirectRDMA無効コード
n 05_no_gdr コード

openacc_mpi_basic/ 05_no_gdr

ü main.c 自体は 04_cuda_aware と同じ。run.sh のみ変更。

n 通信方法ü MPI_Send, MPI_Recv に対してhost_data 指示文でデバイス側のアドレスを渡している。CUDA-aware MPI であれば正しく動作する。
ü 04_cuda_aware と比較し、2つのGPUはノードが異なるため、MPI通信は InfiniBand を通過する。デフォルトでは、ホストメモ
リ経由のノード間GPU間通信となる。
n コンパイルと実行$ make$ qsub ./run.sh$ cat run.sh.o??????
num of GPUs = 2
Rank 0: hostname = a086, GPU num = 0
Rank 1: hostname = a087, GPU num = 1
mean = 30.00
Time =
0.034 [sec]143GPUDirectRDMA有効コード（１）n 06_gdr コード

openacc_mpi_basic/ 06_gdr

n main.c 自体は 04_cuda_aware と同じ。run.sh のみ変更。

n GPUDirectRDMA の有効化
ü OpenMPI では、デフォルトで GPUDirectRDMA は無効である。

これを有効にするには、実行時に btl_openib_want_cuda_gdrを1とする。
mpirun -np 2 -mca btl_openib_want_cuda_gdr 1 ./runn 通信方法ü 05_no_gdr と比較し、 GPUDirectRDMA を有効とする。ホストメ

モリを経由することなく、直接データ転送 (RDMA) をする。

144

GPUDirectRDMA有効コード（２）n コンパイルと実行

openacc_mpi_basic/ 06_gdr$ make$ qsub ./run.sh$ cat run.sh.o??????
num of GPUs = 2
Rank 1: hostname = a087, GPU num = 1
Rank 0: hostname = a086, GPU num = 0
mean = 30.00
Time =
0.051 [sec]145通信方法によるバンド幅比較
n GPU-GPU間通信（一方向通信）

IPC

GDR

No GDRü openmpi-gdr/1.10.7/pgi による測定

No IPC146

OpenACC+MPI化のステップのまとめ
n OpenACC+MPI化
ü acc_set_device_numによるMPIランクとGPUの対応付け
ü kernels, loop, data 指示文を用いてOpenACC化
ü host_data指示文を用い、MPIにデバイス側のアドレスを渡す
double sum = 0.0;
#pragma acc data create(a[0:n], b[0:n]){#pragma acc kernels copyout(a[0:n], b[0:n])#pragma acc loop independentfor (unsigned int i=0; i<n; i++) {a[i] = 3.0 * rank * ny;b[i] = 0.0;}
const int dst_rank = (rank + 1) % nprocs;
const int tag= 0;
if (rank == 0) {MPI_Status status;#pragma acc host_data use_device(b)MPI_Recv(b, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD, &status);} else {#pragma acc host_data use_device(a)MPI_Send(a, w * nx, MPI_FLOAT, dst_rank, tag, MPI_COMM_WORLD);}#pragma acc kernels copyin(b[0:n])#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += b[i];}}openacc_mpi_basic/06_gdr

147

OPENACCとMPIによるマルチGPUプログラミング実習
148

実習
n MPI並列化された3次元拡散方程式のOpenACC化
n サンプルコード： openacc_mpi_diffusion/01_original
ü 3次元拡散方程式のCPUコードにOpenACC のacc_set_device_num

でGPUを割り当て、 kernels, data, loop, host_data 指示文を追加し、複数GPUで高性能で実行しましょう。
for(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+mgn) + 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 = (rank == nprocs - 1 && k == nz - 1) ? ix : ix + nx*ny;
const int km = (rank == 0&& 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_mpi_diffusion/01_original149マルチGPUのための領域分割
n 並列計算するために計算領域を分割する
ü z 方向で分割する。各計算領域はデータ交換のための袖領域（halo）を持つ。
袖領域
厚さ mgn =1

rank0,
GPU0
袖領域
交換
(MPI)
rank1,
GPU1

nz0

袖領域
交換
(MPI)
nz

lnz

nx0
nx(=nx0)

rank2,
GPU2
150

時間発展
n diffusion3d 関数の前に袖領域のデータ交換のために

MPI_Send/MPI_Recv が実行される。
for (; icnt<nt && time + 0.5*dt < 0.1; icnt++) {if (rank == 0 && icnt % 100 == 0)
fprintf(stdout, "time(%4d) = %7.5f¥n", icnt, time);

main.c, main 関数内

const int tag = 0;
MPI_Status status;
MPI_Send(&f[nx*ny*nz]
MPI_Recv(&f[0], nx*ny, MPI_FLOAT, rank_up , tag, MPI_COMM_WORLD);, nx*ny, MPI_FLOAT, rank_down, tag, MPI_COMM_WORLD, &status);

MPI_Send(&f[nx*ny*mgn], nx*ny, MPI_FLOAT, rank_down, tag, MPI_COMM_WORLD);MPI_Recv(&f[nx*ny*(nz+mgn)], nx*ny, MPI_FLOAT, rank_up , tag, MPI_COMM_WORLD, &status);
flop += diffusion3d(nprocs, rank, nx, ny, nz, mgn, dx, dy, dz, dt, kappa, f, fn);swap(&f, &fn);time += dt;}

openacc_mpi_diffusion/01_original
nx*ny*(nz+mgn)
nx*ny*nznz
nx*ny*mgn
0

nx
mgn

151

CPUコードn CPUコードのコンパイルと実行$ cd openacc_mpi_diffusion/01_original$ make$ qsub ./run.sh# cat run.sh.o??????
nprocs = 4
Rank 1: hostname = a086
Rank 0: hostname = a086
Rank 2: hostname = a087
Rank 3: hostname = a087
time(0) = 0.00000
time( 100) = 0.00610...
time(1500) = 0.09155
time(1600) = 0.09766
Time =
6.208 [sec]Performance=
1.80 [GFlops/nprocs]
7.19 [GFlops]
Error[128][128][128] = 4.556413e-06

実行性能（プロセスあたり）
実行性能（合計）
解析解との誤差n OpenACCコードでは、どのくらいの実行性能が達成できるでしょうか？
152

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

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

153

OpenACC化(1): GPU割り当て
n ヘッダー追加とmain関数に acc_set_device_num等を追加...
#include <openacc.h>...

main.c, main 関数内

int main(int argc, char *argv[]){...
const int ngpus = acc_get_num_devices(acc_device_nvidia);
if (rank == 0) {fprintf(stdout, "num of GPUs = %d¥n", ngpus);}
const int gpuid = ngpus > 0 ? rank % ngpus : -1;
if (gpuid >= 0) {acc_set_device_num(gpuid, acc_device_nvidia);}
if (rank == 0) {fprintf(stdout, "OMPI_MCA_btl_smcuda_use_cuda_ipc = %s¥n",getenv("OMPI_MCA_btl_smcuda_use_cuda_ipc"));fprintf(stdout, "OMPI_MCA_btl_openib_want_cuda_gdr = %s¥n",getenv("OMPI_MCA_btl_openib_want_cuda_gdr"));}...

OMPI_MCA_btl_smcuda_use_cuda_ipc環境変数などを確認することで CUDA IPC や GDRが有効/無効か確認できる

154

OpenACC化(2): kernels, loop
n diffusion3d関数に kernels, loopを追加。present 指定。

diffusion.c, diffusion3d 関数内#pragma acc kernels present(f, fn)#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+mgn) + 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 = (rank == nprocs - 1 && k == nz - 1) ? ix : ix + nx*ny;
const int km = (rank == 0&& 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;}

155OpenACC化(3): data, host_datan main関数で data, host_data を追加#pragma acc data copy(f[0:ln]) create(fn[0:ln]){start_timer();main.c, main 関数内

for (; icnt<nt && time + 0.5*dt < 0.1; icnt++) {if (rank == 0 && icnt % 100 == 0)
fprintf(stdout, "time(%4d) = %7.5f¥n", icnt, time);
const int tag = 0;
MPI_Status status;#pragma acc host_data use_device(f){MPI_Send(&f[nx*ny*nz]
MPI_Recv(&f[0], nx*ny, MPI_FLOAT, rank_up , tag, MPI_COMM_WORLD);, nx*ny, MPI_FLOAT, rank_down, tag, MPI_COMM_WORLD, &status);

MPI_Send(&f[nx*ny*mgn], nx*ny, MPI_FLOAT, rank_down, tag, MPI_COMM_WORLD);MPI_Recv(&f[nx*ny*(nz+mgn)], nx*ny, MPI_FLOAT, rank_up , tag, MPI_COMM_WORLD, &status);}
flop += diffusion3d(nprocs, rank, nx, ny, nz, mgn, dx, dy, dz, dt, kappa, f, fn);swap(&f, &fn);time += dt;}
MPI_Barrier(MPI_COMM_WORLD);
elapsed_time = get_elapsed_time();}156

OpenACC化(4): 性能測定
n まずはそのまま make して実行してみましょう。
n 条件を変更して性能測定してみましょう。
ü 計算格子サイズを 1283 から 5123 へ変更
ü GPUでは計算領域が大きいほうが性能が出やすい。
ü 2 GPU x 2 node から 1 GPU x 2 node または 2 GPU x 1 node へ変更
ü CUDA IPC や GDR の有効/無効の変更

n PGI_ACC_TIME によるOpenACC 実行の確認

OpenACC化の例は、openacc_mpi_diffusion/02_openacc

157

複数GPUを用いた
FDTD法による電磁波伝搬計算
158

電磁波の方程式n 真空での電場Eと磁場Hの時間発展

Maxwell 方程式の一部
1
⇥E
=
H
⇥t
(

: 誘電率)

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.

159

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

13
12
11
10

inside̲global.length[1] ( = 7 )98プログラムでは、配列は1次元的に確保します。
whole_global.length[0]= inside_global[0] +2*mgn + 1

であることに注意

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)160

FDTD法（１）n 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
161

FDTD法（２）n 例えば、

162

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

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

openacc_mpi_fdtd/02_openacc1

calc_ex_ey, pml_boundary_ex, pml_boundary_ey,
がOpenACC。

openacc_mpi_fdtd/03_openacc2

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

openacc_mpi_fdtd/04_openacc3

初期化を含め OpenACC。

openacc_mpi_fdtd/05_openacc4

データ移動の最適化。

163

ソースコード（２）n それぞれのファイルの内容
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” のソー
スコードを追記・修正していきます。

164

計算条件
n 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

165

マルチGPU化のための領域分割
mgn

GPU0 Subdomain 0
分割領域の中心領域mgn

nyy 方向で全領域を分割
各分割領域に袖領域を持たせます。

mgndomain分割領域の全体領域

mgn

GPU1

ny

ny0

Subdomain 1

mgn

全領域の中心領域

mgn

mgn

GPU2

mgn

nx0

mgn

ny

Subdomain 2

全領域の全体領域mgn

全体領域：袖領域（PML、青格子）を含む領域
中心領域：袖領域を含まない領域

mgn

nx

166

mgn

コード全体の流れ（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 関数
計算結果出力167
MPI_Gather, write_bmp関数

計算領域の設定（１）n 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
} };

168

mgn

計算領域の設定（２）

mgn

mgn

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

mgn

ny

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

ny

const struct Range whole

分割領域の
中心領域

mgn

ny0

const struct Range insidemgn

struct Range {int length[2];int begin [2];};座標(inside.begin[0], inside.begin[1])座標(whole.begin[0], whole.begin[1])

mgn

mgn

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

nyn Range 構造体

mgn

nx

mgn

169

配列の確保n 物理変数配列は 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);n 多くの配列は whole.length[0] * whole.length[1]n ex_global, ey_global, hz_global はファイル出力に使うため、
whole_global.length[0] * whole_global.length[1]170計算結果出力n 各ランクの ex を ex_global へMPI_Gatherで収集
ü ey, hz も同様

mgn

ny

mgn

// main.cconst int rank_root = 0;
const int sendnelems = whole.length[0] * inside.length[1];const int src= whole.length[0] * (inside.begin[1] - whole.begin[1]);const int dst= whole.length[0] * (inside.begin[1] - whole.begin[1]);

ny

ny0

mgn

mgn

MPI_Gather(&ex[src], sendnelems, MPI_FLOAT_T, &ex_global[dst],
sendnelems, MPI_FLOAT_T, rank_root, MPI_COMM_WORLD);MPI_Gather(&ey[src], sendnelems, MPI_FLOAT_T, &ey_global[dst],
sendnelems, MPI_FLOAT_T, rank_root, MPI_COMM_WORLD);MPI_Gather(&hz[src], sendnelems, MPI_FLOAT_T, &hz_global[dst],
sendnelems, MPI_FLOAT_T, rank_root, MPI_COMM_WORLD);

mgn

mgn

mgn

ny

mgn

袖領域を落とすため、y方向には中心領域のみ切り出して、
nx0
MPI̲Gather する

mgn

mgn

if (rank == rank_root) {write_bmp(icnt, time, whole_global.length, dx, dy, ex_global, ey_global, hz_global);}nx
src = whole.length[0] * (inside.begin[1] - whole.begin[1])
mgn

171

mgn

時間発展（１）n 前半
ü hz の袖領域更新（MPI_Send/MPI_Recv)、
ü 電場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;
(後半へ）

172

時間発展（２）n 後半
ü ex の袖領域更新（MPI_Send/MPI_Recv)、
ü 磁場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++;(出力など)}

173

境界領域のデータ交換（１）hz の境界領域のデータ交換

データアクセスが非対称のため、ここでは hz のみ通信すればよい。
lnx * (jb - jbw - 1)

= src_hz

Datasize = lnx * 1= nhalomgn

subdomains[i]

ny

mgn

n

(ib, jb)lnx * (jb + ny - jbw - 1)

= dst_hz

mgn

subdomains[i-1]

ny

mgn

(ibw, jbw)

(ibw, jbw)

mgn(ib, jb)nx

mgn

nx = inside.length[0]ny = inside.length[1]ib = inside.begin[0]jb = inside.begin[1]lnx = whole.length[0]lny = whole.length[1]ibw = whole.begin[0]jbw = whole.begin[1]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);174
MPI_Recv(&hz[dst_hz], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD, &status);

境界領域のデータ交換（２）ex の境界領域のデータ交換

データアクセスが非対称のため、ここでは ex のみ通信すればよい。

mgnn= src_exnylnx * (jb - jbw)

mgn

subdomains[i+1](ib, jb)

Datasize = lnx * 1= nhalolnx * (jb + ny - jbw)

= dst_ex

mgn

subdomains[i]

ny

mgn

(ibw, jbw)

(ibw, jbw)

mgn(ib, jb)nx

mgn

nx = inside.length[0]ny = inside.length[1]ib = inside.begin[0]jb = inside.begin[1]lnx = whole.length[0]lny = whole.length[1]ibw = whole.begin[0]jbw = whole.begin[1]const int nhalo= whole.length[0];const int inside_end1 = inside.begin[1] + inside.length[1];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);175
MPI_Recv(&ex[dst_ex], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD, &status);

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

プログラムのコンパイルと実行（１）n 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_mpi_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 すると出力なしで実行する。性能測定用。
177

プログラムのコンパイルと実行（２）n プログラムの実行時オプション$ cat run.sh
openacc_mpi_fdtd/01_original#! /bin/sh#PBS -q h-tutorial#PBS -l select=2:mpiprocs=2:ompthreads=0（省略）
mkdir -p sim_run
cd sim_run
nprocs=4
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 を指定すると出力しない。
178

計算結果の表示
n 計算結果は sim_run に BMP として出力される$ cd sim_run/

openacc_mpi_fdtd/01_originaln 計算結果の表示ü 1枚のBMPを見る$ display e05000.bmp
ü 複数のBMPファイルをアニメーションで表示$ animate *.bmp

なおssh -Y txxxxx@reedbush.cc.u-tokyo.ac.jpと –Y をつけていないと表示されない。うまく表示できない場合
は画像を手元にコピーして表示してください。
179計算結果の例
n 出力されたBMPファイルの一例

ny0

mgn

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

y

mgn

入射光

mgn

z

x

nx0

mgn

180

実習1
n calc_ex_ey, pml_boundary_ex, pml_boundary_ey をOpenACC化しましょう。

n Makefileü コンパイルオプションの修正

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

n main.cü OpenACCヘッダーの追加
ü GPUの割り当て
ü data 指示文の追加
ü MPI_Send/MPI_Recvに対するhost_data指示文の追加
n fdtd2d.cü kernels 指示文、loop 指示文の追加
n run.shü GPUDirectRDMA の有効化

解答例は、 openacc_mpi_fdtd/02_openacc1

181

data , host_data指示文
n main関数のwhile 内で data, host_dataを追加#pragma acc data ¥copy(ex[0:nelems], ey[0:nelems], hz[0:nelems])¥
copyin(cexly[0:nelems], ceylx[0:nelems], chzlx[0:nelems], chzly[0:nelems]) ¥copy(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]){#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);}
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);} // acc data
const int j_in = 0;
plane_wave_incidence(&whole, &inside, time, j_in, wavelength, ex, ey);time += 0.5*dt;
182

kernels, loop指示文
n 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]);
}
}
(省略）
}

183

実習2
n main 関数内の while 内をすべて OpenACCにしましょう。
n main.cüüüdata 指示文の移動と copyin などの最適化
MPI_Send/MPI_Recvに対するhost_data指示文の追加
MPI_Gatherに対する host_data 指示文の追加

n fdtd2d.cü

残りの関数にkernels 指示文、loop 指示文の追加

n fdtd2d_sources.cükernels 指示文、loop 指示文の追加

解答例は、 openacc_mpi_fdtd/03_openacc2

184

data 指示文
n 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);}
185

host_data 指示文
n MPI_Gather に対する host_data の追加
const int rank_root = 0;
const int sendnelems = whole.length[0] * inside.length[1];const int src= whole.length[0] * (inside.begin[1] - whole.begin[1]);const int dst= whole.length[0] * (inside.begin[1] - whole.begin[1]);#pragma acc host_data use_device(ex)MPI_Gather(&ex[src], sendnelems, MPI_FLOAT_T, &ex_global[dst],
sendnelems, MPI_FLOAT_T, rank_root, MPI_COMM_WORLD);#pragma acc host_data use_device(ey)MPI_Gather(&ey[src], sendnelems, MPI_FLOAT_T, &ey_global[dst],
sendnelems, MPI_FLOAT_T, rank_root, MPI_COMM_WORLD);#pragma acc host_data use_device(hz)
MPI_Gather(&hz[src], sendnelems, MPI_FLOAT_T, &hz_global[dst],
sendnelems, MPI_FLOAT_T, rank_root, MPI_COMM_WORLD);
if (rank == rank_root) {write_bmp(icnt, time, whole_global.length, dx, dy, ex_global, ey_global,hz_global);
}

186

実習3
n 初期化を含めて全てOpenACCにします。ただし、

set_object_er がCPU上のユーザ定義関数のため、これ以降の初期化関数をOpenACCにします。
n main.cüüdata 指示文の移動と最適化（多くが create になるはずです）
初期化後のMPI_Gatherに対する host_data 指示文の追加

n setup.cükernels 指示文、loop 指示文の追加

解答例は、 openacc_mpi_fdtd/04_openacc3

187

実習4
n 計算領域のサイズやGPU数も変更して性能測定してみましょう。
ü

計算格子サイズを 4096 x 4096 など大きくしてみましょう。

n 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_mpi_fdtd/05_openacc4

188