ドメイン特化による高性能・高生産性プログラミング
CC2:4C
9A4p) (. / *) (.主な共同研究者

SLWPO GLSTM
2:4Ckrq

8C:4q

2:4C

GPUアクセラレータ•)•/• ) (.

ʼF:5:2 8D)• ( D77•

8AEAC CA 0 877• ( , D7 AC CA 5A• 93 ) ., 83 ]AC 5A

• 4E52 (• 644 “ʻ A4:•
D] MLWP ( )• 2 5 Z P X (. DP]VL C( . 8AEp + 8D) 8AE] ZP M L O• /. D7V Z] ,- S TX 7LVV f 0 D A,F:5:2 AL]NLV 8AE•••

8AEE 2pt

黎明期のアプリケーション事例•JC4e(q2CE42r3P] C OPX ALZP 7TXLVT] K• 7LXy 4E52w• ++ 87 CA aT S XP DP]VL C( . 8AE• 4AE(• (, D7V Z] aT S ,() 8AE]
• DCE32 6p

•
•pà

u

Physisステンシルフレームワーク
JC4e((K S Z1 RT S M N W XL cLW ZSc]T]

2 D B )2
• 5PNVL L T P• A LMVP• 8V MLV TPa• 4 ML]POvoid diffusion(int x, int y, int z,
PSGrid3DFloat g1, PSGrid3DFloat g2) {float v = PSGridGet(g1,x,y,z)+PSGridGet(g1,x-1,y,z)+PSGridGet(g1,x+1,y,z)+PSGridGet(g1,x,y-1,z)+PSGridGet(g1,x,y+1,z)+PSGridGet(g1,x,y,z-1)+PSGridGet(g1,x,y,z+1);
PSGridEmit(g2,v/7.0);})2 (ECF B G• DL RP ]ZPNT TN
N OP RPXP L T X
LXO Z TWTdL T X]
• 2 WL TN
ZL LVVPVTdL T X

PhysisCC+MPI
CUDA
CUDA+MPI

Physis DSL概要
•4••
•p©
p

2A:

struct point {double x, y, z;};DeclareGrid3D(struct point, Point);
void kernel(const int i, const int j, const int k,PSGrid3DPoint g1, PSGrid3DDouble g2) {double v = PSGridGet(g1,i,j,k).x+PSGridGet(g1,i-1,j,k).x + PSGridGet(g1,i+1,j,k).x+PSGridGet(g1,i,j-1,k).y + PSGridGet(g1,i,j+1,k).y+PSGridGet(g1,i,j,k-1).z + PSGridGet(g1,i,j,k+1).z;
PSGridEmit(g2,v/7);}

最適化例: データレイアウト変換•p©
••à 2 Lc] C N
à C N P]p©

•P] 2 C2 Lc] C 2pp

p

2 Cstruct Point {float u;float v;float w;
float x[19];};DeclareGrid3D(Point, struct Point);
PSGrid3DPoint g;...
PSGridGet(g.u, i, j, k);
PSGridGet(g.x[l], i, j+1, k);

8AE

C 2struct Point {float *u;float *v;
float *w;
float *x;};...
// PSGridGet(g.u, i, j, k);g->u[i+j*nx+k*nx*ny]// PSGridGet(g.x[2], i, j+1, k);g->x[i+j*nx+k*nx*ny+l*nx*ny*nz];最適化例: 通信コスト隠蔽( 4 Zc M XOL TP]
W 8AE
4AE
X X XT ] TOP NL]P]:X P TZ TX ]) 4 WZ P]
TX P T Z TX ]+ 4 WZ P] M
3XOL c* 3 XOL c
PbNSLXRP] aT S
XPTRSM ]XOL TP]

DTWP

生成されるコード（ホスト）

4 WZ TXR :X P T

A TX ]for (i = 0; i < iter; ++i) {__PSStencilRun_kernel_interior<<<s0_grid_dim,block_dim,0, stream_interior>>>(__PSGetLocalOffset(0),__PSGetLocalOffset(1),__PSDomainShrink(&s0 -> dom,1),*((__PSGrid3DFloatDev *)(__PSGridGetDev(s0 -> g))), *((__PSGrid3DFloatDev *)(__PSGridGetDev(s0 int fw_width[3] = {1L, 1L, 1L}; int bw_width[3] = {1L, 1L, 1L};
3 XOL c 6bNSLXRP__PSLoadNeighbor(s0 -> g,fw_width,bw_width,0,i > 0,1);__PSStencilRun_kernel_boundary_1_bw<<<1,(dim3(1,128,4)),0,stream_boundary_kernel[0]>>>(__PSDomainGetBoundary(&s0 -> dom,0,0,1,5,0),*((__PSGrid3DFloatDev *)(__PSGridGetDev(s0 -> g))), *((__PSGrid3DFloatDev *)(__PSGridGetD__PSStencilRun_kernel_boundary_1_bw<<<1,(dim3(1,128,4)),0,stream_boundary_kernel[1]>>>(__PSDomainGetBoundary(&s0 -> dom,0,0,1,5,1),*((__PSGrid3DFloatDev *)(__PSGridGetDev(s0 -> g))), *((__PSGrid3DFloatDev *)(__PSGridGet…
__PSStencilRun_kernel_boundary_2_fw<<<1,(dim3(128,1,4)),0,stream_boundary_kernel[11]>>>(__PSDomainGetBoundary(&s0 -> dom,1,1,1,1,0),*((__PSGrid3DFloatDev *)(__PSGridGetDev(s0 -> g))), *((__PSGrid3DFloatDev *)(__PSGrid
cudaThreadSynchronize();}cudaThreadSynchronize();}4 WZ TXR 3 XOL c

AVLXP] 4 XN

PX Vc

フレームワーク実装概要ASc]T] p

ppp

p

• 5C
• ASc]T]•1 A: 4E52à8AE p
A: p• B C6pp

•
•p”

pp

• ASc]T]• 5C©
p

GPU性能評価”
p•”
pp]T X X DP]VL) ,

90
80
7060
50
403020100
Hand-tuned

No-opt3Register
Blocking

Offset CSE

X DP]VL

Full Opt) .,100
Elapsed Time (sec)• 5C

Bandwidth (GB/s). Z OT80

Manual60

Physis

40200128192
Grid size (N^3)

256マルチノード並列（ウィークスケーリング）

25Cp,

••p

ʻ

+,+*,*),)(,(,()/ C*GE ,+,,F G DE:

A MVPW ]TdP1 ),-b),-b),- ZP 8AE 7P WT) ,-マルチノード並列（ストロングスケーリング）•pp+( 5, BEF*) 5* 5)(,(/ C GE ,A MVPW ]TdP1 ,()b,()b+ 0-(,ドメイン特化による最適化例: カーネル融合
do i=1, N
do j=1, N
do k=1, N…
end doend doend do• (©
••op

•••p©
ZPX Apzp yp”nnpp

do i=1, N
do j=1, Ndo k=1, N…
end doend doend dodo i=1, N
do j=1, N
do k=1, N…
end doend doend do…

do i=1, N
do j=1, N
do k=1, N…
…
…
end doend doend do…

実際には・・・••

sed inwell, emmprovememt. Forfusionloops
ops ofsions,

pypz

•
•

np of differentpweather applications”
TABLE I: Features
Application

No. of Kernels

No. of Arrays

Reducible GMEM traffic

SCALE-LES [1]
WRF [15]ASUCA [3]MITgcm [16]HOMME [6]COSMO [4]142122115
9443356446583127
2441%
24%
17%
22%
21%
38%

������������

����and adjusting accordingly would relax the dependences byremoving the order-of-execution constraint and changing thenew arrays to be of the strict type. Note that this comes at theexpense of GMEM capacity.
Output-only arrays: That can not be reused. The data de����

����

����

����

����

�������

����������������

����������

����

����

����

������

����

����

����

����

����

������

����������

������

����

����

����

����

������

����

����

����

do i=1, N
do j=1, N
do k=1, N…
end doend doend do
do i=1, N
do j=1, N
do k=1, N…
end doend doend dodo i=1, N
do j=1, N
do k=1, N…
end doend doend doKernel and array dependency graph (SCALE)

����do i=1, N
do j=1, N
do k=1, N…
end doend doend do����������������

����������������������

��������������

�������

do i=1, N
do j=1, N
do k=1, N…
end doend doend dodo i=1, N
do j=1, N
do k=1, N…
end doend doend do組み合わせ最適化による最適化 [SCʼ14]
•WUT
•
•

Wp

WSb1

OPa

1p

XFp

à

6bLWZVPà

T PW] ZLNUPOMTX]

性能予測モデル•E

ECCE:G =GF

M ED G D

Rf etch + RegF ac ⇤ max(OrdShr(x)) + c ⇤ HT H + Radr + 1 < RT  RM ax••

EECCE:

2F

M ED G D(1 + c ⇤ HT H ) ⇤ TB ⇤ BlocksSM X ⇤ Size(ShrLst)  ShSM X

2N D:

BE A D=D C GE

EG:G
G : GG M
2DED E

D C

GE

BSh = TB ⇤ BlocksSM X /((1 + c ⇤ HT H ) ⇤ Size(ShrLst))
••

EG2

2G

F G EGC DE D: : M ,Bef f ⇤ GM EMBWPM emBound =GFLOPS
8FGE: E D: G D C 2E: A GD B LE B:,PN
M
XF lop(x) ⇤ Hal
Tpro = (Fl +) ⇤ 10 9
PM emBound , x 2 S
B ⇤ T hr

G:BE





End-to-End Transformation [HPDCʼ15]
Construct* 











OEG*

$//$data$arrays$

Threads'per'
192'
Block'
Ac2ve'Blocks''
16'per'SMX'Shared'Mem.'
3KB'per'Block'…'…'Threads'Accessing'5'
Data'Element'
FLOP(s)'in'K_1'12'…'…'

Via$Instrumentation$

Global''Memory'BW'214''
GB/S'

DDG*


 



No.'SM'15'…'…'



Constraint*

Only%Automated%Step%%
in%Manual%Approach%Recommend***Best*Fusions*





Construct*




New*
 DDG*

 




 

Generate*





OEG*
NoYes

Via$Device$Query$$





currGen == maxGen

Shared'Memory'
48'KB'
per'SM'

Construct*
Replace Old Individuals

Return overall Best

currGen++*Grouped*Genetic*****Algorithm*

Device**
Metadata*

32'

Shared'Arr.'in'K_1' V'

Via$Static$Analysis$

Ac2ve'Threads''
512'
per'Block'Reg.'per'Thread'Genetic Operators

Operations**
Metadata*

Run2me'K_1' 1.6'μs'

Construct*

Performance*
Projection*
Model*

Evaluate Objective Function

*K__1<<<G,B>>>(…);**K__2<<<G,B>>>(…);*'

Performance**
Metadata*

Metadata*

currGen = 0, InitPopulation

$//$kernels$$

GFLOPS'of'K_1' 54'

Extract*

Constraint**double*V[X][Y[Z];*Original******Code*

currGen = 0, InitPopulation

Performance*

New*
Code*

$//$kernels$$*K__A<<<G,B>>>(…);**K__B<<<G,B>>>(…);*' 

アプリケーション適用事例FFBED)G F EDC42 6[Weather]

Next generation mesoscale weather model[Four years in development]

96[Climate]

Dynamical core of Community Atmospheric Model (CAM)

7V LW[Hydrodynamics]

A fluctuating particle hydrodynamics application based on an hybrid
Eulerian-Lagrangian approach:DRNW[Oceanic]

An oceanic general circulation model relying on a finite volume numericalmethod [18 years in development]2GA 54 8AE[Seismic]

An earthquake wave propagation simulator[ACM Gordon Bell finalist]

3 42[FDTD]

2 *5 75D5 ]TW VL

aSTNS W OPV] SP ZP WT T T c
WL P TLVOT]ZP ]T P

性能評価結果
2.0

2.0
FusionFusion+Fission
Fusion+Tune

Fusion+Fission+Tune
Programmer−Guided Fu.
Programmer−Guided Fu.+Fi.

Programmer−Guided Fu.+Tu.
Programmer−Guided Fu.+Fi.+Tu.
Manual

1.5

Speedup

1.5

1.0

1.0

0.5

0.5

0.0

SCALE−LES

HOMME• DP]VL + 8AE•4E52•(.
•FluampMITgcm

AWP−ODC−GPUB−CALM0.0Adaptive Mesh Refinements (AMR)•• p• ʻyʻkàzzzʻv

D 2 UTl”kp©

GPU実装例q

D 2 UTr

Octreedata structureleaf memory
is managed
on CPUleaf GPU

16x16x16

HilbertSpace-filling Curve

Domain size
Mesh resolution (Δx/d)
Blockage ratio (d/D)
Total mesh numberGPU

D×D×2D
0.001 min.
0.512 max.
0.01529
2.7×108
NVIDIA K20X × 64

格子解像度が異なる格子間の補間

D 2 UT•()+/y•()()

y+•z•m•m% , -./ 0%!"= !"++ , - 0%, -3/ 0%#,()#0%!"= !"+2 , 0%
#$%

#,()

!"

#− !"

#,()!"

#− !"

#,()t
*Filippova Olga and Hnel Dieter. Grid Refinement for Lattice-BGK Models.

Journal of Computational Physics, Vol. 147, No. 1, pp. 219-228, 1998.
23

格子解像度ごとの時間ステップの設定
3s5
Δ7

D 2 UT

Daino: AMR並列フレームワーク
Application
C CodeStencil CodeAMR
Parameter
file

6WT

Translator

Application
LLVM IR
Stencil IR

Generate

Stencil GPU
Kernel(NVVM IR)

Refine Kernel(NVVM IR)
Coarsen Kernel(NVVM IR)
Error Kernel(CUDA)IR PassApplication
LLVM IRAMR Driver
IRCompile

Object Files
Link

Executable

NVPTX

CUDA Driver API CallPTXAPICall•F•à8AE p• 2 B©

Daino Runtime
AMR libraryComm. libraryJC4e(-Kアプリケーション事例
Application

Description

A 2nd order directionally split hyperbolic
Hydrodynamics
schemes to solve Euler equations.
Solver[RTVD scheme modified from GAMER1]
We model shallow water simulations by depthShallow-water
averaging Navier–Stokes equations. [2nd order
SolverRunge-Kutta method]
3D dendritic growth during binary alloy
Phase-fieldsolidification2Simulation[Time integartion by Allen-Chan equation][1] H.-Y. Schive, U.-H. Zhang, and T. Chiueh. Directionally Unsplit Hydrodynamic Schemes with Hybrid MPI/OpenMP/GPU Parallelization in AMR. Int. J. High Perform. Comput. Appl., 26(4):367–377, Nov. 2012[2] T. Shimokawabe et. Al, Peta-scale Phase-Field Simulation for Dendritic Solidification on the TSUBAME 2.0
Supercomputer, SC’11

ウィークスケーリング性能
8AE

PHASE-FIELD

3.0E+02

Runtime (Seconds)

1.0E+03

1.0E+03

5.0E+02

5.0E+02
0.0E+00

48163264128256512996

Number GPUs (Mesh size per GPU: 4,096x512x512)••• (•k”

p0.0E+00

48163264128256512996

Number GPUs (Mesh size per GPU: 4,0963)

pp8AE
p”

2 Bk
ʻ

3.1 x

2.5E+02

4.1 x

Runtime (Seconds)

Auto AMR (GAMER)

2.0E+03

1.5E+03

1.5E+03

Uniform Mesh

Runtime (Seconds)

2.0E+03

SHALLOW-WATERS
Auto AMR
Hand-written AMR

HYDRODYNAMICS
Auto AMR (Daino)

Uniform Mesh
Hand-written AMR

Hand-written AMR

8.7 x

Auto AMR

1.7 x

Uniform Mesh

1.9 x

2.5E+03(9.6 x

DCE32 6) ,2.0E+02
1.5E+02
1.0E+02
5.1E+01
1.0E+00

48163264128256512

Number GPUs (Mesh size per GPU: 8,1923)

”
p2 B
p996

ストロングスケーリング性能
Uniform Mesh
Auto AMR
Hand-written AMR

3.5E+04

1.7 x

3.0E+04

HYDRODYNAMICS

2.1E+03

3.0E+02

9.6 x

Runtime (Seconds)

Runtime (Seconds)

2.5E+02

2.5E+04

SHALLOW-WATER S

Uniform Mesh
Auto AMR (Daino)
Hand-written AMR
Auto AMR (GAMER)

3.5E+02

2.0E+02

2.0E+04

1.5E+02

1.5E+04
1.0E+04

1.5 x

1.0E+02

8.9 x

5.1E+01

5.0E+03
1.0E+00

4••8163264128p

•
•

512

3

Number GPUs (Mesh size 4,096 )p•

256996

1.0E+00

48162 Bop

p

128

256

512

3996

2.5E+02
2.0E+02

4.1 x

1.5E+02
1.0E+02

3.2 x
5.1E+01
1.0E+00

48163264128p

q82 6B ”p

GLSTM L cLWL LXO 2 UT g5LTX 1 2 9TRS P PV 7 LWPa U
2 B X 8AE] i C4() (F G L G:

256

512

Number GPUs (Mesh size per GPU: 8,1923)

p

p

8AE

64

Number GPUs (Mesh size per GPU: 4,096 )ʻ8AE(.
p32Uniform Mesh
Auto AMR
Hand-written AMR

3.0E+02

Runtime (Seconds)

PHASE-FIELD

AL LVVPV LXO 6 TNTPX996

AIMES: Advanced Computation and I/O Methods forEarth-System Simulations• 7 XOPO Mc CD 578 LXO 2 B* cPL ]CAA6H2) ) (- h ) (/• A:1 DS WL] OaTOR 5 BI
L cL
L cLWL B: 6DLULc UT 2 UT D Uc
DPNS DS WL] 5 M ] 6N VP A Vc PNSXT[ P(EBBEG EG全球正20面体気候モデル•:42 q2:4C••••prp©kp©
全球870mメッシュシミュレーション（京の20480ノード、163840計算コアを使用）

à A:àpp•p©
••pqpr• 8AE HP X ASTzp

GridTools••4pppp• 4AEv”pp••

CD62z”

ppk8AE•)pw

p©pq4C4Crx”‒p

NICAM on GridTools
5cXLWTN]

8 TOD
4V]

V TN P
4AE
FP ]T X

ASc]TN]:42
7LX
A:ASc]TN]

8 TOD V]Ac S X

8AE
FP ]T X

ASc]TN]TRTXLV
4 WW
O VP

HP X ASTFP ]T X

おわりに•k••p
• 2 B

•

5Cpp”
•p©

:42p
8 TODV]••SLWPO GLSTM• 4C4C 8 TOD V] p
• CD 4B6CD©

kpp”xv