CUDA・OpenACCによる
GPUコンピューティング
Akira Naruse, 23th June. 2016

GPUコンピューティング

Low latency + High throughput
CPU

GPU

2

アプリケーション実行
アプリケーション・コード

GPU

Do i=1,N

CPU

計算の
重い部分
End do

逐次部分はCPU上で実行

並列部分はGPUで実行3

GPUの構造(TESLA K40)

192 CUDA core/SM

SM
Instruction Cache

大量のCUDAコア
並列性の抽出が鍵

Tesla K40, 15 SM/chip

Warp Scheduler
Dispatch Unit

CUDA Core
Dispatch Port

Dispatch Port

Dispatch Unit

Warp Scheduler
Dispatch Unit

Dispatch Unit

Warp Scheduler
Dispatch Unit

Dispatch Unit

Warp Scheduler
Dispatch Unit

Dispatch Unit

Register File (65,536 x 32-bit)

Operand Collector

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

Core Core Core Core Core Core LD/ST SFU

ALUResult Queue

Uniform Cache
64 KB Shared Memory / L1 Cache
Interconnect Network4

GPUアプリケーション

数百のアプリケーションがGPUに対応
www.nvidia.com/object/gpu-applications.html

5アプリをGPU対応する方法
Application

Library

OpenACC

CUDA

GPU対応ライブラリにチェンジ

既存コードにディレクティブを挿入

主要処理をCUDAで記述

簡単に開始

簡単に加速

高い自由度

6

GPU対応のライブラリ (一部)

NVIDIA cuBLAS

NVIDIA cuSPARSEIMSL Library

NVIDIA cuFFT

NVIDIA cuRAND

NVIDIA cuDNN

NVIDIA AmgXGPU Accelerated
Linear Algebra

Vector Signal
Image Processing

Matrix Algebra on
GPU and Multicore

Sparse Linear
AlgebraC++ STL Features
for CUDA
7CUDA計算ライブラリ
高性能な計算ライブラリを提供








cuDNN
cuBLAS
cuFFT
cuRAND
cuSPARSE
cuSOLVER
Thrustディープラーニング向けライブラリ
BLASライブラリ
Fast Fourier Transformsライブラリ
乱数生成ライブラリ
疎行列ライブラリ
Lapackライブラリの一部
並列アルゴリズム C++ STL

CUDAツールキットに標準搭載 (一部、デベロッパー登録必要)
developer.nvidia.com/cuda-downloads
8

アプリをGPU対応する方法
Application

Library

OpenACC

CUDA

GPU対応ライブラリにチェンジ

既存コードにディレクティブを挿入

主要処理をCUDAで記述

簡単に開始

簡単に加速

高い自由度

9

SAXPY (Y=A*X+Y)CPUvoid saxpy(int n, float a,
float *x, float *y){for (int i = 0; i < n; ++i)y[i] += a*x[i];}...
saxpy(N, 3.0, x, y);
...

CUDA
__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx;
if (i < n)y[i] += a*x[i];}...
size_t size = sizeof(float) * N;
cudaMemcpy(d_x, x, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, size, cudaMemcpyHostToDevice);
saxpy<<< N/128, 128 >>>(N, 3.0, d_x, d_y);cudaMemcpy(y, d_y, size, cudaMemcpyDeviceToHost);...

10

アプリをGPU対応する方法
Application

Library

OpenACC

CUDA

GPU対応ライブラリにチェンジ

既存コードにディレクティブを挿入

主要処理をCUDAで記述

簡単に開始

簡単に加速

高い自由度

11

SAXPY (Y=A*X+Y)OpenMP

OpenACC

void saxpy(int n,float a,float *x,
float *restrict y){#pragma omp parallel forfor (int i = 0; i < n; ++i)y[i] += a*x[i];}void saxpy(int n,float a,float *x,
float *restrict y){#pragma acc parallel copy(y[:n]) copyin(x[:n])for (int i = 0; i < n; ++i)y[i] += a*x[i];}...
saxpy(N, 3.0, x, y);...

...
saxpy(N, 3.0, x, y);...
12

CUDAプログラミング プログラミングモデル
 アーキテクチャ
 性能Tips

13

GPUコンピューティングCPUCPU Memory

GPU Memory

GPU

PCI

 高スループット指向のプロセッサ 分離されたメモリ空間14

GPUプログラム
CPUvoid saxpy(int n, float a,
float *x, float *y){for (int i = 0; i < n; ++i)y[i] += a*x[i];}...
saxpy(N, 3.0, x, y);...

GPU__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx;
if (i < n)y[i] += a*x[i];}...
size_t size = sizeof(float) * N;
cudaMemcpy(d_x, x, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, size, cudaMemcpyHostToDevice);
saxpy<<< N/128, 128 >>>(N, 3.0, d_x, d_y);cudaDeviceSynchronize();
cudaMemcpy(y, d_y, size, cudaMemcpyDeviceToHost);...
15

GPU実行の基本的な流れ
CPUGPU
 GPUは、CPUからの制御
で動作

入力データ
転送
GPUカーネル
投入


GPU上で
演算

同期




入力データ: CPUからGPUに
転送 (H2D)
GPUカーネル: CPUから投入
出力データ: GPUからCPUに
転送 (D2H)

出力データ
転送
16

GPUプログラム
GPU

CPU
void saxpy(int n, float a,
float *x, float *y){for (int i = 0; i < n; ++i)y[i] += a*x[i];}__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx;
if (i < n)y[i] += a*x[i];}入力データ転送
カーネル起動
...
saxpy(N, 3.0, x, y);...

同期
出力データ転送...
size_t size = sizeof(float) * N;
cudaMemcpy(d_x, x, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, size, cudaMemcpyHostToDevice);
saxpy<<< N/128, 128 >>>(N, 3.0, d_x, d_y);cudaDeviceSynchronize();
cudaMemcpy(y, d_y, size, cudaMemcpyDeviceToHost);...
17

GPUプログラム (Unified Memory)GPU

CPU
void saxpy(int n, float a,
float *x, float *y){for (int i = 0; i < n; ++i)y[i] += a*x[i];}__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx;
if (i < n)y[i] += a*x[i];}...

カーネル起動
...
saxpy(N, 3.0, x, y);...

同期

saxpy<<< N/128, 128 >>>(N, 3.0, x, y);cudaDeviceSynchronize();...
18

GPUカーネル
GPU

CPU

}

__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx.x;
if (i < n)y[i] += a*x[i];Global スレッドID}...
saxpy(N, 3.0, x, y);...

...
saxpy<<< N/128, 128 >>>(N, 3.0, d_x, d_y);...

void saxpy(int n, float a,
float *x, float *y){for (int i = 0; i < n; ++i)y[i] += a*x[i]; GPUカーネル: 1つのGPUスレッドの処理内容を記述
 基本: 1つのGPUスレッドが、1つの配列要素を担当

19

Execution Configuration(ブロック数とブロックサイズ)
スレッドID__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx.x;
if (i < n)y[i] += a*x[i];}ブロックサイズ

ブロックID...
saxpy<<< N/128, 128 >>>(N, 3.0, d_x, d_y);
...

ブロック数

ブロックサイズ

ブロック数 x ブロックサイズ = 配列要素数20

スレッド階層y[i] = a*x[i] + y[i] (スレッド、ブロック、グリッド)
x[]

0

127

128

255

256

383

384

y[]

0

127

128

255

256

383

384

0

127

128
0

255
127

256
0

383
127

384
0

スレッド
(global)(local)

ブロック0

ブロック2

ブロック1

グリッド
 ブロックサイズ(スレッド数/ブロック)は、カーネル毎に設定可能
 推奨: 128 or 256 スレッド

21

Execution Configuration(ブロック数とブロックサイズ)__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx.x;
if (i < n)y[i] += a*x[i];}...
saxpy<<< N/
N/128,
128
64, 25664 >>>(N, 3.0, d_x, d_y);
N/256,...

ブロック数

ブロックサイズ

ブロック数 x ブロックサイズ = 配列要素数22

2D配列のGPUカーネル例__global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]){Globalスレッド ID (x)int i = threadIdx.x + blodkDim.x * blockIdx.x;
int j = threadIdx.y + blodkDim.y * blockIdy.y;if ( i < N && j < N )GlobalスレッドID (y)C[i][i] = A[i][j] + B[i][j];}ブロックサイズ (x,y)
...
4
32,
8 );
dim3 sizeBlock( 64,
16, 16
dim3 numBlocks( N/sizeBlock.x, N/sizeBlock.y );MatAdd<<< numBlocks, sizeBlock >>>(A, B, C);...

ブロック数 (x,y)

 ブロックサイズ(ブロック形状)は、1D～3Dで表現可能
23

ブロック・マッピング、スレッド・マッピング
dim3 sizeBlock(16,16)
(0,0)

(0,0)

(1,0)

(0,0)

dim3 sizeBlock(32,8)

(0,2)

(1,1)

(1,2)

(1,0)

(0,1)

(1,1)

(0,2)

(1,2)

(0,3)

(1,3)

(0,4)

(1,4)

(2,0)

(15,15)
(0,1)

(0,0)

(31,7)

(2,1)

(2,2)

ブロックID(blockIdx)

スレッドID(threadIdx)

24

CUDAプログラミング プログラミングモデル
 アーキテクチャ
 性能Tips28

GPUアーキテクチャ概要
 PCI I/F
 ホスト接続インタフェース

 Giga Thread Engine

 SMに処理を割り振るスケジューラ

 DRAM I/F (GDDR5, 384-bit)
 全SM、PCI I/Fからアクセス可能なメモ
リ (デバイスメモリ, フレームバッファ) L2 cache (3MB)

 全SMからアクセス可能なR/Wキャッシュ

Maxwell GM200

 SM (Streaming Multiprocessor)

 「並列」プロセッサ、Maxwell:最多24
29

SM(Stream Multi-Processor)
 CUDAコア
 GPUスレッドはこの上で動作
 Maxwell: 128個

 Other units LD/ST, SFU, etc レジスタ(32bit): 64K個
 共有メモリ: 96KB
 Tex/L1キャッシュ
Maxwell
GM200

30

GPUカーネル実行の流れ CPUが、GPUに、グリッドを投入

 具体的な投入先は、Giga Thread Engine
グリッド
ブロック
スレッド
スレッド

ブロック
31

GPUカーネル実行の流れ
 Giga Thread Engine(GTE)が、SMに、ブロックを投入
 GTEは、ブロックスケジューラ
 グリッドをブロックに分解して、ブロックを、空いているSMに割当てる
グリッド
ブロック
スレッド
スレッド

ブロック
32

ブロックをSMに割り当て
 各ブロックは、互いに独立に実行

 ブロック間では同期しない、実行順序の保証なし 1つのブロックは複数SMにまたがらない 1つのSMに、複数ブロックが割当てられることはある
グリッド
ブロック
ブロック
ブロック
ブロック

33

GPUカーネル実行の流れ
 SM内のスケジューラが、スレッドをCUDAコアに投入

グリッド
ブロック
スレッド
スレッド

ブロック
34

GPUカーネル実行の流れ
 SM内のスケジューラが、ワープをCUDAコアに投入

 ワープ: 32スレッドの塊
 ブロックをワープに分割、実行可能なワープを、空CUDAコアに割当てる
グリッド
ブロック
ワープ
ワープ
ブロック
35

ワープのCUDAコアへの割り当て
 ワープ内の32スレッドは、同じ命令を同期して実行
SIMT(Single Instruction Multiple Threads)
 各ワープは、互いに独立して実行 同じブロック内のワープは、明示的に同期可能(__syncthreads())
グリッド
ブロック
ワープ
ワープ
ブロック

ワープ
36

GPUアーキの変化を問題としないプログラミングモデル

Pascal, CC 6.0
64 cores /SM
Kepler, CC 3.5
192 cores /SM

Maxwell, CC 5.0
128 cores /SM

37

CUDAプログラミング プログラミングモデル
 アーキテクチャ
 性能Tips

42

リソース使用率 (Occupancy)
SMの利用効率を上げる≒ SMに割当て可能なスレッド数を、上限に近づける

 レジスタ使用量(/スレッド)

CUDAコア数: 128
最大スレッド数: 2048
最大ブロック数: 32
共有メモリ: 96KB
レジスタ数(32-bit): 64K個 できる限り減らす
リソース量/SM (GM200)
 DP(64bit)は、2レジスタ消費
 レジスタ割り当て単位は8個
 レジスタ使用量と、割当て可能なスレッド数の関係
 32レジスタ: 2048(100%), 64レジスタ: 1024(50%)
 128レジスタ:512(25%), 256レジスタ: 256(12.5%)43リソース使用率 (Occupancy)
SMの利用効率を上げる≒ SMに割当て可能なスレッド数を、上限に近づける スレッド数(/ブロック) 64以上にする
 64未満だと最大ブロック数がネックになる

CUDAコア数: 128
最大スレッド数: 2048
最大ブロック数: 32
共有メモリ: 96KB
レジスタ数(32-bit): 64K個
リソース量/SM (GM200)

 共有メモリ使用量(/ブロック) できる限り減らす
 共有メモリ使用量と、割当て可能なブロック数の関係
48KB:2ブロック, 12KB:8ブロック, 3KB:32ブロック

44

リソース使用率 (Occupancy)
空き時間を埋める CUDAストリーム (≒キュー)(*) 操作: GPUカーネル, データ転送

 同じCUDAストリームに投入した操作: 投入順に実行
 別のCUDAストリームに投入した操作: 非同期に実行 (オーバラップ実行)[CUDAストリームの効果例]
GPUカーネルとデータ転送がオーバラップして同時に実行されている
45

デバイスメモリへのアクセスは、まとめて
 コアレス・アクセス

 32スレッド(ワープ)のロード/ストアをまとめて、メモリトランザクションを発行
 トランザクションサイズ: 32B, 64B or 128B

 トランザクション数は、少ないほど良い
配列
128B x1

Best
Good
Bad

128B境界0 1

0

128B境界28 29 30 31

0 1 2 314 15 16 17

1

30 31

64B x232B x322
46

分岐を減らす ワープ内のスレッドが別パスを選択すると遅くなる

 ワープ内のスレッドは、命令を共有 (SIMT)
 ワープ内のスレッドが選んだ全パスの命令を実行 あるパスの命令を実行中、そのパスにいないスレッドはinactive状態

 Path divergenceを減らす

 できる限り、同ワープ内のスレッドは同じパスを選択させる

A[id] > 0

true
0

処理 X

false

1
2

3

処理 Y
47

OPENACCプログラミング 概要紹介 プログラムのOpenACC化
 OpenACC化事例

48

アプリをGPU対応する方法
Application
Library

OpenACC

CUDA

GPU対応ライブラリにチェンジ

既存コードにディレクティブを挿入

主要処理をCUDAで記述

簡単に開始

簡単に加速

高い自由度

49

OPENACCCPUGPU

簡単: 既存のコードにコンパイラへのヒントを追加

Program myscience... serial code ...
!$acc kernels
do k = 1,n1
do i = 1,n2... parallel code ...
enddo
enddo
!$acc end kernels... serial code …
End Program myscience

既存のC/Fortranコード

強力: そこそこの労力で、コンパイラがヒントの追加

コードを自動で並列化

オープン: 複数コンパイラベンダが、
複数アクセラレータをサポート
NVIDIA, AMD, Intel(予定)50アプリケーション実行
アプリケーション・コード

GPU

Do i=1,N

CPU

計算の
重い部分
End do

逐次部分はCPU上で実行

並列部分はGPUで実行51

SAXPY (Y=A*X+Y)CPUvoid saxpy(int n, float a,
float *x, float *y){for (int i = 0; i < n; ++i)y[i] += a*x[i];}...
saxpy(N, 3.0, x, y);
...

CUDA
__global__ void saxpy(int n, float a,
float *x, float *y){
int i = threadIdx.x + blodkDim.x * blockIdx;
if (i < n)y[i] += a*x[i];}...
size_t size = sizeof(float) * N;
cudaMemcpy(d_x, x, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, size, cudaMemcpyHostToDevice);
saxpy<<< N/128, 128 >>>(N, 3.0, d_x, d_y);cudaMemcpy(y, d_y, size, cudaMemcpyDeviceToHost);...

52

SAXPY (y=a*x+y)
OpenMP

OpenACC

void saxpy(int n,float a,float *x,
float *restrict y){#pragma omp parallel forfor (int i = 0; i < n; ++i)y[i] += a*x[i];}void saxpy(int n,float a,float *x,
float *restrict y){#pragma acc parallel copy(y[:n]) copyin(x[:n])for (int i = 0; i < n; ++i)y[i] += a*x[i];}...
saxpy(N, 3.0, x, y);...

...
saxpy(N, 3.0, x, y);...

 omp  acc

 データの移動

53

SAXPY (y=a*x+y, FORTRAN)
OpenMP

OpenACC

subroutine saxpy(n, a, X, Y)
real :: a, X(:), Y(:)
integer :: n, i

subroutine saxpy(n, a, X, Y)
real :: a, Y(:), Y(:)
integer :: n, i!$omp parallel dodo i=1,nY(i) = a*X(i)+Y(i)enddo!$omp end parallel doend subroutine saxpy!$acc parallel copy(Y(:)) copyin(X(:))
do i=1,nY(i) = a*X(i)+Y(i)enddo!$acc end parallel
end subroutine saxpy...
call saxpy(N, 3.0, x, y)...

...
call saxpy(N, 3.0, x, y)...

 FORTRANも同様
54

簡単にコンパイル
OpenMP / OpenACC
void saxpy(int n, float a,
float *x,
float *restrict y){ -acc saxpy.c$ pgcc -Minfo#pragma acc parallel copy(y[:n]) copyin(x[:n])saxpy:#pragmaomp parallel for
16, Generatingpresent_or_copy(y[:n])for present_or_copyin(x[:n])(int i = 0; i < n; ++i)
Generating+= codea*x[i];
Generatingy[i]Tesla} parallelizable
19, Loop is
Accelerator kernel generated...
19, #pragma
acc loop gang, vector(128) /* blockIdx.x threadIdx.x */saxpy(N, 3.0, x, y);...
55

簡単に実行OpenMP / OpenACC
void saxpy(int n, float a,
float *x,
float *restrict y){ -acc saxpy.c$ pgcc -Minfo$ nvprof ./a.out#pragma acc parallel copy(y[:n]) copyin(x[:n])saxpy:==10302== NVPROF
is profilingprocess#pragmaomp parallelfor 10302, command: ./a.out
16, Generatingpresent_or_copy(y[:n])==10302== Profiling
application:
for present_or_copyin(x[:n])(inti = 0; i ./a.out< n; ++i)
Generating==10302== Profiling
result:+= codea*x[i];
Generatingy[i]TeslaTime(%)
Time
Calls
Avg
Min
Max Name} parallelizable
19, Loop is
62.95% 3.0358ms
2 1.5179ms 1.5172ms 1.5186ms [CUDA memcpy HtoD]
Accelerator kernel generated
31.48% 1.5181ms
1 1.5181ms 1.5181ms 1.5181ms [CUDA memcpy DtoH]...
19, #pragma
acc loop gang, vector(128) /* blockIdx.x threadIdx.x */5.56% 268.31us268.31us 268.31us saxpy_19_gpu
saxpy(N, 3.0,1 x,268.31us
y);
...
56

OPENACCプログラミング 概要紹介 プログラムのOpenACC化
 OpenACC化事例

57

事例: ヤコビ反復法
A(i,j+1)

while ( error > tol ) {error = 0.0;
for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]));}}A(i-1,j)

A(i,j) A(i+1,j)
A(i,j-1)

for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];
}
}
}
58

並列領域 (kernels construct)while ( error > tol ) {error = 0.0;
#pragma acc kernelsfor (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}#pragma acc kernelsfor (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];}}} 並列領域の指定
 Parallels とKernels Parallels
 OpenMPと親和性
 開発者主体 Kernels

 複数kernelの生成 コンパイラ主体
59

[PGI tips] コンパイラメッセージ$ pgcc –acc –Minfo=accel jacobi.cjacobi:44, Generating copyout(Anew[1:4094][1:4094])
Generating copyin(A[:][:])
Generating Tesla code
45, Loop is parallelizable
46, Loop is parallelizable
Accelerator kernel generated
45, #pragma acc loop gang /* blockIdx.y */46, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */49, Max reduction generated for error

60

並列領域 (kernels construct)while ( error > tol ) {error = 0.0;
#pragma acc kernelsfor (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;$ pgcc-acc jacobi.cerror -Minfo=acc= max(error, abs(Anew[j][i]- A[j][i]);}jacobi:}

 並列領域の指定
 Parallels とKernels Parallels
 OpenMPと親和性

59, Generating present_or_copyout(Anew[1:4094][1:4094])
 開発者主体
Generatingpresent_or_copyin(A[:][:])#pragma acc
kernelsfor (int jGenerating= 1; j < N-1;j++) code{Tesla
for 61,(intLoopi = 1;< M-1; i++) {is iparallelizable
A[j][i] = Anew[j][i]; 複数kernelの生成} 63, Loop is parallelizable
Accelerator kernel generated}
 コンパイラ主体}
61, #pragma acc loop gang /* blockIdx.y */63, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */61
Max reduction generated for error

 Kernels

データ転送(data clause)
while ( error > tol ) {error = 0.0;
#pragma acc kernelsfor (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;$ pgcc-acc jacobi.cerror -Minfo=acc= max(error, abs(Anew[j][i]- A[j][i]);}jacobi:}

 並列領域の指定
 Parallels とKernels Parallels
 OpenMPと親和性

59, Generating present_or_copyout(Anew[1:4094][1:4094])
 開発者主体
Generatingpresent_or_copyin(A[:][:])#pragma acc
kernelsfor (int jGenerating= 1; j < N-1;j++) code{Tesla
for 61,(intLoopi = 1;< M-1; i++) {is iparallelizable
A[j][i] = Anew[j][i]; 複数kernelの生成} 63, Loop is parallelizable
Accelerator kernel generated}
 コンパイラ主体}
61, #pragma acc loop gang /* blockIdx.y */63, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */62
Max reduction generated for error

 Kernels

データ転送 (DATA CLAUSE)
while ( error > tol ) {error = 0.0;
#pragma acc kernels \pcopyout(Anew[1:4094][1:4094]) pcopyin(A[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}#pragma acc kernels \pcopyout(A[1:4094][1:4094]) pcopyin(Anew[1:4094][1:4094])

for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];
}
}







copyin (HostGPU)
copyout (HostGPU)
copy
createpresent

pcopyinpcopyout
pcopy
pcreate}63

配列形状の指定
 配列は、全要素でなく、一部だけ指定して転送することも可能 注意: C/C++ と Fortranでは指定方法が異なる

 C/C++: array[ start : size ]float Anew[4096][4096]
pcopyout( Anew[1:4094][1:4094]) pcopyin( A[:][:]) ) Fortran: array( start : end )real Anew(4096,4096)
pcopyout( Anew(2:4095, 2:4095) ) pcopyin( A(:,:) )64

データ転送 (data clause)
while ( error > tol ) {error = 0.0;
#pragma acc kernels \pcopy(Anew[:][:]) pcopyin(A[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}#pragma acc kernels \pcopy(A[:][:]) pcopyin(Anew[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];
}
}







copyin (HostGPU)
copyout (HostGPU)
copy
createpresent

pcopyinpcopyout
pcopy
pcreate}65

[PGI tips] PGI_ACC_TIME$ PGI_ACC_TIME=1 ./a.out
Accelerator Kernel Timing data

統計データ

/home/anaruse/src/OpenACC/jacobi/C/task1-solution/jacobi.cjacobiNVIDIAdevicenum=0time(us): 649,886
44: data region reached 200 times
44: data copyin transfers: 800
device time(us): total=14,048 max=41 min=15 avg=1753: data copyout transfers: 800
device time(us): total=11,731 max=43 min=6 avg=14
44: compute region reached 200 times
46: kernel launched 200 times
grid: [32x4094]

block: [128]device time(us): total=382,798 max=1,918 min=1,911 avg=1,913
elapsed time(us): total=391,408 max=1,972 min=1,953 avg=1,957
46: reduction kernel launched 200 times
grid: [1]block: [256]device time(us): total=48,235 max=242 min=241 avg=241
elapsed time(us): total=53,510 max=280 min=266 avg=267

66

[PGI tips] PGI_ACC_NOTIFY$ PGI_ACC_NOTIFY=3 ./a.out

トレースデータ...
upload CUDA data file=/home/anaruse/src/OpenACC/jacobi/C/task1solution/jacobi.c function=jacobi line=44 device=0 variable=A bytes=16777216...
launch CUDA kernel file=/home/anaruse/src/OpenACC/jacobi/C/task1solution/jacobi.c function=jacobi line=46 device=0 num_gangs=131008
num_workers=1 vector_length=128 grid=32x4094 block=128 shared memory=1024...
download CUDA data file=/home/anaruse/src/OpenACC/jacobi/C/task1solution/jacobi.c function=jacobi line=53 device=0 variable=Anew bytes=16736272...

67

NVIDIA VISUAL PROFILER (NVVP)

68

NVVPによる解析: データ転送がボトルネック

利用率:低い1 cycleGPU
kernelGPUkernel69

過剰なデータ転送
while ( error > tol ) {error = 0.0;
#pragma acc kernels \pcopy(Anew[:][:]) pcopyin(A[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}#pragma acc kernels \pcopy(A[:][:]) pcopyin(Anew[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];
}
}
}

70

Host

過剰なデータ転送GPU

while ( error > tol ) {error = 0.0;}

#pragma acc kernels \pcopy(Anew[:][:]) \pcopyin(A[:][:]){

copyin}copyout#pragma acc kernels \pcopy(A[:][:]) \pcopyin(Anew[:][:]){

copyin}copyout#pragma acc loop reduction(max:error)for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];}}71データ領域 (data construct)#pragma acc data pcopy(A, Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}#pragma acc kernels pcopy(A[:][:]) pcopyin(Anew[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];
}
}







copyin (CPUGPU)
copyout (CPUGPU)
copy
createpresent

pcopyinpcopyout
pcopy
pcreate}72

データ領域 (data construct)#pragma acc data pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}#pragma acc kernels pcopy(A[:][:]) pcopyin(Anew[:][:])for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];
}
}







copyin (CPUGPU)
copyout (CPUGPU)
copy
createpresent

pcopyinpcopyout
pcopy
pcreate

}
73

Host
#pragma acc data \pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;

適正なデータ転送GPUcopyin#pragma acc kernels \pcopy(Anew[:][:]) \pcopyin(A[:][:]){for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);}}}#pragma acc kernels \pcopy(A[:][:]) \pcopyin(Anew[:][:]){for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {A[j][i] = Anew[j][i];
}
}

}
}

copyout74

データ転送が減少 (NVVP)
利用率:高い1 cycle75

2つの処理CPU Memory

データ転送

GPU Memory

PCI

計算オフロード

計算オフロード、データ転送、両方を考慮する必要がある

76

その他のデータ管理方法
float *array;
Init( ) {...
array = (float*)malloc( … );input_array( array );#pragma enter data copyin(array)...
}Fin( ) {...
#pragma exit data copyout(array)
output_array( array );free( array );...
} Enter data Copyin Create Pcopyin Pcreate Exit data Copyout
 Delete77その他のデータ管理方法#pragma acc data pcopy(A,B)
for (k=0; k<LOOP; k++) {#pragma acc kernels present(A,B)
for (i=0; i<N; i++) {A[i] = subA(i,A,B);}#pragma acc update self(A[0:1])

 Update self
CPU  GPU Update device
CPU  GPUoutput[k] = A[0];A[N-1] = input[k];#pragma acc update device(A[N-1:1])#pragma acc kernels present(A,B)
for (i=0; i<N; i++) {B[i] = subB(i,A,B);}}
78

リダクション(縮約計算)
while ( error > tol ) {error = 0.0;
#pragma acc kernelsfor (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);
}
}
...
}

79

リダクション(縮約計算)
while ( error > tol ) { 演算の種類+和#pragma acc kernels* for (int積j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Max Anew[j][i]最大 = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;Min
最小$ pgcc -Minfo=acc-acc
error = max(error, jacobi.cabs(Anew[j][i] - A[j][i]);}|jacobi:
ビット和} 59, Generating present_or_copyout(Anew[1:4094][1:4094])& #pragmaGenerating
ビット積
present_or_copyin(A[:][:])
acc kernels
Generating
code{
(int
j = 1; j <Tesla
N-1; j++)
^ for
XOR
61,
Loop iis
parallelizablefor (int= 1;
i < M-1; i++) {= Anew[j][i];Loopisparallelizable|| 63,A[j][i]論理和}
Accelerator kernel generated&& }論理積
61, #pragma acc loop gang /* blockIdx.y */error = 0.0;}

63, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */Max reduction generated for error
80リダクション (REDUCTION CLAUSE)
while ( error > tol ) {error = 0.0;
#pragma acc kernels#pragma acc loop reduction(max:error)for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);
}
}
...
}

 演算種類(C/C++)
+
*
max
min
|
&
^
||
&&

和
積
最大
最小
ビット和
ビット積
XOR
論理和
論理積
81

並列方法の指示#pragma acc data pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])#pragma acc loop reduction(max:error)for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);
}
}
...
}

82

並列方法の指示#pragma acc data pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])#pragma acc loop reduction(max:error)for (int j = 1; j < N-1; j++) {$ pgcc-Minfo=acc -acc jacobi.c#pragma acc loop reduction(max:error)
jacobi:for (int i = 1; i < M-1; i++) {59,
Generatingpresent_or_copyout(Anew[1:4094][1:4094])
Anew[j][i]= (A[j][i+1]+ A[j][i-1] ++ A[j+1][i]) * 0.25;
GeneratingA[j-1][i]present_or_copyin(A[:][:])
error = max(error, abs(Anew[j][i] - A[j][i]);Generating Tesla code}
} 61, Loop is parallelizable
63, Loop is parallelizable...
}
Accelerator kernel generated

61, #pragma acc loop gang /* blockIdx.y */63, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */Max reduction generated for error
83並列方法の指示 (loop construct)#pragma acc data pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])#pragma acc loop gang vector(1) reduction(max:error)for (int j = 1; j < N-1; j++) {#pragma acc loop gang vector(128)for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);
}
}





Gang
WorkerVector … SIMD幅







Collapse
Independent
Seq
Cache
Tile

...
}

84

実行条件設定 (gang, vector)#pragma acc loop gang vector(4)for (j = 0; j < 16; j++) {#pragma acc loop gang vector(16)for (i = 0; i < 16; i++) {...

#pragma acc loop gang vector(8)for (j = 1; j < 16; j++) {#pragma acc loop gang vector(8)for (i = 0; i < 16; i++) {
...

i

i

4 x 164 x 16j

4 x 16

j

8x8

8x8

8x8

8x8

85

ループを融合 (collapse)#pragma acc data pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])#pragma acc loop reduction(max:error) \collapse(2) gang vector(128)for (int j = 1; j < N-1; j++) {for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);
}
}
...
}





Gang
WorkerVector … SIMD幅








Collapse
Independent
Seq
Cache
Tile
...

86

ループを融合 (collapse)#pragma acc data pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])#pragma acc loop reduction(max:error) gang vector(128)for (int ji = 0; ji < (N-2)*(M-2); ji++) {j = (ji / (M-2)) + 1;
i = (ji % (M-2)) + 1;
Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);
}
...
}





Gang
WorkerVector … SIMD幅








Collapse
Independent
Seq
Cache
Tile
...

87

並列実行可能(independent)#pragma acc data pcopy(A) create(Anew)
while ( error > tol ) {error = 0.0;
#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])#pragma acc loop reduction(max:error) independentfor (int jj = 1; jj < NN-1; jj++) {int j = list_j[jj];for (int i = 1; i < M-1; i++) {Anew[j][i] = (A[j][i+1] + A[j][i-1] +A[j-1][i] + A[j+1][i]) * 0.25;
error = max(error, abs(Anew[j][i] - A[j][i]);
}
}
...
}





Gang
WorkerVector … SIMD幅








Collapse
Independent
Seq
Cache
Tile
...

88

逐次に実行 (seq)#pragma acc kernels pcopy(Anew[:][:]) pcopyin(A[:][:])#pragma acc loop seq
for (int k = 3; k < NK-3; k++) {#pragma acc loopfor (int j = 0; j < NJ; j++) {#pragma acc loopfor (int i = 0; i < NI; i++) {Anew[k][j][i] = func(
A[k-1][j][i], A[k-2][j][i], A[k-3][j][i],A[k+1][j][i], A[k+2][j][i], A[k+3][j][i], ...
);
}
}
}





Gang
WorkerVector … SIMD幅








Collapse
Independent
Seq
Cache
Tile
...

89

OPENACCプログラミング 概要紹介 プログラムのOpenACC化

 OpenACC化事例

90

NICAM 気象・気候モデル by 理研AICS/東大
 膨大なコード (数十万行)
 ホットスポットがない (パレートの法則)

 特性の異なる2種類の処理

 力学系 … メモリバンド幅ネック
 物理系 … 演算ネック

91

NICAM: 力学系(NICAM-DC)
Tsubame 2.5 (GPU:K20X)

 OpenACCによるGPU化 良好なスケーラビリティ

 Tsubame 2.5, 最大2560
GPUs
 Scaling factor: 0.8

1.E+05
Performance (GFLOPS)

 主要サブルーチンは、全てGPU上で動作(50以上)
 MPI対応済み
 2週間

K computer(*) weak scaling

Tsubame 2.5 (CPU:WSM)

1.E+04

1.E+03

1.E+02

1.E+01
1.E+00

Weak scaling

1.E+01
1.E+02
1.E+03
Number of CPUs or GPUs

1.E+0492

Courtesy of Dr. Yashiro from RIKEN AICSNICAM: 力学系(NICAM-DC)
Measured Performance(GFLOPS)

1.E+05

Tsubame 2.5 (GPU:K20X)
K computer

1.E+04

Tsubame 2.5 (CPU:WSM)

1.E+03
1.E+02
1.E+01
1.E+02

1.E+03
1.E+04
1.E+05Aggregate Peak Memory Bandwidth (GB/s)

1.E+0693

Courtesy of Dr. Yashiro from RIKEN AICSNICAM: 物理系(SCALE-LES) Atmospheric radiation transfer
 物理系の中で、最も重い計算
 OpenACCによるGPU対応、完了
Speedup
vs. CPU 1-core

Better

160
140
120
100
80
60
40
20
0

151

76.0
37.8
1.00

1.99

3.88

8.51

1 core2 core4 core10 coreXeon E5-2690v2(3.0GHz,10-core)

1 GPU

2 GPUs

4 GPUsTesla K40

94(*) PCIデータ転送時間込み, グリッドサイズ:1256x32x32

SEISM3D
 地震シミュレーション by 古村教授(東大地震研)
 主要サブルーチンのGPU対応が完了

 メモリバンド幅ネック、 3次元モデル(2次元分割)、隣接プロセス間通信
SEISM3D (480x480x1024, 1K steps)
600

Time (sec)

500

Better

605
459

400

3.4x speedup(アプリ全体)300
200

GPUの実行時間内訳
140
Others (CPU, MPI and so on)120[CUDA memcpy DtoH]

100[CUDA memcpy HtoD]

80(other subroutines)
update_vel_pml

60

100

134

0
K: 8x SPARC64
VIIIfxCPU: 8x Xeon
E5-2690v2

GPU: 8x Tesla
K40

update_vel

40

update_stress_pml
update_stress20

diff3d_*0
GPU: 8x Tesla K40

95

FFR/BCM (仮称)次世代CFDコード by 坪倉准教授(理研AICS/北大)MUSCL_bench: MUSCLスキームに基づくFlux計算 (とても複雑な計算)
 CFD計算の主要部分 (60-70%)
 OpenACCによるGPU対応、完了
BetterSpeedupvs. 1 CPU core




35
30
25
20
15
10
5
0

33.21

1.00

1.93

1 core2 core

4.55
5 coreXeon E5-2690v2(3.0GHz,10-core)

8.30

10 core1 GPU
Tesla K40

96(*) PCIデータ転送時間込み、サイズ:80x32x32x32

CM-RCM IC-CG (PRELIMINARY) IC-CG法のベンチマークコード by 中島教授(東大) CM-RCM法(Cyclic Multi-coloring Reverse Cuthill-Mckee)を使用

 メインループ内のサブルーチンを全てOpenACCでGPU化

Time (second)4

Better

CM-RCM ICCG (100x100x100)
3.36

3
2

1.58

1.53

1.39

1

0.655

0
OpenMP

OpenMP

OpenMP

OpenMP

OpenACC

Opteron 6386SE(2.8GHz,16core)

SPARC64 Ixfx
(1.85GHz,16core)

Xeon E5-2680v2(2.8GHz,10core)

Xeon-Phi 5110P

Tesla K40
97

Courtesy of Dr. Ohshima from U-TokyoCCS-QCD QCDコード by 石川准教授(広島大) BiCGStab計算を全てOpenACCでGPU化
 データレイアウトを変更: AoS  SoA

CCS-QCD: BiCGStab Total FLOPSGFLOPS

Better

90
80
70
60
50
40
30
20
10
0

Xeon E5-2690v2(3.0GHz,10core,OpenMP)

63.4

60.8

57.9

53.3

Tesla K40(OpenACC)

42.5
32.3

8x8x8x32

24.3

8x8x8x64

22.1

8x8x16x64
Problem Size

20.9

8x16x16x64

19.7

16x16x16x64
98