1. No category

KeplerアーキテクチャGPUにおける高速なSGEMVの実装
椋木 大地†，今村 俊幸†，高橋 大介‡（†理化学研究所計算科学研究機構，‡筑波大学システム情報系）
GEneral Matrix Vector multplication (GEMV）は密行列とベクトルの積：y=αAx+βyを計算する
Basic Linear Algebra Subprograms（BLAS）ルーチンである．本研究ではNVIDIA Keplerアーキテ
クチャGPUを対象に高速なSGEMV（Single-precision GEMV）の実装を行った．行列はColumnmajorで，転置されていない状態とする．高い性能を得るためにブロッキングやベクトルロード命令の
活用を行ったほか，問題サイズによって性能が上下することを防ぐために，マルチプロセッサ数や問題
サイズに応じて最適なスレッド数を自動決定する手法を検討した．これによりスペックの異なる複数種
類のGPUにおいて，CUBLAS等の既存の実装と比べて高速かつ安定した性能を実現した．
2. 実装
2.1. 基本実装
!  基本実装を下図に示す．1つのスレッドブロックは2次元（TxB, TyB）のスレッドを持ち，グリッド内
でスレッドブロックはx次元方向に1次元で起動する．
!  ベクトルxは同じデータへの複数回アクセスが発生するため共有メモリに格納し，ブロッキングを
行った．ベクトルxのデータはスレッドブロック内の全スレッドを使用して共有メモリに読み込む．
!  グローバルメモリ上のデータへのアクセスにはfloat4型によって128bitのベクトルロード命令を使用
した．これによりスレッドレベルでは4x4のレジスタブロッキングがなされている．
!  行列Aへのアクセスにはリードオンリーデータキャッシュを適用した．
!  x次元方向スレッド数はTxB=32とし，y次元方向スレッド数TyBは性能を安定化させるために，GPU
のマルチプロセッサ数と問題サイズによって可変させる手法を検討した（2.2参照）．
1スレッドはfloat4で
4要素ずつアクセス
…
#(1,1) #(2,1)
TxB×TyB×4
…
#(TxB,1) #(1,2)
#(TxB,TyB)
Vector x
テキスト
ベクトルxをスレッドブロッ
ク内のすべてのスレッド
（TxB TyB個）を使って共有
メモリに格納
Thread#(1,2)
ThreadBlock#1
Thread#(1,TyB)
起動するスレッドブロック
はBD=ceil(M/(4 TxB))個
Thread
Block#BD
Matrix A
N
1スレッドブロック
あたり (TxB,TyB)個
のスレッドで構成
Thread
#(TxB,1)
TyB×4
TxB×4
Thread
#(TxB,TyB)
4
4
備考
計算する問題サイズ（行列の行数）
N
問題依存
ワープあたりのスレッド数
TW = 32
マルチプロセッサあたりのレジスタ数
RM = 65536
レジスタ割当単位
Runit = 256
共有メモリ割当単位
Sunit = 256
ワープ割当粒度
Wgran = 4
マルチプロセッサあたりの最大ワープ
数（物理的制約）
WM-pmax = 64
マルチプロセッサあたりの最大スレッ
ドブロック数（物理的制約）
BM-pmax = 16
マルチプロセッサあたりの共有メモリ
量 [Bytes]
SM = 49152
デバイスあたりのマルチプロセッサ数
MD
スレッドあたりの使用レジスタ数
RT = 39
ベクトルxの読込に用いるロード命令
のベクトル長
NxT = 4
行列Aの読込に用いるロード命令のベ
クトル長
NaT = 4
1要素の大きさ[Bytes]
None = 4
スレッドブロックあたりのx次元方向
スレッド数（blockDim.x）
TxB = 32
スレッドブロックあたりのy次元方向
スレッド数（blockDim.y）
TyB = {4~16}
ブロックあたりのワープ数
WB = TxB × TyB / TW
ワープあたりの使用レジスタ数
RW = ceil (RT × TW, Runit)
マルチプロセッサあたりの最大ワープ
数（レジスタに制約）
WM-rmax = min (floor (RM /
RW, Wgran), WM-pmax)
マルチプロセッサあたりの最大スレッ
ドブロック数（レジスタに制約）
BM-rmax
= floor (WM / WB, 1)
スレッドブロックあたりの共有メモリ
使用量 [Bytes]
SB = ceil(TxB × TyB × NxT ×
None × NaT/NxT, Sunit)
マルチプロセッサあたりの最大スレッ
ドブロック数（共有メモリに制約）
BM-smax = floor (SM / SB, 1)
マルチプロセッサあたりの最大スレッ
ドブロック数
BM-max = min
(BM-rmax, BM-smax, BM-pmax)
デバイスあたりの最大スレッドブロッ
ク数
BD-max = BM-max × MD
計算に必要なデバイスあたりの総起動
スレッドブロック数
BD
= ceil (N / (TxB × NaT), 1)
スレッドブロック占有率 [%]
OB = BD / ceil (BD, BD-max)
1
0.8
0.6
%
0.4
0.2
0
5000
10000
15000
20000
25000
0
30000
1
0.8
0.6
%
GFlops
78
76
74
72
70
68
66
64
0.4
0.2
0
5000
10000
15000
20000
25000
GK110
GK110
GK110B
GPU Clock (GPU Boost)
706 MHz
732 MHz
837 (876) MHz
889 (980) MHz
Multiprocessors (SMX)
13
14
14
15
CUDA Cores
2496 (192x13)
2688 (192x14)
2688 (192x14)
2880 (192x15)
Single-Precision Flops
3.52 TFlops
3.95 TFlops
4.5 TFlops
5.12 TFlops
Memory
320-bit GDDR5
384-bit GDDR5
384-bit GDDR5
384-bit GDDR5
Memory Capacity
5120 MB (*)
6144 MB (*)
6144 MB
6144 MB
Memory Clock
5200 MHz
5200 MHz
6008 MHz
7000 MHz
Memory Bandwidth
208 GB/s (*)
250 GB/s (*)
288.4 GB/s
336 GB/s
ECC Support
Yes
Yes
N/A
N/A
CPU
Xeon E5-2609 2.4GHz
Xeon E5-2609 2.4GHz
Xeon E5-2650 v2 2.6GHz
Core i7 920 2.67GHz
CUDA Driver
331.62
331.67
331.67
331.49
OS
CentOS 6.4
CentOS 6.5
CentOS 6.5
CentOS 5.10
SGEMV on Tesla K20Xm
100
50
0.8
0.6
0.4
0.2
5000
10000
15000
20000
25000
30000
40
30
20
MUBLAS
0.8!0.8
MUBLAS
CUBLAS
6.0"6.0
CUBLAS
KBLAS
1.0! 1.0
KBLAS
MAGMA
MAGMA1.4.1!
1.4.1
10
0
1
0
80
60
TyB=5
TyB=5
TyB=6
TyB=6
TyB=7
TyB=7
TyB=8
TyB=8
TyB=9
TyB=9
TyB=10
TyB=10
TyB=11
TyB=11
TyB=12
TyB=12
0
0
5000
10000
15000
20000
25000
"!
"!
"!
"!
60
40
MUBLAS
0.8!
MUBLAS
0.8
CUBLAS
6.0"
CUBLAS 6.0
KBLAS
1.0! 1.0
KBLAS
MAGMA
MAGMA1.4.1!
1.4.1
20
0
30000
0
5000
10000
15000
20000
25000
30000
N
N
SGEMV on GeForce GTX Titan
SGEMV on GeForce GTX Titan Black
140
160
120
140
"!
"!
"!
"!
35000
1
0.8
0.6
0.4
0.2
0
5000
10000
15000
20000
25000
30000
TyB=13
TyB=13
TyB=14
TyB=14
TyB=15
TyB=15
TyB=16
TyB=16
0
0.6
%
GFlops
0.8
0.4
0.2
5000
10000
15000
80
60
40
1
0
100
20000
25000
30000
TyB=17
TyB=17
TyB=18
TyB=18
TyB=19
TyB=19
TyB=20
TyB=20
MUBLAS
0.8!0.8
MUBLAS
CUBLAS
6.0"6.0
CUBLAS
KBLAS
1.0! 1.0
KBLAS
MAGMA
MAGMA1.4.1!
1.4.1
20
0
0
5000
0.8
0.4
0.2
TyB=21
TyB=21
TyB=22
TyB=22
TyB=23
TyB=23
TyB=24
TyB=24
0
5000
10000
15000
20000
25000
30000
1
0.8
0.6
%
GFlops
78
76
74
72
70
68
66
64
0.4
0.2
TyB=25
TyB=25
TyB=26
TyB=26
TyB=27
TyB=27
TyB=28
TyB=28
0
0
5000
10000
15000
20000
25000
30000
1
0.8
0.6
%
GFlops
N
78
76
74
72
70
68
66
64
0.4
0.2
TyB=29
TyB=29
TyB=30
TyB=30
TyB=31
TyB=31
TyB=32
TyB=32
0
0
5000
10000
15000
20000
25000
30000
N
実装依存の定数
図2. TyBを変化させた時の性能（実線）とOB（破線）（Tesla K20m）
TyBを4∼16の中から選択
可変
0.8
25
0.6
20
15
0.4
10
0.2
0
OB
TyB
30
TyB
%
1
5
0
5000
10000
15000
20000
25000
30000
0
N
OBを最大化 するTyBを選択
図3. 選択されたTyBとOB（Tesla K20m）
78
76
74
72
70
68
66
64
Best
Obtained
0
5000
10000
15000
20000
25000
40
MUBLAS
0.8!
MUBLAS
0.8
CUBLAS
6.0"
CUBLAS 6.0
KBLAS
1.0! 1.0
KBLAS
MAGMA
MAGMA1.4.1!
1.4.1
20
0
0
N
1
0.6
"!
"!
"!
"!
10000 15000 20000 25000 30000 35000
0
N
78
76
74
72
70
68
66
64
80
60
N
78
76
74
72
70
68
66
64
120
100
GFlops
78
76
74
72
70
68
66
64
%
GFlops
N
N
GPU依存の定数
GK110
0
30000
%
GFlops
78
76
74
72
70
68
66
64
0
アーキテクチャ
依存の定数
Code Name
70
TyB=1
TyB=1
TyB=2
TyB=2
TyB=3
TyB=3
TyB=4
TyB=4
GFlops
GFlops
78
76
74
72
70
68
66
64
%
値・計算式
GeForce GTX Titan Black
80
N
GFlops
項目
表1. TyBの選択に必要なパラメータと計算式
GeForce GTX Titan
SGEMV on Tesla K20m
N
GFlops
!  上記の実装を任意のTyBで起動すると問題サイズ
に対して性能が周期的に変動する．この原因の
一つは，起動するスレッドブロック数が，GPU
のスレッドブロックの処理単位で割り切れない
数となるからである [1]．
!  GPUのスレッドブロックの処理単位は，マルチ
プロセッサの個数とマルチプロセッサあたりの
最大スレッドブロック数から決まる．後者はTyB
によって変化する．
!  本研究ではGEMVの性能を，全マルチプロセッ
サに対するスレッドブロックの占有率という形
でモデル化し，問題サイズに対する性能の周期的
変動を再現した．そして問題サイズに応じてス
レッドブロックの占有率が最大となるようなTyB
を選択することで，性能の変動を最小限にする
方法をとった．
!  実験の結果，TyB>16ではN>24000近辺から性
能低下が見られた．またTyBが小さすぎると物理
コア数に対してスレッド数が不足する．そこで
TyBを4 16（1スレッドブロックあたり
128 512スレッド）から選択するようにした．
Tesla K20Xm
!  4種類のGPUにおいて，既存の実装を上回る性能と，問題サイズに対する性能の安定性を実現した．
!  N=20000におけるメモリスループットと理論メモリバンド幅に対する効率は，Tesla K20m：
146.9GB/s（71%），Tesla K20Xm：179.0GB/s（72%），GeForce GTX Titan：255.1GB/
s（88%），GeForce GTX Titan Black：2.97GB/s（88%）であった．（TeslaシリーズGPUの理
論メモリバンド幅はECC無効時の値を用いて計算．ECCにより理論メモリバンド幅が12.5%低下する
と仮定すると，それぞれ理論値の81%，82%の効率となる）
図1. 実装の概要
2.2. 性能の安定化
Tesla K20m
3.2. 結果
Thread#(1,1)
float4
…
M
…
Thread#(2,1)
GPU
(*) ECC無効時の値．ECC有効時はメモリを12.5%使用する
1スレッドはfloat4で
4要素ずつアクセス
Thread#(1,1)
…
Vector y
ベクトルyの1要素を計算する
ためにTyB個のスレッドで総
和を計算（共有メモリでス
レッド間の値を交換）
…
Thread
Block#1
Thread
Block#2
表2. 評価環境
GFlops
float4
!  CUDA6.0環境において我々の実装（MUBLAS 0.8）とCUBLAS 6.0[2]，KBLAS 1.0[3][4]，MAGMA
1.4.1[5]の性能を比較した．（MAGMAの最新版1.5.0beta2のSGEMV性能は1.4.1と同一であった）
!  性能評価に用いたGPUは，Tesla K20m, Tesla K20Xm, GeForce GTX Titan, GeForce GTX Titan
Blackの4種類．このうちTeslaシリーズのGPUではECCを有効化している．
!  nvccのコンパイルオプションは -arch=sm_35 -O3 とした．
!  y=αAx+βyにおいて行列は大きさNxNの正方行列，入力データはα・βも含め乱数で初期化．Nは1
ずつ変化させ，GPUのグローバルメモリを使い切るサイズまで測定した．
!  性能はFlopsで表し，10回実行したうちの最も良い値を採用した．実行前にデバイスのマルチプロセッ
サ数を取得する初期化関数を1度だけ実行するが，この実行時間は性能には含めていない．
!  行列を格納するGPU上のグローバルメモリの確保にはcudaMallocPitchを用い，lda（leading
dimension of matrix A）が128bit単位になるようにした．（MUBLASの実装はベクトルロード命令
float4を使用するためにldaが128bit単位であることを要求するが，CUDAにおいて2次元配列を確保
する際にはcudaMallocPitchの使用が推奨されている）
GFlops
Thread#(1,1)
3.1. 評価手法
Host
4
3. 性能評価
GPU
1. 概要
30000
N
図4. TyB=1 32における最高性能（Best）と，選択された
TyBによって得られた性能（Obtained）（Tesla K20m）
5000
"!
"!
"!
"!
10000 15000 20000 25000 30000 35000
N
図4. 性能評価の結果
4. まとめと今後の課題
!  複数のKeplerアーキテクチャGPUにおいて，既存の実装と比べて高いピーク性能，および問題サイズ
に対する性能の変動の少ない安定した性能を実現した．
!  問題サイズに対する周期的な性能の変動を抑えるために，問題サイズによってスレッド数を可変させ
る方法を検討した．理論的な性能モデルに基づきスレッド数を自動で決定するため，実性能の測定結
果をフィードバックさせる手法[6][7]と比べ，チューニングに時間を要さず，個々のGPUごとのチュー
ニングも必要としない．
!  性能モデルにより性能の周期性を再現することができたが，性能の上下を完全に再現するには至って
おらず，性能モデルには改善の余地がある．
!  他のLevel-2ルーチンや，異なる精度，Kepler以外のアーキテクチャ（Fermi，Maxwellなど）への適
用が今後の課題である．
● 参考文献
[1] 今村俊幸: CUDA環境下でのDGEMV関数の性能安定化・自動チューニングに関する考察 , 情報処理学会論文誌コンピューティングシス
テム, Vol.4, No.4, pp.158-168 (2011).
[2] NVIDIA Corporation: The NVIDIA CUDA Basic Linear Algebra Subroutines (cuBLAS) , https://developer.nvidia.com/cublas.
[3] Abelfatteh, A.: KAUST BLAS (KBLAS) , http://cec.kaust.edu.sa/Pages/Abelfatteh.aspx.
[4] Abelfatteh, A., Keyes, D., Ltaief, H.: Systematic approach in optimizing numerical memory-bound kernels on GPU ,
Proceedings of the 18th international conference on Parallel processing workshops (Euro-Par 12), LNCS 7640, pp. 207-216
(2013).
[5] University of Tennessee: Matrix Algebra on GPU and Multicore Architectures (MAGMA) , http://icl.cs.utk.edu/magma.
[6] 今村俊幸, 林熙龍, 内海貴弘, 山田進, 町田昌彦: Fermi,Kepler複数世代GPUに対するSYMVカーネルの性能チューニング , 情報処理学会
研究報告, 2013-HPC-138(7), pp. 1-7 (2013).
[7] Sørensen, H. H. B.: Auto-tuning Dense Vector and Matrix-Vector Operations for Fermi GPUs , Parallel Processing and
Applied Mathematics (PPAM2011), LNCS Vol. 7203, pp. 619-629 (2012).
● 謝辞
本研究の一部は，JSPS特別研究員奨励費（課題番号 251290）ならびにJST CREST「進化的アプローチによる超並列複合システム向け開
発環境の創出」の助成によるものである．
2014年7月16日 GTC Japan 2014