1 / 27

格子 QCD における GPU 計算

格子 QCD における GPU 計算. 広大 理　尾崎裕介 共同研究者　石川健一. 1. Introduction. 格子 QCD 計算では HMC(Hybrid Monte-Carlo) 法がよく用いられる この手法において、差分方程式を解く計算が最も時間のかかる計算 この計算を GPU に計算させることによって加速 solver のアルゴリズムは Bi- CGStab 法 この際に作成した CUDA のコードに対して行った工夫等について紹介. Outline. Introduction GPU を使って倍精度の結果を得るために

penelope-lyons
Télécharger la présentation

格子 QCD における GPU 計算

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

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

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

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

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

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

  6. 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 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の単精度高速計算！

  7. 2.4 収束の様子

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

  9. 3.2 の計算の様子 quark : 3×4ベクトル hopping : 12×12行列

  10. 3.3 Basic strategy NvidiaCuda Programming Guide より、 • なるべく並列度を上げる • 1 thread あたりの計算を 1 格子点の計算に assign • (12×12行列) × (12次元ベクトル) • メモリアクセスの最適化 • global memory に対するアクセス速度が遅い • 400～600 clock • shared memory、texture cache 等の有効利用 • 4 clock • 計算によって生成できるデータはGPU上で計算 • memory access が遅いので計算した方が速い場合がある

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

  12. 3.5 効率的なメモリアクセス　データ量とロード回数 • 計2回のロード • (6×2 float = 48Byte) × 4set • block 内ではshared memory を使ってthread 間のmemory 共有ができる • shared memory : 16KByte • できるだけ多くの格子点を共有したい • ロード回数の多い格子点をshared • 43×2=128の格子点 = 12.6KByte • リンクは texture を用いた • 計8回のロード • 12×2 float = 96Byte

  13. 3.6 solverの性能 • この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の性能

  14. 3.7 coalesced access これまでのデータ構造は coalesced access の条件を満たしていなかった 96Byte 96Byte 96Byte 96Byte 96Byte 16B 16B 16B このようにデータを並び替えてcoalesced accessの条件を達成

  15. 3.8 divergent branch • 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 !! • ある格子点(K)を計算する際、その隣の格子点(S)にあるデータが必要 • shared memory を用いる場合

  16. 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 • 他の計算（内積等）が足を引っ張っている

  17. 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冪以外に対応するように改造

  18. 4.2 reduction : kernel 7 template <unsigned intblockSize> __global__ void reduce6(int *g_idata, int *g_odata, unsigned int n) { extern __shared__ intsdata[]; unsigned inttid = threadIdx.x; unsigned inti = blockIdx.x* (blockSize*2) + tid; unsigned intgridSize = 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冪以外にも 対応可能

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

  20. 4.4 各成分の計算部分 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]; = 1thread 当りの計算 a*b = c このやり方だと約20GFlops

  21. 4.5 各成分の計算部分の改良 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]; = 1thread 当りの計算 = 1thread 当りの計算 a*b = c • 32GFlops 強引にfloatアクセスに

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

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

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

  25. // s_bicgstab_gpu.cu cuhgvfielddA; 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); ⋮ } 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

  26. 5.3 Fortran+cuda (Makefile) Intel Fortran + nvcc # GPUSolverLIB/Makefile all :: libgpusolver.a s_bicgstab_gpu.o : s_bicgstab_gpu.cu nvcc –c $< –o $@ libgpusolver.a : s_bicgstab_gpu.o arcr $@ $< # 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) : (cdGPUSolverLIB; make) cutil.hをincludeする場合 cudaをcallする場合

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

More Related