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

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

講習会スケジュール
n 開催日時ü 10月5日（水） 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の登録メールの配送に小一時間かかることがあります4

Zoom: 「手を挙げる」方法
1.

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

2.

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

5

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

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

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

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

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

2.

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

7

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

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

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

メッセージを送信する

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

9

サンプルコードのダウンロード
n 講習会で使うサンプルコードはGitHubにありますn https://github.com/hoshino-UTokyo/lecture_openacc

n ダウンロード手順1.

$ ssh tXXXXX@wisteria.cc.u-tokyo.ac.jp
n wisteriaへのssh。tXXXXXはアカウント名

2.
3.
4.
5.

$ cd /work/gt00/tXXXXX$ git !"#$%&'(()*+,,-.('/01!#2,'#*'.$#345#67#,"%!(/8%9#)%$:!!1-.($ cd lecture_openacc/$ cd C or F
n

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

10

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

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

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

54.9 TFLOPS, 11.2 TB

Oakforest-PACS
Fujitsu, Intel KNL

T2K Tokyo

140TF, 31.3TB

Big Data &Extreme Computing

25PFLOPS, 919.3TB

Oakleaf-FX: Fujitsu PRIMEHPC
FX10, SPARC64 IXfx

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

12

既存システム：利用者2,000+，学外50+%n Reedbush (SGI, Intel BDW + NVIDIA P100 (Pascal))n データ解析・シミュレーション融合スーパーコンピュータ
n 3.36 PF, 2016年7月〜 2021年11月末
n 東大ITC初のGPUシステム (2017年3月より), DDN IME (Burst Buffer)n Oakforest-PACS (OFP) (富士通, Intel Xeon Phi (KNL))
n 2016年10月〜2022年3月末
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を搭載
13

Wisteria/BDEC-01

BDEC:「計算・データ・学習（S+D+L）」
融合のためのプラットフォーム（Big Data & Extreme Computing）

n 2021年5月14日運用開始
n 東京大学柏Ⅱキャンパス
n 33.1 PF, 8.38 PB/sec. ，富士通製

Platform for Integration of (S+D+L)
Big Data & Extreme Computingn ~4.5 MVA（空調込み）, ~360m2

n Hierarchical, Hybrid, Heterogeneous (h3)n 2種類のノード群
n シミュレーションノード群 (S, SIM)：
Odysseyn 従来のスパコン
n Fujitsu PRIMEHPC FX1000 (A64FX), 25.9 PF
n7,680ノード（368,640 コア），20ラック，Tofu-Dn データ・学習ノード群（D/L, DL）：Aquarius
n データ解析，機械学習
n Intel Xeon Ice Lake + NVIDIA A100, 7.2 PF
n 45ノード（Ice Lake：90基，A100：360基）, IB-HDRn 一部は外部リソース（ストレージ，サーバー，センサーネットワーク他）に直接接続
n ファイルシステム：共有（大容量）＋高速

シミュレーションノード群

Odyssey
Fujitsu/Arm A64FX
25.9PF, 7.8PB/s

Shared File
System

2.0TB/s

Fast File
System

共有ファイル
システム（SFS）

データ・学習ノード群
Aquarius

高速ファイルシステム（FFS）

25.8 PB, 500GB/s

Intel Ice Lake + NVIDIA A100

7.20PF, 578.2TB/s

1PB, 1.0TB/s

800 Gbps

External
Resources
External Network
外部ネットワーク

外部リソーススパコン料金表（2021年4月時点 ）n 最小セット料金表
トークン※1

料金（大学・
公共機関）

ストレージ
容量

利用期間

Wisteria/BDEC-01

720 ※25,000円

2TB

2021/8/2以降～年度末まで

Reedbush

720 ※36,300円

1TB

2021/11/30まで

Oakforest-PACS

720

4,200円

1TB

年度末まで (2021年度で運用終了予定）

Oakbride-CX

720

8,400円

4TB

年度末まで※1 トークン≒ノード時間 ※2, 3。720トークンなら、1ノードを720時間利用できる。
※2 Odyssey(CPUノード)基準。Aquarius(GPUノード)を使う場合、利用するGPU数x3.0倍のトークンを消費する。
つまり1ノード(8GPU)を利用する場合、720/24=30時間しか利用できない。
※3 Reedbush-H,L はそれぞれ2.5, 4.0倍のトークンを消費する。つまり、Hの場合は720/2.5=288時間, Lの場合は720/4=180時間しか利用できない

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 cuda ompi-cuda
NVIDIA コンパイラ
言語
GNUコンパイラ Intelコンパイラ(旧PGI)CUDAコンパイラ
C

gcc

icc

nvc (pgcc)

C++

g++

icpc

nvc++(pgc++)

gfortran

ifort

Fortran
OpenACC

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

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

22

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

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

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

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

並列計算
n 実行時間 T の逐次処理のプログラムを p 台の計算機で並列計算することで、実行時間を T / p にする。

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

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

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

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

仕事4

仕事5

時間
25

データ並列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];26GPUにおけるループ並列化
n GPU における高速化は通常、プログラム中の重たいループ構造を並列化することで達成する
n 今回学ぶOpenACC は特定のループ構造を簡単に並列化できる
n 全てのループ構造を並列化できるわ

けではないn どのようなループなら並列化可能か知る必要がある

世の中のすべてのループ

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

OpenACCで高速化可能なループ

27

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 を駆使する必要がある)

28

簡単に並列化できるループ 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]29

簡単に並列化できるループ 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]30

簡単に並列化できるループ 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]31

簡単に並列化できるループ 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]32

簡単に並列化できるループ 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]33

簡単に並列化できるループ 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]34

簡単に並列化できるループ 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]35

簡単に並列化できるループ 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]36

簡単に並列化できるループ 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]37

簡単に並列化できるループ 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]38

簡単に並列化できないループ
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]39

簡単に並列化できないループ
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]40

簡単に並列化できないループ 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]41

簡単に並列化できないループ 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]42

簡単に並列化できないループ 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]43

簡単に並列化できないループ 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]44

簡単に並列化できないループ 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]45

簡単に並列化できないループ 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]46

簡単に並列化できないループ 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]47

簡単に並列化できないループ 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]48

簡単に並列化できないループ 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]49

簡単に並列化できないループ 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]50

簡単に並列化できないループ 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]51

簡単に並列化できないループ 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]52

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

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

※6です

４や！
5やな！

４！

メモリ

4

6

1

３ + 1;
A[0] = A[2]CPUで実行される足し算は、
1. データの読み込み
2. 足し算3. データの書き込みの3パートからなる。
スレッドは各々独立に１〜３を実行するため、タイミング
によって結果が変わる！
（この例の場合は4,5,6のいず
れかになる）

A[0] A[1] A[2]53

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

配列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

54

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

配列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

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

55

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

配列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)
という）

56

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

配列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)から繰り返し
57

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

配列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倍の高速化
メモリなどを介してスレッドの間でデータのやり取りをすることをスレッド間通信というスレッドの同期・通信が入ると途端に難しくなる！

58

GPU入門59

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

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

60

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

61

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

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

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

n コアレスドアクセス
62

NVIDIA A100 Tensor Core GPU (1/2)
NVIDIA A100 Tensor コア GPU アーキテクチャの詳細

n 108 SM (Streaming Multiprocessor)

図4.

出典:
NVIDIA
128 個の SM を搭載した GA100 フル GPU (A100 Tensor
コア
GPUA100
は Tensor core GPUアーキテクチャ

63

NVIDIA A100 Tensor Core GPU (2/2)

NVIDIA A100 Tensor コア GPU アーキテクチャの詳細

n 倍精度にも対応したTensor Coreを搭載
n 19.5 TF @ FP64, FP32
n 156 TF @ TF32 (実質19bit)
n 312 TF @ FP16, BF16n 624 TF @ INT8n 1248 TF @ INT4

n メモリ HBM2 40GB 搭載
n 1.555 TB/s

図 5.
GA100 ストリーミング マルチプロセッサ (SM)出典:
NVIDIA
A100 Tensor core GPUアーキテクチャ

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

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

OSは存在しない

２．計算を行う

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

メインメモリ

デバイス
メモリ
n
n

計算はOSのあるCPUから始まる
物理的に独立のデバイスメモリとデータのやり取り必須

物理的に独立65

どんなアプリならGPUで高速化できる？
n 原則：GPUに一度送ったデータを使い回せるアプリケーションn 最低でも100回は使いまわしたい
n 例：データ量 N に対して計算量 O(N2) 以上の計算（行列積、多体問題など）や、反復法など

n 思考実験
n 次のプログラムを、右の表のコンピュータの(1) CPUを使った時の実行時間は？
(2) GPUを使った時の実行時間は？
double precision :: A(1:N), B(1:N)
if(GPU) BをCPUからGPUにコピー
do i = 1, N
A(i) = B(i)end doif(GPU) AをGPUからCPUにコピーあるコンピュータの性能
データ転送性能

CPUのメモリ

100 GB/sec

GPUのメモリ

1000 GB/sec

CPU-GPU間のバス

20 GB/sec(1) 配列A・Bのbyte数 / CPUのメモリ性能= N * 2 * 8 / 100(2) 配列A・Bのbyte数 / GPUのメモリ性能+ 配列A・Bのbyte数 / CPU-GPU間バスのメモリ性能= N * 2 * 8 / 1000 + N * 2 * 8 / 20
N = 109 (1G) なら？

(1) 0.16 sec (2) 0.816 sec

66

どんなアプリならGPUで高速化できる？
n 原則：GPUに一度送ったデータを使い回せるアプリケーションn 最低でも100回は使いまわしたい
n 例：データ量 N に対して計算量 O(N2) 以上の計算（行列積、多体問題など）や、反復法など

n 思考実験
n 次のプログラムを、右の表のコンピュータの(1) CPUを使った時の実行時間は？
(2) GPUを使った時の実行時間は？

double precision :: A(1:N), B(1:N)
if(GPU) BをCPUからGPUにコピー
do t = 1, 100
do i = 1, N
100回使いA(i) = B(i)
回してみる
end doend doif(GPU) AをGPUからCPUにコピーあるコンピュータの性能
データ転送性能

CPUのメモリ

100 GB/sec

GPUのメモリ

1000 GB/sec

CPU-GPU間のバス

20 GB/sec(1) 100 * 配列A・Bのbyte数 / CPUのメモリ性能= 100 * N * 2 * 8 / 100(2) 100 * 配列A・Bのbyte数 / GPUのメモリ性能+ 配列A・Bのbyte数 / CPU-GPU間バスのメモリ性能= 100 * N * 2 * 8 / 1000 + N * 2 * 8 / 20
N = 109 (1G) なら？

(1) 16 sec (2) 2.4 sec

67

OPENACC入門
68

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

簡単

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

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

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

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

ンターフェースü OpenMP のようなディレクティブ（指示文）・ベース
ü C 言語/C++, Fortran に対応
ü 2011年秋に OpenACC1.0、最新は 3.1
ü コンパイラ：［商用］PGI → NVIDIA HPC SDK,
Cray,［フリー］GCC（NVIDIA HPC SDKは無料版あり）
ü 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）用コード、複数のアクセラレー

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

OpenACCでできること
n OpenACC は特定のループ構造を簡

単に並列化できるü 全てのループ構造を並列化できるわけではない

世の中のすべてのループ

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

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

ループが、データ独立か、リダクション
か、それ以外か71

OpenACCでできないこと
n CUDAならshared memoryなど使って頑張れば並列化できる、

データ依存性のあるループの並列化
n 例外：atomic演算で解決可能な書き込み競合を含むループ

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

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

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

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

73

OpenACCを推奨する理由
CPU
main

GPU

実行時間
12

subA

10

subA

8
6

subB
subC

subB

4
2

subC

0
CPU

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

74

OpenACCを推奨する理由
CPU

GPU

main

実行時間
12

subA

10
8

subA
subB
subC

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

6

subB

4
2

subC

0
CPU

75

OpenACCを推奨する理由
CPU

GPU

実行時間

main

12
10
8

subA
subB

6

データ転送

subB

データ転送

4
2

subC

アレ？！

0
CPU

CPU+GPU

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

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

77

OpenACCを推奨する理由
CPU

GPU

main

subA

12

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

subC

実行時間
10
8
6

subB

4
2

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

CPU

CPU+GPU

78

OpenACCを推奨する理由
CPU

GPU

main

実行時間
12

データ転送

10

subA

8
6

subB
subC

4
2
0

データ転送

CPU

CPU+GPU

結局全部CUDA化した…

79

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

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

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

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

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

以上により、CUDA化と遜色ない性能を少ない工数で達成できる81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
82

OpenACCの実行イメージ
1スレッド

OpenMP

OpenACC・CUDA

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

CPU

CPU

CPU

デバイス
(GPU)
83

はじめての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から
copyout84

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

85

OpenACCの主な指示文
n アクセラレータ（GPU）実行領域指定指示文（必須）ü kernels, paralleln ループ最適化指示文（オプションだがほぼ必須）
ü loopn データ移動指示文（オプションだがほぼ必須）
ü data, enter data, exit data, updaten その他
ü host_data, atomic, routine, declare

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

86アクセラレータ実行領域の指定
n kernels 指示文（必須）ü 囲まれた領域がアクセラレータで実行されるカーネルにü 複数のループネストを囲んだ時、一般にはそれぞれのループネストが別々のカーネルに
ü 右の例ではカーネルが２つ生成されると思われるが、コンパイラの実装次第であるため、２つに分ける必要があるならkernels指示文を２つ使うべきü 推奨：基本的には、ループネスト一つにつき一つのkernels指示文
ü 注意点： kernels 指示文終了時に暗黙の同期(GPU内のスレッド）が取られる。
ü 似た指示文として、領域内が一つのカーネ
ルとして生成される parallel 指示文もある

int main() {#pragma acc kernels{for (int i=0; i<n;kerneli++) {
A[i] = 0;
ループネストが独立なら、まと}めて囲んでも大丈夫。
for (int i=0; i<n;どのように実行されるかはコンi++) {
パイラ次第。
B[i] = 0;}}}
int main() {推奨#pragma acc kernelsfor (int i=0; i<n; i++) {kernel1
A[i] = 0;
ここで同期。つまり}
kernel1 の終了が保#pragma acc kernels
証される。
for (int i=0; i<n; i++) {B[i] = A[i];kernel2}}
kernel2 が kernel1 に依存して
いる

87

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 化
1.

GPUで実行したいループを
kernelsで囲む#pragma acc kernels

kernels直後の{}で囲まれる領域がfor (int i=0; i<n; i++) {b[i] = a[i] + c; GPU上で実行される}

ループはベストエフォートで並列化される（C言語ではほとんどされない）

double sum = 0;
必要なデータ転送はベストエフォートで行われる（C言語ではよく失敗する）
for (int i=0; i<n; i++) {sum += b[i];}fprintf(stdout, "%f¥n", sum/n);free(a); free(b);return 0;
}

88

F

CPUコードのOpenACC化
openacc_hello/01_hello_acc
program main
implicit none
! 変数宣言allocate(a(n),b(n))
c = 2.0n ループのOpenACC 化
1.

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

do i = 1, na(i) = 10.0
end do!$acc kernels
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 main

Fortranの場合、kernels ~ end kernels の間がGPUで実行される
ループはベストエフォートで並列化される（Fortranでは概ね成功する）
必要なデータ転送はベストエフォートで行われる（Fortranでは概ね成功する）

89

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

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

91

ループ最適化指示文
n loop 指示文（オプションだがほぼ必須）
ü ループの並列化の可否を教える
Ø データ独立なループ(independent)Ø リダクションループ (reduction)
Ø 並列化すべきでないループ (seq)ü ループマッピングのパラメータの調整Ø 難しいので、最初は考える必要はない•

コンパイラがある程度最適な値を決定してくれるので任せていいØ gang, worker, vector を用いて指定する

gang: CUDA で言う thread block 数の指定。グルー
プ単位での処理の分散を行う際に用いる。よほどの玄人以外はgangの数まで指定すべきではない。
• worker: GPU では使わない• vector: CUDA で言う thread block 内の thread 数の指定。グループ内での処理の分散を行う際に用いる。
数を指定するなら、1024以下の32の倍数が良い。
•

ループ指示文指定例#pragma acc kernels#pragma acc loop independentfor (int i=0; i<n; i++) {A[i] = 0;
データ独立ループ}
double sum = 0;
#pragma acc kernels#pragma acc loop reduction(+:sum)for (int i=0; i<n; i++) {sum += A[i];リダクションルー}
プ
double sum = 0;
#pragma acc kernels#pragma acc loop independent gangfor (int j=0; j<n; j++) {#pragma acc loop independentvector(64)for (int i=0; i<n; i++) {多重ループへのsum += A[i];gang, vector適
}
用
}
92

データの独立性
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];}

93

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

4.

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

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

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 化
1.
2.

GPUで実行したいループを
kernelsで囲む
loop independent でループが並列化可能であることを教える#pragma acc kernels#pragma acc loop independent
ループが並列化可能と見なされる（並列化可能でないループに
for (int i=0; i<n; i++) {independent を付けると結果が間違う）
b[i] = a[i] + c;}

必要なデータ転送はベストエフォートで行われる（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;
}

96

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 indepnedent
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 化
1.
2.

GPUで実行したいループを
kernelsで囲む
loop independent でループが並列化可能であることを教える

ループが並列化可能と見なされる（並列化可能でないループに
independent を付けると結果が間違う）
必要なデータ転送はベストエフォートで行われる（Fortranでは概ね成功する）

97

Kernels指示文の自動データ転送
n kernels 構文に差し掛かると、
ü OpenACCコンパイラは実行に必要なデータを自動で転送する。
Ø 往々にして失敗するため、後述のdata指示文、またはGPUのUnified memory機能を利用すべきü 配列はGPUのメモリに確保され、shared 変数として扱われる。
Ø デバイスメモリに動的に確保され、スレッド間で共有。
Ø デバイスからホストへコピーすることが可能。
Ø C言語の場合特に、配列のサイズがわからないなどで失敗する。
Ø 各スレッドでprivateに扱うべき小さな配列は、acc kernels private(配列名）とする。
ü スカラ変数は firstprivate または private 変数として扱われる。
Ø ホストからデバイスへコピーが渡され初期化。ホストに戻せない。
Ø スカラ変数に関しては、自動転送に任せていいü 構文に差し掛かるたびに転送を行う。data 指示文で制御できる。
98

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が行われる

99

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などの組み込み関数をサポートしている)、
基本的にサイズを書く必要がない。

100

data 指示文の指示節
n

copy
ü

n

allocate, memcpy(H->D), memcpy(D->H), deallocate

copyinallocate, memcpy(H->D), deallocateü 解放前にホストへデータをコピーしないü

ncopyout
allocate, memcpy(D->H), deallocate
ü 確保後にホストからデータをコピーしないüncreateallocate, deallocateü コピーしないünpresentün

何もしない。既にデバイス上で確保済みであることを伝える。

copy/copyin/copyout/create は既にデバイス上確保されているデータに対しては何もしない。
present として振る舞う。（OpenACC2.5以降）
101

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);}102}データの移動範囲の指定
n ホストとデバイス間でコピーする範囲を指定
ü 部分配列の転送が可能
ü Fortran と C言語で指定方法が異なるので注意

n 二次元配列A転送する例
ü Fortran: 下限と上限を指定!$acc data copy(A(lower1:upper1, lower2:upper2) )...
!$acc end data

ü C言語: 始点とサイズを指定#pragma acc data copy(A[begin1:length1][begin2:length2])...
103

Unified Memoryn Unified Memory とは…
ü 物理的に別物のCPUとGPUのメモリをあたかも一つのメモリのように扱う

機能
ü NVIDIA A100 GPUではハードウェアサポート
Ø ページフォルトが起こると勝手にマイグレーションしてくれるn OpenACC と Unified Memory
ü OpenACCの仕様にUnified Memoryを直接使う機能はない
Ø nvidia コンパイラではオプションを与えることで使える
Ø nvfortran –acc –ta=tesla,managedü 使うとデータ指示文が無視され、代わりにUnified Memoryを使う
Ø ハイエンドのNVIDIA GPU + NVIDIA compilerの環境が揃いさえすれば、データ転送を考える必要がなく非常に楽 (Wisteria環境では利用推奨）
Ø data指示文が間違ったOpenACCコードでも正しく動いてしまうØ Unified memory を使うと、GPU direct というMPIの直接通信機能が使えない
Ø CPU側のメモリアロケーションを全部監視してるので、遅くなるケースがある

104

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

n ループのOpenACC 化
1.
2.
3.

GPUで実行したいループを
kernelsで囲む
loop independent でループが並列化可能であることを教える
data 指示文でデータ転送を行う
n このケースではあまり data 指示文の意味はない。後の最適化で本領発揮。

105

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 data copyin(a) copyout(b)!$acc kernels!$acc loop indepnedent
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 mainn ループのOpenACC 化
1.
2.
3.

GPUで実行したいループを
kernelsで囲む
loop independent でループが並列化可能であることを教える
data 指示文でデータ転送を行う
n

このケースではあまり data 指示文の意味はない。後の最適化で本領発揮。

106参考：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);
...
}

107

OpenACCコードのコンパイルn

NVIDIAコンパイラによるコンパイル
ü WisteriaではOpenACCはNVIDIAコンパイラで利用できます。
$ module load nvidia$ nvc -O3 -acc -Minfo=accel –gpu=cc80 -c main.c-acc: OpenACCコードであることを指示-Minfo=accel:

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

GPUの種類を指定する。compute capability 8.0 (cc80) のコードを生成する。

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

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

110

配列のインデックス計算

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) = c(i,j) + a(i,j) + b(i,j)
end doend doend subroutine calc

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

111

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

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

112

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

C

F

C
F

openacc_basic/02_kernels

void calc(unsigned int nx, unsigned int ny, const float *a, const float *b, float *c){const unsigned int n = nx * ny;#pragma acc 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];}}}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) = c(i,j) + a(i,j) + b(i,j)
end doを shared 変数として自動で転送しend doてくれるはずだが…
!$acc end kernels
end subroutine calc113

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

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

CFC$ makenvc -O3 -acc -Minfo=accel -ta=tesla,cc80 -c main.c
NVC++-S-0155-Compiler failed to translate accelerator region (see -Minfo messages): Could
not find allocated-variable index for symbol - b (main.c: 11)
calc:
14, Complex loop carried dependence of a->,c->,b-> prevents parallelization
Accelerator serial 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
NVC++-F-0704-Compilation aborted due to previous errors. (main.c)
NVC++/x86-64 Linux 21.3-0: compilation aborted
make: *** [Makefile:33: main.o] エラー 2

F$ makenvfortran -O3 -mp -acc -ta=tesla,cc80 -Minfo=accel -c main.f90
calc:
13, Generating implicit copyin(b(:nx,:ny)) [if not already present]
Generating implicit copy(c(:nx,:ny)) [if not already present]
Generating implicit copyin(a(:nx,:ny)) [if not already present]
14, Loop is parallelizable
15, Loop is parallelizable
Generating Tesla code
14, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
15,! blockidx%x threadidx%x auto-collapsed

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

114

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

115

簡単なOpenACC: kernels 指示文（４）n 03_kernels_copyコード:初期化
ü CPUコードにkernels 指示文の追加

CF

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 はスカラ変数のため自動的に各スレッドへコピーが渡される。
116

簡単なOpenACC: kernels 指示文（５）Cn コンパイル
ü データの独立性がコンパイラにはわからず、並列化されない。
$ makenvc -O3 -acc -Minfo=accel -ta=tesla,cc80 -c main.c
calc:
11, Generating copy(a[:n],c[:n],b[:n]) [if not already present]
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 serial 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:
44, Generating copyout(c[:16777216],b[:16777216]) [if not already present]
45, Loop is parallelizable
Generating Tesla code
45, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */48, Loop is parallelizable
Generating Tesla code
48, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */117

簡単なOpenACC: kernels 指示文（５）Fn コンパイル
ü データの独立性を見切り、並列化。
nvfortran -O3 -mp -acc -ta=tesla,cc80 -Minfo=accel -c main.f90
calc:
13, Generating copyin(a(:,:)) [if not already present]
Generating copyout(c(:,:)) [if not already present]
Generating copyin(b(:,:)) [if not already present]
14, Loop is parallelizable
15, Loop is parallelizable
Generating Tesla code
14, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
15,! blockidx%x threadidx%x auto-collapsed
main:
61, Generating copyout(b(:,:),c(:,:)) [if not already present]
62, Loop is parallelizable
63, Loop is parallelizable
Generating Tesla code
62, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
63,! blockidx%x threadidx%x auto-collapsed
68, Loop is parallelizable
Generating Tesla code68,
! blockidx%x threadidx%x auto-collapsed!$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x

118

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

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

n 不明瞭な書き込み参照先
ü インデックス計算
Ø 計算結果がループ変数に対して独立かどうかわからないØ Fortranでも、多次元配列を一次元化すると起こるØ 逆に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];}ü 間接参照119

簡単なOpenACC: loop 指示文（１）Cn 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;}}120

簡単な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 do各次元についてloop指示文を指定する（並列サイズなどを指定したいなど）場合、do文で書き下す必要がある。

c(:,:) = 0.0!$acc end kernels
121簡単なOpenACC: loop 指示文（２）C

openacc_basic/04_loop

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

$ makenvc -O3 -acc -Minfo=accel -ta=tesla,cc80 -c main.c
calc:
11, Generating copy(a[:n],c[:n],b[:n]) [if not already present]
15, Loop is parallelizable
17, Loop is parallelizable
Generating Tesla code
15, #pragma acc loop gang, vector(128) collapse(2) /* blockIdx.x threadIdx.x */17,/* blockIdx.x threadIdx.x auto-collapsed */main:
46, Generating copyout(c[:16777216],b[:16777216]) [if not already present]
48, Loop is parallelizable
Generating Tesla code
48, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */52, Loop is parallelizable
Generating Tesla code
52, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */※Fortran版は既に並列化されていたため省略。loop
independent をつける事による挙動の変化はない。
122

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

openacc_basic/04_loop

ü 答えは正しいが、実行時間が大変長い。
$ pjsub ./run.sh$ cat run.sh.??????.out
mean = 3000.00
Time =
42.990 [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
}

123

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

124

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

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

openacc_basic/05_data$ pjsub ./run.sh$ cat run.sh.??????.out
mean = 3000.00
Time =
0.376 [sec]126簡単な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) = c(i,j) + a(i,j) + b(i,j)
end doend do!$acc end kernels
end subroutine

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

簡単なOpenACC: reduction指示節（１）C

n 07_reductionコード
ü 06_presentコードで reductionを使用

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コード
ü リダクションコードが生成された。
$ makenvc -O3 -acc -Minfo=accel -ta=tesla,cc80 -c main.c
(省略)
main:
(省略)
67, Loop is parallelizable
Generating Tesla code
67, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */Generating reduction(+:sum)128簡単なOpenACC: reduction指示節（１）F

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

openacc_basic/07_reduction

並列ループ毎にreduction

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

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

openacc_basic/07_reduction

ü 配列 c の転送が削減されたこと、リダクションがGPU上で行われることによる性能向上。
$ pjsub ./run.sh$ cat run.sh.??????.out
mean = 3000.00
Time =
0.353 [sec]130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

131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_reduction132

OPENACC入門実習

133

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
134

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

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

136

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

t = n+1

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

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

137

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

平均後の
自分自身の値138

上下左右の値

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

自分自身の値の４倍

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

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

138

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

139

CPUコードn CPUコードのコンパイルと実行$ cd openacc_diffusion/01_original$ make$ pjsub ./run.sh# cat run.sh.??????.out
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 =
8.677 [sec]Performance=
5.15 [GFlops]
Error[128][128][128] = 4.556413e-06

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

OpenACC化(0): Makefile の修正
n Makefile に OpenACC をコンパイルするよう –acc などを追加しましょう

C

CC
= nvc
CXX = nvc++
GCC = gcc
RM = rm -f
MAKEDEPEND = makedepend
CFLAGS
GFLAGSCXXFLAGSLDFLAGS...

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

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

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

141

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 して実行してみましょう。
142

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 して実行してみましょう。
143

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 してジョブ投入 pjsub ./run.shしてみましょう。遅いですが実行できます。

144

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 してジョブ投入 pjsub ./run.shしてみましょう。遅いですが実行できます。

145

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

146

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

147

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

148

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

149

NVCOMPILER_ACC_TIME によるOpenACC 実行の確認
nn

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

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

150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

カーネル
実行時間

151

Unified MemoryによるOpenACC 実行の確認
n

NVIDIA compiler で Unified memory を利用する場合、-ta オプションに managed を付ける
だけ！
data指示文が無視されるようになるü サンプルコードでは、data指示文をコメントアウトしているü 性能はどのくらい違うだろうか？
ü

C

CC
= nvc
CXX = nvc++
GCC = gcc
RM = rm -f
MAKEDEPEND = makedepend
CFLAGS...

F= -O3 -acc -Minfo=accel -ta=tesla,cc80,managed

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

= -O3 -mp -acc -ta=tesla,cc80,managed -Minfo=accel

サンプルコードは、 openacc_diffusion/04_openacc_managed

152

チャレンジ課題：GPUを用いた
FDTD法による電磁波伝搬計算(C言語版のみ）
153

光の屈折と回折（１）
n

屈折
ü
ü

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

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

154

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

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

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

E電場

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

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

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

Maxwell 方程式の一部

1
⇥E
=
⇥t
(

: 誘電率)

H

H
=
t

1
µ

E

( µ : 透磁率)

この方程を、2次元FDTD法（Finite-difference time-domain 法）*を用いて解いて行きます。

* K. S. Yee, "Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in
Isotropic Media,'' IEEE Trans. on Antennas and Propagat., vol. 14, pp. 302-307, May 1966.
157

FDTD法（1）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
158

FDTD法（２）n 例えば、159

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

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

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

データ移動の最適化。

161

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

162

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

163

コード全体の流れ（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関数※MPIはOFFにしてあります

FDTD計算繰り返し時間発展

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

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

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

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

165

n Range 構造体

mgn

ny

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

mgn

mgn

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

const struct Range whole

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

分割領域の
中心領域

mgn

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

ny

ny0

const struct Range insidemgn

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

mgn

mgn

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

ny

mgn

計算領域の設定（２）

mgn

nx

mgn

166

配列の確保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]167時間発展（１）n 前半
ü 電場Eの時間発展（calc_ex_ey）、境界条件(pml_boundary_...)
ü 入射光（plane_wave_incidence）※MPIはOFFにしてあります

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;
(後半へ）
168

時間発展（２）※MPIはOFFにしてあります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++;(出力など)}

169

チャレンジ課題:GPUを用いた
FDTD法による電磁波伝搬計算の実習（C言語版のみ）
170

プログラムのコンパイルと実行（１）n CPUコードのコンパイルと実行$ module load nvidia cuda/11.2 ompi-cuda$ 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_fdtd/01_original
MPIのmoduleが必要? の数字はジョブごとに変わります。
利用したノード（省略）
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 すると出力なしで実行する。性能測定用。

171

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

openacc_fdtd/01_originalmodule load cuda/11.2 #要書き換え（省略）
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 を指定すると出力しない。
172

計算結果の表示
n 計算結果は sim_run に BMP として出力される$ cd sim_run/
openacc_fdtd/01_originaln 計算結果の表示ü 1枚のBMPを見る$ display e05000.bmp
ü 複数のBMPファイルをアニメーションで表示$ animate *.bmp

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

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

ny0

mgn

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

y
z

mgn

入射光

mgn

x

nx0

mgn

174

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

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

n main.cüü

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

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

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

解答例は、 openacc_fdtd/02_openacc1
175

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

176

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

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コードをさらに最適化しましょう。
üü

NVCOMPILER_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

実習5
n 実はUnified Memoryを利用するとずっと簡単に実装できます。
Ø

実習2・3で行ったdata 指示文に関する実装する必要がないため

n -ta=tesla,cc80,managed として性能を比較してみましょうØ

性能がだいぶ違うと思います。改善するならどうすべきでしょうか。
ØØmanaged memory はページ単位（結構大きい）でデータ転送を行います
05_openacc4 で使っている update 指示文は、data 指示文で確保済みのデータを CPU-GPU間でコピーするための指示文です。

#pragma acc update host(ex[src:sendnelems],ey[src:sendnelems],hz[src:sendnelems])
for(i = 0;i < sendnelems;i++){ex_global[dst+i] = ex[src+i];
src から始まりsendnelems
ey_global[dst+i] = ey[src+i];
個の要素をhost(CPU)にコ
hz_global[dst+i] = hz[src+i];}

ピーする

性能改善済みの解答例は、openacc_fdtd/06_openacc5

181

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

182

補足スライド
183

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

メモリread終了

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

階層的スレッド管理とコミュニケーション
スレッドブロック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

185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！
186

Warp内分岐
n Divergent Branch

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

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

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

187コアレスドアクセス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

Thread0

188

ストライドアクセスがあるとどうなるか
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
189