格子QCDにおけるGPU計算 広大理　尾崎裕介 共同研究者　石川健一

格子QCD計算ではHMC(Hybrid Monte-Carlo)法がよく用いられる 1. Introduction 格子QCD計算ではHMC(Hybrid Monte-Carlo)法がよく用いられる この手法において、差分方程式を解く計算が最も時間のかかる計算 この計算をGPUに計算させることによって加速 solver のアルゴリズムはBi-CGStab法 この際に作成したCUDAのコードに対して行った工夫等について紹介 2009/6/24

Outline Introduction GPUを使って倍精度の結果を得るために ホッピングの計算（行列×ベクトル） 内積の計算 Fortran + Cuda おわりに 2009/6/24

単精度演算を行いながら、得られる結果が倍精度であるような手法があれば理想的 反復改良の手法によって、実は可能 2.1 GPUで倍精度 差分方程式の解は倍精度の精度が必要 しかしGPUは単精度計算が高速 単精度：約1TFlops 倍精度：86.4GFlops 単精度演算を行いながら、得られる結果が倍精度であるような手法があれば理想的 反復改良の手法によって、実は可能 2009/6/24

2.2 単精度倍精度混合Solver 差分方程式(連立1次方程式) 前処理付き反復法による解法（倍精度） ：残差ベクトル M：前処理行列 2009/6/24

2.3 単精度倍精度混合Solver GPUの単精度高速計算！ subroutine usual_solver(...) ⋮ do i = 0,1,2,... p = M * r ! Precondition q = A * p x = x + αp r = r – αq enddo end subroutine 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 GPUの単精度高速計算！ 2009/6/24

2.4 収束の様子 2009/6/24

単精度 Solver のアルゴリズムは Bi-CGStab このような反復法による計算では のような行列とベクトルの積の計算が支配的 3.1 ホッピングの計算 単精度 Solver のアルゴリズムは Bi-CGStab Wilson quark の場合によく用いられる。 このような反復法による計算では　　　　　　のような行列とベクトルの積の計算が支配的 まずはこの計算を高速化 2009/6/24

3.2 の計算の様子 quark : 3×4ベクトル hopping : 12×12行列 2009/6/24

Nvidia Cuda Programming Guide より、 なるべく並列度を上げる 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

3.4 データアクセス量の削減 格子の辺に対応する行列は12×12 しかし、12×12の行列を計算機上に保存しておく必要はない 各μごとにコーディングを行えば3×3行列を4組保存しておけば十分 さらにUはリー群SU(3)の元であり、3×3も必要ない 厳密には 8float 今回は少し楽をして 6complex 2009/6/24

3.5 効率的なメモリアクセス データ量とロード回数 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

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

3.7 coalesced access これまでのデータ構造は coalesced access の条件を満たしていなかった block 0 block 1 … site ０ site １ site ２ site 128 site 0 96Byte 96Byte 96Byte 96Byte 96Byte block 0 block 1 … 1 128 16B 16B 16B このようにデータを並び替えてcoalesced accessの条件を達成 2009/6/24

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 multiprocessor あたりの texture cache 8KB shared memory は16KB No divergent branch !! 2009/6/24

coalesced access and using No texture 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

reduction のコードをkernel 5 → kernel 7 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

4.2 reduction : kernel 7 赤を消すと2冪以外にも 対応可能 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] += sdata[tid + 32]; if (blockSize >= 32) sdata[tid] += sdata[tid + 16]; if (blockSize >= 16) sdata[tid] += sdata[tid + 8]; if (blockSize >= 8) sdata[tid] += sdata[tid + 4]; if (blockSize >= 4) sdata[tid] += sdata[tid + 2]; if (blockSize >= 2) sdata[tid] += sdata[tid + 1]; } if (tid == 0) g_odata[blockIdx.x] = sdata[0]; 赤を消すと2冪以外にも 対応可能 2009/6/24

4.3 内積の高速化 reduction の部分は kernel 7 のコードを用いることにして、各成分同士の複素数×複素数の計算部分をコーディング kernel 5とkernel 7は速度が倍違う → 倍速くなるだろう 2009/6/24

4.4 各成分の計算部分 = このやり方だと約20GFlops a*b = c x y z w 実 虚 a x y z w 実 虚 b x block 0 block 1 … 1 float4 float4 float4 実部 虚部 float4 ar,ai,br,bi; ar = (*a).b[blk].ri[0].c[site]; ai = (*a).b[blk].ri[4].c[site]; br = (*b).b[blk].ri[0].c[site]; bi = (*b).b[blk].ri[4].c[site]; x y z w 実 虚 a x y z w 実 虚 b x y z w 実 虚 c = このやり方だと約20GFlops 1thread 当りの計算 a*b = c 2009/6/24

4.5 各成分の計算部分の改良 = = a*b = c 強引にfloatアクセスに 32GFlops x a x b x c y a y b block 0 … float4 float4 float4 実部 虚部 float* xr,xi,yr,yi; float ar,ai,br,bi; xr = &((*a).b[0].ri[0].c[0].x); xi = &((*a).b[0].ri[4].c[0].x); yr = &((*b).b[0].ri[0].c[0].x); yi = &((*b).b[0].ri[4].c[0].x); ar = xr[i]; ai = xi[i]; br = yr[i]; bi = yi[i]; x a x b = x c 1thread 当りの計算 y a y b = y c 1thread 当りの計算 a*b = c 強引にfloatアクセスに 32GFlops 2009/6/24

内積を計算する際の複素数×複素数の計算部分をfloat4からfloatに変更 4.6 内積の高速化 内積を計算する際の複素数×複素数の計算部分をfloat4からfloatに変更 見方を変えると並列度を上げたことに対応 結果、20GFlopsから32GFlopsにspeed up!! 元の17GFlops から約2倍 全体のsolverも50GFlops → 55GFlops なぜ速くなったかはまだよくわかっていない 単純にfloat4アクセスよりfloatアクセスの方が速いのか？ 並列度が上がった影響なのか？ この結果は内積部分での話。他の演算ではどうか？ texture経由でのfloat4アクセスはどうか？ 2009/6/24

5. Fortran + Cuda 今回のプログラムはFortanからcudaのコードを呼ぶように書いている この書き方には標準がなく、プラットフォームやコンパイラによって異なる 今回はLinux上のIntelコンパイラを使った時の話を紹介 2009/6/24

5.1 Fortran + cuda 基本的には Intel Fortran と C 間の場合と同じ cuda 側の関数宣言時、先頭に「extern “C”」をつける 引数を受け取る際、cuda側ではポインタとして受け取る 配列の順序 コンパイル時にcudart等へのリンク　　　　　　　　　　　　　　等 ただし、Fortran上でGPU上のメモリを扱うことが難しかった (プログラムの先頭でメモリ確保し、後の呼び出しで引数として渡す) 今回はグローバル変数を用いて、Fortran側にデバイスメモリが現れないようにコーディングを行った 2009/6/24

5.2 Fortran+cuda // 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); void s_bicgstab_gpu_(hqfield *rs, hqfield *ps, .. ) { cuhqfield *dr; cudaMalloc((void**)&dr),...); ⋮ cudaMemcpy(ps,dp,...); cudaFree((void *)dr); 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

5.3 Fortran+cuda (Makefile) Intel Fortran + nvcc # Makefile all :: solver_main GPULIB = GPUSolverLIB/libgpusolver.a 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) # GPUSolverLIB/Makefile all :: libgpusolver.a s_bicgstab_gpu.o : s_bicgstab_gpu.cu nvcc –c $< –o $@ libgpusolver.a : s_bicgstab_gpu.o ar cr $@ $< cutil.h をincludeする場合 cudaをcallする場合 2009/6/24

格子QCD計算の分野では大規模連立1次方程式を解く計算が大変 この計算をGPUによって加速 6. おわりに 格子QCD計算の分野では大規模連立1次方程式を解く計算が大変 この計算をGPUによって加速 GPUで高速計算を実現するためにはメモリアクセスにかなり気を配らないといけない Coalessed Access Divergent Warp Bank conflict float access ? 2009/6/24