格子ＱＣＤ計算のＧＰＵを用いた加速の可能性について 石川健一(広島大理) @CCSコロキウム,　筑波大学, ２００９年６月１９日

1. 目次 ２．格子ＱＣＤについて ３．ハイブリッドモンテカルロ法 ４．ＧＰＵを用いたソルバーの加速 単精度前処理 ベンチマーク結果 ５．ＧＰＵを用いたソルバーの並列化への展望 ６．まとめ

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

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

２. 格子ＱＣＤについて 分配関数（ユークリッド化された経路積分） 物理量の期待値 統計力学で用いられてきた方法 作用 有効作用 有効作用を重みとする多重積分　＝＞　モンテカルロ積分 統計力学で用いられてきた方法

２. 格子ＱＣＤについて 特徴、 自由度がとても大きい 16^4 ～ 32^4 格子 系統誤差 自由度がとても大きい　　16^4 ～　32^4　格子 グルーオン：８ｘ４ｘ（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 どのようにしてモンテカルロ積分のための配位を生成するか？

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

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

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

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

ＧＰＵ例 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

先行研究(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. Latice2008 045 [hep-lat/0810.5365]. M.Clark et al. Work shop@CCS, 10-12 March, 2009. AMD　製　GPU　例（つい最近） V.Demchik and A.Strelchenko, hep-lat/0903.3053v2

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

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

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

Mixed precision / Inner-Outer solver ４．ＧＰＵを用いたソルバーの加速 Mixed precision / Inner-Outer solver Flexible Preconditioner 可能な反復法で計算 Richardson 反復（Iterative refinement）法を基に作る より複雑な反復法の中に組み入れることも可能 ＧＣＲ法、ＦＧＭＲＥＳ法、ＣＧ法、BiCGStab法 今回の計算では　Nested-BiCGStab (櫻井-多田野)を用いた

４．ＧＰＵを用いたソルバーの加速 Nvidia CUDA を用いたプログラミング 基本は超並列プログラム 基礎言語はＣ言語 情報源：[http://www.nvidia.co.jp/object/cuda_home_jp.html] 基本は超並列プログラム 基礎言語はＣ言語 ベクトル並列計算機向プログラムを移植するのが簡単 ベクトルループ => CUDA スレッド並列 ＭＰＩ並列 => CUDA ブロック並列 ＧＰＵメモリからデータをＧＰＵへ読み出すときはメモリアドレスとスレッド番号が連続したとき最も高速（プログラムがベクトルが出来ていればＯＫ） 作成するプログラムはベクトルループ内の本体部分のみ作成することに対応。 ループ本体部分をたくさん同時に動かす感じ。 ループ本体部分=>　ＣＵＤＡカーネル ホストから CUDAカーネルをＧＰＵに並列度を指定して投げる感じ ＧＰＵで計算するため（した結果）のデータはあらかじめ（後から）転送する必要がある。 ＣＵＤＡ　ＡＰＩによる　GPU管理、データ転送、ＧＰＵ上のメモリ管理

４．ＧＰＵを用いたソルバーの加速 Nvidia CUDA を用いたプログラミング例

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

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

４．ＧＰＵを用いたソルバーの加速 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　をスレッドに残りをブロックに 　メモリアクセスがスレッド番号に対して連続になるようにデータをＣＰＵで並べ替えてからＧＰＵで計算 　倍精度単精度変換もＣＰＵで行う

メモリアクセスを連続に！(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

４．ＧＰＵを用いたソルバーの加速 Nvidia CUDAへの QCD Mult の実装 …… …… …… スレッド長Ｎ スレッド１ Float x 4 が基本となる　　： 構造体 float4 { float x, y, z, w}; 単位 スピノールデータのレイアウト（complex 12 成分= real 24 成分） （カラー，スピン） （１，１）実 （２，１）実 （３，１）実 （１，２）実 …… （２，２）実 （３，２）実 （１，３）実 （２，３）実 …… （３，３）実 （１，４）実 （２，４）実 （３，４）実 …… スレッド長Ｎ スレッド１ スレッド２ スレッドＮ

４．ＧＰＵを用いたソルバーの加速 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％ほどさらに高速化

４．ＧＰＵを用いたソルバーの加速 必要なデータはあらかじめ送っておく。 ＧＰＵでの計算は、 リンク場　: U, 右辺ベクトル: b すべての作業ベクトルはＧＰＵ上にメモリを確保 ＧＰＵでの計算は、 Mult　行列ベクトル積, フェルミオン場ベクトル同士の線形代数計算 　　の作業単位のＧＰＵカーネルを作成。 単精度ソルバーはホスト側から上記カーネルをアルゴリズムに沿って順番に呼び出す。 残差ノルムはＧＰＵからホストへ適宜転送し収束判定する 収束したら結果ベクトルをＧＰＵからホストへ転送 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); …. }

４．ＧＰＵを用いたソルバーの加速 結果： 性能の体積依存性（石川-尾崎） CPU: Core2Duo@3GHz 結果：　性能の体積依存性（石川-尾崎） CPU: Core2Duo@3GHz 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　のアルゴリズムに則った順番で呼び出す。

４．ＧＰＵを用いたソルバーの加速 結果： 性能の体積依存性（石川-尾崎） 結果：　性能の体積依存性（石川-尾崎） CPU(DP)のみの計算時間と CPU(DP)+GPU(SP)の計算時間

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

４．ＧＰＵを用いたソルバーの加速 結果：　性能の体積依存性 CPU(DP)のFlops値 キャシュが効かないところで大体 2 GFlops

４．ＧＰＵを用いたソルバーの加速 結果： 性能の体積依存性 Mult 部分だけは60-100 GFlops 結果：　性能の体積依存性 CPU(DP)+GPU(SP)のGPU(SP)部分のFlops値 単精度ソルバー全体では性能は半分。ＢＬＡＳ２相当の計算が足を引っ張る。 Mult 部分だけは60-100 GFlops

ランダムではないちゃんとしたゲージ配位での比較 ４．ＧＰＵを用いたソルバーの加速 ランダムではないちゃんとしたゲージ配位での比較 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倍以上速くなるのは大きい

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

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

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

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

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

５．ＧＰＵを用いたソルバーの並列化への展望 通信の隠蔽：　領域をオーバーラップさせる、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

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

５．ＧＰＵを用いたソルバーの並列化への展望 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

５．ＧＰＵを用いたソルバーの並列化への展望 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割は単精度計算 ブロックソルバーの時間は薄い部分（ＧＰＵ加速可能部分） 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で比較する予定

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

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