Presentation is loading. Please wait.

Presentation is loading. Please wait.

格子QCD計算のGPUを用いた加速の可能性について

Similar presentations


Presentation on theme: "格子QCD計算のGPUを用いた加速の可能性について"— Presentation transcript:

1 格子QCD計算のGPUを用いた加速の可能性について
石川健一(広島大理) @CCSコロキウム, 筑波大学, 2009年6月19日

2 1. 目次 2.格子QCDについて 3.ハイブリッドモンテカルロ法 4.GPUを用いたソルバーの加速
単精度前処理 ベンチマーク結果 5.GPUを用いたソルバーの並列化への展望 6.まとめ

3 2. 格子QCDについて 陽子、中性子、π中間子などに代表される粒子(ハドロン)の性質を第一原理から理解したい。
量子色力学(Quantum ChromoDynamics[QCD]) 基礎方程式は電磁気学Maxwell方程式の拡張          => Yang-Mills 方程式

4 時空を格子化した有限自由度の格子QCDという方法を用いる。
2. 格子QCDについて QCDを解くためには 時空を格子化した有限自由度の格子QCDという方法を用いる。 K.G.Wilson (1974) 連続時空上の場の変数 → 格子上の場の変数 格子化した時空 連続時空

5 2. 格子QCDについて 分配関数(ユークリッド化された経路積分) 物理量の期待値 統計力学で用いられてきた方法 作用 有効作用
有効作用を重みとする多重積分 => モンテカルロ積分 統計力学で用いられてきた方法

6 2. 格子QCDについて 特徴、 自由度がとても大きい 16^4 ~ 32^4 格子 系統誤差
自由度がとても大きい  16^4 ~ 32^4 格子 グルーオン:8x4x(16^4~32^4)=200万~3000万自由度(実数勘定) クォーク:3x4x (16^4~32^4) =160万~1500万自由度(実数勘定) 系統誤差 格子化に伴う誤差:格子間隔 a = 0.1 fm ~ 0.06 fm → 0 fm 有限体積による誤差: 核子が大体収まる大きさ L= 3fm~4fm →∞ クォークの質量が現実世界と異なるための誤差       m_q=40MeV~100MeV →  m_ud = 2MeV~ 10MeV どのようにしてモンテカルロ積分のための配位を生成するか?

7 3. ハイブリッドモンテカルロ(HMC)法 格子QCD分配関数 ユークリッド化されているので、統計力学の分配関数と同等の計算になる。

8 クォーク部分の計算が大変 3.ハイブリッドモンテカルロ(HMC)法
真空からクォーク-反クオーク対が対生成、対消滅している効果を表している。真空偏極。

9 Hybrid Monte Carlo(HMC)法
分子動力学法+メトロポリステスト   U を重み                    で生成 計算が困難な点 分子動力学法の部分で の大規模連立1次方程式を何回も解く必要がある。 幸いD[U]の要素はほとんどゼロ(大規模疎行列)であるので、反復法 を用いて連立方程式を解く。 HMC法で最も時間がかかる部分 クォークの質量が軽いとD[U]の条件は悪化する K(D)∝1/mq 全計算時間の80-90%を占めるようになる。 クォークソルバー クォークソルバーの加速が重要

10 格子QCD計算で最も時間のかかるクォークソルバーの加速をGPUを用いて出来ないか考えた。 GPU:グラフィックプロセッシングユニット
3.ハイブリッドモンテカルロ(HMC)法 格子QCD計算で最も時間のかかるクォークソルバーの加速をGPUを用いて出来ないか考えた。 GPU:グラフィックプロセッシングユニット パソコンについているグラフィックカードのプロセッサ メーカー : AMD と Nvidia  2強? 近年とても速くなっている 単精度で 1 TFlops, 倍精度でも 100GFlops 程度 グラフィック計算以外の一般的な計算も出来るようになってきた 開発言語: Nvidia CUDA, AMD CAL Brook+.. Nvidia CUDAが先行 コモディティなので安いやつは安い(ハイエンドゲームマシン用) HPC向けのも売り出されている(ゲーム用と比べると割高)

11 GPU例 NVIDIA GT200 (Tesla 10series) 240 SP (SP cores), 30 DP cores
~1,000(or 600)Glops(SP), ~90GFlops(DP) High Memory bandwidth 120 GB/s

12 先行研究(QCD) GPUの強力な計算力をQCDに使えないか? Nvidia CUDA以前 Nvidia CUDA後
G.I.Egri, et al.,”Lattice QCD as a video game” Comput. Phys. Commun. 177(2007)631. Nvidia CUDA後 K.Barros et al. Latice [hep-lat/ ]. M.Clark et al. Work March, 2009. AMD 製 GPU 例(つい最近) V.Demchik and A.Strelchenko, hep-lat/ v2

13 今回の話 Nvidia製GPU (GeForce GTX 280)を使って Lattice QCD のクォークソルバーの加速を試みた。その結果の紹介 PC一台での計算(並列化はまだ) 並列化(複数PCや複数GPUカード)に向けての試行錯誤の紹介 今後の展望

14 4.GPUを用いたソルバーの加速 解きたい連立一次方程式: Dx = b 一階の差分方程式
行列サイズ:格子点数を N^4 とすると 12xN^4 大規模疎行列 反復法で解 x を求める

15 4.GPUを用いたソルバーの加速 GPUは単精度が圧倒的に速い (peak 1TFlops) 単精度を利用して解く
しかし解の精度は倍精度に保つ 反復改良法/Richardson反復のテクニックで可能 一般の反復法にも組み入れることが出来る(可変前処理) 前処理として単精度ソルバーを使う 残差は倍精度で正しくなるように組む FGMRES, GCR, CG, BiCGStab….. (with flexible prec.)

16 Mixed precision / Inner-Outer solver
4.GPUを用いたソルバーの加速 Mixed precision / Inner-Outer solver Flexible Preconditioner 可能な反復法で計算 Richardson 反復(Iterative refinement)法を基に作る より複雑な反復法の中に組み入れることも可能 GCR法、FGMRES法、CG法、BiCGStab法 今回の計算では Nested-BiCGStab (櫻井-多田野)を用いた

17 4.GPUを用いたソルバーの加速 Nvidia CUDA を用いたプログラミング 基本は超並列プログラム 基礎言語はC言語
情報源:[ 基本は超並列プログラム 基礎言語はC言語 ベクトル並列計算機向プログラムを移植するのが簡単 ベクトルループ => CUDA スレッド並列 MPI並列 => CUDA ブロック並列 GPUメモリからデータをGPUへ読み出すときはメモリアドレスとスレッド番号が連続したとき最も高速(プログラムがベクトルが出来ていればOK) 作成するプログラムはベクトルループ内の本体部分のみ作成することに対応。 ループ本体部分をたくさん同時に動かす感じ。 ループ本体部分=> CUDAカーネル ホストから CUDAカーネルをGPUに並列度を指定して投げる感じ GPUで計算するため(した結果)のデータはあらかじめ(後から)転送する必要がある。 CUDA APIによる GPU管理、データ転送、GPU上のメモリ管理

18 4.GPUを用いたソルバーの加速 Nvidia CUDA を用いたプログラミング例

19 Nvidia CUDA を用いたプログラミング例
L次元ベクトルの和の計算例 (L=N×M)    c = a + b //=== host から呼び出される GPU code ==== _global_ void vadd_kernel(float *x, float *y, float *a) { int i = threadIdx.x+blockIdx.x*blockDim.x; x[idx] = a*x[i] + y[idx]; } //==== host 側 code === void main() { …… // GPU上にメモリ確保 cudaMalloc((void**)&a,….); ….. // c = a+b カーネルをGPUへ投げる // thread数/block=N, block数=M で並列実行 vadd_kernel<<<M,N>>>(a,b,c); 高い並列度をうまく利用する必要がある thread 1 thread 2 thread 3 thread 4 thread N block 1 block 2 block 3 block 4 block M block grid thread : 最小の実行単位 (max 512/block) thread block : 同一のmultiprocessor上で 実行されるthreadの集まり (max 65535) grid :thread blockの集まり      並列化されたカーネルの全体 19 19

20 Theread番号、 block番号は多次元化も可能 常に最大のthread長を使えるとは限らない
(カーネルが使用するレジスタの数)×(スレッド長) <= (ハードウェアの持つ実レジスタ数) 8つの実プロセッサが4回同じ命令を異なるデータに対して実行するため、効率のよいスレッド長の最小は32thread 20 20

21 4.GPUを用いたソルバーの加速 Nvidia CUDAへの QCD Mult の実装
Block for 1x1x16x16 lattice Grid for 16^3x16 lattice Thread for single site  1格子点を 1スレッドに割り当てる T,Z,Y,Xの順に連続、T or T-Z をスレッドに残りをブロックに  メモリアクセスがスレッド番号に対して連続になるようにデータをCPUで並べ替えてからGPUで計算  倍精度単精度変換もCPUで行う

22 メモリアクセスを連続に!(Coalesced Access)
Nvidia CUDA を用いたプログラミング例 メモリアクセスを連続に!(Coalesced Access) 4,8,or 16Byte = float, floatx2 or floatx4 格子点 0 格子点 1 格子点 2 thread 0 thread 1 thread 2 格子点 0 格子点 1 格子点 2 thread 0 thread 1 thread 2 22 22

23 4.GPUを用いたソルバーの加速 Nvidia CUDAへの QCD Mult の実装 …… …… …… スレッド長N スレッド1
Float x 4 が基本となる  : 構造体 float4 { float x, y, z, w}; 単位 スピノールデータのレイアウト(complex 12 成分= real 24 成分) (カラー,スピン) (1,1)実 (2,1)実 (3,1)実 (1,2)実 …… (2,2)実 (3,2)実 (1,3)実 (2,3)実 …… (3,3)実 (1,4)実 (2,4)実 (3,4)実 …… スレッド長N スレッド1 スレッド2 スレッドN

24 4.GPUを用いたソルバーの加速 Nvidia CUDAへの QCD Mult の実装 普通のメモリアクセスはキャシュされない
読み書き可 Texture fetching を使うとキャッシュされる 読み出し専用 管理構造体:  texture<float4,1,cudaReadModeElementType> tex_u; 元データ:    gauge_field *u; 元データをTex.fetch. できるように結び付ける:               cudaBindTexture(0,tex_u,u,sizeof(gauge_field)); データをTex.fetch.でロード:               float4 uu = tex1Dfetch(tex_u,isite); 再利用性があるデータには Texture fetching を使う スピノールは1格子点のデータが周りの格子点から8回参照される リンク変数は2回参照される 30-40%ほどさらに高速化

25 4.GPUを用いたソルバーの加速 必要なデータはあらかじめ送っておく。 GPUでの計算は、
リンク場 : U, 右辺ベクトル: b すべての作業ベクトルはGPU上にメモリを確保 GPUでの計算は、 Mult 行列ベクトル積, フェルミオン場ベクトル同士の線形代数計算   の作業単位のGPUカーネルを作成。 単精度ソルバーはホスト側から上記カーネルをアルゴリズムに沿って順番に呼び出す。 残差ノルムはGPUからホストへ適宜転送し収束判定する 収束したら結果ベクトルをGPUからホストへ転送 Allocate working vectors on GPU (cudaMalloc) Send U, b, kappa (cudaMemcpy) ….. while (iter < maxiter) { // qd = D(kappa,U) pd cuMult_kernel<<<Nth,Nbk>>>(kappad,pd,qd,ud); // ctmpd = <rtd|qd> cuProd<<<Nth,Nbk>>>(rtd,qd,ctmpd); // alpd = rho0/ctmpd cuComplexDiv<<<1,1>>>(rho0d,ctmpd,alpd); // ctmpd = - alpd cuNegative<<<1,11>>>(alpd,ctmpd); // xd = xd + alp * pd cuUpdateVec<<<Nth,Nbk>>>(alpd,pd,xd); // rd = rd – alp * qd cuUpdateVec<<<Nth,Nbk>>>(ctmpd,qd,rd); …. }

26 4.GPUを用いたソルバーの加速 結果: 性能の体積依存性(石川-尾崎) CPU: Core2Duo@3GHz
結果: 性能の体積依存性(石川-尾崎) CPU: GPU: GeForce GTX 280 O(a)-improved Wilson-Dirac Red/black site prec’d, Z-T 平面で128格子点を thread並列, 残りをblock並列 Nested-BiCGStab テストゲージ配位:exp(0.001*i*ランダムP) ホストコード Fortran90, BiCGStab, 反復で2回単精度ソルバーを前処理として呼び出す。 単精度ソルバー CUDA/C, BiCGStab, Mult やベクトル計算のカーネルを BiCGStab のアルゴリズムに則った順番で呼び出す。

27 4.GPUを用いたソルバーの加速 結果: 性能の体積依存性(石川-尾崎)
結果: 性能の体積依存性(石川-尾崎) CPU(DP)のみの計算時間と CPU(DP)+GPU(SP)の計算時間

28 4.GPUを用いたソルバーの加速 結果: 性能の体積依存性 軽く10倍以上速くなる
結果: 性能の体積依存性 CPU(DP)のみの計算時間と CPU(DP)+GPU(SP)の計算時間の比 軽く10倍以上速くなる 体積がある程度大きい必要がある 並列度を稼ぐ

29 4.GPUを用いたソルバーの加速 結果: 性能の体積依存性 CPU(DP)のFlops値 キャシュが効かないところで大体 2 GFlops

30 4.GPUを用いたソルバーの加速 結果: 性能の体積依存性 Mult 部分だけは60-100 GFlops
結果: 性能の体積依存性 CPU(DP)+GPU(SP)のGPU(SP)部分のFlops値 単精度ソルバー全体では性能は半分。BLAS2相当の計算が足を引っ張る。 Mult 部分だけは GFlops

31 ランダムではないちゃんとしたゲージ配位での比較
4.GPUを用いたソルバーの加速 ランダムではないちゃんとしたゲージ配位での比較 Time CPU only: 184sec, CPU: 1.9GFlops CPU+GPU: 8.6 sec, CPU: 1.7GFlops GPU: 58GFlops GPU D mult.: 102GFlops この場合 184/8.6 = 21倍の高速化! 10倍以上速くなるのは大きい

32 4.GPUを用いたソルバーの加速 更なる高速化
M.Clark et al. Work March, 2009 らによる更なる改良例 半精度で2割性能上昇 TeX.Fetch. で自動的に  short(16bit) => float(32bit) 変換できる。 そのほかの工夫 ゲージ固定: U_4=I SU(3)行列を実8パラメータで保存 私も試してみました 半精度をリンクに対して試したところ確かに2割Flops値は増加 100GFlops -> 120GFlops しかし反復回数増加のため実時間では遅くなりました。

33 まとめ 4.GPUを用いたソルバーの加速 さらに大きい格子の計算は1台のPC+GPUでは無理。やはり並列化を考える必要がある。
PC1台に1枚の高性能GPUカードをさし場合についてベンチマークテストした(20^4格子まで) 大きめの仕事をGPUにさせないと効率が悪い。10^4以上の格子 単精度加速により10倍以上の性能向上を見た(同一アルゴリズム同士の比較) 10倍速いといっても大きな格子サイズの計算なので絶対時間はそれなりに遅い。とても大きい格子はメモリの制限、CUDA制限で動かない。 コスト? GPUの価格: GeForce GTX285 : 3-5万円 Tesla C1060 : 20-30万円? HostPCの価格:C2D,C2Q一式: 10-20万円? ちゃんとしたサーバーだと: 50万円? PC本体で10倍の性能向上は10台以上並べて並列化? さらに大きい格子の計算は1台のPC+GPUでは無理。やはり並列化を考える必要がある。

34 5.GPUを用いたソルバーの並列化への展望
格子QCD計算の並列計算は基本的に時空を領域分割 Mult y=D x 計算では隣のノードとの通信が発生する。 Mult D 計算を領域分割して計算する方法では通信が2段階必要になる GPUCPUCPUGPU 担当領域が小さいと効率が悪い。 GPUに大きめの仕事を割り当てつつ、通信回数を削減したい GPUで1/Dを解くのではなく 前処理として 1/D の何らかの近似を計算させることに専念 領域分割前処理(Domain-Decomposition) まだGPUを使った並列計算は完成していません。 CPUのみの計算ですが、D.D. を使ってどの部分が 加速できるかを見てみました

35 5.GPUを用いたソルバーの並列化への展望
領域分割前処理 (Domain–decomposition) 問題の微分方程式離散化 問題空間をいくつかに分割(重複も可) 流れ(基本的に反復改良法) 1.初期解と、初期残差 2.残差に対して問題の式を分割された空間で解く 3.解いた結果を使って近似解を構成 4.残差を計算=>2.にもどる 解の更新の順番とか、境界の扱い方とか、部分空間に制限する方法とか、重複部分をどう扱うかとか、、、、 さまざまなバリエーション この方法だけでは収束しないのでこの反復をKrylov部分空間法の前処理として採用する。 Ω1 Ω2 Γ2 Γ1

36 5.GPUを用いたソルバーの並列化への展望
Luescher の導入した領域分割前処理 : Schwartz alternating method (SAP): 2 no-overlapping domain = block Schur complement (Luscher) = Multiplicative Schwarz Method 小領域ソルバーをアクセラレータで加速 並列度、ノードに2色入れる必要がある。Overlap 無しのためノードの担当する格子点数は少ない C.f. SAP=Multiplicative Schwarz Solve in Even domain Solve in Odd domain

37 5.GPUを用いたソルバーの並列化への展望
通信の隠蔽: 領域をオーバーラップさせる、Overlapped D.D. Multiplicative Schwarz (MS) vs Additve Schwarz (AS) MS : generalized Block Gauss-Seidel AS: generalized Block Jacobi             (MS > AS, factor 2) Restricted (Overlapping) Additive Schwarz (RAS) method   Projection on a fermion field [Cai & Sarkis, SIAM J.S.C.21(1999)] Overlapped region Depth 2 Solve in i-th domain D_i x_i =r_i Residual vector field r

38 5.GPUを用いたソルバーの並列化への展望
RAS: 領域を重ね、領域間の依存性を無視、重なり領域は戻さない: 並列度が上がる。(1ノード1領域) 小領域を大きくしてアクセラレータの効率を上げる。 領域の重なりにより前処理性能を上げる。 Return only original region Overlapped region Depth 2 Solution vector field x x = x + \sum_i x_i

39 5.GPUを用いたソルバーの並列化への展望
Restricted (Overlapping) Additive Schwarz (RAS) method Iterate Overlapping improves performance. But task increased. Prolongation is not overlapping (Restricted). This also improves the performance. RASのみでは収束が遅い、MG やDeflation が必要 Each blocks are independent. solved in each block in parallel. GPGPU! C.f. SAP=Multiplicative Schwarz Solve in Even domain Solve in Odd domain

40 5.GPUを用いたソルバーの並列化への展望
Test on small lattice (16^3x32), timing comparison. PC Cluster: 16 nodes. Block size: SAP: 8^4 RAS: (8+2d)^3x(16+2d) Deflate10 small eigenvalues. Best case comparison. 計算時間の9割は単精度計算 ブロックソルバーの時間は薄い部分(GPU加速可能部分) SAP+Defl RAS(d=1)+Defl Fast Slow SAP with Deflation is the best. RAS(d=1) approaches SAP w/o deflation. RAS(d=2,3) reduce iteration count by 1/2-1/3. But the task in each node is rapidly increasing by overlapping reagion. AS は MS にはやはり勝てない? 加速不可能な部分が残る=>全体で2倍速くなれば良いほうか? GPUで比較する予定

41 5.GPUを用いたソルバーの並列化への展望
前処理性能:  加法的シュワルツ法より積的シュワルツ法のほうが良い 積的シュワルツ法での領域オーバーラップ版の性能評価をする必要がある ブロックソルバーをGPUで加速する: 領域のオーバーラップが大きいほど収束が良いが計算量が増えるため時間は増大 全体計算時間のすべてを加速できない 加速できなかった部分のプロパティを調べる

42 6.まとめ GPUを使った格子QCD計算の加速の試みについて紹介した。 1台のPCで並列化を考えなければ非常に有効である。
格子サイズ 10^4 ~ 16^4 程度で効率が良い Schrodinger functional の計算がこのサイズ パラメータ並列でたくさんのPCに計算をばら撒けばよい 大規模計算のためには並列化を考えなければならない。 32^4 以上 GPUの効率を考えるとGPUに与える格子サイズは 10^4~16^4 全体格子サイズが大きくなりすぎないように、GPU間の並列度をあげるのは困難 たくさんのGPUが利用可能なときは、領域分割をオーバーラップさせるのがひとつの手段 シュワルツ法による前処理にGPU利用を提案


Download ppt "格子QCD計算のGPUを用いた加速の可能性について"

Similar presentations


Ads by Google