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

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

講習会スケジュール
n 開催日時ü 12月15日（火） 10:00 – 17:00

n プログラムü 10:00 - 10:50 スパコンの使い方など
ü 11:00 - 11:50 GPUとOpenACC基礎（座学）
ü （昼休み）ü 13:30 – 14:20 OpenACC演習Ⅰ
ü 14:30 – 15:20 OpenACC演習Ⅱ
ü 15:30 – 16:20 OpenACC演習Ⅲ
ü 16:30 – 17:00 質問など

2

講習会についてn 本講習会はü GPUに関する基礎知識
ü OpenACCを用いたGPUプログラミングの基礎を中心に扱います。
n その他の講習会

https://www.cc.u-tokyo.ac.jp/events/lectures/

Youtubeにて過去の講習会を配信中！
https://www.youtube.com/channel/UC
2CHaGp1AO-vqRlV7wmU0w/videos?view=0&sort=p&flow=gridn スパコンイベント情報メール配信サービス

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

講習会の進め方
n Zoomを利用したオンライン講習会ですü この講義は録画されていますü 質問があるとき以外はミュートでお願いしますü ビデオもオフを推奨しますnslackを使って質問に対応しますü slackはリンクを知っている人は誰でも使える設定になっていますü slackのリンクをzoomのチャットに貼るので、未登録の場合は今のうちに登録

お願いしますü slackの登録メールの配送に小一時間かかることがありますü スクリーンショットなどで画像を共有することで、質問対応しますü Windows：Alt + PrtScn で作業中ウィンドウのスクショがクリップボードにコピーされます。slackのチャット部分で貼り付け(Ctrl + V)することで画像をアップロードできますü Mac : command + shift + control + 4 の同時押し、その後撮りたいウィンドウ上でspaceを押すことで、スクリーンショットがクリップボードにコピーされます。slackのチャット部
4分で貼り付け(command + V)することで画像をアップロードできます東大情報基盤センターのスパコン概要5東大センターのスパコン
FY
11

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

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

54.9 TFLOPS, 11.2 TB

Oakforest-PACS
Fujitsu, Intel KNL

T2K Tokyo

140TF, 31.3TB

Big Data &Extreme Computing

25PFLOPS, 919.3TB

Oakleaf-FX: Fujitsu PRIMEHPC
FX10, SPARC64 IXfx

Wisteria/BDEC-01
33 PFLOPS

1.13 PFLOPS, 150 TB

Oakbridge-CXIntel Xeon CLX

Oakbridge-FX

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

6.61 PFLOPS

Reedbush, HPE
Broadwell + Pascal

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

1.93 PFLOPS

Reedbush-L HPE
1.43 PFLOPS

6

3システム：利用者2,000+，学外50+%n Reedbush (SGI, Intel BDW + NVIDIA P100 (Pascal))n データ解析・シミュレーション融合スーパーコンピュータn 3.36 PF, 2016年7月〜 2021年3月末（予定）n 東大ITC初のGPUシステム (2017年3月より), DDN IME (Burst Buffer)n Oakforest-PACS (OFP) (富士通, Intel Xeon Phi (KNL))n JCAHPC (筑波大CCS＆東大ITC)
n 25 PF, TOP 500で6位 (2016年11月) (日本1位)（初登場時）n Omni-Path アーキテクチャ, DDN IME (Burst Buffer)n Oakbridge-CX (富士通, Intel Xeon Platinum 8280)
n 大規模超並列スーパーコンピュータシステム
n 6.61 PF, 2019年7月 〜 2023年6月
n 全1,368ノードの内128ノードにSSDを搭載
7

Reedbushn GPUクラスタ
NVIDIA Tesla P100 SXM2

SGI C1102-GP8

8

Oakforest-PACS(OFP)n Top500: 世界第22位
n HPCG: 世界第16位（2020年11月現在）

Fujitsu
PRIMERGY CX1640 M1

Fujitsu PRIMERGY CX600 M1
シャーシ当たりCX1640 M1 ×8搭載

9

Oakbridge-CX

計算ノード
Chassis：Node:

PRIMERGY CX400 M4 x342 <4node / Chassis>
PRIMERGY CX2550 M5 x1,240, CX2560 M5 x128

x1,368 node

n 世界第68位（2020年11月現在）

Fujitsu PRIMERGY CX2550 M5

全体性能

理論演算性能：主記憶容量：
メモリバンド幅：
ラック数：
SSD搭載：

6.61PF
256.5TiB
385.1TB/s
21ラック
128ノード理論演算性能：4.8384 TF
手記憶容量： 192GiB
メモリバンド幅： 281.6GB/s

計算ノード間ネットワーク (Omni-Path Architecture)
通信性能 100Gbps
管理サーバ群

ログインノードx10

Fujitsu PRIMERGY CX400 M1
シャーシ当たりCX2550 M5×4搭載

ノード単体

FUJITSU Server
PRIMERGY CX2560 M5 x 10

並列ファイルシステム

x15
FUJITSU Server
PRIMERGY RX2530 M4 x 15
(ジョブ、運用、認証、Web、
セキュリティログ保存)

ストレージ装置

12.4PB

ストレージ装置：DDN ES18KE x2セット
ファイルシステム：DDN ExaScaler
10(Lustreベースファイルシステム)

Wisteria/BDEC-01
n 「シミュレーション・データ・

学習（S+D+L）」融合，2種類のノード群

詳細はQRコードまたは下記のurlよりhttps://www.cc.utokyo.ac.jp/public/pr/pr-wisteria.php

n シミュレーションノード群（Odyssey）
n 従来のスパコン
n A64FX（富岳と同じ） CPU x 7,680 noden 25.9PF, 7.8PB/secn データ・学習ノード群（Aquarius）
n データ解析，機械学習n (Intel Xeon CPU x2 + NVIDIA A100 GPU x8)x 45 node
n 7.2PF, 578.2TB/secn データ・学習ノード群の一部は外部リソース（ストレージ，サーバー，センサーネット
ワーク他）に直接接続

n 2021年5月稼働開始予定

11

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

12スパコン料金表（2020年4月時点 ）n 最小セット料金表
トークン

料金

ストレージ容量

利用期間

Reedbush

720*

6,300円

1TB

年度末まで

Oakforest-PACS

720

4,200円

1TB

年度末まで

Oakbride-CX

720

8,400円

4TB

年度末まで

n 「トークン制」で運営
n トークン≒ノード時間。720トークンなら、1ノードを720時間利用できる。
n ただし、Reedbush-H,Lにはそれぞれ2.5, 4.0の消費係数が設定されていて、消費係数倍のトークンを消費する。つまり、Hの場合は720/2.5=288時間, Lの場合は720/4=180時間しか利用できないn トークン、容量共に料金積み増しで増量可能。詳しくは以下

http://www.cc.u-tokyo.ac.jp/supercomputer/ofp/service/
http://www.cc.u-tokyo.ac.jp/supercomputer/reedbush/service/
http://www.cc.u-tokyo.ac.jp/supercomputer/obcx/service/13

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

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

14

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

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

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

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

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

15

Reedbushでのプログラムの実行
n ジョブスクリプト(〇〇.sh)を作成し、ジョブとして投入、実行する。
$ qsub ./〇〇.shn 投入されたジョブを確認する。（qstatではないので注意）$ rbstatn 実行が終了すると、以下のファイルが生成される。
〇〇.sh.o??????
〇〇.sh.e??????

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

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

GPUプログラミングを始める前に！
1. 並列プログラミングって？

17

GPUプログラミングを始める前に！
n GPUは並列計算機です！よって本講習会で学ぶのは並列プログラミングになります！
n 並列プログラミングの例： MPI, OpenMP など

n 並列プログラミングは、プログラムを高速化するために行います！

並列プログラミングについての解説動画はこちら
https://www.youtube.com/channel/UC2CHaGp
1AO-vqRlV7wmU0w/videos?view=0&sort=p&flow=grid

並列プログラミング・
高性能計算についての事前知識があると
有利！
18

並列計算
n 実行時間 T の逐次処理のプログラムを p 台の計算機で並列計算することで、実行時間を T / p にする。
Tp台
T/pn 実際にできるかどうかは、処理内容（アルゴリズム）による。アルゴリズ

ムによって難易度は異なる。
ü 並列化できないアルゴリズム、通信のオーバーヘッド
ü 部分的にでも並列化できないアルゴリズムがあると、どれだけ並列数を上げても、その時間は短縮されない。
n 並列処理（計算）の種類
ü 「タスク並列」と「データ並列」
19

タスク並列n タスク（仕事）を分割することで並列化する。

GPUは苦手n タスク並列の例：カレーを作る
ü 仕事１：野菜を切る
ü 仕事２：肉を切る
ü 仕事３：水を沸騰させる
ü 仕事４：野菜と肉を入れて煮込む
ü 仕事５：カレーのルウを入れる

n 並列化
仕事1
仕事2
仕事3

仕事4

仕事5

時間
20

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

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

2+1=
3 + 19 =
4 + (-6) =-8 + 10 =10 + 3 =12 + (–88)=-20 + 29 =4 + 10 =-32 + 12 =-5 + 5 =GPUの並列計算はこれが原則。
プログラムでは普通、配列とループで記述する
for (i = 0; i < N; i++) C[i] = A[i] + B[i];21GPUにおけるループ並列化
n GPU での高速化は通常、プログラ

ム中の重たいループ構造を並列化することで達成する
n 今回学ぶOpenACC は特定のループ構造を簡単に並列化できる
n 全てのループ構造を並列化できるわ

けではないn どのようなループなら並列化可能か

世の中のすべてのループ

GPUで(≒CUDAで)高速化可能なループ

OpenACCで高速化可能なループ

知る必要がある22

OpenACCで並列化できるループ
データ独立なループの例

リダクションの必要なループの例

データ依存のあるルー
プの例

for ( i = 0; i < N; i++)C[i] = A[i] + B[i];sum = 0;for ( i = 0; i < N; i++)sum += A[i];B[0] = 0;
for ( i = 0; i < N; i++)B[i+1] = B[i] + A[i];OpenACCで簡単に並列化できる

CUDAでも比較的簡
単に並列化できる※OpenACCでも、GPUで正しく動く
コードを書くことはできる。しかし遅いので意味がない

CUDAで高速実装可能だが難しい (shared memory や
warp shuffle を駆使する必要がある)

23

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

A[0] = A[0] + 1;

A[1] = A[1] + 1;

A[2] = A[2] + 1;

分担で計算を行う
スレッドさん

メモリ

3

6

1

A[0] A[1] A[2]24

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

A[0] = A[0] + 1;

A[1] = A[1] + 1;

A[0]て
なんや

メモリ

A[1]て
なんや

3

6

A[2] = A[2] + 1;
A[2]て
なんや

1

A[0] A[1] A[2]25

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

A[0] = A[0] + 1;

メモリ

A[1] = A[1] + 1;

3

6

A[2] = A[2] + 1;

1

A[0] A[1] A[2]26

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

A[0] = A[0] + 1;

A[1] = A[1] + 1;

６やな

A[2] = A[2] + 1;

１やて

３や

メモリ

3

6

1

A[0] A[1] A[2]27

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

A[0] = A[0] + 1;

A[1] = A[1] + 1;

６

A[2] = A[2] + 1;

１

３

メモリ

3

6

1

A[0] A[1] A[2]28

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

３ + 1;
A[0] = A[0]メモリ

６ + 1;
A[1] = A[1]3

6

１ + 1;
A[2] = A[2]1

A[0] A[1] A[2]29

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

３ + 1;
A[0] = A[0]４や

６ + 1;
A[1] = A[1]

２やわ

７やで

メモリ

3

１ + 1;
A[2] = A[2]61

A[0] A[1] A[2]30

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

３ + 1;
A[0] = A[0]６ + 1;
A[1] = A[1]
７

２

４

メモリ

3

１ + 1;
A[2] = A[2]61

A[0] A[1] A[2]31

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

３ + 1;
A[0] = A[0]メモリ

６ + 1;
A[1] = A[1]47

１ + 1;
A[2] = A[2]2

A[0] A[1] A[2]32

簡単に並列化できるループ for ( i = 0; i < 3; i++)A[i] = A[i] + 1;

３ + 1;
A[0] = A[0]６ + 1;
A[1] = A[1]このようなデータ並列を簡単に適用できるループを、
n データ独立(independent)なループ
n 依存性のないループn 自明な並列性を持つループ
などと呼ぶ

メモリ

１ + 1;
A[2] = A[2]

成功！

4

7

2

A[0] A[1] A[2]33

簡単に並列化できないループ
A[0] = A[0] + 1;

for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0]に3回1を足してるだけなので、最終結果は 3 + 1 + 1 + 1 = 6。
足し算なのでどんな順番で足しても結果は変わらないはずだが…

メモリ

3

6

1

A[0] A[1] A[2]34

簡単に並列化できないループ
A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0]て
なんや

メモリ

for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

A[0]て
なんや

3

6

A[0] = A[0] + 1;
A[0]て
なんや

1

A[0] A[1] A[2]35

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0] = A[0] + 1;
少し休ん
でからで
ええか

メモリ

3

6

1

A[0] A[1] A[2]36

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0] = A[0] + 1;
少し休ん
でからで
ええか

３やな
３や

メモリ

3

6

1

A[0] A[1] A[2]37

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0] = A[0] + 1;

A[0] = A[0] + 1;

３
３

メモリ

3

6

1

A[0] A[1] A[2]38

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]メモリ

３ + 1;
A[0] = A[1]3

6

A[0] = A[0] + 1;

1

A[0] A[1] A[2]39

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]４や３ + 1;
A[0] = A[1]A[0] = A[0] + 1;

４やで

メモリ

3

6

1

A[0] A[1] A[2]40

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]A[0] = A[0] + 1;

４
そろそろ
行くか

４

メモリ

3

6

1

A[0] A[1] A[2]41

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]A[0] = A[0] + 1;

４やで
４やろ？

メモリ

4

6

1

A[0] A[1] A[2]42

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]A[0] = A[0] + 1;

４やで
4やて
４やろ？

メモリ

4

6

1

A[0] A[1] A[2]43

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]A[0] = A[0] + 1;

４やで
4
４やろ？

メモリ

4

6

1

A[0] A[1] A[2]44

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]4 + 1;
A[0] = A[0]

４やで
４やろ？

メモリ

4

6

1

A[0] A[1] A[2]45

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]4 + 1;
A[0] = A[0]
5やわ

４やで
４やろ？

メモリ

4

6

1

A[0] A[1] A[2]46

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]4 + 1;
A[0] = A[0]
5

４やで
４やろ？

メモリ

4

6

1

A[0] A[1] A[2]47

簡単に並列化できないループ for ( i = 0; i < 3; i++)A[0] = A[0] + 1;

３ + 1;
A[0] = A[0]３ + 1;
A[0] = A[1]

※6です

メモリ

CPUで実行される足し算は、
1. データの読み込み
2. 足し算3. データの書き込みの3パートからなる。

４や！
5やな！

は？

4

6

３ + 1;
A[0] = A[2]1

スレッドは各々独立に１〜３を実行するため、タイミング
によって結果が変わる！
（この例の場合は4,5,6のいず
れかになる）

A[0] A[1] A[2]48

どうやって並列化するか？

配列A

1

2

3

4

5

6

7

8

9

例えば以下を8スレッドで並列化

sum = 0;for ( i = 0; i < 16; i++)sum = sum + A[i];10 11 12 13 14 15 16

配列B

49

どうやって並列化するか？

配列A

1

配列B

2

13

1.

3

4

37

5

6

11
5

7

8

15
7

9

例えば以下を8スレッドで並列化

sum = 0;for ( i = 0; i < 16; i++)sum = sum + A[i];10 11 12 13 14 15 16

19
9

27
13

31
15

各々自分の担当領域で足し算（結果を別の場所に保存）

50

どうやって並列化するか？

配列A

1

配列B

2

3

1.
2.

3

4

7

5

6

11

7

8

15

9

例えば以下を8スレッドで並列化

sum = 0;for ( i = 0; i < 16; i++)sum = sum + A[i];10 11 12 13 14 15 16

19

23

27

31

各々自分の担当領域で足し算（結果を別の場所に保存）遅れているスレッドを待つ！（これを同期(thread synchronization)という）

51

どうやって並列化するか？

配列A

1

配列B

2

3

3

4

7

5

6

11

7

8

15

9

例えば以下を8スレッドで並列化

sum = 0;for ( i = 0; i < 16; i++)sum = sum + A[i];10 11 12 13 14 15 16

19

23

27

31

配列C
1.
2.
3.

各々自分の担当領域で足し算（結果を別の場所に保存）遅れているスレッドを待つ！（これを同期(thread synchronization)という）
一部のスレッドを寝かせて、起きてるスレッドで(1)から繰り返し
52

どうやって並列化するか？

配列A

1

配列B
配列C

2

3

3

4
7

6

7

11

10

配列D

5

8
15

26
36

9

例えば以下を8スレッドで並列化

sum = 0;for ( i = 0; i < 16; i++)sum = sum + A[i];10 11 12 13 14 15 16
19

23

27

42

31
58

100

sum 136nnn

これは一般的にリダクションと呼ばれる演算パターン
一番働くスレッドが4回の足し算。逐次の場合と比較して4倍の高速化
メモリなどを介してスレッドの間でデータのやり取りをすることをスレッド間通信というスレッドの同期・通信が入ると途端に難しくなる！

53

GPU入門54

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

3D GameComputer Graphics

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

55

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

56

抑えておくべきGPUの特徴
n 最低限知っておくべきこと
ü 超並列計算が必須！
n 物理コア数が1000以上、論理コア数（スレッド）は数十万以上
n プログラムの並列性（スレッド分割可能数）が小さいと速くならないü CPU と GPUの間でのデータ転送が必須！
n GPU は CPU の指示なしでは動けない
n CPU と GPU は独立に動く
ü CPUとGPUの同期を行い、データの一貫性を保つ必要がある

n さらなる高速化のためには
n 階層的スレッド管理と同期・通信n Warp 単位の実行

これらはプログラミング言語が CUDA
か OpenACCに関わらず、GPUプログ
ラミングでは考慮する必要がある。

n コアレスドアクセス
57

NVIDIA Tesla P100

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

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

バス (PCIe など)
CPU

GPU

~32GB/s

OSが動いている

OSは存在しない

２．計算を行う

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

メインメモリn

計算はOSのあるCPUから始まる

n

物理的に独立のデバイスメモリとデータのやり取り必須

デバイス
メモリ

物理的に独立59

どんなアプリならGPUで高速化できる？
n 原則：GPUに一度送ったデータを使い回せるアプリケーションn 最低でも100回は使いまわしたい
n 例：データ量 N に対して計算量 O(N2) 以上の計算（行列積、多体問題など）や、反復法など
n 思考実験
n 次のプログラムを、右の表のコンピュータの (1) CPUを使った時 (2)GPUを使った時の実行時間は？
if(GPU) BをCPUからGPUにコピー
do i = 1, N
A(i) = B(i)end doif(GPU) AをGPUからCPUにコピーあるコンピュータの性能
データ転送性能

CPUのメモリ

65 GB/sec

GPUのメモリ

540 GB/sec

CPU-GPU間のバス

10 GB/sec(1) 2 * 8 * N / 65(2) 2 * 8 * N / 540+ 2 * 8 * N / 10
N = 109 (1G) なら？

(1) 0.246 sec (2) 1.629 sec

60

どんなアプリならGPUで高速化できる？
n 原則：GPUに一度送ったデータを使い回せるアプリケーションn 最低でも100回は使いまわしたい
n 例：データ量 N に対して計算量 O(N2) 以上の計算（行列積、多体問題など）や、反復法など
n 思考実験
n 次のプログラムを、右の表のコンピュータの (1) CPUを使った時 (2)GPUを使った時の実行時間は？
if(GPU) BをCPUからGPUにコピー
do time_step = 1, 100
do i = 1, N
A(i) = B(i)end doend doif(GPU) AをGPUからCPUにコピーあるコンピュータの性能
データ転送性能

CPUのメモリ

65 GB/sec

GPUのメモリ

540 GB/sec

CPU-GPU間のバス

10 GB/sec(1) 2 * 8 * N * 100 / 65(2) 2 * 8 * N * 100 / 540+ 2 * 8 * N / 10
N = 109 (1G) なら？

(1) 24.6 sec (2) 4.56 sec

61

OPENACC入門
62

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

簡単

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

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

難しいü プログラミングにはGPGPU用言語を使用する必要あり。
63

OpenACCn OpenACCとはü 新しいアクセラレータ（≒GPU）用プログラミ

ングインターフェースü OpenMP のようなディレクティブ（指示文）・
ベース
ü C 言語/C++, Fortran に対応
ü 2011年秋に OpenACC1.0、最新は 3.1
ü コンパイラ：［商用］PGI, Cray,［フリー］GCC（PGIは無料版あり）
ü WEBサイト：http://www.openacc.org/
n 指示文ベースの利点

C言語#pragma acc directive-name [clause, …]{// C code}

Fortran!$acc directive-name [clause, …]! Fortran code!$acc end directive-nameü 指示文：コンパイラへのヒントü アプリケーションの開発や移植が比較的簡単
ü ホスト（CPU）用コード、複数のアクセラレー

タ用コードを単一コードとして記述。メンテナ
ンスが容易。高生産性。

64

OpenACCでできること
n OpenACC は特定のループ構造を簡単に並列化できるü 全てのループ構造を並列化できるわけではない

世の中のすべてのループ

GPUで(≒CUDAで)高速化可能なループn 主に以下の３つを記述できる
ü どこでデータを移動するか

OpenACCで高速化可能なループü どこを GPU で実行するか
ü （GPUで実行する領域ないに出てくる）

ループが、データ独立か、リダクション
か、それ以外か65OpenACCでできないこと
n CUDAならshared memoryなど使って頑張れば並列化できる、データ依存性のあるループの並列化
n 例外：atomic演算で解決可能な書き込み競合を含むループ

n shared memoryなど使った性能限界を目指す最適化

こういうのはアプリの一部分であることが多いので、ここだけCUDAやライブラリを使えば良い。

OpenACC と CUDAやライブラリの併用など、上級者は楽するためにOpenACCを使う

OpenACCを推奨する理由
n CPUプログラムの一般的なGPU化手順
1. プログラムのプロファイリング（重い部分を特定する）
2. 重い部分を並列化し、GPU上で実行する

67

OpenACCを推奨する理由
CPU
main

GPU

実行時間
12

subA

10

subA

8
6

subB
subC

subB

4
2

subC

0
CPU

とあるプログラムの構造と実行時間を調べた結果

68

OpenACCを推奨する理由
CPU

GPU

main

実行時間
12

subA

10
8

subA
subB
subC

subB だけ
CUDA化したら良さそう！

6

subB

4
2

subC

0
CPU

69

OpenACCを推奨する理由
CPU

GPU

実行時間

main

12
10
8

subA
subB

データ転送
データ転送

6

subB

4
2

subC
アレ？！

0
CPU

CPU+GPU

データ転送分遅くなった！

70

OpenACCを推奨する理由
n CPUプログラムの一般的なGPU化手順
1. プログラムのプロファイリング（重い部分を特定する）
2. 重い部分を並列化し、GPU上で実行する
3. CPU-GPU間のデータ転送を最小化する

OpenACCであってもCUDAであっても、
結局ここまでが必須！

71

OpenACCを推奨する理由
CPU

GPU

main

subA

12

subAの結果をsubBに、subBの結果をsubCに使っている…
データ転送
データ転送

subC

実行時間
10
8
6

subB

4
2

データ転送をループの外に追い出すた 0
めには…

CPU

CPU+GPU72

OpenACCを推奨する理由
CPU

GPU

main

実行時間
12

データ転送

10

subA

8
6

subB
subC

4
2
0

データ転送

CPU

CPU+GPU

結局全部CUDA化した…

73

OpenACCを推奨する理由
n CPUプログラムの一般的なGPU化手順
1. プログラムのプロファイリング（重い部分を特定する）
2. 重い部分を並列化し、GPU上で実行する
3. CPU-GPU間のデータ転送を最小化する
4. GPU実行部でなお重い場所を最適化する

1,2,3をOpenACCで実装することで、最低限の実装までの工数を減らす。
4の最適化をCUDAで行う。OpenACCにはCUDAと組み合わせるためのインターフェースが用意されている。

74

OpenACCを推奨する理由
n 実アプリをGPU化する場合、データ転送を最小化するためには、

結局大部分をGPU化する必要がある
n しかし実アプリ全体をCUDA化するのは非常に工数が掛かるため、
まずはOpenACCで全体をGPU化する
n この時点で性能が十分であれば、GPU化を終了する

n OpenACCで並列化できないループや、OpenACCでは性能が十分ではないループに関して、CUDA化を行う
n 多くの場合このようなループは、アプリケーションの一部に限られる

以上により、CUDA化と遜色ない性能を少ない工数で達成できる75

OpenACC と CUDA の組み合わせn host_data指示文を使う：data指示文でCPU・GPUでペアで確保されたデータの、GPU側のアドレスをゲットできる → 後はやりた
い放題
n GPU側のアドレスを使いたい例
ü 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
76OpenACCの実行イメージ
1スレッド

OpenMP

OpenACC・CUDA

CPU

デバイス
(GPU)

int main(){...
#このループを並列実行
for (i=0; i<n; i++) {
...
}
...
}

CPU

CPU

77

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

CPU

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

a b#pragma acc data copyin(a[0:n]), copyout(b[0:n])#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;
}

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

解放
double sum = 0;
for (int i=0; i<n; i++) {sum += b[i];}fprintf(stdout, "%f¥n", sum/n);free(a); free(b);return 0;}

GPUから
copyout78

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

CPU

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

a b#pragma acc data copyin(a[0:n]), copyout(b[0:n])#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {b[i] = a[i] + c;
}

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

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

GPUから
copyout

79

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

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

80

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

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という仕様になっているが、混乱の元なので、一つのループネストに一つのkernels 指示文を推奨！

kernel 2}}

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

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

kernel 1という仕様になっているが、混乱の元なので、一つのループネストに一つのkernels 指示文を推奨！

kernel 2

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

C

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 化
ü GPUで実行したいループを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;
}

84

C

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 化
ü GPUで実行したいループを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;
}

85

C

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

ü GPUで実行したいループをkernels

で囲むü 加えて、ループが並列化可能かどうかを示すために、loop 指示
文を追加
直後のループがデータ独立(independentなルー
プであることを表している

kernels直後の{}で囲まれる領域がGPU上で実行される

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

n ループのOpenACC 化

CPU-GPU間のデータ転送は、一定の規則に基づき
ベストエフォートで行われる

86

F

CPUコードのOpenACC化
openacc_hello/01_hello_acc
program main
implicit none
! 変数宣言allocate(a(n),b(n))
c = 2.0
do i = 1, na(i) = 10.0
end do!$acc kernels!$acc loop independent
do i = 1, nb(i) = a(i) + c
end do!$acc end kernels
sum = 0.d0
do i = 1, n
sum = sum + b(i)end doprint *, sum/ndeallocate(a,b)
end program mainn ループのOpenACC 化
ü GPUで実行したいループをkernels

で囲むü 加えて、ループが並列化可能かどうかを示すために、loop 指示
文を追加

Fortranの場合は、kernels ~ end kernelsの間がGPUで実行される
CPU-GPU間のデータ転送は、一定の規則に基づき
ベストエフォートで行われる

87

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

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

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

OpenACCのデフォルトの変数の扱い
n スカラ変数
ü firstprivate または private

つまりは気にしなくていいü ホストからデバイスへコピーが渡され初期化。ホストに戻せない。

n 配列
ü shared

配列をどう扱うかが重要ü デバイスメモリに動的に確保され、スレッド間で共有。
ü デバイスからホストへコピーすることが可能。

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

C

データ管理・移動
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;}#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;}

openacc_hello/01_hello_acc

直後の { のタイミングで、mallocとCPU -> GPUのデータコピーが行われる
直後の } の終了タイミングで、GPU -> CPUのデータコピーとFreeが行われる

90

F

データ管理・移動
n data 指示文

ü デバイス(GPU)メモリの確保と解放、ホスト(CPU)とデバイス(GPU)間のデータ転送を制御

kernels指示文では、データ転送は自動的に行われる。data指示文でこれを制御することで、不要な転送を避け、性能向上できるü CUDA で言うところの cudaMalloc, cudaMemcpy に相当
program main
implicit none
integer,parameter :: n = 1000real(KIND=4),allocatable,dimension(:) :: a,breal(KIND=4) :: cinteger :: ireal(KIND=8) :: sumallocate(a(n),b(n))
c = 2.0
do i = 1, na(i) = 10.0
end do!$acc data copyin(a) copyout(b)!$acc kernels!$acc loop independent
do i = 1, nb(i) = a(i) + c
end do!$acc end kernels!$acc end data

openacc_hello/01_hello_acc

Fortranでは配列サイズ情報が変数に付随するため(lbound,ubound,sizeな
どの組み込み関数をサポートしている)、基本的にサイズを書く必要がない。

91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 は既にデバイス上確保されているデータに対しては何もしない。
92
present として振る舞う。（OpenACC2.5以降）

data 指示文の指示節
xxxの選択肢はcopy
copyincopyout
createpresent#pragma acc data XXX(a[0:N]){/* C コード */}if(配列aのペア、a_GPUがGPU上にまだない) {if(XXX == copy, copyin, copyout, create){
a_GPU を GPU上に確保}
if(XXX == copy, copyin){a_GPU[0:N] = a[0:N];}
if(XXX == present){
print(エラー! a はGPU上にありません！);}}{/* C コード */}if(上のif文がtrueだった時) {if(XXX == copy, copyout){a[0:N] = a_GPU[0:N];}
if(XXX == copy, copyin, copyout, create){
free(a_GPU);}93}データの移動範囲の指定
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])
...
94

C

ループ指示文による並列化
openacc_hello/01_hello_acc

CPU

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

a b#pragma acc data copyin(a[0:n]), copyout(b[0:n])#pragma acc kernels#pragma acc loop independent
loop指示文for (int i=0; i<n; i++) {b[i] = a[i] + c;
}

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

解放
double sum = 0;
for (int i=0; i<n; i++) {sum += b[i];}fprintf(stdout, "%f¥n", sum/n);free(a); free(b);return 0;}

GPUから
copyout95

F

ループ指示文による並列化
openacc_hello/01_hello_acc

CPUprogram main
implicit none
! 変数宣言allocate(a(n),b(n))
c = 2.0a b

do i = 1, na(i) = 10.0
end do!$acc data copyin(a) copyout(b)!$acc kernels!$acc loop independent
do i = 1, nb(i) = a(i) + c
end do!$acc end kernels!$acc end data
sum = 0.d0
do i = 1, n
sum = sum + b(i)end doprint *, sum/ndeallocate(a,b)
end program main

GPU

メモリ確保

GPUへ
copyin a

b

メモリ確保

カーネル
実行

loop指示文解放GPUから
copyout

96

loop 指示文
n ループの種類を指定
ü データ独立なループ(independent)ü リダクションループ (reduction)
ü 並列化すべきでないループ (seq)n ループマッピングのパラメータの調整ü 難しいので、最初は考える必要はない
ü gang, worker, vector を用いて指定する。大まかに以下のように考えると良い。
- gang: CUDA の thread block 数の指定
- vector: CUDA の block 内の threads 数の指定
97

loop 指示文
n ループの種類を指定
ü データ独立なループ(independent)ü リダクションループ (reduction)

これを見極めるのはプログラマの役目！

ü 並列化すべきでないループ (seq)n ループマッピングのパラメータの調整ü 難しいので、最初は考える必要はない
ü gang, worker, vector を用いて指定する。大まかに以下のように考えると良い。
- gang: CUDA の thread block 数の指定
- vector: CUDA の block 内の threads 数の指定
98

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

ü

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

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

99

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

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

}
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(){...

ü kernels 指示文でGPUでの実行領域を指定。
ü loop 指示文でループの並列化
ü data 指示文でデータ転送を制御。

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

100

OpenACCコードのコンパイル
n

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

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

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

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

簡単なOpenACCコードn サンプルコード: openacc_basic/
ü OpenACC指示文 kernels, data, loop を利用したコード
ü 計算内容は簡単な四則演算

CF

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];}}do j = 1,ny
do i = 1,nxc(i,j) = c(i,j) + a(i,j) + b(i,j)
end doend doü ソースコードopenacc_basic/01_original

CPUコード。

openacc_basic/02_kernels

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

openacc_basic/03_kernels_copy

OpenACCコード。上にcopy指示節追加。

openacc_basic/04_loop

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

openacc_basic/05_data

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

openacc_basic/06_present

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

openacc_basic/07_reduction

OpenACCコード。上にreduction指示節を使用。

102

C

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

103

配列のインデックス計算

Fn サンプルコード: openacc_basic/
ü OpenACC指示文 kernels, data, loop を利用したコード
ü 計算内容は簡単な四則演算
subroutine calc(nx, ny, a, b, c)
implicit none
integer,intent(in) :: nx,nyreal(KIND=4),dimension(:,:),intent(in) :: a,breal(KIND=4),dimension(:,:),intent(out) :: cinteger :: i,j
do j = 1,ny
do i = 1,nxc(i,j) = a(i,j) + b(i,j)
end doend doend subroutine calc

Fortran版では多次元配列を利用

104

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

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

105簡単なOpenACC: kernels 指示文（１）CF

n 02_kernelsコード: calc関数
ü CPUコードにkernels 指示文の追加

openacc_basic/02_kernels

Cvoid calc(unsigned int nx, unsigned int ny, const float *a, const float *b, float *c){const unsigned int n = nx * ny;#pragma acc kernelsfor (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];
}
}
}

F

subroutine calc(nx, ny, a, b, c)
implicit none
integer,intent(in) :: nx,nyreal(KIND=4),dimension(:,:),intent(in) :: a,breal(KIND=4),dimension(:,:),intent(out) :: cinteger :: i,j!$acc kernels
do j = 1,ny
do i = 1,nx
OpenACC コンパイラは配列 (a, b, c)c(i,j) = a(i,j) + b(i,j)
を shared 変数として自動で転送しend doてくれるはずだが…
end do!$acc end kernels
end subroutine calc106

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

データサイズがわからずコンパイルエラー
C言語では配列サイズの指定がほぼ必須！

CFC$ makepgcc -O3 -acc -Minfo=accel -ta=tesla,cc60 -c main.c
PGC-S-0155-Compiler failed to translate accelerator region (see -Minfo messages): Could
not find allocated-variable index for symbol (main.c: 13)
calc:
14, Complex loop carried dependence of a->,c->,b-> prevents parallelization
Accelerator serial kernel generated
Accelerator kernel generated
Generating Tesla code
14, #pragma acc loop seq
15, #pragma acc loop seq
15, Accelerator restriction: size of the GPU copy of c,b,a is unknown
Complex loop carried dependence of a->,c->,b-> prevents parallelization
PGC-F-0704-Compilation aborted due to previous errors. (main.c)
PGC/x86-64 Linux 18.4-0: compilation aborted
make: *** [main.o] Error 2F$ makepgfortran -O3 -mp -acc -ta=tesla,cc60 -Minfo=accel -c main.f90
calc:
13, Generating implicit copyin(b(:nx,:ny))
Generating implicit copyout(c(:nx,:ny))
Generating implicit copyin(a(:nx,:ny))
14, Loop is parallelizable
15, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
14, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
15, !$acc loop gang, vector(32) ! blockidx%x threadidx%x

データサイズを検知して自動転送
Fotranではサイズ情報が配列に付
随するため
107

簡単なOpenACC: kernels 指示文（３）Cn 03_kernels_copyコード: calc関数
ü 配列サイズを明示的に指定

openacc_basic/03_kernels_copy

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

108

簡単なOpenACC: kernels 指示文（４）CF

n 03_kernels_copyコード:初期化
ü CPUコードにkernels 指示文の追加

openacc_basic/03_kernels_copy

CF

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

program main...
!$acc kernels copyout(b,c)
do j = 1,ny
do i = 1,nxb(i,j) = b0
end doend doc(:,:) = 0.0!$acc end kernels...
end program...
}

Fortranの配列代入形式も使える
b0 はスカラ変数のため自動的に各スレッドへコピーが渡される。
109

簡単なOpenACC: kernels 指示文（５）Cn コンパイル
ü データの独立性がコンパイラにはわからず、並列化されない。
$ makepgcc -O3 -acc -Minfo=accel -ta=tesla,cc60 -c main.c
calc:
13, Generating copy(a[:n],c[:n],b[:n])14, Complex loop carried dependence of a-> prevents parallelization
Loop carried dependence due to exposed use of c[:n] prevents parallelization
Complex loop carried dependence of c->,b-> prevents parallelization
Accelerator scalar kernel generated
Accelerator kernel generated
Generating Tesla code
14, #pragma acc loop seq
15, #pragma acc loop seq
15, Complex loop carried dependence of a->,c->,b-> prevents parallelization
Loop carried dependence due to exposed use of c[:i1+n] prevents parallelization
main:
43, Generating copyout(c[:n],b[:n])45, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
45, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */48, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
48, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */pgcc -O3 -acc -Minfo=accel -ta=tesla,cc60
main.o-o run110簡単なOpenACC: kernels 指示文（５）Fn コンパイル
ü データの独立性を見切り、並列化。
pgfortran -O3 -mp -acc -ta=tesla,cc60 -Minfo=accel -c main.f90
calc:
13, Generating copyin(a(:,:))
Generating copyout(c(:,:))
Generating copyin(b(:,:))
14, Loop is parallelizable
15, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
14, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
15, !$acc loop gang, vector(32) ! blockidx%x threadidx%x
main:
61, Generating copyout(c(:,:),b(:,:))
62, Loop is parallelizable
63, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
62, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
63, !$acc loop gang, vector(32) ! blockidx%x threadidx%x
68, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
68, !$acc loop gang, vector(4) ! blockidx%y threadidx%y!$acc loop gang, vector(32) ! blockidx%x threadidx%x
pgfortran -O3 -mp -acc -ta=tesla,cc60 -Minfo=accel main.o -o run
111

Tips: なぜデータの独立性を見切れないか
n エイリアス （変数の別名）
n 主にポインタの利用
n 右は一見データ独立でも…
n foo(&a[0],&a[1]) のような呼び出し

これってデータ独立？
void foo(float *a, float *b){for (int i=0; i<N; i++)b[i] = a[i];}をすればデータ独立でない！

n 不明瞭な書き込み参照先
n インデックス計算
n 計算結果がループ変数に対して独立かどうかわからない
n Fortranでも、多次元配列を一次元化すると起こるn 逆にCでも、多次元配列を使えば独立性を見切れる

インデックス計算
for (int i=0; i<N; i++){j = i % 10;
b[j] = a[i];}間接参照
for (int i=0; i<N; i++){b[idx[i]] = a[i];}n 間接参照112

C

簡単なOpenACC: loop 指示文（１）n 04_loopコード
ü 03_kernelsコードにloop independent の追加

openacc_basic/04_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;
}
}

113

簡単なOpenACC: loop 指示文（１）F

n 04_loopコード
ü 03_kernelsコードにloop independent の追加

openacc_basic/04_loop

subroutine calc(nx, ny, a, b, c)...
!$acc kernels copyin(a,b) copyout(c)!$acc loop independent
do j = 1,ny!$acc loop independent
do i = 1,nxc(i,j) = a(i,j) + b(i,j)
end doend do!$acc end kernels
end subroutine! main 関数内!$acc kernels copyout(b,c)!$acc loop independent
do j = 1,ny!$acc loop independent
do i = 1,nxb(i,j) = b0
end doend doc(:,:) = 0.0!$acc end kernels

各次元についてloop指示文を指定する（並列サイズなどを指定したいなど）場合、do文で書き下す必要がある。

114

簡単なOpenACC: loop 指示文（２）C

openacc_basic/04_loop

n コンパイルü ループが並列化され、カーネルが生成された。

$ 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 */※Fortran版は既に並列化されていたため省略。loop
independent をつける事による挙動の変化はない。（少なくとも
PGI compiler ver. 18.7 では）

115

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

openacc_basic/04_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
}

116

C

簡単なOpenACC: data指示文（１）n 05_dataコード
openacc_basic/05_data

ü 04_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 は利用用途に合わせた指示節を指定。
117

F

簡単なOpenACC: data指示文（１）n 05_dataコード
ü 04_loopにdata指示文追加! main関数内!$acc data copyin(a) create(b) copyout(c)!$acc kernels copyout(b,c)!$acc loop independent
do j = 1,ny!$acc loop independent
do i = 1,nxb(i,j) = b0
end doend doopenacc_basic/05_dataa: allocate, H -> D
b: allocate
c: allocate

present として振舞う。

c(:,:) = 0.0!$acc end kernels
do icnt = 1,nt
call calc(nx, ny, a, b, c)
end do!$acc end data

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

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

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

openacc_basic/05_data$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
1.174 [sec]119簡単なOpenACC: present指示節
n 06_presentコード

Cü 05_dataコードで present 指示節を使用

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

Fsubroutine calc(nx, ny, a, b, c)...
!$acc kernels present(a, b, c)!$acc loop independent
do j = 1,ny!$acc loop independent
do i = 1,nxc(i,j) = a(i,j) + b(i,j)
end doend do!$acc end kernels
end subroutine

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

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

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

リダクション計算（２）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
122簡単なOpenACC: reduction指示節（１）n 07_reductionコード
ü 06_presentコードで reductionを使用

C

openacc_basic/07_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 07_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)123簡単なOpenACC: reduction指示節（１）n 07_reductionコード
ü 06_presentコードで reductionを使用
sum = 0!$acc kernels present(c)!$acc loop reduction(+:sum)
do j = 1,ny!$acc loop reduction(+:sum)do i = 1,nx
sum = sum + c(i,j)
end doend do!$acc end kernels

F

openacc_basic/07_reduction

並列ループ毎にreduction

ü data 指示文で c を create に変更。

n 07_reductionコード
ü リダクションコードが生成された。
$ makepgfortran -O3 -mp -acc -ta=tesla,cc60 -Minfo=accel -c main.f90
(省略)
main:
(省略)
86, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
84, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
86, !$acc loop gang, vector(32) ! blockidx%x threadidx%x
Generating reduction(+:sum)124簡単なOpenACC: reduction指示節（２）n 07_reductionコードの実行ü 答えは正しく、速度が上がった。

openacc_basic/07_reduction

ü 配列 c の転送が削減されたこと、リダクションがGPU上で行われることによる性能向上。
$ qsub ./run.sh$ cat run.sh.o??????
mean = 3000.00
Time =
1.089 [sec]125OpenACC化のステップのまとめ

Cn 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/07_reduction

126OpenACC化のステップのまとめFn OpenACC化のための3つの指示文の適用
ü kernels 指示文を用いてGPUで実行する領域を指定
ü data 指示文を用い、ホスト-デバイス 間の通信を最適化ü loop 指示文を用い、並列処理の指定!$acc data copyin(a) create(b,c)!$acc kernels present(b,c)!$acc loop independent
do j = 1,ny!$acc loop independent
do i = 1,nxb(i,j) = b0
end doend doc(:,:) = 0.0!$acc end kernels
do icnt = 1,nt
call calc(nx, ny, a, b, c)
end do!続き
sum = 0!$acc kernels present(c)!$acc loop reduction(+:sum)
do j = 1,ny!$acc loop reduction(+:sum)do i = 1,nx
sum = sum + c(i,j)
end doend do!$acc end kernels!$acc end data

!続く

openacc_basic/07_reduction
127

OPENACC入門実習

128

C

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

F

実習
n 3次元拡散方程式のOpenACC化ü サンプルコード： openacc_diffusion/01_originaln 3次元拡散方程式のCPUコードにOpenACC の kernels, data, loop 指示文を追加し、GPUで高性能で実行しましょう。
do k = 1, nz
do j = 1, ny
do i = 1, nx

diffusion.f90, diffusion3d 関数内

w = -1; e = 1; n = -1; s = 1; b = -1; t = 1;
if(i == 1) w = 0
if(i == nx) e = 0
if(j == 1) n = 0
if(j == ny) s = 0
if(k == 1) b = 0
if(k == nz) t = 0
fn(i,j,k) = cc * f(i,j,k) + cw * f(i+w,j,k) &+ ce * f(i+e,j,k) + cs * f(i,j+s,k) + cn * f(i,j+n,k) &+ cb * f(i,j,k+b) + ct * f(i,j,k+t)
end doend doend doopenacc_diffusion/01_original130

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

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

131

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

t=n

t = n+1

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

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

132

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

平均後の
自分自身の値133

上下左右の値

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

自分自身の値の４倍

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

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

133

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

134

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

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

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

F90 = pgfortran
RM = rm -f
FFLAGS...

= -O3 -mp -acc -ta=tesla,cc60 -Minfo=accel

136

C

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

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

F

OpenACC化(1): kernels
n diffusion3d関数に kernelsを追加しましょう!$acc kernels copyin(f) copyout(fn)do k = 1, nz
do j = 1, ny
do i = 1, nx
w = -1; e = 1; n = -1; s = 1; b = -1; t = 1;
if(i == 1) w = 0
if(i == nx) e = 0
if(j == 1) n = 0
if(j == ny) s = 0
if(k == 1) b = 0
if(k == nz) t = 0
fn(i,j,k) = cc * f(i,j,k) + cw * f(i+w,j,k) &+ ce * f(i+e,j,k) + cs * f(i,j+s,k) + cn * f(i,j+n,k) &+ cb * f(i,j,k+b) + ct * f(i,j,k+t)
end doend doend do!$acc end kernels

diffusion.f90, diffusion3d 関数内make して実行してみましょう。
138

C

OpenACC化(2): loop
n diffusion3d関数に loopを追加しましょう#pragma acc kernels copyin(f[0:nx*ny*nz]) copyout(fn[0:nx*ny*nz])#pragma acc loop independentfor(int k = 0; k < nz; k++) {#pragma acc loop independentfor (int j = 0; j < ny; j++) {#pragma acc loop independentfor (int i = 0; i < nx; i++) {const int ix = nx*ny*k + nx*j + i;const int ip = i == nx - 1 ? ix : ix + 1;
const int im = i == 0
? ix : ix - 1;
const int jp = j == ny - 1 ? ix : ix + nx;
const int jm = j == 0
? ix : ix - nx;
高速化よりも、まずは正し
const int kp = k == nz - 1 ? ix : ix + nx*ny;
const int km = k == 0
? ix : ix - nx*ny; い計算を行うコードを保つ
fn[ix] = cc*f[ix]+ ce*f[ip] + cw*f[im]+ cn*f[jp] + cs*f[jm]+ ct*f[kp] + cb*f[km];ことが大事です。
末端の関数から修正を進め
ます。

}
}
}
return (double)(nx*ny*nz)*13.0;}

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

139

F

OpenACC化(2): loop
n diffusion3d関数に loopを追加しましょう!$acc kernels copyin(f) copyout(fn)!$acc loop independent
do k = 1, nz!$acc loop independent
do j = 1, ny!$acc loop independent
do i = 1, nx

高速化よりも、まずは正しい計
算を行うコードを保つことが大事です。
末端の関数から修正を進めます。

w = -1; e = 1; n = -1; s = 1; b = -1; t = 1;
if(i == 1) w = 0
if(i == nx) e = 0
if(j == 1) n = 0
if(j == ny) s = 0
if(k == 1) b = 0
if(k == nz) t = 0
fn(i,j,k) = cc * f(i,j,k) + cw * f(i+w,j,k) &+ ce * f(i+e,j,k) + cs * f(i,j+s,k) + cn * f(i,j+n,k) &+ cb * f(i,j,k+b) + ct * f(i,j,k+t)
end doend doend do!$acc end kernels

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

140

OpenACC化(3): データ転送の最適化(1)C

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

141

OpenACC化(3): データ転送の最適化(1)F

n diffusion3d関数で present とし、main関数で data を追加!$acc kernels copyin(f) copyout(fn)!$acc loop independent
do k = 1, nz!$acc loop independent
do j = 1, ny!$acc loop independent
do i = 1, nx

diffusion.f90, diffusion3d 関数内

w = -1; e = 1; n = -1; s = 1; b = -1; t = 1;
if(i == 1) w = 0
if(i == nx) e = 0
if(j == 1) n = 0
if(j == ny) s = 0
if(k == 1) b = 0
if(k == nz) t = 0
fn(i,j,k) = cc * f(i,j,k) + cw * f(i+w,j,k) &+ ce * f(i+e,j,k) + cs * f(i,j+s,k) + cn * f(i,j+n,k) &+ cb * f(i,j,k+b) + ct * f(i,j,k+t)
end doend doend do!$acc end kernels

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

142

OpenACC化(4): データ転送の最適化(2)C

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

143

OpenACC化(4): データ転送の最適化(2)F

n diffusion3d関数で present とし、main関数で data を追加!$acc data copy(f) create(fn)call start_timer()main.f90, main 関数内

do icnt = 0, nt-1
if(mod(icnt,100) == 0) write (*,"(A5,I4,A4,F7.5)"), "time(",icnt,") = ",timeflop = flop + diffusion3d(nx, ny, nz, dx, dy, dz, dt, kappa, f, fn)call swap(f, fn)
time = time + dt
if(time + 0.5*dt >= 0.1) exit
end doelapsed_time = get_elapsed_time()!$acc end data

copy/create など適切なものを選びます。
make して実行してみましょう。 どのくらいの実行性能が出ましたか？
OpenACC化の例は、openacc_diffusion/02_openacc

144

PGI_ACC_TIME によるOpenACC 実行の確認
nn

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

Reedbush でジョブに環境変数PGI_ACC_TIME を設定する場合は、ジョブスクリプト中に記載する。
$ cat run.sh...
. /etc/profile.d/modules.shmodule load pgi/18.7
export PGI_ACC_TIME=1./run

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

145

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

カーネル
実行時間

146

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

光の屈折と回折（１）n屈折ü 光が異なる媒質の境界で進行方向を変えること
ü 波の進む速度（位相速度）が媒質によってことなるため

http://www.kagakunavi.jp
http://www2e.biglobe.ne.jp/~shinzo/jikken/pikatto/pikatto.html

148

光の屈折と回折（２）n 回折ü 光の進路に障害物があるとき、その障害物の陰など、一見すると幾何学的には到達できない領域に回り込んで伝わる現象

http://rikanet2.jst.go.jp
http://web.me.com/takeda_t/atelier_bonryu/PH_Salon_R.4.html

149電磁波の伝播n 光は電磁波の一種
n 電場と磁場と電磁波の進行方向

E電場

H磁場
電磁波の進行方向ü 電磁波は、空間の電場と磁場がお互いの電磁誘導によって相互に発生して、空間を横波となって伝播する150

http://www.tecnoao-asia.com/about_abc.html

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

Maxwell 方程式の一部

1
⇥E
=
⇥t

H

H
=
t

1
µ

E

( : 誘電率)( µ :法）*を用いて解いて行きます。
透磁率)
この方程を、2次元FDTD法（Finite-differencetime-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.
151

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
152

FDTD法（２）n 例えば、153

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

13
12
11

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

9

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

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

7
6
5
4
3

mgn ( = 4 )

2
1
0

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

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

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

openacc_fdtd/02_openacc1

calc_ex_ey, pml_boundary_ex, pml_boundary_ey,
がOpenACC。

openacc_fdtd/03_openacc2

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

openacc_fdtd/04_openacc3

初期化を含め OpenACC。

openacc_fdtd/05_openacc4

データ移動の最適化。

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

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

156

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

157

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

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

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

mgndomain分割領域の全体領域

mgn

GPU1

ny

ny0

Subdomain 1

mgn

全領域の中心領域

mgn

mgn

GPU2

mgn

nx0

mgn

ny

Subdomain 2

全領域の全体領域mgn

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

mgn

nx

mgn

158

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

159

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

160

mgn

計算領域の設定（２）

mgn

ny

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

mgn

mgn

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

const struct Range whole

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

分割領域の
中心領域

mgn

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

ny

ny0

const struct Range insidemgn

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

mgn

mgn

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

nyn Range 構造体

mgn

nx

mgn

161

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

ny

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

mgn

mgn

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

163

時間発展（１）n 前半
ü 電場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;
(後半へ）
164

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

165

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

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

subdomains[i]

ny

mgn

n

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

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]= dst_hz

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

166

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

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

mgnn= src_exnylnx * (jb - jbw)

mgn

subdomains[i+1](ib, jb)

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

= dst_ex

mgn

subdomains[i]

ny

mgn

(ibw, jbw)

(ibw, jbw)

mgn(ib, jb)nx

mgn

nx = inside.length[0]ny = inside.length[1]ib = inside.begin[0]jb = inside.begin[1]lnx = whole.length[0]lny = whole.length[1]ibw = whole.begin[0]jbw = whole.begin[1]const int nhalo= whole.length[0];const int inside_end1 = inside.begin[1] + inside.length[1];const int src_ex
const int dst_ex= whole.length[0] * (inside.begin[1] - whole.begin[1]);= whole.length[0] * (inside_end1- whole.begin[1]);MPI_Send(&ex[src_ex], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD);MPI_Recv(&ex[dst_ex], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD, &status);

167

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

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

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

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

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

計算結果の表示
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 をつけていないと表示されない。うまく表示できない場合
は画像を手元にコピーして表示してください。
171計算結果の例
n 出力されたBMPファイルの一例

ny0

mgn

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

y

mgn

入射光

mgn

z

x

nx0

mgn

172

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

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

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

n main.cüü

OpenACCヘッダーの追加
data 指示文の追加

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

解答例は、 openacc_fdtd/02_openacc1
173

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

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

175

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

n fdtd2d.cü

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

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

解答例は、 openacc_fdtd/03_openacc2
176

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

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

178

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

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

解答例は、 openacc_fdtd/04_openacc3
179

実習4
n 計算領域のサイズなどを変更して性能測定してみましょう。
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_fdtd/05_openacc4

180

Q&A
n アカウントは1ヶ月有効です。
n 資料のPDF版はWEBページに掲載します。
n https://www.cc.u-tokyo.ac.jp/events/lectures/133/
n アンケートへの協力をお願いします。
n とくにオンラインの講習会については初めてなので、不便な点等お知らせいただけ

ますと幸いです。

181

補足スライド
182

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

メモリread終了

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

階層的スレッド管理とコミュニケーション
スレッドブロックn 階層的なコア/スレッド管理
n P100は56 SMを持ち、1 SMは64 CUDA

coreを持つ。トータル3584 CUDA coren 1 SMが複数のスレッドブロックを担当
し、1 CUDA core が複数スレッドを担当
n スレッド間のコミュニケーション
n 同一スレッドブロック内のスレッドは高速コミュニケーション可能
n 異なるスレッドブロックに属するスレッ

ド間はコミュニケーションが低速
n いったんメモリに書き出したり、CPUに処理を戻さなくてはならない

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

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

1

2

3

…

31 32

スレッド

1

2

3

…

31 32

配列 A

4

3

5

…

8

配列 A

4

3

5

…

8

0

÷ × ＋ …

−

×

2

1

9

配列 B

0

× × × …

× ×

2

1

3

1

…

OK！

9

配列 B

3

1

…

NG！
185

Warp内分岐
n Divergent Branch

CUDA 8 以前のバージョン
CUDA 9 以上では多少マシになるが、ペナルティがあることに変わりはない

n Warp 内で分岐すること。Warp単位の分岐ならOK。
Divergent branchなし::
if ( TRUE ) {::} else {::}::else 部分は実行せずジャンプ

Divergent branchあり::
if ( 奇数スレッド ) {::} else {::}::一部スレッドを眠らせて全分岐を実行
最悪ケースでは32倍のコスト

186コアレスドアクセスn

同じWarp内のスレッド（連続するスレッド）は近いメモリア
ドレスへアクセスすると効率的ü コアレスドアクセス（coalesced access）と呼ぶ
ü メモリアクセスは128 Byte 単位で行われる。128 Byte に収まれば1回のアクセス、超えれば128 Byte アクセスをその分繰り返す。

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

160

192

224

256

...

15 16

...

31

160

192

224

256

...

15 16

...

31

Thread0

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

Thread0187

ストライドアクセスがあるとどうなるか
n GPUはストライドア

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

クセスに弱い！

Memory Bandwidth [GB/sec]

600
500

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

400
300
200
100
0
0

20

40

60

80
Stride

100

120

140
188

OpenACC と Unified Memoryn Unified Memory とは…
n 物理的に別物のCPUとGPUのメモリをあたかも一つのメモリのように扱う機能
n Pascal GPUではハードウェアサポートn ページフォルトが起こると勝手にマイグレーションしてくれるn OpenACC と Unified Memoryn OpenACCにUnified Memoryを直接使う機能はない
n PGIコンパイラではオプションを与えることで使える
n pgfortran –acc –ta=tesla,managed

n 使うとデータ指示文が無視され、代わりにUnified Memoryを使う
n ハイエンドのNVIDIA GPU + PGI compilerの環境が揃いさえすれば、データ転送を考える必要がなく非常に楽n OpenACCのコードとして間違ったコードでも動いてしまう189