An	  Introduction	
  to

CUDA	
  ProgrammingT.	
  Aoyama	  HPC-­‐Phys	
  workshop	
  at	
  Hiroshima	
  Univ.,	
  March	
  18th,	  2019

Table	  of	  Contents

CF• Notes	  on	
  Hands-­‐on• Introduction	
  to	
  CUDA	  • Himeno	
  benchmark	
  • Programming	  tutorial	
  • Kernels	  • Memory	
  management	  • Launch	
  kernel	  • Compile	
  and	  Run
• Advanced	
  topics2

Hands-­‐on• 実習用ワークステーション	  • sim01,	
  sim02	
  • NVIDIA	  GeForce	
  GTX	
  680M	  (Kepler	
  architecture)	
  • CUDA	  Toolkit	
  10.0	  • PGI	  Compiler	
  Community	
  Edition	
  18.4	  • sim03	  • NVIDIA	  GeForce	
  GTX	  1080	
  (Pascal	
  architecture)	
  • CUDA	  Toolkit	
  10.1	
  • PGI	  Compiler	
  Community	
  Edition	
  18.103

Hands-­‐on• 実習用プログラム:	  	
  /home/aoym/workshop/	  • himeno_C/	
  	
  	
  	
  	
  	
  himenoBMTxpa_mod2.c
himeno_F90/	  	
  himenoBMTxp_mod.f90	
  • Host用コード。これをベースに書き換え。	
  • 各自ログイン後にコピー、コンパイル・実行してみる	  • オリジナルから変更:	
  	
  • offload	
  用に	
  jacobi2()	
  を用意	  • i,j,k	  ループ内の処理を変更:	
  差分計算・更新・誤差評価	
  • himeno_C_CUDA/
himeno_F90_CUDA/

	
  ——	
  C	
  /	
  Fortran90	
  での実装例	
  • himeno_C_OpenACC/
himeno_F77_OpenACC/

	
  	
  	
  ——OpenACC	
  での実装例(松古氏による)

4

CUDA	
  Tutorial
5GPU	
  Computing• GPUシステムの構成
ホスト(CPU)CPU	  CoreCPU

Host	
  Memory

デバイス(GPU)

Streaming
Multiprocessor
PCI	
  Express	
  
or	
  
NVLINK

Shared	  MemoryCUDA	  CoreO(103)

Device	  Memory

6

GPU	
  Computing• NVIDIA	  GPUs
GTX680M

GTX1080

P100

V100

Generation

KeplerPascalPascal

Volta

#SM

7

20

56

80

FP32	
  core

1344

2560

3584

5120

FP64	
  core

56

80

1792

2560

Peak	
  FP32

2.03TFlops

8.9TFlops

10.6TFlops

15.7TFlopsPeak	
  FP640.083TFlops

0.26TFlops

5.3TFlops

7.8TFlops

SP/DP

24:1

32:1

2:1

2:1

Mem	
  B/W

115.2GB/s

320GB/s

732GB/s

900GB/s

CC

3.0

6.1

6.0

7.0

Compute	  Capability

7

GPU	
  Computing• GPUを効果的に使う	
  program• 計算時間が集中する箇所(Hot	  Spot)
その部分をGPUにオフロードする	
  • 多数のスレッドによる並列実行	  • 独立な計算に分割可能な処理	
  • データをホスト(CPU)とデバイス(GPU)Hot	
  SpotGPU
offloadの間でデータ転送する必要がある	
  • オーバーヘッドを考慮
8GPU	
  Computing• GPGPU	
  =	
  General-­‐Purpose	
  computing	  on	
  Graphic	
  Processing	  Unit	  • GPUを汎用の計算に利用	
  • Programming	  environments:	  • Directive-­‐based	
  • OpenACC,	
  OpenMP	
  (4.5-­‐),	  …	
  • 既存コードからincrementalにGPU化	  • API-­‐based	
  • CUDA,	
  OpenCL,	
  …	
  • GPUでの処理を細かく制御	  • Libraries	  • GPU向けに移植されたLibraryやPackageを利用	
  • cuBLAS,	
  cuFFT,	
  magma,	
  …	
  

9

⚫

device ( list ) host → device のデータ転送⚫

host ( list ) device → host のデータ転送。self clause は host と同じ。

⚫

async ( expression ) parallel 構文の場合と同様、データ転送を非同期にできる。Parallel 構文によるデバイスカーネルの実行とも非同期に実行可能。wait 構文で同期をとる。

GPU	
  Computing:	
  OpenACC
C 言語コードによる例 (based on Himeno benchmark)#pragma acc data copy(pp[0:size]), ¥copyin(size, mrows, mcols, mdeps, omega, pwrk1[0:size], ¥pbnd[0:size], pa[0:4*size], pb[0:3*size], pc[0:3*size]), ¥create(pwrk2[0:size])// 転送する配列データは起点とサイズで指定する{
for(n=0 ; n<nn ; n++){gosa = 0.0;

GPUにoffloadする
部分をdirectiveで指示// 構造化ブロックを囲む#pragma acc parallel num_workers(1) vector_length(256){// 構造化ブロックを囲む#pragma acc loop gang worker vector reduction(+ : gosa)
for(int site = 0 ; site < size; site++){// 3-level のループを 1 次元化済み
int k = site % mdeps;if(k == 0 || k == mdeps-1) continue;… (略)float s0 = pa[IDX(0,i,j,k)]*pp[IDX(0,i+1,j, k)]+ pa[IDX(1,i,j,k)]*pp[IDX(0,i, j+1,k)]… (略)float ss = …
gosa += ss*ss;
pwrk2[IDX(0,i,j,k)] = pp[IDX(0,i,j,k)] + omega * ss;}} //acc parallel// thread 間の同期
#pragma acc parallel num_workers(2) vector_length(32){#pragma acc loop gang worker vector
for(int site = 0 ; site < size; site++){… (略)}} // acc parallel// thread 間の同期
} /* end n loop */} // acc data// copy, copyout のデータがホストに送られる。メモリ領域解放

10

GPU	
  Computing
CUDA	
  =	
  Compute	
  Unified	
  Device	
  Architecture• NVIDIAが開発・提供する	
  GPGPUプラットフォーム
およびプログラミングモデル	
  • C/C++	  	
  (NVIDIA	
  CUDA	  SDK)	
  • C/C++	
  の拡張として導入。CUDA	
  C	
  と	
  API	
  が提供される。	
  • Fortran	
  (PGI)	
  • PGI	
  compiler	
  でサポート。Fortran90	
  の拡張として導入。
11

GPU	
  Computing• 一般的な処理の流れ	  1. デバイスを選択・初期化	  API呼び出しmalloc,	
  memcpy	
  に相当

2. デバイス上にメモリ領域を確保	
  3. データをホストからデバイスに転送	
  カーネル(GPUスレッドでの

4. カーネルを起動して計算を行う	
  処理)を記述5. 計算結果をデバイスからホストへ転送	
  拡張構文を用いて

6. デバイス上のメモリ領域を解放	
  カーネルを起動

7. デバイスを解放

12

CUDA	
  Tutorial• 今回は	
  Himeno	
  benchmark	
  コードを題材にお手軽にCUDAを導入してみる。	
  • シンプルにカーネル化	  • shared	  memory	
  などの	
  optimize	
  は後回し	  • ライブラリの利用	  • 集約処理(reduction)に	
  cuBLAS	
  を使う
姫野ベンチマーク	  • 理研情報基盤センター・姫野龍太郎氏による	  • http://accc.riken.jp/supercom/documents/himenobmt/

13

Himeno	  benchmark• 非圧縮性流体などに現れるPoisson方程式をPoint-­‐Jacobi法で解く	
  @@22pp @@22pp @@22pp@@22pp@@22pp@@2 p2 p
+
+
+↵
+
+
=⇢
2 + @y 2 + @z 2 + ↵@x@y + @x@z + @y@z = ⇢
@x
2
2
2
@x
@y
@z
@x@y
@x@z
@y@z

2次中心差分で離散化
A0 pi+1,j,k + A1 pi,j+1,k + A2 pi,j,k+1
A0 pi+1,j,k + A1 pi,j+1,k + A2 pi,j,k+1+B0 (pi+1,j+1,k pi 1,j+1,k pi+1,j 1,k + pi 1,j 1,k )+B0 (pi+1,j+1,k pi 1,j+1,k pi+1,j 1,k + pi 1,j 1,k )+B1 (pi+1,j,k+1 pi+1,j,k 1 pi 1,j,k+1 + pi 1,j,k 1 )+B1 (pi+1,j,k+1 pi+1,j,k 1 pi 1,j,k+1 + pi 1,j,k 1 )+B2 (pi,j+1,k+1 pi,j+1,k 1 pi,j 1,k+1 + pi,j 1,k 1 )+B2 (pi,j+1,k+1 pi,j+1,k 1 pi,j 1,k+1 + pi,j 1,k 1 )+C0 pi 1,j,k + C1 pi,j 1,k + C2 pi,j,k 1 Dpi,j,k = W+C0 pi 1,j,k + C1 pi,j 1,k + C2 pi,j,k 1 Dpi,j,k

=W

• p	
  は従属変数、係数A,B,C,D,Wは定数行列(空間座標i,j,kに依存)

14

Himeno	  benchmark• 19点のステンシル	
  • x,y,z	
  および	
  xy,	
  xz,	  yz	  方向の隣接格子点を参照	
  • 問題サイズ	  S

64x64x128

M

128x128x256

L

256x256x256

XL

512x512x1024

• 演算	
  
• 疎行列ベクトル積	
  • 誤差	
  =各格子点での残差の2乗和

15Code	
  flow

himeno_C/himenoBMTxpa_mod2.c

src_body.cMar 16, 19 11:32
1
2
3
4
5
6

P

jacobi2(…)

float
jacobi2(int nn, Matrix* a, Matrix* b, Matrix* c,
Matrix* p, Matrix* bnd, Matrix* wrk1, Matrix* wrk2){
for(n=0 ; n<nn ; n++){gosa = 0.0;

iteration	
  loopsite	
  loop7

for(int site = 0 ; site < size; site++){int k = site % mdeps;
int j = (site/mdeps) % mcols;
int i = site/(mdeps * mcols);

8
9
10
11

差分を計算,	
  wrk2	
  にストア12

float s0;s0 = pa[IDX(0,i,j,k)]*pp[IDX(0,i+1,j, k)]+ pa[IDX(1,i,j,k)]*pp[IDX(0,i, j+1,k)]+ pa[IDX(2,i,j,k)]*pp[IDX(0,i, j, k+1)]// ...
+ pwrk1[IDX(0,i,j,k)];

13
14
15
16
17
18
19

ss = (s0 * pa[IDX(3,i,j,k)] − pp[IDX(0,i,j,k)]) * pbnd[IDX(0,i,j,k)];20
21

pwrk2[IDX(0,i,j,k)] = ss;

22

}

23

誤差は差分の2乗和(norm)	
  24

for(int site = 0; site < size; site++) {float v = pwrk2[site];gosa += v * v;
}

25
26
27
28

で計算される29

for(int site = 0 ; site < size; site++){int k = site % mdeps;
int j = (site/mdeps) % mcols;
int i = site/(mdeps * mcols);

30
31
32
33
34

pp[IDX(0,i,j,k)] += omega * pwrk2[IDX(0,i,j,k)];

35

}

36
37

} /* end n loop */38
39

return(gosa);

40
41

}

p	
  を更新するaxpy	
  の形になっている

16

312 float
313 jacobi(int nn, Matrix* a,Matrix* b,Matrix* c,
314
Matrix* p,Matrix* bnd,Matrix* wrk1,Matrix* wrk2)
315 {
316
int
i,j,k,n,imax,jmax,kmax;
317
float gosa,s0,ss;318• 格子点をまわるループ	
  319
imax= p->mrows-1;
320
jmax=
p->mcols-1;• ループ内に依存性はない	
  321
kmax= p->mdeps-1;
322• 1格子点の演算	
  ↔	
  1スレッドで実行 としてモデル化、
323
for(n=0 ; n<nn ; n++){格子点数のスレッドを	
  GPU	
  上で並列処理	
  
324
gosa
= 0.0;
325
326
for(i=1 ; i<imax; i++)327
for(j=1 ; j<jmax ; j++)スレッドに分散328
for(k=1 ; k<kmax ; k++){329
s0= MR(a,0,i,j,k)*MR(p,0,i+1,j, k)
格子点での演算330+//	  MR(a,1,i,j,k)*MR(p,0,i,
j+1,k)“カーネル”331+ MR(a,2,i,j,k)*MR(p,0,i, j, k+1)
332
+ MR(b,0,i,j,k)
333*( MR(p,0,i+1,j+1,k) - MR(p,0,i+1,j-1,k)334- MR(p,0,i-1,j+1,k) + MR(p,0,i-1,j-1,k) )ループカウンタ(i,j,k)	
  ↔	
  スレッドを識別する	  Index•335+ MR(b,1,i,j,k)336*( MR(p,0,i,j+1,k+1) - MR(p,0,i,j-1,k+1)
337- MR(p,0,i,j+1,k-1) + MR(p,0,i,j-1,k-1) )338+ MR(b,2,i,j,k)Strategy	  for	  parallelization17

Thread,	
  Block,	
  Grid• スレッドは階層的にグループ化されて実行される	  • スレッドの集まり	
  =	
  ブロック	  • ブロックの集まり	
  =	
  グリッド	  cf.	
  OpenACC	
  では	
  Vector/Worker/Gang	
  に抽象化されている• 各スレッドは一意的な	
  ID	
  で識別される	  • 組み込み変数としてカーネル内から参照可	  • threadIdx	
  	
  	
  	
  block内のthreadのID
blockIdx	
  	
  	
  	
  	
  	
  grid内のblockのID	
  • ブロック・グリッドの形状・サイズは	
  blockDim,	
  gridDim	
  
• ID	
  は	
  1,2,3次元の形状をとる	  • これらの変数は	
  dim3	
  型

18

Thread	  Index• スレッドのインデックスを	
  threadIdx,	
  blockIdx	
  から計算	  idx = threadIdx.x + blockDim.x * blockIdx.x;

CF

idx = threadIdx%x + blockDim%x * (blockIdx%x - 1)

• Fortranの場合、index	
  は	
  1から	
  であることに注意。
19

CUDA-­‐nizing	  code	  flowsrc_body.cMar 16, 19 11:32
1
2
3
4
5
6

float
jacobi2(int nn, Matrix* a, Matrix* b, Matrix* c,
Matrix* p, Matrix* bnd, Matrix* wrk1, Matrix* wrk2){
for(n=0 ; n<nn ; n++){gosa = 0.0;

メモリの確保とデータ転送

cudaMalloc(…)

Page 1/1

cudaMemcpy(…)

7

for(int site = 0 ; site < size; site++){int k = site % mdeps;
int j = (site/mdeps) % mcols;
int i = site/(mdeps * mcols);

8
9
10
11
12

float s0;s0 = pa[IDX(0,i,j,k)]*pp[IDX(0,i+1,j, k)]+ pa[IDX(1,i,j,k)]*pp[IDX(0,i, j+1,k)]+ pa[IDX(2,i,j,k)]*pp[IDX(0,i, j, k+1)]// ...
+ pwrk1[IDX(0,i,j,k)];

13
14
15
16
17
18

カーネル実行
jacobi_step<<<grid,block>>>(…)

19

ss = (s0 * pa[IDX(3,i,j,k)] − pp[IDX(0,i,j,k)]) * pbnd[IDX(0,i,j,k)];20
21

pwrk2[IDX(0,i,j,k)] = ss;

22

}

23
24
25
26
27
28

for(int site = 0; site < size; site++) {float v = pwrk2[site];gosa += v * v;}

(総和計算	
  →ライブラリ利用)

for(int site = 0 ; site < size; site++){int k = site % mdeps;
int j = (site/mdeps) % mcols;
int i = site/(mdeps * mcols);

カーネル実行
jacobi_update<<<grid,block>>>(…)

29
30
31
32
33
34

pp[IDX(0,i,j,k)] += omega * pwrk2[IDX(0,i,j,k)];

35

}

36
37

20

} /* end n loop */38
39

return(gosa);

40
41

}

計算結果をホストに転送

cudaMemcpy(…)

C

Kernel• スレッドで実行する処理を関数の形式で定義	
  __global__ void kernel (int n, float *dA){// kernel contents}__global__

カーネルの記述__device__

デバイスコード用の関数	
  カーネルや他のデバイスコードから呼ばれる__host__ホストコード向けの関数

• 引数	
  
• デバイス上の配列を指すポインタ	  • スカラー値21

C

Writing	
  Kernels__global__
void jacobi_step_dev(float *pa, floag *pb, float *pc, float *pp,
float *pbnd, float *pwrk1, float *pwrk2, int mrows, int mcols, int mdeps){
int index = threadIdx.x + blockDim.x * blockIdx.x;

1

インデックス計算

int k = index % mdeps;
int j = (index / mdeps) % mcols;
int i = (index / (mdeps * mcols)) % mrows;
if (i == 0 || i == mrows-1 || j == 0 || j == mcols-1 || k == 0 || k == mdeps-1){
pwrk2[IDX(0,i,j,k)] = 0.0f;} else {float s0 = 0.0f;
s0 = pa[IDX(0,i,j,k)]*pp[IDX(0,i+1,j, k)]+ pa[IDX(1,i,j,k)]*pp[IDX(0,i, j+1,k)]+ pa[IDX(2,i,j,k)]*pp[IDX(0,i, j, k+1)]+ pb[IDX(0,i,j,k)]*( pp[IDX(0,i+1,j+1,k)] - pp[IDX(0,i+1,j-1,k)]- pp[IDX(0,i-1,j+1,k)] + pp[IDX(0,i-1,j-1,k)] )+ pb[IDX(1,i,j,k)]*( pp[IDX(0,i,j+1,k+1)] - pp[IDX(0,i,j-1,k+1)]- pp[IDX(0,i,j+1,k-1)] + pp[IDX(0,i,j-1,k-1)] )+ pb[IDX(2,i,j,k)]*( pp[IDX(0,i+1,j,k+1)] - pp[IDX(0,i-1,j,k+1)]- pp[IDX(0,i+1,j,k-1)] + pp[IDX(0,i-1,j,k-1)] )+ pc[IDX(0,i,j,k)] * pp[IDX(0,i-1,j, k)]+ pc[IDX(1,i,j,k)] * pp[IDX(0,i, j-1,k)]+ pc[IDX(2,i,j,k)] * pp[IDX(0,i, j, k-1)]+ pwrk1[IDX(0,i,j,k)];float ss = (s0 * pa[IDX(3,i,j,k)] - pp[IDX(0,i,j,k)]) * pbnd[IDX(0,i,j,k)];loop内の処理	
  (Cプログラムから切り出す)

pwrk2[IDX(0,i,j,k)] = ss;}}22

KernelF• global属性付きサブルーチンの形式で定義	
  • カーネルはモジュールの中に置く必要がある	
  module kernel_module
contains
attributes(global) subroutine kernel (n,dA)
integer,value :: nreal,dimension(n),device :: dA!…
end module kernel_module

• 引数	
  
• デバイス配列は	
  device	
  を付ける	  (e.g.	
  dA)	  • スカラー値の受け渡しは	
  value	
  を付ける	  (e.g.	
  n)23

Writing	
  Kernelstmp_f90.cuf
カーネルはmodule内に配置

Nov 29, 18 12:56
1
2
3
4
5
6
7
8
9
10
11
12
13

16
17

デバイス配列は	
  device
スカラー引数は	  value	  で修飾

i = threadIdx%x + blockDim%x * (blockIdx%x−1)
j = threadIdx%y + blockDim%y * (blockIdx%y−1)
k = threadIdx%z + blockDim%z * (blockIdx%z−1)

18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

Page 1/1

module cuda_kernel
use cudafor
contains
attributes(global) subroutine jacobi_step(p,a,b,c,bnd,wrk1,wrk2,&
imax,jmax,kmax)
implicit none
real,dimension(:,:,:,:),device :: a,b,c
real,dimension(:,:,:),device :: preal,dimension(:,:,:),device :: bnd,wrk1,wrk2
integer,value :: imax,jmax,kmax
integer :: index
integer :: i,j,kreal :: s0,ss

14
15

F

if ( (i.lt.1).or.(i.gt.imax).or.&(j.lt.1).or.(j.gt.jmax).or.&(k.lt.1).or.(k.gt.kmax) ) then! index out of range.
returnelse if ( (i.eq.1).or.(i.eq.imax).or.&(j.eq.1).or.(j.eq.jmax).or.&(k.eq.1).or.(k.eq.kmax) ) then! edge points
ss=0.0else
s0 = a(I,J,K,1)*p(I+1,J,K) &+a(I,J,K,2)*p(I,J+1,K) &+a(I,J,K,3)*p(I,J,K+1) &+b(I,J,K,1)*(p(I+1,J+1,K)−p(I+1,J−1,K) &−p(I−1,J+1,K)+p(I−1,J−1,K)) &+b(I,J,K,2)*(p(I,J+1,K+1)−p(I,J−1,K+1) &−p(I,J+1,K−1)+p(I,J−1,K−1)) &+b(I,J,K,3)*(p(I+1,J,K+1)−p(I−1,J,K+1) &−p(I+1,J,K−1)+p(I−1,J,K−1)) &+c(I,J,K,1)*p(I−1,J,K) &+c(I,J,K,2)*p(I,J−1,K) &+c(I,J,K,3)*p(I,J,K−1)+wrk1(I,J,K)
ss=(s0*a(I,J,K,4)−p(I,J,K))*bnd(I,J,K)
end ifwrk2(i,j,k) = ss
end subroutine jacobi_step
end module cuda_kernel

3次元のindexを使用	
  i,j,k	
  に対応させる

loop内の処理(Fortranコードから切り出し)

24

MemoryC• メモリ管理	
  • APIを通じてホスト側から制御	
  • 基本的な流れ	  • デバイス上にメモリを確保	
  (global	
  memory)	
  cudaMalloc(&devPtr, size)• devPtr	
  はデバイス上のメモリへのポインタ
ホスト側からは単なる識別子(直接読み書きはできない)	
  • データ転送	  	  	
  (次ページ)	
  • 確保したメモリを解放	
  cudaFree(devPtr)

25

C

Memory• 基本的な流れ	  • データ転送	  cudaMemcpy(dstPtr, srcPtr, size, kind)• kindは転送の方向	
  cudaMemcpyHostToDevice

ホスト→デバイス

cudaMemcpyDeviceToHost

デバイス→ホスト

cudaMemcpyDeviceToDevice

デバイス上のメモリコピー• dstPtr,	
  srcPtr	
  は方向に応じてデバイス上またはホスト上のポインタを指定
間違えると実行時エラー
26

MemoryF• デバイス配列は	
  device	
  属性をつけて宣言する
allocate文で動的に確保できる	  • real,device,allocatable,dimension(:,:) :: dA
allocate(dA(M,N))
MxN	
  の2次元配列をデバイス上に確保	
  • 配列の代入文でデータ転送	
  • dA = hA• ホスト上の配列	
  hA	
  からデバイス配列	
  dA	
  へ転送	
  
• 部分配列も可	  • deallocate文で解放

27

MemoryF• APIを利用することも可	
  • cudafor	
  モジュールを	
  use	
  する	
  
• 例	
  
• real,device,allocatable :: dB(:)
status = cudaMalloc(dB, 1024)• 多次元配列は扱いにくい

28

09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

}
//===============================================================
float
jacobi(int nn, Matrix* a,Matrix* b,Matrix* c,
Matrix* p,Matrix* bnd,Matrix* wrk1,Matrix* wrk2){
int
i,j,k,n,imax,jmax,kmax;• 関数呼び出しと同じ形式	
  float gosa,s0,ss;

Kernel	  Launch• ブロック・グリッドのサイズと形状を特別な構文で指定	
  imax= p->mrows-1;• kernel<<<grid,block>>>(args,…)
jmax=
p->mcols-1;
kmax=
p->mdeps-1;• callkernel<<<grid,block>>>(args,…)

block	  は	
  int型(1次元)	
  または	
  dim3型(2,3次元)
for(n=0• ;grid,	  n<nn; n++){
gosa = 0.0;
for(i=1 ; i<imax; i++)
並列度の情報for(j=1 ; j<jmax ; j++)for(k=1 ; k<kmax ; k++){s0= MR(a,0,i,j,k)*MR(p,0,i+1,j,kernel<<<grid,block>>>(args)k)+//	
  MR(a,1,i,j,k)*MR(p,0,i,
j+1,k)
格子点での演算+ MR(a,2,i,j,k)*MR(p,0,i, j, k+1)
+ MR(b,0,i,j,k)*( MR(p,0,i+1,j+1,k) - MR(p,0,i+1,j-1,k)- MR(p,0,i-1,j+1,k) + MR(p,0,i-1,j-1,k) )+ MR(b,1,i,j,k)*( MR(p,0,i,j+1,k+1) - MR(p,0,i,j-1,k+1)

29

Kernel	
  Launch• カーネル実行は非同期	
  • カーネル起動後、すぐにホスト側に処理が復帰する	
  • カーネルの終了を待つには	
  • cudaDeviceSynchronize()30

Initialize	  Devices• デバイスの明示的な初期化は不要	
  • 最初のAPI呼び出しで自動的に初期化が行われる	
  • (GPUが複数あるとき)	
  使うデバイスを	
  id	
  で指定	  • cudaSetDevice(devid)• デバイスの機能や設定を調べるには	  • cudaGetDeviceProperties(&prop, devid)• prop	
  は	
  cudaDeviceProp型の構造体	
  • compute	
  capability	
  やデバイスメモリのサイズなど取得できる	
  • デバイスの解放は	  • cudaDeviceReset()• 現在activeなデバイスのcontextを破棄する

31

cuBLAS	
  Library• BLAS	  (Basic	
  Linear	  Algebra	
  System)	
  • 線形代数演算のセット	
  • 使い方	  • #include <cublas_v2.h>use cublas(C)(Fortran)• cublasHandle_t型のハンドルを作成	
  • cublasCreate(&handle)• BLASルーチンの呼び出し例	
  • 2乗ノルムの平方根	
  (単精度)	
  • cublasSnrm2(handle,	
  size,	
  x,	
  xstride,	
  &val)• axpy	  • cublasSaxpy(handle,	
  size,	
  a,	
  x,	
  xstride,	
  y,	
  ystride)32

Compile	
  and	  Run

C• nvcc	  [options]	
  source.cu	  • コンパイラ	
  nvcc	  	  (NVIDIA	
  CUDA	  SDK)	
  • ソースファイル名の	  convention	
  は	
  〜.cu	  • オプション	  	
  (詳しくは	
  nvcc	
  —help)	
  • -arch sm_30 (Kepler)	
  compute	
  capability	
  の指定-arch sm_60 (Pascal)	
  • —ptxas-options=-v• kernel	
  の使用リソースなど詳細情報を出力	  • カーネルコードとホストコードを分離してそれぞれ
コンパイルし、まとめて一つの実行形式にして出力する

33

Compile	
  and	  Run

F

• pgf90	
  [options]	
  source.cuf	
  • コンパイラは	
  pgf90	
  (PGI	  Fortranコンパイラ)	
  • ソースファイル名の	  convention	
  は	
  〜.cuf	  • オプション	  	
  (詳しくは	
  man	
  pgf90)	
  • -Mcuda=cc30,ptxinfo• cc30	
  :	
  (Kepler)	
  compute	
  capability	  • cc60	
  :	
  (Pascal)	
  • ptxinfo	
  :	
  kernel	
  関数の情報を表示	  • -Mfixed	  	
  	
  固定形式のソースファイル
34

Profiling	
  and	
  Debug• 実行時性能を評価	
  • 各関数の経過時間・実行回数をカウント	  • ロード・ストアの効率などをモニタ	  • nvprof	
  [type]	
  ./a.out	  [args]	
  • type	  にprofiileのモードやカウンタを指定
無指定では	  summary	  mode	  • GUI	  の	
  nvvp	  なども利用可能	
  • デバッグ	  • cuda-­‐gdb	  • コンパイル時に	  -­‐G	
  -­‐lineinfo	
  をつけてコンパイル	  • 使い方は	
  gdb	
  とほぼ同じ。特定のスレッドを調査できる	  • printf	
  デバッグ	  • カーネル内で	  printf	
  が使える

35Performance• 演算量	  
• 格子点あたり	
  34演算	  • データ量	
  • 12個の係数・1個の一時変数・1個の変数	
  =	
  14変数	  • 隣接格子点を毎回ロード→+	
  18変数	  • Byte/Flop	
  =	
  4byte	
  x	
  14	
  or	
  32	
  /	
  34	
  ops	  • メモリバンド幅が律速とすると、期待される性能は	
  • 再利用あり	  →	
  115.2	
  GB/s	  /	
  (	
  4byte	
  x	
  14	
  /	
  34	
  ops)	
  	
  =	
  70	
  GFlops	
  • 再利用なし	
  →	
  115.2	
  GB/s	  /	
  (	
  4byte	
  x	
  32	
  /	
  34	
  ops)	
  	
  =	
  30	
  GFlops	
  36

Performance• 実習環境	  	
  GeForce	
  GTX680M	
  	
  (Kepler)	
  /	
  GTX	  1080	
  (Pascal)
GTX680M

GTX1080

P100

V100

Generation

KeplerPascalPascal

Volta

#SM

7

20

56

80

FP32	
  core

1344

2560

3584

5120

FP64	
  core

56

80

1792

2560

Peak	
  FP32

2.03TFlops

8.9TFlops

10.6TFlops

15.7TFlopsPeak	
  FP640.083TFlops

0.26TFlops

5.3TFlops

7.8TFlops

SP/DP

24:1

32:1

2:1

2:1

Mem	
  B/W

115.2GB/s

320GB/s

732GB/s

900GB/s

CC

3.0

6.1

6.0

7.0

Compute	  Capability

37

Performance• Himeno	
  Benchmark: 実習環境で測定	
  (単位は	  GFlops)
GTX	
  680MGTX	
  1080

C

0.56

0.98

C	
  with	  OpenACC

20.22

103.35

CUDA	
  C

33.25

103.25CUDA	
  C	
  (optimized)

43.46

192.72Fortran90

1.07

3.39

F90	
  with	  OpenACC

23.71

88.75

CUDA	
  Fortran

20.35

94.69

• Mサイズ(128x128x256)

38

実習(1)C/F90版をそのまま	
  compile・実行してみる(2)CUDA化する	  a) kernel	
  作成	
  (差分計算・更新)、kernel	
  呼び出しで置き換え	  b) デバイスメモリの確保と解放処理を追加	
  c) データ転送を追加	
  d) 総和計算をcuBLASルーチン呼び出しに置き換え	  (3)パラメータを変えてみる、など	  • ブロックサイズ・グリッドサイズを変更

39

Advanced	
  Topics
1. Shared	
  memory	
  の利用	  • block内のthreadで共有	
  →	
  thread間のcommunication	  • Global	
  memory	
  より高速	  →	
  Programmable	
  な	
  cache	  2. Block	
  size/Grid	
  size	  の最適化	  • block,	
  grid	  のサイズを調整し、SMの占有率を上げる	
  • なるべく多数のthreadをSMに割り当て	  • SMのリソース(レジスタ・共有メモリ)は有限40

Advanced	
  Topics
3. Memory	  layout	
  の検討	  • threadが協調的にglobal	
  memoryにアクセスすると効率がよい→	
  coalesced	
  access	  • 今の場合、係数行列の読み出しがメイン	
  (再利用されない)→	  データストリームを最適化	  4. 複数GPUの利用	
  • single	
  node/multi-­‐node上の複数のGPUを並列に使う	
  • 通信と演算のoverlapによりデータ転送遅延を隠蔽
41

References• CUDA	
  C	
  プロフェッショナル	
  プログラミング,	
  John	
  Cheng,	
  Max	
  Grossman,	
  Ty	
  McKercher	
  著,	
  株式会社クイープ	
  訳	
  • CUDA	
  C	
  PROGRAMMING	
  GUIDE,	
  v10.0,	
  NVIDIA• PGI	  Compilers	
  &	
  Tools	  	
  CUDA	  FORTRAN	
  PROGRAMMING	
  GUIDE	
  AND	
  
REFERENCE,	
  Version	  2018,	
  NVIDIA

42

Backup43

Optimization(a) shared	
  memory	
  の利用	  
• 圧力変数	
  p(i,j,k)	
  は複数の隣接サイトから参照される→	
  shared	  memory	
  を介して隣接サイト間で再利用(b) stream	
  optimization	
  • 係数行列	
  A,	
  B,	
  C,	
  bnd,	
  wrk1	
  は1回のみ	
  load	
  される→別々の配列をまとめることで	
  load	  を効率化
データレイアウトを最適化

44

Optimization(a) shared	
  memory	
  の利用	  • 格子を	
  tile	  に分割、1	
  tile	  を	
  1	
  thread	  block	
  で計算する	  • thread	  block	
  内で	
  p(i,j,k)	
  を共有し、差分計算で再利用する(shared	  memory	
  へのアクセスは	
  global	
  memory	
  より高速)• パラメータ	  • tile	  のサイズ・形状	
  	  —	
  e.g.	
  1	
  x	
  8	
  x	
  128	  • thread	  block	
  のサイズ	
  —	
  e.g.	
  256	
  • Note:	
  1	
  thread	
  は	
  1	
  x	
  8	
  x	
  128	  /	  256	  =	  4	  site	  の演算を担当
45

Optimization• 計算手順	  • shared	  memory	
  上に袖領域を含めて	
  array	
  を確保	  • tile	
  部分	
  (+袖領域)	
  の	
  p(i,j,k)	
  を	
  shared	  memory	
  にコピー	  • shared	  memory	
  を参照して差分を計算• あまり大きな効果はない	  • L1	
  cache	
  が有効

46

Optimization(b) stream	
  optimization	
  • 係数行列	
  A,	
  B,	
  C,	
  bnd,	
  wrk1	
  を一つの配列に	
  pack	  する:	  • A[0..3][Nx][Ny][Nz]	
  	
  →	
  	
  Coeff[0..3][Nx][Ny][Nz]	
  • B[0..2][Nx][Ny][Nz]	
  	
  →	
  	
  Coeff[4..6][Nx][Ny][Nz]	
  • C[0..2][Nx][Ny][Nz]	
  	
  →	
  	
  Coeff[7..9][Nx][Ny][Nz]	
  • bnd[Nx][Ny][Nz]	
  	
  	
  	
  	
  	
  	
  →	
  	
  Coeff[10][Nx][Ny][Nz]	
  • wrk1[Nx][Ny][Nz]	
  	
  	
  	
  	
  →	
  	
  Coeff[11][Nx][Ny][Nz]• 配列の	  index	
  を入れ替え:	
  • Coeff[12][Nx][Ny][Nz]	
  →	
  Coeff[Nx][12][Ny][Nz]	
  • tile	  と格子のサイズを考慮
47• かなり効果がある	
  • memory	
  access	
  の効率化が重要References• CUDAによる定常反復Poissonベンチマークの高速化
小川慧・青木尊之 情報処理学会研究報告	
  2008-­‐HPC-­‐115• GPGPU上での流体アプリケーションの高速化手法
成瀬彰・住元真司・久門耕一
情報処理学会研究報告	
  2008-­‐HPC-­‐117• 姫野ベンチマークのマルチノード実行における通信と演算のオーバー
ラップによる高速化
加藤季広・青木尊之・額田彰・遠藤敏夫・松岡聡・長谷川篤史
情報処理学会研究報告	
  2009-­‐HPC-­‐120
48

Unified	
  Memory• ホストとデバイスのメモリ空間を統合	
  • メモリ管理を単純化	  • 別々のメモリ空間によるポインタの重複を排除	
  • 明示的なデータ転送が不要	
  • 必要に応じてドライバが裏で転送を行う	
  • (Pascal以降)	
  on-­‐demand	
  でpagingする.
GPUメモリ容量を超える	
  oversubscribe	  可能.	  • 使い方	  • cudaMallocManaged(&ptr, size)• ホスト側の	
  malloc	
  呼び出しを置き換え	  • データを書き込み、そのままカーネルの引数に渡せる

49

newMat(Matrix* Mat, int mnums,int mrows, int mcols, int mdeps){
Mat->mnums= mnums;
Mat->mrows= mrows;
Mat->mcols= mcols;
Mat->mdeps= mdeps;
Mat->m= NULL;
Mat->m= (float*)malloc(mnums * mrows * mcols * mdeps * sizeof(float));

Unified	
  MemoryC• Himeno	
  benchmark	  では	
   return(Mat->m != NULL) ? 1:0;}

• Matrix	
  型の初期化	
  newMat()	
  の中で行列のデータを確保して

int
newMat(Matrix* Mat, int mnums,int mrows, int mcols, int mdeps){
Mat->mnums= mnums;
int
Mat->mrows= mrows;
newMat(Matrix* Mat, int mnums,int Mat->mcols=
mrows, int mcols,
mcols; int mdeps){
Mat->mdeps= mdeps;
Mat->mnums= mnums;
Mat->m= NULL;
Mat->mrows= mrows;
int vol = mnums * mrows * mcols * mdeps;
Mat->mcols= mcols;
cudaMallocManaged((void**)&(Mat->m), sizeof(float) * vol);Mat->mdeps= mdeps;
Mat->m= NULL;
return(Mat->m != NULL) ? 1:0;
Mat->m= (float*)}malloc(mnums * mrows * mcols * mdeps * sizeof(float));いる箇所を置き換える1

return(Mat->m != NULL) ? 1:0;}

int
newMat(Matrix* Mat, int mnums,int mrows, int mcols, int mdeps){
Mat->mnums= mnums;
Mat->mrows= mrows;

50

Unified	
  MemoryF• Fortran	
  では	
  managed	
  属性を付ける	  • real,managed,allocatable :: A(:)
allocate(A(N))• API	
  を通じて確保もできる	
  • stat =cudaMallocManaged(A,N,cudaMemAttachHost)

51