格子QCDにおけるGPU計算
広大理 尾崎裕介
共同研究者 石川健一
1. Introduction
 格子QCD計算ではHMC(Hybrid Monte-Carlo)法がよ
く用いられる
 この手法において、差分方程式を解く計算が最も時
間のかかる計算
 この計算をGPUに計算させることによって加速
solver のアルゴリズムはBi-CGStab法
 この際に作成したCUDAのコードに対して行った工夫
等について紹介
2009/6/24
2
Outline
1. Introduction
2. GPUを使って倍精度の結果を得るために
3. ホッピングの計算（行列×ベクトル）
4. 内積の計算
5. Fortran + Cuda
6. おわりに
2009/6/24
3
2.1 GPUで倍精度
 差分方程式の解は倍精度の精度が必要
 しかしGPUは単精度計算が高速
単精度：約1TFlops
倍精度：86.4GFlops
 単精度演算を行いながら、得られる結果が倍精度で
あるような手法があれば理想的
 反復改良の手法によって、実は可能
2009/6/24
4
2.2 単精度倍精度混合Solver
 差分方程式(連立1次方程式)
 前処理付き反復法による解法（倍精度）
：残差ベクトル
M：前処理行列
 単精度倍精度混合Solver
2009/6/24
（単精度）
5
2.3 単精度倍精度混合Solver
subroutine usual_solver(...)
⋮
do i = 0,1,2,...
⋮
p = M * r
! Precondition
q = A * p
x = x + αp
r = r – αq
⋮
enddo
⋮
end subroutine
GPUの単精度高速計算！
2009/6/24
subroutine mixed_solver(...)
⋮
As = A
! 倍→単 変換
⋮
do i = 0,1,2,...
⋮
rs = r
! 倍→単 変換
ps = (As)-1 * rs ! 単精度計算
p = ps
! 単→倍 変換
q = A * p
x = x + αp
r = r – αq
⋮
enddo
⋮
end subroutine
6
2.4 収束の様子
2009/6/24
7
3.1 ホッピングの計算
 単精度 Solver のアルゴリズムは Bi-CGStab
Wilson quark の場合によく用いられる。
 このような反復法による計算では
行列とベクトルの積の計算が支配的

のような
 まずはこの計算を高速化
2009/6/24
8
3.2
の計算の様子
quark : 3×4ベクトル
hopping : 12×12行列
2009/6/24
9
3.3 Basic strategy
Nvidia Cuda Programming Guide より、
 なるべく並列度を上げる
1 thread あたりの計算を 1 格子点の計算に assign
 (12×12行列) × (12次元ベクトル)
 メモリアクセスの最適化
global memory に対するアクセス速度が遅い
 400～600 clock
shared memory、texture cache 等の有効利用
 4 clock
計算によって生成できるデータはGPU上で計算
 memory access が遅いので計算した方が速い場合がある
2009/6/24
10
3.4 データアクセス量の削減
 格子の辺に対応する行列は12×12
 しかし、12×12の行列を計算機上に保存しておく必要はない
 各μごとにコーディングを行えば3×3行列を4組保存しておけば十分
 さらにUはリー群SU(3)の元であり、3×3も必要ない
 厳密には 8float
 今回は少し楽をして 6complex
2009/6/24
11
3.5 効率的なメモリアクセス
データ量とロード回数
 計8回のロード
12×2 float = 96Byte
 計2回のロード
(6×2 float = 48Byte) × 4set
 block 内ではshared memory を使っ
てthread 間のmemory 共有ができる
 shared memory : 16KByte
 できるだけ多くの格子点を共有したい
 ロード回数の多い格子点をshared
 43×2=128の格子点 = 12.6KByte
 リンクは texture を用いた
2009/6/24
12
3.6 solverの性能
 その他Bi-CGStab法で必要な演算（複素数の内積等）をSDK
等を参考に作成
 これらを組み合わせてGPUを使った単精度solverを作成
 even/odd preconditioning ← 格子サイズを半分にする手法
 clover term の計算 ← 格子間隔による誤差を減少させる
 さらに前処理付きBi-CGStab倍精度solverとマージ
 このようにしてできたsolverを実際に走らせてみた
 GPU : NVIDIA GeForce GTX 280
 CPU : Intel Core i7 2.67GHz
 結果約20GFlopsの性能
 このsolverを以下のようにして高速化した
 これらの手法を紹介
2009/6/24
13
3.7 coalesced access
これまでのデータ構造は coalesced access の条件を満たしていなかった
block 0
site ０
site １
site ２
96Byte
96Byte
96Byte
block 1
…
site 128
site 0
96Byte
96Byte
…
block 0
0
1
16B
16B
…
128
0
1
…
…
block 1
…
0
1
…
0
1
…
…
…
16B
このようにデータを並び替えてcoalesced accessの条件を達成
2009/6/24
14
3.8 divergent branch
 ある格子点(K)を計算する際、その隣の格子点(S)にあるデー
タが必要
 shared memory を用いる場合
 if (SはKと同じブロック) → shared memory load
else → global memory load
divergent branch !!
 代わりに texture fetch を使う
 任意の格子点 → texture fetch load
No divergent branch !!
 multiprocessor あたりの texture cache 8KB
 shared memory は16KB
2009/6/24
15
3.9 ここまでの結果
 coalesced access and using No texture
 40GFlops
 coalesced access and using texture
 50GFlops
 coalesced access にすることで 速度は倍以上!
 texture fetch により、さらに速度up!
ただし、データはcacheに乗りきらない
 格子点12.3KB + リンク 24.6KB > 8KB
それでも一部のデータはcacheの恩恵を受けたため？
 特にホッピングの計算部分の性能は107GFlops
他の計算（内積等）が足を引っ張っている
2009/6/24
16
4.1 内積の高速化
 内積について高速化を行った
 現在の内積は17GFlops
bandwidthTest → 114GByte/s (GeForce GTX 280)
内積計算 : 2Byte/Flop
理論性能 : 57GFlops
 reduction のコードをkernel 5 → kernel 7
NVIDIA_CUDA_SDK/projects/reduction/doc/reduction.pdf
要素数2冪以外に対応するように改造
2009/6/24
17
4.2 reduction : kernel 7
template <unsigned int blockSize>
__global__ void reduce6(int *g_idata, int *g_odata, unsigned int n)
{
extern __shared__ int sdata[];
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x* (blockSize*2) + tid;
unsigned int gridSize = blockSize*2*gridDim.x;
sdata[tid] = 0;
while (i < n) { sdata[tid] += g_idata[i] + g_idata[i+ blockSize]; i += gridSize; }
__syncthreads();
if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); }
if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); }
if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); }
if (tid < 32) {
if (blockSize >= 64) sdata[tid]
if (blockSize >= 32) sdata[tid]
if (blockSize >= 16) sdata[tid]
if (blockSize >= 8) sdata[tid]
if (blockSize >= 4) sdata[tid]
if (blockSize >= 2) sdata[tid]
}
if (tid == 0) g_odata[blockIdx.x]
}
2009/6/24
+=
+=
+=
+=
+=
+=
sdata[tid
sdata[tid
sdata[tid
sdata[tid
sdata[tid
sdata[tid
+ 32];
+ 16];
+ 8];
+ 4];
+ 2];
+ 1];
赤を消すと2冪以外にも
対応可能
= sdata[0];
18
4.3 内積の高速化
 reduction の部分は kernel 7 のコードを用いることにして、各
成分同士の複素数×複素数の計算部分をコーディング
 kernel 5とkernel 7は速度が倍違う → 倍速くなるだろう
2009/6/24
19
4.4 各成分の計算部分
block 0
0
1
float4 float4
…
0
…
block 1
0
…
0
…
0
1
…
…
…
float4
実部
虚部
x x
x x
x x
float4 ar,ai,br,bi;
y y
y y
y y
z z
z z
w w
w w
w w
ar
ai
br
bi
実虚
実虚
実虚
a
b
c
=
z z
=
=
=
=
(*a).b[blk].ri[0].c[site];
(*a).b[blk].ri[4].c[site];
(*b).b[blk].ri[0].c[site];
(*b).b[blk].ri[4].c[site];
 このやり方だと約20GFlops
1thread 当りの計算
a*b = c
2009/6/24
20
4.5 各成分の計算部分の改良
block 0
…
float4
float4
…
float4
実部
x x
x x
a
b
虚部
=
x x
c
1thread 当りの計算
y y
y y
a
b
=
y y
c
1thread 当りの計算
a*b = c
2009/6/24
float* xr,xi,yr,yi;
float ar,ai,br,bi;
xr
xi
yr
yi
=
=
=
=
&((*a).b[0].ri[0].c[0].x);
&((*a).b[0].ri[4].c[0].x);
&((*b).b[0].ri[0].c[0].x);
&((*b).b[0].ri[4].c[0].x);
ar
ai
br
bi
=
=
=
=
xr[i];
xi[i];
yr[i];
yi[i];
 強引にfloatアクセスに
 32GFlops
21
4.6 内積の高速化
 内積を計算する際の複素数×複素数の計算部分を
float4からfloatに変更
見方を変えると並列度を上げたことに対応
 結果、20GFlopsから32GFlopsにspeed up!!
元の17GFlops から約2倍
 全体のsolverも50GFlops → 55GFlops
 なぜ速くなったかはまだよくわかっていない
単純にfloat4アクセスよりfloatアクセスの方が速いのか？
並列度が上がった影響なのか？
この結果は内積部分での話。他の演算ではどうか？
texture経由でのfloat4アクセスはどうか？
2009/6/24
22
5. Fortran + Cuda
 今回のプログラムはFortanからcudaのコードを呼ぶように書
いている
 この書き方には標準がなく、プラットフォームやコンパイラに
よって異なる
 今回はLinux上のIntelコンパイラを使った時の話を紹介
2009/6/24
23
5.1 Fortran + cuda
 基本的には Intel Fortran と C 間の場合と同じ
 cuda 側の関数名の後ろにアンダースコア「＿」を１つ付け加える
 cuda 側の関数宣言時、先頭に「extern “C”」をつける
 引数を受け取る際、cuda側ではポインタとして受け取る
 配列の順序
 コンパイル時にcudart等へのリンク
等
 ただし、Fortran上でGPU上のメモリを扱うことが難しかった
 (プログラムの先頭でメモリ確保し、後の呼び出しで引数として渡す)
 今回はグローバル変数を用いて、Fortran側にデバイスメモリ
が現れないようにコーディングを行った
2009/6/24
24
5.2 Fortran+cuda
subroutine use_gpu_solver(...)
⋮
As = A
! 倍→単 変換
call new_gpu_solver(As,dA)
⋮
do i = 0,1,2,...
⋮
rs = r
! 倍→単 変換
call s_bicgstab_gpu(rs,ps,dA,...)
p = ps
! 単→倍 変換
q = A * p
x = x + αp
r = r – αq
⋮
enddo
call delete_gpu_solver(dA)
⋮
end subroutine
2009/6/24
// s_bicgstab_gpu.cu
cuhgvfield dA;
extern “C”
void new_gpu_solver_(cuhgvfield *As) {
cudaMalloc((void**)&dA,...);
cudaMemcpy(dA,As,...);
}
extern ”C”
void delete_gpu_solver_() {
cudaFree((void *)dA);
}
extern “C”
void s_bicgstab_gpu_(hqfield *rs,
hqfield *ps,
..
)
{
cuhqfield *dr;
cudaMalloc((void**)&dr),...);
⋮
cudaMemcpy(ps,dp,...);
cudaFree((void *)dr);
⋮
}
25
5.3 Fortran+cuda (Makefile)
Intel Fortran + nvcc
# Makefile
all :: solver_main
# GPUSolverLIB/Makefile
all :: libgpusolver.a
GPULIB = GPUSolverLIB/libgpusolver.a
s_bicgstab_gpu.o : s_bicgstab_gpu.cu
nvcc –c $< –o $@
LIBDIR = -L$(CUDADIR)/lib \
–L$(CUDASDKDIR)/lib \
-L$(CUDASDKDIR)/common/lib
LIBS = $(LIBDIR) –lcudart –lcutil -lm
solver_main : $(OBJ) $(GPULIB)
ifort $(LIBS) $(OBJ) $(GPULIB) \
–o $@
$(GPULIB) :
(cd GPUSolverLIB; make)
2009/6/24
libgpusolver.a : s_bicgstab_gpu.o
ar cr $@ $<
cutil.h をincludeする場合
cudaをcallする場合
26
6. おわりに
 格子QCD計算の分野では大規模連立1次方程式を
解く計算が大変
 この計算をGPUによって加速
 GPUで高速計算を実現するためにはメモリアクセス
にかなり気を配らないといけない
Coalessed Access
Divergent Warp
Bank conflict
float access ?
2009/6/24
27

TID News! - 株式会社ティ・アイ・ディ

解答編 完全な解答ではありません。 各自途中の文章

外注管理の課題解決セミナー 【外注QCDマネジメントの

解答 編 完全な解答ではありません。

ISO9001内部品質監査員養成コース

設計・開発部門“見える化”改革

CUDA並列処理 プログラミング初体験

第1章 なんと25,000円! 手のひら スパコンJetson TK1初体験!

RIC-Nays

第2章 最新GPUのメカニズム

VM-FMS