（様式 20）成果報告書

TSUBAME 共同利用 平成 27 年度 学術利用 成果報告書

利用課題名 GPGPU における核融合プラズマコードの性能評価
英文：Performance analysis of fusion plasma code on GPGPU
井戸村泰宏
Yasuhiro Idomura
日本原子力研究開発機構
Japan Atomic Energy Agency
https://www.jaea.go.jp
邦文抄録（300 字程度）
セミラグランジアン法および差分法に基づく 2 つの核融合プラズマ乱流コードの高次元ステンシル演算に対するGPU 向け最適化技術を開発した。セミラグランジアン法におけるメモリのリストアクセスをテクスチャメモリの利用によって大幅に高速化した。差分法においては、演算子の対称性を活用したレジスタ再利用によってメモリアクセス数を削減した。これらの最適化により、Sandy Bridge で最適化されたコードに対して、GPU を用いてセミラグランジアン法で 7.6 倍、差分法で 4 倍の高速化に成功した。
英文抄録（100 words 程度）
We developed optimization techniques for high dimensional stencil computations on GPUs, which are
extracted from two fusion turbulence codes based on Semi-Lagrangian and Finite-Difference schemes.
The indirect memory access of the Semi-Lagrangian scheme is dramatically accelerated by using the
texture memory. In the Finite-Difference scheme, the reuse of registers by taking account of the physicalsymmetry of the operator reduces the amount of memory accesses. Through these optimizations, weachieve acceleration of 7.6x for the Semi-Lagrangian scheme and of 4x for the Finite- Different schemeon GPUs compared with the fully optimized codes on Sandy Bridge.

Keywords: Plasma turbulence, performance analysis, GPGPU, Semi-Lagrangian, Finite-Difference
背景と目的

大規模科学技術計算となっており、特に、次世代の核磁場閉じ込め核融合炉は約１億度の燃料プラズマを融合実験炉 ITER の解析を行う上でエクサスケール計

トーラス形状の磁場で閉じ込めて核融合反応を持続さ算を必要としている。このような、エクサスケール計算

せて、核融合エネルギーを取り出す。炉心プラズマのを実現する計算手法を確立することは、将来の核融合

閉じ込め性能はプラズマ乱流による燃料粒子やエネル

炉心解析コード開発に必要不可欠であり、核融合エネ

ギーの輸送で決まっており、プラズマ乱流の第一原理

ルギー開発に大きく貢献する。

計算が炉心性能評価に必須のツールとなっている。核エクサスケール計算機では、GPU や Xeon Phi を始融合炉心プラズマは粒子間の衝突が弱く（クヌーセン数

めとするアクセラレータが主要な役割を果たすと考えらが大きく）、粒子軌道やランダウ共鳴といった運動論効れる。しかしながら、これらアクセラレータ上で既存の果が顕著となるため、その第一原理モデルはサイクロト

CPU 向け計算手法や最適化手法が有効とは限らない。

ロン運動を平均化した 5 次元のボルツマン方程式（ジャそこで本課題では、異なるアルゴリズムを用いる第一イロ運動論方程式）で与えられる。炉心プラズマを構成

原理的核融合プラズマ乱流コードからホットスポットをするイオンと電子はその質量比のため、熱運動や粒子抽出したのち、それらについて GPU 上で最適化し、性

軌道半径の時空間スケールが大きく乖離し、イオンに能比較を行う。これにより、アルゴリズムごとの GPU に

ついては燃料の重水素、三重水素、および、核融合反対する適正を明らかにするとともに最適化手法を確立

応で生成するヘリウムを含む多種イオン系となる。このすることを目的とする。

ように、核融合プラズマ乱流コードは、5 次元問題、マルチスケール問題、多相問題という極めてユニークなstep t = tn as illustrated in Fig. 1 (b). Giv
φ. P1 and 3.2.2 Finite-Difference Scheme
chemesis row major (C style), it can be seen thaters, g is thesolve Eq.
with 4th-order
Finite-Difference
method,（様式
is continuous
accessin20）成果報告書“j ” direction anm electronn,
we brieflyIf we
describe
the(3)numericalit (Seeis discretizedi” 上で高い性能を発揮するための最適化技術の開
direction. Thus, the Finite-Differencein the kernels
Ref. [25] asfor details), “GPU概要
n+1 schemes.
nn the combinationn
nof continuous aanFinite-Difference
Forn involves
viaanda 4D発を行った。
本課題では、核融合プラズマ乱流シミュレーションコfi,j− fi,j−f
i+2,j + 8fi+1,j − 8fi−1,j + fi−2,j+ vi but
sider
1D ード
examplestheseからホットスポットに対応
schemesor a given
GT5D およびof
GYSELA
∆t
12∆x
d
to to
extendthem to higher
するカーネルコード（それぞれ差分法カーネルおよびセnds
then dimensionalnnn−fi,j+2 + 8fi,j+1 − 8fi,j−1 + fi,j−2
上で最適
e simple
1D
equation(the advecor
Eq. (2)ミラグランジアン法カーネル）を抽出し、GPU
is Vlasov
= 0.
(6)
+ a
i,j
12∆v
化技術の開発を行った。
hase
thisspace)study.
これらのコードは 4 次元の移流問題に対応する同一
s. [17], [18].
The values of f at next time step t = tn+1 can be∂f の方程式を異なるアルゴリズムで解く。アクセラレータ∂fby0,
the Eq. (6) usingt)+ a (x,calculatedv, t)=(3) the values of f at current time
においては、メモリアクセスがボトルネックとなりやすい∂x∂vtn as illustrated in Fig. 1 (b). Given the f [i][j] arraystep t =ことが知られているため、各アルゴリズムにおけるメモis row(C style),
bution function,v is major
the velocity
anditacan
is be seen that the memory access
リアクセスパターンに着目した。各スキームにおけるメ
ishow
continuous
accessin“j ” direction and strided access inWenumerical
will explain
Eq.2 次元（空間(3) canbe solved
モリアクセスパターンを
x、速度
v）の移流
i” direction. Thus,
the Finite-Difference method generally
or details),
grangianand “Finite-Difference
方程式に基づき説明する。スカラー量
f schemes,の移流方程式
involves the combination of continuous and strided access.
hemes. For
は、
hemes but
Fig. 1. Numerical schemes. (a) Semi-Lagrangia
values of f at t = tn+1 at grid points are estima
imensional
ngian Schemeing backward characteristics s (tn ) (foot point)
the advec6], [27], we
will introduce some basics
on Difference scheme where the values of f at tで表される。セミラグランジアン法では、fが特性曲線
cheme. In（粒子の軌道を意味する）に沿って一定となる性質を利the Semi-Lagrangian approach, are estimated by using the values of f at t = t
points (the number of required points can vary on
f value along characed as propagating
用し、あるタイムステップにおける各格子点上のf の値(3)を粒子の出発点の値から評価する。通常、出発点は格= s (t), whereΓ = (x, v) is an Eulerian
子点上の点とならないため、出発点近傍のデータの補

ty and a間によってisf を評価する（図１(a)参照）。粒子軌道は空n be solved間の各点によって異なるため、出発点の座標も空間の
e schemes,
各点によって異なる。これは、データ補間のために、出発点近傍の格子にアクセスする際、リストアクセスが生じることを意味する。ここで、GYSELAにおいては 4(a)次 Semi-Lagrangian
図１ (a)セミラグランジアン法カーネル、および、(b)差Fig. 1. Numerical schemes.
scheme where thevalues off at t 次元+1
= tn+1at grid points分法カーネルにおけるデータアクセスパターン。
are estimated by the correspond元移流の演算子を
2 次元+1次元の演算子とing backward characteristics s (tn ) (foot point) at t = tn , (b) Finite分割することでデータ局所性を向上させている。今回用
Difference scheme where the values of f at t = tn+1 at grid points
いるのは最も計算コストの高い 2 次元の補間カーネル結果および考察are estimated by using the values of f at t = tn on the adjacent gridである。
本課題では、上述した GPU(Tesla K20X)におけるカe basics on
n approach,points (the number of required points can vary on the scheme.).
ng charac一方、差分法においては、移流方程式を差分化して、ーネル最適化による性能向上の程度を評価するため、n Eulerian
次のタイムステップの f の値を評価する。図１(b)に示す
通 常 比 較 対 象 と な る 汎 用 メ ニ ー コ ア CPU(Sandy
ように、差分法においては隣接格子の参照が必要とな

Bridge)に 加 え 、同 じく汎 用 メニ ー コ ア CPU で あ るる。メモリアクセスパターンを考えると、最内の差分におFX100(SPARC64XIfx) 、 ア ク セ ラ レ ー タ Xeon

いては連続アクセスとなるが、それ以外の方向においPhi(Xeon Phi 5110P)において最適化されたコードとては、ストライドアクセスが生じる。

の演算性能比較を行った。FX100 や Xeon Phi は

本課題では、リストアクセスおよびストライドアクセスと

1TFops 級と GPU と同程度の演算性能を持つため、最いうアルゴリズム特有のメモリアクセスパターンに着目

適化効果を確認する上での比較対象としてより適切でし、これらのメモリアクセスパターンを持つカーネルが

ある。

5
theta
AoS[i] = feet[l][k][j][2∗4.3(ii + i)O
Code 2. Access to the foot point arrayin AoSor star
AoSoA
6phi star AoS[i] = feet[l][k][j][2∗(ii+ i)fir+1#define VSIZE 16
We（様式 20）成果報告書
7 {2
for(int ii = 0; ii < Nr; ii += VSIZE)
that inw// Load the coordinate of foot points
3
for(int i = 0; i < VSIZE; i++){ 8
クスチャマッピングと同様である。そこで、GPU が有す64∗(w
9 points
thetastarstyle
AoSoA[i] = feet[l][k][j][2ii +4// Load the coordinate of footin AoS
る画像処理パイプラインを有効活用するため、リストア
セミラグランジアンカーネルの最適化
10∗ii + i35
theta star AoS[i] = feet[l][k][j][2∗クセスの参照先配列をテクスチャメモリに配置した。最(ii phi+ i)];star AoSoA[i] = feet[l][k][j][2size,
前述したように、セミラグランジアン法では粒子軌道を
11 ∗(ii }6
phi star AoS[i] = feet[l][k][j][2+ i) + 1];
終的にこれらの最適化によりオリジナルの 1.79 倍の性 thread
たどることで、分布関数の時間発展を計算する。プラズ
12}7
能向上を得た。これは、最適化された Sandy Bridge loop s
マ乱流中においては、空間の各点で粒子が受ける力が// Load the coordinate of foot points in AoSoA style
8
版の約 7.55 倍に対応する。
異なるため、粒子軌道は各格子点で異なったものとな
best co9
theta star AoSoA[i] = feet[l][k][j][2∗ii + i];る。従って、粒子軌道をたどり、その出発点を追跡する
10
phi star AoSoA[i] = feet[l][k][j][2∗ii + i + VSIZE];(64, 8,際、格子点ごとに参照する格子点は異なる。このメモリ
11}
registe12}アクセスパターンをリストアクセスと呼ぶ。

オリジナルの CPU 版カーネルでは、各格子点において出発点の座標データのペア(θ,φ)を保持する Arrayof Structure（AoS）型のデータ構造を取っていた(図and “
storedFig. 2. The memory layout storing foot point
The foot points are stored in the Array of Stru
ネルを以降ではオリジナルと呼ぶ。このデータ構造は、 the strided access is required. (b) The
Code
foot3
キャッシュ局所性を向上させるため、従来の CPU では Array of Structure of Array (AoSoA) 1 layout
i = thw
有効な最適化である。しかしながら、この配列へのメモ is sequential.
2
j = thr
リアクセスはストライドアクセスとなるため、GPU には
出発点配列の AoS 形式(a)と AoSoA 形式(b)によ3
k = blFig.
2. The memory layout storing foot図２
pointsin (θ, ϕ) dimension. (a)適さない。そこで、アクセスパターンが連続アクセスとなる格納。
Thefoot points are stored in the Array ofStructure (AoS) layout where 4 do l =theるように、配列を構成する構造体を座標データのペア
strided access is required. (b) The foot points are stored in the 5 ! Load6
thet(θ,φ)でなく、SIMD
幅単位のθ座標列とφ座標列の差分法カーネル最適化
load
both
array
Array
of Structureof Array (AoSoA)SIMD
layoutwhereinstruction
the memorysinceaccessphi
ペアからなる構造体へと変更した。図２(b)に示されるよ
図１(b)に示すように差分法では、隣接格子点に対す7
is sequential.
corresponding
pseudo code is given! Corr
うに、この配列へのメモリアクセスは連続アクセスとなる参照を行う。GT5D では、4 次元配列 f(l,k,j,i)に対し8With
thememorylayoutoptimization,
thet
る。ここで、図２(b)では便宜的に SIMD 幅 4 の場合をて差分演算を行うが、この際最内の l 方向については910phiof 1.10×.
示しているが、実コードでは Warp 数 32 の整数倍とし up連続アクセスとなるものの、それ以外の
k,j,i 方向につた 。 こ のload
ような配
列 の デ ー タ 構 造 since
は Arraybothof
いてはストライドアクセスとなる。
SIMD
instruction
arrays are contiguous. The 11 ! Com
12
call
Structure of Array（AoSoA）型と呼ばれ、Xeon Phi
先行研究[Fujita et al., PDSEC2014]によって、2(a)参照)。このデータ構造のまま GPU に移植したカーbalance
optimization
corresponding pseudo code is4.2.2
given Loadin Code2 (line9-10). 13などの SIMD 幅が大きいアーキテクチャで有効な最適GT5D の GPU 移植が行われていたため、本課題ではWith
the memory layout optimization,
we measurea speed- we
14
call
further
optimization,chan
化技術として知られる。これにより出発点の座標計算 For
移植済みのカーネルをもとに最適化を行った。このカー
15
upにおけるメモリアクセスを
of 1.10×. Coalescing load とした。これ scheduling
ネルを以降ではオリジナルと呼ぶ。差分演算の係数行
from static to dynamic.
The
16
!
Load
によってオリジナルに対し、1.24 倍の性能向上を得た。
列は一般的には 4 次元の依存性を持つが、ITER を始
is helpful
for improving the load
balan
17
vda
加えて、計算された出発点の座標を元に参照される
めとするトカマク型装置における軸対称磁場配位の場
4.2.2 Load balance optimization
vda
threads.
A possible candidate 18for the
配列へのリストアクセスについても最適化を行った。プ合には係数が k 方向に依存しないという特徴を持つ（ト19vdac
indirect
access thepart where
the loadingForfurther optimization, wechangedOpen-MPラズマ乱流中における粒子軌道は空間位置によって変ロイダル対称性）。オリジナルのカーネルでは、最内220...
化するものの、電場や磁場といった軌道を決める物理
by次元の(l,k)方向にスレッドを割り当てていた上、係数のthread.
As a result,
we obtainedaspschedulingfrom static to dynamic.
The dynamicscheduling
21
vda
計算結果を
へ格納にするために、ifrespect
tosharedtheoriginal
performance.
is 量が空間内で連続的に変化することから、近接格子点
helpful for improving the loadbalanceby memory
reducing
idling 22 ! Cub
における軌道も連続的に変化し、近接格子点における文を用いた分岐によって係数の計算を k 方向の index
23
fval
threads.
A possible candidate forが１のスレッドのみで計算していた。このthe load imbalance Warp
is the
参照点もまた空間的に近接したものとなる。これは、リ
分岐を
24
(vd
4.2.3
Resultsindirect
access part where the loading
costs can differ thread 25 + vda
ストアクセスにおけるメモリアクセスパターンにある程度解消するために、スレッド割り当て方向を(l,j)方向へと変更し、係数の計算結果はレジスタへ格納するように
by空間局所性があることを意味する。このアクセスパター
thread. As a result, we obtained
aimpact
speed-upof 1.15× withon
26 Xeon+ (vPThe
of optimization
ンは、GPU による画像処理で非常によく用いられるテ変更した。この変更によってオリジナルの 5.43 倍の性
27+ vda
respect to the original performance.
Table 3. Firstly, the memory layout
o
28
+ (v
modified from the original AoS
style
29
+
vda
4.2.3 Results
and the speed-up of 1.10× is30obtaine+ (vbe the sameeffective.
on all theenvironments,
total
be the attainsameexcept
on allthat
the theenvironments,ex
Bothkernels
reasonably
goodperformanc
duetoFinite-Difference
kernel,
where
the
mem
alesced)effective.
Bothkernels attain reasonably
good performances
flop of theSemi-Lagrangiankernelon Semi-Lagrangian
FX100
by trueof thekernelonwith minimum flop
optimization.
Thisisisreduced
especially
forFX
t
The
mo
with
minimum
optimization.
Thisisespecially
true
for
the
is
well
suitedfor
multi-coreCPUs.
We
about 10%.Since
the
FX100
has
commonfeatures
as
Sandy
about
10%.Since
the
FX100
has
commonFinite-Difference kernel, where the memory access patte
usagefFinite-Difference
kernel,where
memoryaccesspattern
6.1.2 w
performanceontheaccelerators.
Bridge (large
cacheandtheSIMD
widthequals4), discuss
Bridge(largeshared
cache
and
SIMD
is shared
well suitedformulti-core
CPUs.
We to
will
t
is well suited for multi-core CPUs. We will discuss the in the S
performancestrategy
onatheaccelerators.
a similar optimization
as in optimization
Secs. 4.1 andstrategy
5.1 was as As
similarin weSec
access t
performance on the accelerators.
alescedeffective. Both kernels attaineffective.
reasonably
goodperformancesBothkernels
attain reasonablygThe shaThemwith minimum optimization.
This
is especiallytrue for thewithminimum
optimization.
Thisbeishardesp
usage
Finite-Difference kernel, where
the memory kernel,
access pattern
Finite-Difference
where themem
Despitein theis well suited for multi-coreCPUs.
We will
discuss theCPUs.
is wellsuitedfor multi-coreWekerneloaccessperformance on the accelerators.
performance on the accelerators. with
re
The
sh
the divebe harlarge
imDespiteffective
kernel4.03 倍となった。この差分法カーネルの最適化に関すFig. 10. The obtained
performancetheSemi-Lagrangian
Difference
kernelsofon
Sandy
Bridge, Xeon andPhi, FinGP
Fig. 10. The obtained
performanceofthe Semi-Lagrangian
andFinite-and FX100.
withrT
2 T ESTBED DESCRIPTION
2.2 GPGPU
Difference
kernelson
Sandy
Bridge,
Xeon
Phi,
GPGPU
left plot shows the speed-up relative to Sandy
Bri
る成果を GTCJapan2015 において口頭発表した[1]。
6.2 FuDifference
kernelson
Sandy
Bridge,
Xeon
Phi,
GPGPU
and
FX100.
The
described
above, GPGPU
hasto
theSandy
SIMT architecture.
In order to understand the difference inleft
the plotarchitectures,showsAsthespeed-up
Bridge,andtheoretical
the
right
p
the
div
shows
therelative
achieved
percentage
ofplottheA thread percentage
group
consistingBridge,of 32
threads,
calledleft
plot showsthespeed-up
relative
to Sandy
andwhich
the isright
we employ typicalmulti-coreCPUs andaccelerators.
shows
the Asachievedofthetheoreticalpeakperformance.
For
acce
works
with the same
instruction.
A divergent
the conventional multi-corechoose the Intel
Sandy as “warp”,
largeishowsCPU,the weachieved
percentageof thetheoreticalpeakperformance.
ベンチマーク結果
Bridge-E5 2680 [20] (Called “Sandy Bridge” in this work). In branch can be problematic, because all paths in the branchportedeffectiv
addition, as the latest multi-core CPU, we use the Fujitsu are executed by unit of warp, and a part of threads are
図３に示すように最適化されたカーネルコードの性能

complex memory access via large shared caches, but the Instruction, Multiple Thread (SIMT)” manner. In contrast,
accelerators would not be. Therefore, it is worth evaluating Xeon Phi allocates the fine-grained tasks to each thread. It
whether these kernels can sustain good performance on the is a kind of “Same Instruction, Multiple Data (SIMD)” archidiscussed
accelerators, since the most efficient memory access patterns tecture. The difference in the memory hierarchy（様式 is20）成果報告書
below.
may be different on each architecture.
In this paper, we present the optimization of the SemiLagrangian and Finite-Difference kernels for the 4D advecParallel
能向上を得た。
2.1 Xeon
Phi and Distributed Systems 紙へ論文投稿を
tion term on the Intel Xeon Phi coprocessor and the Nvidia
Phi is characterized by multiple SIMD cores. The numk GPGPU.
方向の差分時に参照される隣接点のデータを再
Tesla
The strategy to extract high performance Xeon 行った[3]。
on each architecture is investigated through the detailed ber of cores up to 60 (Xeon Phi 5110 P) is larger than conven利用するため、カーネルの最内ループを
k 方向としたcomparison of different memory access patterns,
the indi- tional CPUs. In this study, we keep the number of threads
rect access in the Semi-Lagrangian scheme and the strided as 240 (4 Hyper Threads in each core) and “KMP AFFINITY
上、k 方向の隣接点データをレジスタに格納し、サイク = compact” is specified for thread allocation style. Anotheraccesses in the Finite-Difference scheme.
This paper is organized as follows. Section 2 describes
リックに再利用した。これによって、オリジナルの7.44 important aspect is a L2 cache, which consists of private L2
the characteristics of accelerators used in this paper. In Sec. caches on each core connected with the “Ring Interconnect”
倍の性能向上を得た。i方向 technology. A remarkable characteristics is the coherency in
3, the physical model is 方向の差分については、ifirstly provided. Then the twodifferent numerical schemes, the Semi-Lagrangian and the the L2 cache, and thus it can be considered as the convenにアンロールすることにより、従来メモリアクセスが生じ
Finite-Difference schemes are explained focusing on their tional shared L2 cache in the multi-core CPUs. The peakmemory access patterns. The optimization approaches for bandwidth is not large enough to get a sustained computaていた部分の一部をレジスタへのアクセスで置き換え
the Semi-Lagrangian scheme is shown in Sec. 4 and for the tional performance about 1 TFlops. Although three different
Finite-Difference scheme in Sec. 8.53
5. In Sec.
6, the optimization programming models are available on Xeon Phi, we use theた。これによりオリジナルの倍の性能向上を得た。
results are discussed with the cross-platform benchmarks. so-called “native” approach for simplicity. Another benefit
The obtained results are summarized in Sandy
Sec. 7. Bridge の of this approach is to avoid the data transfer between host
最終的に得られたコードの性能はCPUs andXeon
Phi Theco-processors.
Fig.
10.
obtained performance of the Semi-L

because(e. g. F
so-called “shared memory”. The usage of “shared memory”
for accelerators, we employ the Intel Xeon
Phi coprocessor,
6.1.1
XeonthePhi
Fig.
obtained
performanceof the Semi-L
Fig.
The GPGPU,
obtainedof10.
theTheSemi-Lagrangian
andpatFiniteflexible
programmingof various memory
accessXeon Phi 5110P [22]and 10.
the Nvidia
K20X allowsのベンチマーク結果。性能向上比（上）および演算性能化されたコードと比較した。各ハードウェアの性能につ
6.1.1
Xeon
Phi Tesla performance
GT5D
cFl
6.2
Firstly,
itPhi,
turned
out
that
the
high
cache
terns,
which
is notallowed
on
Xeon
Phi.
In
addition,
there The
[23]. The hardwareDifference
features of thesekernels
architectures
given
Difference
kernelson
Sandy
Bridge,Xeon
Phi, GP
onareSandy
Bridge,
Xeon
GPGPU
and
FX100.
Firstly,
it turned
out
that
the
high
cache
localityisessentいては、表１に示す。
の対ピーク比（下）。
somehighcached on-chip
memories
such asis
theessential“constantin Table 1.
Firstly,
it turnedout thatarerelative
the
cache
locality
porting
left
plot
shows
the
speed-up
relative
to Sandy
Brleft plot shows
the speed-up
togoodSandy
Bridge,andtheThe
right
plot
for
performance
on
Phi.
Th
memory”
and “texture
memory”,
which canprovidehigherXeon
for
good
performance
on
Xeon
Phi.
performance
For acc
Processor
Sandy Bridge
Xeon Phi GPGPU
FX100
shows
the
achieved
percentage
of usage
the
for
good
performance
on
Xeon
Phi.
The
performanceis theoreticalshows
the achieved
percentageofthe theoreticalpeakperformance.
aswithhost-aeffective
bandwidth
in somesituations.
Number of cores
8
60 improved896
32+2 by effective
improvedbyeffective
cache
cache
usage
with
appropriate
thre
portedShared Cache [MB]20
0.5 × 60
1.5
24
まとめ、今後の課題
improvedbyeffectiveout
for G
Memory [GB]648632 cache usage with appropriate thread
allocationin OFtheFinite-Differencekern
Finite-Difference
kernel, and thehiPeak performance
172.8
1010allocation
1310
1000 in3 TtheHE CHARACTERISTICS
THE KERNELS
becaus
本課題では、核融合プラズマ乱流シミュレーションコ
allocation in the Finite-Difference
kernel, and the high it was s[GFlops]
SPARC64 XIfx [21] (“FX100”), which is considered as a kindidling during when they wait for the working threads. The

of
GPGPU is the direct access to the local cache,of SandyreferenceBridge、Xeon
architecture in thebenchmark in Sec. 6.
As advantage
をPhi、および、FX100
で最適
図３
セミラグランジアン法カーネルと差分法カーネル
6.1.1
Xeon Phi(e. g. F
6.1.1 Xeon Phi
6.1.1 Xeon Phi
GT5DFirstly,
it turnedout that
the high cacheFirstly, it turned out that thehigh cache
locality
is essential
porting
goodon XeonTh
for good performance on for
Xeon
Phi.performance
The performanceis Phi.
as hosby appropriate
effective cache
usageout
withimproved by effective cacheimproved
usage withthreadfora
allocationkernel,
in the andFinite-Difference
kern
allocation in the Finite-Difference
the high it was

Let us describe the kernels targeted in this work. The first
GT5D および GYSELAからホットスポットに対応
one isード
the Semi-Lagrangiankernel extracted
from GYSELA
code, and the other is the Finite-Difference kernel extracted
from するカーネルコード（それぞれ差分法カーネルおよびセGT5D code. Both codes solve the same equation withdifferent numerical schemes. The most time-consuming part
0.3
0.3
0.19
0.5
ミラグランジアン法カーネル）を抽出し、GPU
上で最適
TABLE 1is the 4D advection part solved with the Semi-Lagrangian
Hardware description. Thermal Design Power (TDP) is extracted fromschemeand Finite-Difference scheme, respectively. Since化技術の開発を行った。本課題では特に差分法カーネ
表１
ベンチマークで用いたハードウェアの仕様。
venders
data-sheets [20], [22], [23]. As a reference, power efficiency
these two codes are highly parallelized and based on a hydata are derived from the system data of Helios (for Sandy Bridge),
Helios-MIC node (for Xeon Phi), TSUBAME2.5 (for TeslaK20X) and
ルに存在するストライドアクセス、セミラグランジアン法
brid MPI/OpenMP
programming model, appropriate parPlasma simulator (for FX100) in the Green 500 list Nov. 2015 [24].
allel distribution of memory accesses and computations is
カーネルにおけるリストアクセスといった複雑なメモリア
GPU 版のセミラグランジアン法カーネルは、最速とな the key
issue.
Peak B/W [GB/s]
SIMD width
TDP [W]
Power efficiency[GFlops/W]
B/F ratio

51.2
256 bit (AVX)
130
562

320
512 bit
225
1501

250
235
2973

480
256 bit−
1910っており、差分法カーネルは FX100 についで速い。セミクセスパターンに着目し、それぞれレジスタの有効活用、

ラグランジアン法カーネルでは、演算密度が高い上、テテクスチャメモリの活用といった最適化を行った。開発

クスチャメモリによるリストアクセスの高速化が有効でした最適化技術により、GT5D の差分法カーネルおよあったため、GPU で高い性能が得られたと考えられる。

び GYSELA コードのセミラグランジアン法カーネルは、

一方、差分法カーネルは、演算密度が低くメモリバンド

Sandy Bridge において最適化されたものと比べそれ幅がボトルネックとなる上、FX100 において有効であっ

ぞれ 7.55 倍、4.03 倍という高い性能向上を得た。

た容量の大きい共有キャッシュを活用した最適化[2,3]

今後の課題としては、残るコード全体の移植、メモリ

が GPU では有効ではないため、GPU より大きなメモリ

アクセスの削減、およびエクサスケール計算に向けた

バンド幅を有する FX100 で高い性能が得られたと考え

通信コストの削減などが挙げられる。コードの移植につられる。

いては、まず OpenACC による移植を行って、性能測

このベンチマークで得られた成果については、世界定を行ったのち、ボトルネックについては改めて CUDA

最大のスーパーコンピュータ関連学会である SC15 でで移植する予定である。メモリアクセスの削減について

ポスター発表を行い[2]、また IEEE Transactions on

は、今まで演算結果を高次元配列へ格納しメモリアク（様式 20）成果報告書
セスで演算を行っていた箇所を低次元配列への演算で置き換えるなどの対応を行う。通信コストの削減については、アルゴリズムレベルでの変更を検討する必要がある。以下に、GT5D に関する具体例を述べる。GT5D
ではクリロフ部分空間解法によって差分・陰解法行列演算を行うが、そこで発生する通信は、縮約演算を含む。縮約演算については、通信と演算のオーバーラッ
プによる通信コストの削減は困難である。そのため、先行研究によって提案されている省通信クリロフ部分空間解法の適用により、通信回数自体を削減し、通信コストを削減するといったアプローチの検討を開始した。

成果リスト
[1]

朝比祐一, GPU における核融合プラズマ乱流コードの最適化,GTCJapan2015[2]

Y.Asahi et al, Optimization of stencil-based
fusion kernels on Tera-flops many-corearchitectures,

ACM/IEEE

conferenceonSupercomputing (SC15).
[3]

Y. Asahi et al., Optimization of fusion
kernels on accelerators with indirect or
strided memory access patterns, submitted
to IEEE Transactions on Parallel andDistributed Systems.