Transcript PPTX file
格子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 0
site 1
site 2
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 側の関数名の後ろにアンダースコア「_」を1つ付け加える
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