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

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

1

講習会スケジュールn 開催日時ü 12月4日（月） 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/support/kosyu/ü 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 pgin CUDA開発環境を使う場合$ module load cudan Intelコンパイラを使う場合$ module load inteln MPIを使う場合$ module load openmpi-gdr/2.1.1/{gnu,intel,pgi}$ module load mvapich2-gdr/2.2/{gnu,intel,pgi}
ü PGIなどのコンパイラに追加して loadn モジュールはジョブ実行時にもコンパイル時と同じもの

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

8

サンプルコードのコンパイルn

Reedbush へのログイン$ 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.gzn

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

モジュールをロードする。
$ module load pgin

コンパイルする。
$ cd openacc_hello/01_hello_acc$ maken

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

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

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

10

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

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

2基の大型システム，6年サイクル
09 10 11 12 13 14 15 16 17 18 19 20 21 22

Hitachi	  SR11K/J2IBM	
  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	
  IXfxBDEC	
  system50+	
  PFLOPS	
  (?)

1.13	
  PFLOPS,	
  150 TB

Oakbridge-­‐FX

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

Big Data &Extreme Computing

Reedbush,	
  SGI
Broadwell	
  +	
  Pascal1.93	
  PFLOPS
Reedbush-­‐L
長時間ジョブ実行用演算加速装置
Pascal	  	
  1.4+PF付き並列スーパーコンピュータ

Peta

K

K	
  computer

Post-­‐K	  ?

12

4システム運用中
n Oakleaf-FX (富士通 PRIMEHPC FX10)ü 1.135 PF, 京コンピュータ商用版, 2012年4月 〜 2018年3月

n Oakbridge-FX (富士通 PRIMEHPC FX10)
ü 136.2 TF, 長時間実行用（168時間）, 2014年4月 〜 2018年3月

n Reedbush (SGI, Intel BDW + NVIDIA P100 (Pascal))ü データ解析・シミュレーション融合スーパーコンピュータü 3.361 PF, 2016年7月〜 2020年6月
ü 東大ITC初のGPUシステム (2017年3月より), DDN IME (Burst Buffer)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/3）
Fujitsu  PRIMEHPC  FX10 （FX10スーパーコンピュータシステム）
Total  Peak  performance                  :  1.13  PFLOPS
Total  number  of  nodes                      :  4,800
Total  memory                                                    :  150TB
Peak  performance  per  node    :  236.5  GFLOPSMain  memory  per  node                  :  32  GB
Disk  capacity                                                      :  2.1  PB
SPARC64  IXfx 1.848GHz

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

Oakbridge-­‐FX：長時間ジョブ用のFX10
ノード数：24～576
制限時間：最大168時間
（１週間）
14

東京大学情報基盤センター スパコン（2/3）
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
15

16

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

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

16

FX10TOFU
Network



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

Core
#1

Core
#2

Core
#3

L1

L1

L1

L1





…
: L1





20GB/
Core
#12

Core
#13

Core
#14

Core
#15

32KB L1

L1

L1

L1

ICC

L2 (16  12MB)

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

MemoryMemoryMemoryMemory4GB ×2

4GB ×2

4GB ×2

4GB ×28GB×432GB



DDR3 DIMM

17

FX10の通信網
TOFU



1TOFU


6
5GB/


18

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

トーラスになっています19

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

20

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

  !16GB

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

Memory8"128GB

Memory16GB ×2
16GB ×2
76.8 GB/=(8Byte 2400MHz 4 channel)
21

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
22

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
23

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 Phi
of Memoryon package: innovativeLanding)memoryn Intel Xeon
(Knights
architecture for high bandwidth and high capacity
Integration of Fabric on packagen 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

MCDRAM24

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





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

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

100,000円

までn コース２： 200,000円
まで：

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

16ノード(基準)、最大64ノードn グループコースn

400,000円 (企業 480,000円) ： 1口 8ノード（基準）、最大128ノードまで

n 以上は、「トークン制」で運営
n 申し込みノード数

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

東大情報基盤センター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 以上は、「トークン制」で運営
n 申し込みノード数

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

東大情報基盤センター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 以上は、「トークン制」で運営
n 申し込みノード数

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

28

29

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, Oakforest-PACS
システムが使える「無償トライアルユース」および「有償トライアル
ユース」制度があります。
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 双方のコースともに、簡易な利用報告書の提出が必要
30

スーパーコンピュータシステムの詳細
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/
31

GPU入門2

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

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

3

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

n
n

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

4

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の持つ機能がほとんどないか、全くない。

5

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

7

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

n 並列化GPUは苦手
仕事1
仕事2
仕事3

仕事4

仕事5

時間
8

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

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

n 高い性能を得るには、スレッド数 >> コア数（P100 では

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

10

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

スレッド

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

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

12

Thread
0

...

15 16

...

31

OPENACC入門

13

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

簡単

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

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

難しい

14

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

Fortran!$acc directive-name [clause, …]! Fortran code!$acc end directive-name15

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

OpenMP

OpenACC

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

CPU

CPU

CPU

デバイス
(GPU)
16

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

CPU

GPU

a b
メモリ確保

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

17

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

CPU

GPU

a b
メモリ確保

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

18

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

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

19

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 指示文もある
20

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

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

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

23

OpenACCのデフォルトの変数の扱いn スカラ変数
ü firstprivate または privateü ホストからデバイスへコピーが渡され初期化。ホストに戻せな

い。
n 配列
ü sharedü デバイスメモリに動的に確保され、スレッド間で共有。
ü デバイスからホストへコピーすることが可能。

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

24

データ管理・移動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;}
れ、プライベート変数となる。
25

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以降）
26

データの移動範囲の指定
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])...

27

並列処理の指定
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 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;
}

CPU

GPU

a b
メモリ確保

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から
copyout28

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

ク数を指定）
gang, worker, vector を用いて指定する。大まかに以下のように考えると良い。
- gang: CUDA の thread block 数の指定
- vector: CUDA の block 内の threads 数の指定
ü ループがデータ独立であることを指定 (independent clause)
データ独立でないと並列化できない。C言語ではコンパイラがしばしばデータ独立であることを判断できないので、その場合はこれを指定。
ü リダクション処理 (reduction clause)ü 逐次処理 (seq clause)29データの独立性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];}30

参考：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){= blockIdx.x * blockDim.x + threadIdx.x;
const int ikernelint 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);
...
}

31

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) のコードを生成する。

n

Makefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load pgi/17.1$ make32簡単な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指示節を使用。
33

配列のインデックス計算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

34

簡単な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_originaln 計算内容ü 配列 a、b、cをそれぞれ 1.0, 2.0, 0.0 で初期化ü calc関数内で c += a *b を nt(=1000)回実行。
ü この実行時間を測定

35簡単な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間のメモリ転送が行われる。

36

簡単なOpenACC: kernels 指示文（２）n 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;}}

openacc_basic/02_kernels

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

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

37

簡単な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 */main.o-o runpgcc -O3 -acc -Minfo=accel -ta=tesla,cc60
38簡単な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;}}39

簡単な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 */40簡単な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}}

41簡単な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以降）
ü 配列 a, b, c は利用用途に合わせた指示節を指定。

42

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

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

44

リダクション計算（１）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

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

リダクション計算(２）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, 初期値: largest
46簡単な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)47簡単な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]48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_reduction

49

OPENACC入門実習50

実習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_original51

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

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

52

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

t=n

t = n+1

t = n+2
t = n+3
時間ステップを進めるn サンプルコードは、
ü 計算領域: nx * ny * nz （3次元）
ü 最大タイムステップ: ntとなっている。

53

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

平均後の
自分自身の値54

上下左右の値

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

自分自身の値の４倍

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

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

55

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コードでは、どのくらいの実行性能が達成できるでしょうか？
56

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

57

OpenACC化(1): kernelsn 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 して実行してみましょう。
58

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

59

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 にしなくても期待通りに動作します。

60

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

61

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

$ 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

62

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

カーネル
実行時間

63

MPI復習2

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

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

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

データを独立に持つスレッド

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

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

5

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 で指定
6

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
をプログラムの先頭と最後で行う。

変数は各プロセスで独立

7

MPIコードのコンパイルn

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

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

n

Makefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load pgi/17.1 openmpi/1.10.7/pgi$ make8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/1.10.7/pgi
mpirun -np 2./runselect=1:

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

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

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_Waitn 集団通信ü MPI_Bcast, MPI_Scatter, MPI_Gather, MPI_Allgather, MPI_Alltoall,

MPI_Reduce, MPI_Allreducen 時間計測ü MPI_Wtime

10

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

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言語インターフェースで説明します。
12

基礎的な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
整数
エラーコードが入る。
recvbuf13

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

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

sendbuf

n

count

n

datatype

n

dst

n

tag

n

comm

n

err

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

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

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

14

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_Recvb
15

基礎的な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 直後から送受信と関係ない作業を再開できる。
16

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

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

sendbuf

n

count

n

datatype

n

dst

n

tag

n

comm

n

request

n

err

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

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

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

17

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

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

requestn

statusnerr

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

18

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

rank 3

MPI_Irecv
直ちに再開

ノンブロックキング
MPI_Wait
MPI_Waitb

待つ19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がバッファリングされる場合には実行されてしまいますが、されるかどうかは実装依存であるため、それを前提にしてはいけません。
20

基礎的な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

21

集団通信：収集
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 個入れる。
n

sendbufn

sendcount

任意整数

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

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

n
n

n
n
n
n
n

23

集団通信：分配
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
整数
エラーコードが入る。
sendbuf25

集団通信：リダクション
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
整数
エラーコードが入る。
buffer27

集団通信：ブロードキャスト
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を指定すると、何も起こらず（バッファを書き変えずに）、関数が終了する。端領域などで通信をしない場合に、利用できる。

29

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

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

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

3

Reedbush-H の計算ノード

6.8	
  GB/sec16	
  GB/sec16	
  GB/sec20	
  GB/sec	
  x	
  2
4

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

NVLinkPCI-e

DeviceMemoryHostMemory

GPU 0

CPU
Chip
set

DeviceMemoryGPU 1

NVLinkDeviceMemoryHostMemory

GPU 0

CPU

PCI-e

Network
Card

Chip
set
Network
Cardn GPU-GPU通信にはいくつかの種類がある
5

ノード内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
NetworkCard6ノード内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
Network
Card
7

ノード間の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
NetworkCard
8

ノード間の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
9

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)が利用できる。
10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) のコードを生成する。

n

Makefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load pgi/17.1 openmpi-gdr/1.10.7/pgi$ make11簡単な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;= 0;
const int tag
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.012簡単な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 * ny13

簡単な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使用。
14

簡単なOpenACCとMPIコード: CPUコード (1)n CPUコードのコンパイルと実行ü 配列の平均値と実行時間が出力されています。
$ cd openacc_mpi_basic/01_original$ make$ qsub ./run.sh$ cat run.sh.o??????
mean = 30.00
0.078 [sec]Time =? の数字はジョブごとに変わります。
答えは常に30.0
openacc_mpi_basic/01_original

15

簡単な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/1.10.7/pgi
mpirun -np 2 ./runselect=1:

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

各ノードに割り当てるMPIプロセス数。この場合は2プロセス。
mpirun -np 2-np に全並列数（全MPIプロセス数）を指定。 select * mpiprocs となること16が多い。./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
17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;= 0;
const int tag
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
18kernels 指示文追加コード（２）n 通信方法ü MPI前後で配列の全領域をデバイスメモリとホストメモリ間でコピーü ホストメモリ経由のノード内GPU間通信n コンパイルと実行$ make$ qsub ./run.sh$ cat run.sh.o??????
num of GPUs = 2
mean = 30.00
Time =
0.233 [sec]19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

20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;= 0;
const int tag
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

21

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]22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

23

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;= 0;
const int tag
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

24

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 ./run25

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プロセス。

26

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]27GPUDirectRDMA有効コード（１）n 06_gdr コード

openacc_mpi_basic/ 06_gdrn 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) をする。

28

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]29通信方法によるバンド幅比較
n GPU-GPU間通信（一方向通信）

IPC

GDR

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

No IPC

30

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

31

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

実習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_original

33

マルチGPUのための領域分割n 並列計算するために計算領域を分割するü z 方向で分割する。各計算領域はデータ交換のための袖領域（halo）を持つ。
袖領域
厚さ mgn =1

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

nz0

袖領域
交換
(MPI)
nz

lnz

nx0
nx(=nx0)

rank2, GPU2
34

時間発展
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

35

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コードでは、どのくらいの実行性能が達成できるでしょうか？
36

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

37

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が有効/無効か確認できる

38

OpenACC化(2): kernels, loopn 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;}

39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);, nx*ny, MPI_FLOAT, rank_down, tag, MPI_COMM_WORLD);MPI_Send(&f[nx*ny*mgn]
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();}40

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

41

複数GPUを用いた
FDTD法による電磁波伝搬計算
2電磁波の方程式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.
6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)7

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

En
H

1
En 1
=
t

n+ 12

H

n

1
2

t

H
=

1
2

n

1
µ

En

変形して、
n

n 1

+

=H

n

E =E
H

n+ 21

t
1
2

H
t
µ

n

1
2

En
8

FDTD法（２）n 例えば、9ソースコード（１）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

データ移動の最適化。

10

ソースコード（２）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” のソー
スコードを追記・修正していきます。

11

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

12

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

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

mgn

mgn

ny

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

分割領域の全体領域

GPU1

mgn

domainny

ny0

Subdomain 1

mgn

全領域の中心領域

mgn

mgn

GPU2

mgn

nx0

mgn

ny

Subdomain 2

全領域の全体領域mgn

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

mgn

nx

mgn

13

コード全体の流れ（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関数

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

15

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

nyconst 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

16

配列の確保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]17計算結果出力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

mgn

18

時間発展（１）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.begin[1] - 1);= whole.length[0] * (inside_end1= 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;
(後半へ）

19

時間発展（２）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.begin[1]);= whole.length[0] * (inside_end1

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++;(出力など)}

20

境界領域のデータ交換（１）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.begin[1] - 1);= whole.length[0] * (inside_end1= 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);21
MPI_Recv(&hz[dst_hz], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD, &status);

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

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

mgn

n

= 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.begin[1]);= whole.length[0] * (inside_end1

MPI_Send(&ex[src_ex], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD);22
MPI_Recv(&ex[dst_ex], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD, &status);

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

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

プログラムのコンパイルと実行（２）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 を指定すると出力しない。
25

計算結果の表示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 をつけていないと表示されない。うまく表示できない場合
は画像を手元にコピーして表示してください。
26

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

ny0

mgn

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

y

mgn

入射光

mgn

nx0

mgn

x
z

27

実習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
28

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

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]);
}
}
(省略）
}

30

実習2
n main 関数内の while 内をすべて OpenACCにしましょう。
n main.cdata 指示文の移動と 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
31

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

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

33

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

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

解答例は、 openacc_mpi_fdtd/04_openacc334

実習4n 計算領域のサイズや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

35