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

1

講習会スケジュール
n 開催日時ü 6月30日（金） 13:00 – 17:00

n プログラムü 13:00 – 13:30 スパコンの使い方など
ü 13:30 – 14:30 OpenACC + MPI混合プログラミング（座学）
ü 14:45 – 17:00 OpenACC + MPI演習Ⅱ

2

講習会についてn 本講習会はü OpenACCとMPIを使った複数GPUプログラミング

入門を中心に扱います。
n その他の講習会

https://www.cc.u-tokyo.ac.jp/events/lectures/n スパコンイベント情報メール配信サービス

Youtubeにて過去の講習会を配信中！
https://www.youtube.com/channel/UC
2CHaGp1AO-vqRlV7wmU0w/videos?view=0&sort=p&flow=grid

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

講習会の進め方
n Zoomを利用したオンライン講習会ですü この講義は録画されていますü 質問があるとき以外はミュートでお願いしますü ビデオもオフを推奨します

n slackを使って質問に対応しますü slackはリンクを知っている人は誰でも使える設定になっていますü slackのリンクをzoomのチャットに貼るので、未登録の場合は今のうちに登録お願いしますü slackの登録メールの配送に小一時間かかることがありますü スクリーンショットなどで画像を共有することで、質問対応しますü Windows：Alt + PrtScn で作業中ウィンドウのスクショがクリップボー

ドにコピーされます。slackのチャット部分で貼り付け(Ctrl + V)することで画像をアップロードできますü Mac : command + shift + control + 4 の同時押し、その後撮りたいウィ
ンドウ上でspaceを押すことで、スクリーンショットがクリップボードにコピーされます。slackのチャット部分で貼り付け(command + V)することで画像をアップロードできます4Zoom: 「手を挙げる」方法
1.

Zoomメニュー中の「リアクション」をクリック

2.

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

5

Zoom: 手が挙がっていることの確認方法
1.

Zoomメニュー中の「参加者」をクリックして，参加者一覧を表
示2.

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

Zoom: 「手を降ろす」方法
1.

Zoomメニュー中の「リアクション」をクリック

2.

ポップアップで表示された「手を降ろす」をクリック

7

Slack: 質疑応答チャンネルへの移動
n 左側のメニューバーのチャンネル一覧内に「第209回-マルチgpuプログラミング入門」があるので，クリックn 表示されていない場合
1.
2.
3.

「チャンネルを追加する」をクリック「チャンネル一覧を確認する」をクリック「第209回-マルチgpuプログラミング入門」があるので，「参加する」をクリック

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

メッセージを送信する

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

9

サンプルコードのダウンロード
n 講習会で使うサンプルコードはGitHubにありますn https://github.com/hoshino-UTokyo/lecture_openacc_mpi.git
n 2023年4月にシステムの環境変更があったため、動作しない場合はお知らせくださいn ダウンロード手順1.

$ ssh tXXXXX@wisteria.cc.u-tokyo.ac.jp
n wisteriaへのssh。tXXXXXはアカウント名$ cd /work/gt00/tXXXXX$ git !"#$%&'(()*+,,-.('/01!#2,'#*'.$#345#67#,"%!(/8%9#)%$:!!92).1-.(
4. $ cd lecture_openacc_mpi/
5. $ cd C or F
2.
3.

n

CまたはFortran好きな方を利用してください。

10

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

2001-2005
Hitachi SR8000
1,024 GF

Hitachi
SR2201307.2GF

2006-2010
Hitachi SR11000
J1, J25.35 TF, 18.8 TF

2011-20152016-20202026-2030

Hitachi SR16K/M1
Yayoi
54.9 TF

Hitachi
SR8000/MPP

OBCX(Fujitsu)

2,073.6 GF

6.61 PF

OakforestPACS (Fujitsu)

Hitachi HA8000
T2K Todai
Fujitsu FX10
Oakleaf-FX

東京大学情報基盤センターのスパコン

OFP-II

150+ PF

25.0 PF

140 TF

利用者2,600+名
55％は学外2021-2025

Fujitsu
33.1 PF

1.13 PFBDEC02250+ PF

ReedbushU/H/L (SGI-HPE)
3.36 PF

Ipomoea-01 25PB

Ipomoea03

Ipomoea-02
12

2001-2005
Hitachi SR8000
1,024 GF

SR80002006-2010
Hitachi SR11000
J1, J2
5.35
18.8 TF
IBMTF,
Power5+

Hitachi
SR2201

Hitachi
SR8000/MPP

HARP-1E

SR8000307.2GF

2011-20152016-2020

Hitachi SR16K/M1
Yayoi
54.9
TF
IBM
Power7

6.61 PF

OakforestPACS (Fujitsu)

Hitachi HA8000
T2K Todai
140 TFAMD Opteron

汎用CPU
加速装置付

2026-2030

Intel CLX
OBCX(Fujitsu)

2,073.6 GF

疑似ベクトル2021-2025

OFP-II

150+ PF

25.0 PF

AcceleratorsIntel Xeon Phi
Fujitsu FX10
Oakleaf-FX

SPACR64 IXfx

BDEC02

Fujitsu
33.1 PF

1.13 PF

ReedbushU/H/L (SGI-HPE)

A64FX,
Intel Icelake+
NVIDIA A100

250+ PF

Accelerators

3.36 PFIntel BDW +NVIDIA P100

Ipomoea-01 25PB

Ipomoea03

Ipomoea-02
13

14

柏キャンパス
OFP，OBCX

柏IIキャンパス
Wisteria/BDEC-01，mdx
ABCI/ABCI-II

TX「柏の葉キャンパス」駅
14Platform for Integration of (S+D+L)
Big Data & Extreme Computing

Reedbush (HPE, Intel BDW + NVIDIA P100 (Pascal))

Simulation Nodes:
Odyssey
Fujitsu/Arm A64FX
25.9PF, 7.8 PB/s

Shared File
System
（SFS）
25.8 PB, 500 GB/s

2.0 TB/s

Data/Learning
Nodes: AquariusIntel Ice Lake + NVIDIA A100

Fast File
System（FFS）
1 PB, 1.0 TB/s

7.20 PF, 578.2 TB/s

Oakforest-PACS (OFP) (Fujitsu, Intel Xeon Phi (KNL))• JCAHPC （筑波大CCS・東大ITC），2016年10月～2022年3月末（予定）• 25 PF, #39 in 58th TOP 500 (November 2021)

Oakbridge-CX (OBCX) (Fujitsu, Intel Xeon CLX)

800 Gbps

External
Resources

External
Resources
External Network

Simulation Nodes（Odyssey）• データ解析・シミュレーション融合スーパーコンピュータ• 2016年7月～2021年11月末（引退）• 東大ITC初のGPUクラスタ，ピーク性能3.36 PF

Data/Learning Nodes（Aquarius）• 2019年7月～2023年6月末（予定）• 6.61 PF, #110 in 58th TOP500-June 2023 (Plan)

Wisteria/BDEC-01（Fujitsu）• シミュレーションノード群（Odyssey）：A64FX (#17)• データ・学習ノード群（Aquarius）：Intel Icelake+NVIDIA A100) (#106)• 33.1 PF, #13 in 57th TOP 500，2021年5月14日運用開始• 「計算・データ・学習（S+D+L）」融合のためのプラットフォーム• 革新的ソフトウェア基盤「h3-Open-BDEC」（科研費基盤（S） 2019年度～2023年度）

Reedbush

Oakforest-PACS

Oakbridge-CX 15

Aquariusの構成
n Intel Xeon Platinum 8360Y (36c 2.4GHz) x 2ソケット, 512GBメモリ
n ノード当たり8基のNVIDIA A100 GPU
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

x16

IBHDR

IB HDR:
200 Gbps x 4 link
(ノード当たり)x16

PCIe Gen4
x16

A100

UPI

CPU1
Intel Xeon Platinum
8360Y
(ICX)

CPU2
Intel Xeon Platinum
8360Y
(ICX)

x16

x16

PLX
x16

A100

IBHDR

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

x16

PLX

PLX

x16

メモリバンド幅：
409.6 GB/sec

PLX

x16

x16

x16

x16

x16

A100

A100

A100

A100

IBHDR

x16

A100

PCI Express:
Gen4 x16 = 16 GHz x 16 (128b130b)= 31.5 GB/s (片方向)
x16

IBHDR

PCIe Gen4
x16

A100

NVLink3:
GPU当たり 50 GHz x4x2 x6= 300 GB/s (片方向)NVSwitch
NVSwitch
NVSwitch
NVSwitch
NVSwitch
NVSwitch
16

Wisteria 利用上の注意（１）n ディレクトリについて（home と work)ü ログイン時のディレクトリ（/home/gt00/txxxxx）にはログ

イン時に必要なファイルのみを置く
ü プログラム作成や実行などに必要なファイルは /work 以下のディレクトリ（/work/gt00/txxxxx）に置く
ü /home は計算ノードからは参照できない

17

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

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

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

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

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

Wisteriaでのプログラムの実行
n ジョブスクリプト(〇〇.sh)を作成し、ジョブとして投入、実行する。
$ pjsub ./〇〇.shn 投入されたジョブを確認する。（qstatではないので注意）$ pjstatn 実行が終了すると、以下のファイルが生成される。
〇〇.sh.??????.out〇〇.sh.??????.err（?????? はジョブID）n 上記の標準出力ファイルの中身を確認する。
$ cat 〇〇.sh.??????.out

n 必要に応じて、上記のエラー出力ファイルの中身を確認する。
$ cat 〇〇.sh.??????.out
19コンパイラの種類と実行(Aquarius)n ログインノードとAquarius計算ノードとでは、CPUの命令セットが（ほぼ）同じ
n ログインノード：命令セットアーキテクチャ Intel CascadeLake＋AVX512, x86_64
n Aquarius計算ノード： 命令セットアーキテクチャ Intel IceLake＋AVX512, x86_64

n 様々なコンパイラが利用可能: GPU向けには gcc+CUDAか NVIDIAを推奨$ module load gcc cuda ompi-cuda または$ module load nvidia nvmpi <- 2023年4月にmodule環境変更
GNUコンパイラ

Intelコンパイラ

NVIDIA コンパイラ(旧PGI)

C

gcc

icc

nvc (pgcc)

C++

g++

icpc

nvc++(pgc++)

gfortran

ifort

言語

Fortran
OpenACC

CUDAコンパイラ
nvcc

nvfortran (pgfortran)○
20

JOBスクリプトサンプルの説明（Aquarius, MPIなし）#!/bin/bash#PJM -L rscgrp=lecture-a#PJM -L gpu=2#PJM -L elapse=00:01:00#PJM -g gt00
module load nvidia./a.out

リソースグループ名：lecture-a
利用GPU数
実行時間制限：１分（講習会では最大10分）
利用グループ名：gt00
21

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

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

n 複数GPU計算の方法
ü 推奨：複数GPUをMPIで並列化
ü 最大の難点は環境構築の難しさ。

ü 複数GPUをOpenMPで並列化
ü できなくはない。OpenMPを使う場合、OpenACCと同じループに指示文の適用はできないため、omp_get_thread_num()を用いて並列化することになる。
ü 複数GPUを acc_set_device_num (OpenACC) や cudaSetDevice (CUDA) で切り替えながら計算
n 本講習会では、複数ノードに搭載された複数のGPUを活用できる

MPI による並列化を行う。

23

Aquariusの計算ノード
n Intel Xeon Platinum 8360Y (36c 2.4GHz) x 2ソケット, 512GBメモリ
n ノード当たり8基のNVIDIA A100 GPU
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

x16

IBHDR

IB HDR:
200 Gbps x 4 link
(ノード当たり)x16

PCIe Gen4
x16

A100

UPI

CPU1
Intel Xeon Platinum
8360Y
(ICX)

CPU2
Intel Xeon Platinum
8360Y
(ICX)

x16

x16

PLX
x16

A100

IBHDR

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

x16

PLX

PLX

x16

メモリバンド幅：
409.6 GB/sec

PLX

x16

x16

x16

x16

x16

A100

A100

A100

A100

IBHDR

x16

A100

PCI Express:
Gen4 x16 = 16 GHz x 16 (128b130b)= 31.5 GB/s (片方向)
x16

IBHDR

PCIe Gen4
x16

A100

NVLink3:
GPU当たり 50 GHz x4x2 x6= 300 GB/s (片方向)NVSwitch
NVSwitch
NVSwitch
NVSwitch
NVSwitch
NVSwitch24

GPUx2 の計算ノード模式図
DeviceMemoryGPU 1

DeviceMemoryNVLink

PCI-e

GPU 0

Host
Memory

CPU
Chip
set
Network
Card

DeviceMemoryGPU 1

NVLinkDeviceMemoryHostMemory

GPU 0

CPU

PCI-e

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

ノード内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
NetworkCard26

ノード内GPU間通信（２）n NVLink経由のGPU間通信
ü CUDA Inter Process Communication (IPC) を利用し、同じノード内にあるプロセスではホストを経由せずGPU間で直接通信できる。最近のサーバでは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
27

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

ノード間のGPU間通信（２）n GPUDirectRDMAによるノード間GPU間通信
ü ホストメモリを経由することなく、GPUとInfiniBand (Network Card) 間で直接データ転送 (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
29

CUDA-aware MPI
n 従来のMPI
ü MPIの受信領域や送信領域にホストメモリ上のアドレスのみ指定可能。
ü デバイスメモリのデータを転送する際には、一度ホストメモリへコピーが必要。
n CUDA-aware MPI
ü CUDA（NVIDA GPU向けの開発環境）とMPIによる複数GPU計算では、しばしばデバイスメモリの内容を他ノードのGPUへMPIで転送する。CUDAaware MPIでは受信領域や送信領域にデバイスメモリ上のアドレスも指定可能。MPIライブラリ内部ではGPUDirectなどが利用され、MPIライブラ
リを利用するだけでGPU 間の高速通信技術を利用することができる。
ü OpenACC からも利用できる。
ü acc host_data 指示文で、デバイスメモリのアドレスを渡せば良いü Unified Memory 使用時には利用できない

ü Wisteria では、 nvidia module をloadした後に出てくるnvmpi module がCUDA-aware で GPUDirectRDMA(GDR)に対応したOpenMPIです。

30

OpenACCとMPIコードのコンパイルn

OpenACCとMPIコードのコンパイルü ここではnvidiaコンパイラとOpenMPIを利用します。
$ module load nvidia nvmpi$ mpicc -O3 -acc -Minfo=accel -gpu=cc80 -c main.c-acc: OpenACCコードであることを指示-Minfo=accel:

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

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

nMakefileでコンパイル
講習会のサンプルコードには Makefile がついているので、コンパイルするためには、単純に下記を実行すれば良い。
$ module load nvidia nvmpi$ make31簡単な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.032

簡単な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

n = nx * ny
ny

ny
w

w

MPI

nx

nx

Rank 1, array a

Rank 0, array b

赤要素:
3.0 * ny33

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

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

openacc_mpi_basic/01_original35簡単なOpenACCとMPIコード: CPUコード (2)n 1ノード内の2GPUを利用ü ジョブスクリプトを確認します。
$ cat run.sh#!/bin/bash#PJM -L rscgrp=lecture-a#PJM -L gpu=2#PJM --mpi proc=2#PJM -L elapse=00:10:00#PJM -g gt00
module purge
module load nvidia nvmpi
export UCX_MEMTYPE_CACHE=nmpirun -np 2 ./run

gpu=2:

使用するGPU数を指定。Lecture-aキューは、他のユーザと共有の環境になっており、最大
4GPUまで利用できる。講習会では使えないが、ノード占有のキューではnode=2などと書く。
--mpi proc=2:
MPIのプロセス数を書く
mpirun -np 2-np に全並列数（全MPIプロセス数）を指定。 ./run が実行したいプログラム。

36

MPIランクとGPU割り当てn

OpenACC関数を呼ぶためヘッダーを追加
C: #include <openacc.h>
Fortran: use openaccnノード上のGPU数の取得(acc_get_num_devices)
const int ngpus = acc_get_num_devices(acc_device_nvidia);

acc_device_nvidiaを指定することで、NVIDIA GPUを数える。
nあるプロセスに特定のGPUを割り当て(acc_set_device_num)
int rank= 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
acc_set_device_num(rank % ngpus, acc_device_nvidia);

device_num を rank % ngpus とすることで、ノード内のプロセスは異なるGPUを利用することとなる。これを呼ばないと、0番を使いに行く。
openacc_mpi_basic/02_kernels

37API関数の代わりに環境変数を使う方法
n API関数を使わず、CUDA及びMPI処理系の環境変数を使うことで、GPUとMPIプロセスの割付けを行う方法もある。

run.sh#!/bin/bash#PJM -L rscgrp=lecture-a#PJM -L gpu=4#PJM --mpi proc=4#PJM -L elapse=00:10:00#PJM -g gt00

wrapper.sh（chmodで実行権を与える必要あり）#!/bin/bash
export CUDA_VISIBLE_DEVICES=${OMPI_COMM_WORLD_LOCAL_RANK}$1

module purge
module load nvidia cuda ompi-cuda
mpiexec -machinefile $PJM_O_NODEINF -n 4 ./wrapper.sh ./run

CUDA_VISIBLE_DEVICES： デフォルトでは0。APIによるGPU番号のセットがされていない場合、ここにセットされた
番号のGPUを使う。
OMPI_COMM_WORLD_LOCAL_RANK： MPIの処理系としてOpenMPIを使う場合、ノードローカルなランク番号がここに設定される。つまり上の場合、OMPI_COMM_WORLD_LOCAL_RANK番のGPUを使うことになる。
38

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
39kernels 指示文追加コード（２）n 通信方法ü MPI前後で配列の全領域をデバイスメモリとホストメモリ間でコピー
ü ホストメモリ経由のノード内GPU間通信

n コンパイルと実行$ make$ pjsub ./run.sh$ cat run.sh.??????.out
num of GPUs = 2
mean = 30.00
Time =
0.078 [sec]40

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

41

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 present(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 present(b[0:n])#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += b[i];}}deallocate

42

update指示文追加コード（２）n 通信方法ü MPI前後で配列の必要な領域のみを update 指示文でデバイスメモリとホストメモリ間でコピー
ü 02_kernels と比較し、メモリの確保・解放とデータのコピー量が削減されている。
ü ホストメモリ経由のノード内GPU間通信
n コンパイルと実行$ make$ pjsub ./run.sh$ cat run.sh.??????.out
num of GPUs = 2
mean = 30.00
Time =
0.037 [sec]43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
44

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 present(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 present(b[0:n])#pragma acc loop reduction(+:sum)for (unsigned int i=0; i<n; i++) {sum += b[i];}}deallocate

45

CUDA-aware MPIコード（２）n 通信方法ü MPI_Send, MPI_Recv に対してhost_data 指示文でデバイス側のアドレスを渡している。CUDA-aware MPI であれば正しく動作する。
ü 03_update と比較し、CUDA IPC が利用され、ホストメモリを介さない通信となる。NVLink経由でGPU間通信する。
n コンパイルと実行$ make$ pjsub ./run.sh$ cat run.sh.??????.out
num of GPUs = 2
Rank 0: hostname = wa36, GPU num = 0
Rank 1: hostname = wa36, GPU num = 1
mean = 30.00
Time =
0.024 [sec]46（参考）2ノードの2GPUでの実行
n 各ノード1GPUで2ノードを利用

openacc_mpi_basic/ 05_no_gdr

ü ジョブスクリプトを確認します。(講習会では利用できません）$ cat run.sh#!/bin/bash#PJM -L rscgrp=regular-a#PJM -L node=2#PJM --mpi proc=2#PJM -L elapse=00:10:00#PJM -g gz00
module purge
module load nvidia nvmpi
export UCX_MEMTYPE_CACHE=nexport UCX_IB_GPU_DIRECT_RDMA=n
node=2:

mpiexec -machinefile $PJM_O_NODEINF -n $PJM_MPI_PROC -npernode 1 ./run

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

総MPIプロセス数。この場合は2プロセス。
mpiexec -np -machinefile $PJM_O_NODEINF $PJM_MPI_PROC -npernode 1:-n に全並列数（全MPIプロセス数）を指定。$PJM_MPI_PROCが—mpi procの数字になる。-npernodeで47
ノードあたりのプロセス数を指定する。この場合は１。

（参考）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 を通過する。Wisteria-Aquariusのデフォルトでは、GPU Direct RDMA通信となる。
n コンパイルと実行$ make$ pjsub ./run.sh$ cat run.sh.o??????
num of GPUs = 8
Rank 0: hostname = wa07, GPU num = 0
Rank 1: hostname = wa10, GPU num = 1
mean = 30.00
Time =
0.099 [sec]48（参考）GPUDirectRDMA有効コード（１）n 06_gdr コード

openacc_mpi_basic/ 06_gdr

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

n GPUDirectRDMA の有効化
ü WisteriaのOpenMPI では、デフォルトで GPUDirectRDMA は有効である。これを無効にするには、05_no_gdrのスクリプトに出てきた、UCX_IB_GPU_DIRECT_RDMAをnoにする。
export UCX_IB_GPU_DIRECT_RDMA=n
export UCX_IB_GPU_DIRECT_RDMA=y

#無効の時
#有効の時

n 通信方法ü 05_no_gdr と比較し、 GPUDirectRDMA を有効とする。ホストメモリを経由することな

く、直接データ転送 (RDMA) をする。

49

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

openacc_mpi_basic/ 06_gdr$ make$ pjsub ./run.sh$ cat run.sh.??????.out
num of GPUs = 8
Rank 1: hostname = wa10, GPU num = 1
Rank 0: hostname = wa07, GPU num = 0
mean = 30.00
Time =
0.026 [sec]50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

51

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

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

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

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

54

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

t = n+1

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

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

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

平均後の
自分自身の値

上下左右の値

1

0

0

0
1

0

0

2

0

0
2

8
0

0
2

0

j 3

0
1

8
0

32
20
64

8
0

0
1

4

0

0
2

8
0

0
2

0

5

0

0

0
1

0

0

1

2

3

4

5

i

自分自身の値の４倍

１回目の平均後２回目の平均後最初の状態

繰り返し平均化を行うと、インクが拡散します。

56

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

57

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

袖領域
厚さ mgn =1

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

nz0

袖領域
交換
(MPI)
nz

lnz

nx0
nx(=nx0)

rank2,
GPU2
58

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

59

CPUコードn CPUコードのコンパイルと実行$ cd openacc_mpi_diffusion/01_original$ make$ pjsub ./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コードでは、どのくらいの実行性能が達成できるでしょうか？
60OpenACC化(0): Makefile の修正
n Makefile に OpenACC をコンパイルするよう –acc などを追加しましょうCC= mpicc
CXX = mpic++
GCC = gcc
RM = rm -f
MAKEDEPEND = makedepend
CFLAGS
GFLAGSCXXFLAGSLDFLAGS...

= -O3 -acc -Minfo=accel -gpu=cc80= -Wall -O3 -std=c99= $(CFLAGS)=

61

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

62

OpenACC化(2): kernels, loop
n diffusion3d関数に kernels, loopを追加。present 指定。
#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+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;}

63

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

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

65

NVCOMPILER_ACC_TIME によるOpenACC 実行の確認
n ジョブ実行が終わると、標準エラー出力にメッセージが出力される。
$ cat run.sh.??????.err
Accelerator Kernel Timing data
/work/01/gt00/z30108/openacc_samples_test/openacc_samples/C/openacc_dif
fusion/03_openacc_pgi_acc_time/main.c
main NVIDIA devicenum=0time(us): 725
39: data region reached 2 times
データ移動の回数
39: data copyin transfers: 1
device time(us): total=342 max=342 min=342 avg=342
53: data copyout transfers: 1
device time(us): total=383 max=383 min=383 avg=383/work/01/gt00/z30108/openacc_samples_test/openacc_samples/C/openacc_dif
fusion/03_openacc_pgi_acc_time/diffusion.c
diffusion3d NVIDIA devicenum=0time(us): 0
17: compute region reached 1638 times
25: kernel launched 1638 times
起動したスレッド
grid: [16384] block: [128]elapsed time(us): total=67,084 max=54 min=40 avg=40
17: data region reached 3276 times

カーネル
実行時間

66

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

n NVCOMPILER_ACC_TIME によるOpenACC 実行の確認

OpenACC化の例は、openacc_mpi_diffusion/02_openacc

67

Async節, wait指示文による非同期実行
nAsync(整数)節
Kernels指示文などに付与して使う
n CPUスレッドは通常、kernelsの終了を待ってから次のプログラムに移行するが、Asyncが付いている場合は待たずに先へすすむ
n Asyncに与えられた数値が同じkernelは順番に処理されるが、異なる場合は順不同で処理される
nn

Wait(整数)指示文
Async(整数)と同じ数値を持つタスクの終了を待つ
n (整数)を記述しない場合、全てのタスクの終了を待つnfor(i=...){#pragma acc update device(array) async(i)#pragma acc parallel loop async(i)for(...){...}#pragma acc update host(array) async(i)}#pragma acc wait••

ループカウンタ i ごとに更新と計算を組にして処理を発行終了を待たずに次の i の発行を行い、最後に待ち合わせる

MPIの非同期実行と組み合わせることで、うまくいくと高速化できる

68

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

70

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)

71

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
72

FDTD法（２）n 例えば、

73ソースコード（１）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

データ移動の最適化。

74

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

75

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

76

マルチGPU化のための領域分割
GPU0 Subdomain 0

mgn

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

分割領域の中心領域

mgn

domain分割領域の全体領域

mgn

GPU1

ny

ny0

Subdomain 1

mgn

全領域の中心領域

mgn

mgn

GPU2

mgn

nx0

mgn

ny

Subdomain 2

全領域の全体領域mgn

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

mgn

nx

77
mgn

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

電場 E の時間発展
calc_ex_ey 関数
電場 E の境界条件
pml_boundary_ex 関数
pml_boundary_ey 関数
電場の入射（入射光）
plane_wave_incidence 関数

媒質と物体の設定（ユーザ定義）
set_object_er 関数 (main.cで定義)

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

初期化と計算条件の設定
init_vars, set_initial_condition 関数
init_pml_vars 関数
set_pml_initial_condition 関数
set_pml_rer 関数

磁場 H の時間発展
calc_hz 関数

計算結果出力MPI_Gather収集, write_bmp関数

FDTD計算繰り返し時間発展

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

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

磁場 H の境界条件
pml_boundary_hz 関数
計算結果出力MPI_Gather, write_bmp関78

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

79

mgn

計算領域の設定（２）n Range 構造体

ny
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

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

mgn

nx

mgn

80

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

82

時間発展（１）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;
(後半へ）

83

時間発展（２）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++;(出力など)}

84

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

lnx * (jb - jbw - 1)

= src_hz

Datasize = lnx * 1= nhalomgn

subdomains[i]

ny

mgn

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

(ib, jb)

subdomains[i-1]

ny

mgn

(ibw, jbw)

lnx * (jb + ny - jbw - 1)

= dst_hzmgnn(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);MPI_Recv(&hz[dst_hz], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD, &status);85

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

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

mgn

subdomains[i+1](ib, jb)

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

= dst_ex

subdomains[i]

ny

mgn

(ibw, jbw)

mgnn(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);MPI_Recv(&ex[dst_ex], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD, &status);

86

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

プログラムのコンパイルと実行（１）n CPUコードのコンパイルと実行$ cd openacc_mpi_fdtd/01_original$ make$ pjsub ./run.sh$ cat run.sh.??????.out
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]------------------------------

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

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

88

プログラムのコンパイルと実行（２）n プログラムの実行時オプション$ cat run.sh#! /bin/sh#PJM -L rscgrp=lecture-a#PJM -L gpu=2#PJM –mpi proc=2（省略）

openacc_mpi_fdtd/01_original

mkdir -p sim_run
cd sim_run
mpiexec -machinefile $PJM_O_NODEINF -n $PJM_MPI_PROC ../run 512 512 $PJM_MPI_PROC 5000 50

mpirun -np <nprocs> ../run <nx> <ny> <nprocs> <nt> <nout>$PJM_MPI_PROC: 全ランク数（=分割数）（=--mpi procの値）
nx, ny: 計算領域サイズ
nt: 全時間ステップ
nout: 出力を行うタイムステップ数。50 の場合、50ステップに1回出力する。0 を指定すると出力しない。
89

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

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

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

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

ny0

mgn

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

y
z

mgn

入射光

mgn

x

nx0

mgn

91

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

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

n main.cüüüü

OpenACCヘッダーの追加
GPUの割り当て
data 指示文の追加
MPI_Send/MPI_Recvに対するhost_data指示文の追加

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

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

n run.shü

GPUDirectRDMA の有効化
解答例は、 openacc_mpi_fdtd/02_openacc1

92

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

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

94

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

95

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

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

97

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

98

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

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

n OpenACCコードをさらに最適化しましょう。

NVCOMPILER_ACC_TIMEも活用しましょう。
ü 実は単純に fdtd2d.c に kernels と loop を入れても、いくつかの関数で暗黙の copyin が発生します。これも修正していきましょ$ makeう。
calc_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_openacc499