格子ＱＣＤのための線形計算 （連立一次方程式、固有値問題）について

Presentation on theme: "格子ＱＣＤのための線形計算 （連立一次方程式、固有値問題）について"— Presentation transcript:

1 格子ＱＣＤのための線形計算 （連立一次方程式、固有値問題）について
石川健一(広島大理) アルゴリズムによる計算科学の融合と発展 @CCS,　筑波大学, ２００９年４月２２日・２３日

2 1. 目次 ２．格子ＱＣＤについて ３．ハイブリッドモンテカルロ法 ４．格子化クォークのタイプ ５．クォーク伝播関数ソルバー ６．固有値問題
７．まとめ

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

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

5 ２. 格子ＱＣＤについて 特徴、 自由度がとても大きい 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 典型的な計算時間１～２年（論文、予算、、、、） どのようにしてモンテカルロ積分のための配位を生成するか？

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

7 モンテカルロ重み＝行列式＝補助場による分布 D[U]の逆を含む：非局所的
3.ハイブリッドモンテカルロ(HMC)法 クォーク部分の計算が大変 モンテカルロ重み＝行列式＝補助場による分布 D[U]の逆を含む：非局所的 φで積分すると 真空からクォーク－反クオーク対が対生成、対消滅している効果を表している。真空偏極。

8 Hybrid Monte Carlo（HMC）法
Exp(-H) の分布に従うU,P,φ　を HMC法で生成 Uのサンプル{ U(1),U(2),U(3),…..,U(N)} 物理量の期待値(物理量は U のみの関数 O[U])

9 HMC法 3.ハイブリッドモンテカルロ(HMC)法 に対して、マルコフチェーンモンテカルロ法の一種ＨＭＣ法を適用
　　　に対して、マルコフチェーンモンテカルロ法の一種ＨＭＣ法を適用 分子動力学(Molecular Dynamics) (MD) メトロポリス法(Metropolis) 分子動力学ではエネルギーは差分化のため保存しない 分布がほしい分布　　　　　　　　からずれる で補正する 分子動力学法についての発展は今回は省略

10 Hybrid Monte Carlo（HMC）法
計算が困難な点 分子動力学法の部分で の大規模連立１次方程式を何回も解く必要がある。 幸いＤ［Ｕ］の要素はほとんどゼロ（大規模疎行列）であるので、反復法 を用いて連立方程式を解く。 しかし特にクォークの質量が小さいと、D[U]がゼロに近い固有値を持つため連立方程式を解くのが難しくなる。 クォークの質量が小さいと分子動力学の力が大きくなる。分子動力学の時間刻みを小さくする必要がある。

11 問題をまとめると、 3.ハイブリッドモンテカルロ(HMC)法 クォークのタイプによって D[U] がいろいろある
モンテカルロ法による経路積分の評価で困難が生じるのはクォークの寄与を取り入れる部分。 Uの生成には、 D[U] x = b の連立一次方程式を大量に解く必要がある これまでさまざまな前処理方法が開発されてきた 最近の発展について述べたい。 クォークソルバー

12 ４.格子化クォークのタイプ 離散化には任意性がある 良い性質を持つものほど計算コストがかかる クォークの持つ対称性
ゲージ対称性　：　　格子上でも死守すべき対称性 ポアンカレ対称性: 一部破れるけどなるべく保持したい対称性、ユークリッド化後：離散並進、離散回転、反転対称性 カイラル対称性：　　　 　　　　　　　 格子上で厳密に実現することは不可能(ニールセン-二宮の定理) カイラル対称性の保持の仕方による分類 Kogut-Susskind 型：４つの同じ質量をもつクォーク（４フレーバー）を同時に扱う。U(1)xU(1)カイラル対称性を保持 Wilson-Dirac型：１つのクォークから扱える。カイラル対称性は無い Overlap型、Domain-Wall 型：１つのクォークから扱える。カイラル対称性は変形された対称性を保持。計算コストがとても高い(Wilson型に比べて10-100倍)

13 ４.格子化クォークのタイプ Kogut-Susskind(KS) 型：４つの同じ質量をもつクォーク（４フレーバー）を同時に扱う。U(1)xU(1)カイラル対称性を保持 最も簡単なD[U]の構造を持つ。現実世界のクォークは６種類で質量もばらばら。軽いのはアップとダウンの２種類のクォーク　＝　４つ同時では扱えない。スキップします。

14 ４.格子化クォークのタイプ Wilson-Dirac型：１つのクォークから扱える。カイラル対称性は無い。 Ｄ[U]は一階の差分演算子
クォーク質量　m ∝ 1/κ カラー(3),スピン(4)の自由度が各格子点にある ブロックサイズ１２のブロックの帯が並んだ疎行列 M の固有値分布 U=1,4^4 Witzel,Takeda,Wolff, hep-lat[arXiv: ]

15 ４.格子化クォークのタイプ Overlap型、Domain-Wall 型：１つのクォークから扱える。カイラル対称性は変形された対称性を保持。計算コストがとても高い(Wilson型に比べて10-100倍) 質量　m Wilson-Diracクォークの行列 DW のサイン関数を含む Dov は密な行列 対称性：　質量m=0 でGinparg-Wilson関係を満たす はユニタリー行列 Taken from Talk slid by R.Edwards

16 ４.格子化クォークのタイプ まとめ、 Wilson-Dirac型 Dx = b は反復法で解く 係数行列 D の型
１階差分型：さまざまな前処理が適用されてきた Red-brack(2色), ILU,SSOR, ２色ブロック … etc. Dx = b を直接解く: Non-Hermitian sovler, BiCGStab, GMRES などが使われてきた (D†D)x=D†b を解く: positive Hermitian solver, 　　　　　　　　　　CG が使われる。条件数悪し。 γ5Dx=γ5D を解く: indefinite Hermitian solver, MINRESなど、、 前処理つき　BiCGStab, GMRES 系統が良く使われる

17 ４.格子化クォークのタイプ Overlap 型 Dをベクトルにかけるときに sign 関数の計算が必要
Dx=b 自体は Hermite化後 CG などで解く Sign関数の計算 Arnoldi / Lanczos 分解による方法　　 Pade / Rational ：部分分数展開近似による方法 部分分数/連分数展開＋補助場の導入による方法　　　　　　　　　　　　　　　　　　　　　　　　　　=> ５次元化　ソルバーの入れ子を解消 HW や DWの固有値範囲の情報を用い sign関数の近似の精度をコントロールする必要がある。 Overlap 型についてはこれ以上次官の都合上言及しない 日本では　JLQCD　collaboration でくわしく研究されている KEK　松古、金児、橋本…..

18 ５.クォーク伝播関数ソルバー Ｄｘ＝ｂを速く解きたい 格子サイズ: 32^3x64 (現状 O(10)TFlops)
　　　　　　　　　=> 256倍：　128^3x256 (O(2)PFlops) 次元: １２ｘ[格子サイズ]= 2500万　＝＞　64億 (1)前処理：　Wilson型に関しては構造が簡単（Red-Black, ILD, Domain-Decomposition etc.）, Overlap型については難しい⇒5D化で構造が簡単に [省略] (2)固有値固有ベクトルの情報：　　Deflation/Multigrid (3)単精度加速： (4)アクセラレータ：

19 (2) Deflation technique
5.クォーク伝播関数ソルバー (2) Deflation technique LQCD では何回も Dx=b を解く 右辺ベクトルを取り替えながら、係数行列は一定。または、係数行列が少しづつ変化しながら連鎖的に解く必要がある Quark propagator Solver in HMC trajectory ブロック化ソルバー：係数行列固定、右辺ベクトル多数 Deflation technique: 係数行列の固有ベクトルの情報により前処理 Deflation remove/suppress small eigenspace of D.

20 (2) Deflation technique （Solve Ａｘ＝ｂ ・・・(1) ）
5.クォーク伝播関数ソルバー (2) Deflation technique （Solve Ａｘ＝ｂ　・・・(1) ） 行列A ：　小さい固有値に属する　p-次元部分空間を持つとする c,u : この部分空間に属する以下の条件を満たすベクトルとする c から射影演算子を定義 Pを使って問題を前処理 式(1)の解 x は式(2) の解 y を用いて以下のように書かれる “Cp”はどうやって作る？ 反復には　P の掛け算が伴う(overhead)

21 5.クォーク伝播関数ソルバー (2) Deflation technique (cont’d) Many works by
正確な不変部分空間の計算はとてもコストが高い 　 　　低モード密度 ∝ O(V), モード計算コスト O(V) 　　　　　　　　=> 全体で　O(V^2) 　の問題 (a) Dx=b を解きながら、同時に不変部分空間を求め使いまわす GCRO-DR: Parks & Sturler, used by PACS-CS collab. GMRES-DR,GMRES-E..: Wilcox, Morgan & Abdel-Rehim (b)低モード固有ベクトルは滑らかに局在。ブロック分割可能？ Luscher’s Domain decomposed subspace blocking with local coherency. (Used also in HMC algorithm) [Luescher, JHEP07(2007),hep-lat/ ; A.Stathopoulos, K.Orginos, hep-lat/ ; W.Wilcox, PoS(LATTICE2007),hep-lat/ ; A.Abdel-Rehim,R.B.Morgan,W.Wilcox,PoS(LATTICE2007); R.B.Morgan,W.Wilcox,math-ph/ ,math-ph/ ; M.L.Parks, E.De Sturler et al, SIAM J. on Sci.Comp. 28(2006)1651 LATTICE2008: Poster by Abdel-Rehim, Talk by Wilcox] More details see

22 (a)固有空間の計算と連立一次方程式を同時に解く.
5.クォーク伝播関数ソルバー (a)固有空間の計算と連立一次方程式を同時に解く. Very effective for few Near zero modes / negative eigen modes case. Near zero modes case First equation or few equations are solved with GMRES-DR. Once the subspace converged, change solver with GMRES-proj, or Deflated solver. Normal GMRES stagnates [dot-dot-dashed line] Solver with Deflation/Projection converges. [other lines] [Wilcox, LAT2007]

23 (b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する.
5.クォーク伝播関数ソルバー (b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する. 基底ベクトルの構成方法が重要 Low mode ベクトルはブロックベクトルで近似できる [Luscher, JHEP07(2007)081]

24 (b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する.
5.クォーク伝播関数ソルバー [Luscher, JHEP07(2007)081] (b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する. B　:　ブロック化された部分空間でのD[U]の表現（＝サブグリッド上に射影された） 基底ベクトル（C）の構成 ランダムベクトルに1/D を数回かけて低モードを増幅 そのベクトルをサブブロックに分割 分割されたベクトルを正規直交化 Deflation の射影演算子　(B^-1)を含む.

25 (2) MultiGrid Solver 5.クォーク伝播関数ソルバー
MultiGrid solver also removes critical slowing down. Choice of subspace basis is important. (Prolongator) 射影のための基底ベクトルの選び方が重要（Luscher の基底と同じ） QCD, 16^3x32 [Brannick,Brower,Clark,Osborn,Rebbi, PRL100(2008)]

26 (3) 単精度加速 単精度の利点: 有効バンド幅、キャッシュサイズ、（レジスターの数）などが倍になる＝効率（実性能/理論性能）が高い
5.クォーク伝播関数ソルバー (3) 単精度加速 単精度の利点: 有効バンド幅、キャッシュサイズ、（レジスターの数）などが倍になる＝効率（実性能/理論性能）が高い Intel 64/AMD 64; Single prec. > Double prec. Cell PS3/GPGPU; Single >> Double. For Wilson kernel : 3Byte/Flop => 1.5 Byte/Flop 連立一次方程式の解法を加速するのに使えないか？ 反復改良法/Richardson反復のテクニックで可能 一般の反復法にも組み入れることが出来る（可変前処理） 前処理として単精度ソルバーを使う 残差は倍精度で正しくなるように組む FGMRES, GCR, CG, BiCGStab….. (with flexible prec.) [Buttari,Dongarra,Kurzak,Luszczek,Tomov, AMC Trans.Math.Soft]

27 Nested-BiCGStab (櫻井-多田野)
5.クォーク伝播関数ソルバー (3) 単精度加速 Nested-BiCGStab (櫻井-多田野) Outer solver : BiCGStab (D.P.) Inner solver: BiCGStab (S.P.)= Outer solverに対する可変前処理 Intel 64, SSE3 使用+リンク再構築で全て倍精度の時より速度が倍に（演算量は増えているが） [PACS-CS collab.]

28 5.クォーク伝播関数ソルバー (４) アクセラレーター 単精度加速と組み合わせるとさらに有効なもの
GPGPU, CellB.E. は単精度計算が圧倒的に速い GRAPE => QCD? NVIDIA GT200 (Tesla 10series) 240 SP (SP cores), 30 DP cores ～1,000(or 600)Glops(SP), ～90GFlops(DP) We expect > 60 GFlops(SP) for QCD kernel. (assuming 10% efficiency)

29 5.クォーク伝播関数ソルバー (４) アクセラレーター NVIDIA CUDAを使った例（石川-尾崎） CPU: Core2Duo@3GHz
GPU: GeForce GTX 280 O(a)-improved Wilson-Dirac Red/black site prec’d Nested-BiCGStab 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倍以上速くなるのは大きい

30 5.クォーク伝播関数ソルバー (４) アクセラレーター
NVIDIA CUDAを使った例（Clark et al. Work March, 2009） Nvidia GPU でさらに加速する方法 単精度＝＞半精度（16bit）化してさらにバンド幅を稼ぐ 符号付16bit整数をノーコストで[-1.0,1.0] 単精度へ自動変換できる機能を利用（Texture Fetching） 有効バンド幅と、有効　Texture Cache サイズが倍、レジスタは32bit 浮動小数点のものしかないので演算は単精度 メモリとのやり取りはTexture Fetching (ロード)と整数への丸め(ストア)

31 5.クォーク伝播関数ソルバー (４) アクセラレーター NVIDIA CUDAを使った例
[M.Clark et al. Work March, 2009] 半精度(16bit int)でさらに２割性能上昇 そのほかの工夫でさらに加速 ゲージ固定 ＳＵ（３）行列を実８パラメータで保存 GPGPU:演算が圧倒的に速いので、バンド幅、キャッシュの有効利用のほうが重要

32 [H.Baier at al. ,Lat2008, arXive:0801.1559
5.クォーク伝播関数ソルバー [Spary,Hill,Trew hep-lat/ ; S.Motoki & A. Nakamura Lat2007; F.Belletti et al. LAT2007] (４) アクセラレーター Cell B.E.での計算 大規模 Cell B.E.クラスタ： QPACE (QCD Parallel computing on Cell B.E.) project (EU) 2048 PowerXCell8i, + カスタムネットワーク3Dトーラス(FPGA) 2009春完成？ ＧＰＵクラスター: National Taiwa Univ. 16nodes+16TeslaS1070: 64TFlops KEK, 4 nodes + 4 Tesla S1070: 16TFlops ….? 多体系の計算に対する応用がGPGPUでも先行　（GRAPE） AMD/ATI での計算 [V.Demchik and A. Strelchenko arXiv: ] 並列化はどうする？ 領域分割前処理+Deflation/Multigrid ? [H.Baier at al. ,Lat2008, arXive:

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

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

35 Multiplicative Schwarz (MS) vs Additve Schwarz (AS)
5.クォーク伝播関数ソルバー 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 　 [Cai & Sarkis, SIAM J.S.C.21(1999)] Projection on a fermion field Overlapped region Depth 2 Solve in i-th domain D_i x_i =r_i Residual vector field r

36 RAS: 領域を重ね、領域間の依存性を無視、重なり領域は戻さない：
5.クォーク伝播関数ソルバー RAS: 領域を重ね、領域間の依存性を無視、重なり領域は戻さない： 並列度が上がる。（１ノード１領域） 小領域を大きくしてアクセラレータの効率を上げる。 領域の重なりにより前処理性能を上げる。 Return only original region Overlapped region Depth 2 Solution vector field x x = x + \sum_i x_i

37 Restricted (Overlapping) Additive Schwarz (RAS) method
5.クォーク伝播関数ソルバー 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

38 Test on small lattice (16^3x32), timing comparison.
5.クォーク伝播関数ソルバー 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 にはやはり勝てない？ GPUで比較する予定

39 6. 固有値問題 Ｄ［Ｕ］の固有値固有ベクトルが分かるといろいろ便利なことがある。 Low-mode averaging
物理量にはクォークD[U]の小さい固有値固有ベクトルが重要。 並進対称性を利用した統計の向上 Determinant splitting 行列式の評価で、小さい固有値は正確に計算し残りはモンテカルロで評価する Epsilon regime D[U]のゼロ付近の固有値分布とカイラル対称性の破れの関係 Deflation for quark solver 大量のDx=b を解く時に Dの固有空間が分かるとそれを使って前処理して、高速化 Matrix function Overlap 型クォークの sign関数、奇数フレーバーのための sqrt(D),　行列式をTrLog[D]で評価、前処理のための (D†D)(1/n) などなど、、、

40 固有値固有ベクトルを求めるのは Dx=b を一回だけ解くのに比べてはるかに困難 Lattice QCD で使われている手法
6. 固有値問題 固有値固有ベクトルを求めるのは Dx=b を一回だけ解くのに比べてはるかに困難 Lattice QCD で使われている手法 H =γ5D : Hermite 行列 =>　Lanczos アルゴリズムとその系統 D : non-Hermitian 行列 => Arnoldi アルゴリズムとその系統 Lanczos系 Thick-restareted Lanczos [Wu&Simon] PACC-CS collab. で使用中、他のグループでも使われているかも Dに対してはほぼ ARPACK (Implicitly Restarted Arnoldi Method)が使われている。 大体のLQCD論文で ARPACKが言及されている そのほかのアルゴリズムは？ Jacobi-Davidoson 法は？ [Lehoucq, Sorensen, Yang, Maschhoff]

41 6. 固有値問題 Arnoldi/Lanczos系での困難
原点に近い固有値固有ベクトルを求めたいとき、逆反復にする(+shiftもするかも) Arnoldi(Lanczos)反復で問題を小さな問題に変換 1. 適当な初期ベクトルから直交基底ベクトルを構成していく（m反復後） 2. Hm の固有値固有ベクトルから 1/D の大きいほうの固有値固有ベクトルを抽出＝Dの小さい固有値固有ベクトルを抽出 3. ほしいベクトルを何本か残して式の形が変化しないようにユニタリー変換して縮小してリスタート(Implicit restart) 4. ほしいベクトルが収束するまでリスタートを繰り返す 1/Dの計算は反復法で行う 計算コスト: 大体 (m x [リスタート回数]) 回、反復法で Dx=b を解く 得られる固有ベクトルの精度は 1/D の計算精度で決まってしまう。 Dx=b に対する反復法も良い精度で解く必要がある。　

42 6. 固有値問題 Arnoldi/Lanczos系での困難 固有ベクトルの情報を何かの計算の加速に使おうと思うと、
連立一次方程式 Dx=b を精度よく何回も解かないといけないので本末転倒しているような気がする。 私がいろいろ調べた結果 Jacobi-Davidson法はその困難がない？ Deflation付のソルバーを使わないといけない様である。 固有値固有ベクトルペア(Ritz pair)を１個づつ求めていく様である。 ターゲットの Ritz pair の選び方とか、うまくリスタートする実装の仕方がまだ良く理解できていないのでテストしていません。 RESidual Arnoldi Method (RESAM) 部分空間反復法とArnoldi 法の合成みたいな方法 Jacobi-Davidson法にも似ているが Deflation 付ソルバがいらない D論なのでくわしく書いてある 売りは、Dx=b は低精度で解くので十分なので IRAMよりも速い [C.R.Lee & G.W.Stewart, Tech. Report TR ; C.R.Lee, Ph.D. Thesis] IRAM(KSRAM) と　RESAM を　QCD で比較してみました。

43 6. 固有値問題 RESRAM Step. 4 は精度を求めないので単精度でも OK
1. 初期部分空間を Arnoldi法で作っておく 2. その部分空間から Ritz pair を用いてターゲットとなる近似固有値固有ベクトルを求める 3. その近似固有値固有ベクトルの残差を計算する。収束していたら Deflation に加える、次のターゲットを探すべく 2. へ戻る 4. 残差に対して　1/D の近似をかける(Ritz値へshift しても良い) 5. それを以前の部分空間に対して正規直交化して部分空間に付け加えて部分空間を拡張する 6.部分空間の次元がある最大値に達したら Schur型を経由して次元を縮める。その際収束している部分空間が変化しないようにうまくする。 Step. 4 は精度を求めないので単精度でも OK 最終結果に倍精度を求めるなら、それ以外の計算、特に固有方程式の残差の計算は倍精度で行う必要がある。

44 6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
IRAM ARPACK をそのままコンパイル KSRAM IRAMを自分で実装できなかったのでIRAMとほぼ同等な Krylov-Schur Restarted Arnoldi Method [TRLANのnon-Hermitan版] を実装, RESAM 以下を変化させて速度測定 変化させるパラメータ： 線形ソルバー(Dx=b)の精度、 求める固有値固有ベクトルの精度、 線形ソルバーは共通の GCRO-DR法 最大部分空間次元は 40 リスタート時に残す部分空間次元は min((すでに求まった部分空間次元＋20),40) 環境：PC一台 : CPU Core2Duo 原点付近の固有値を２０個求める

45 6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
原点付近の固有値を２０個求める 比較テスト　小問題　[4^4 格子 beta=5.0 クエンチ Matrix market より] 欲しい固有値固有ベクトルの精度よりも線形ソルバーの精度を落としたときの固有値固有ベクトルの残差はどうなるか？ IRAM:線形ソルバーの精度より良くなることはない RESAM:　常に必要な精度が得られる

46 6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
原点付近の固有値を２０個求める 比較テスト　小問題　[4^4 格子 beta=5.0 クエンチ Matrix market より] 倍精度で解を求める場合の計算時間の比較 IRAM,KSRAM:線形ソルバーは倍精度で説き続ける RESAM:線形ソルバーは1/10の精度で求めればよい。 全体の計算時間のほとんどは線形ソルバーの時間に費やされる。 固有値ソルバーの反復回数は IRAM/KSRAM << RESAM 積算された時間は RESRAM << IRAM/KSRAM

47 6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
原点付近の固有値を２０個求める 比較テスト　小問題　[4^4 格子 beta=5.0 クエンチ Matrix market より] RESAMで倍精度で解を求める場合の固有ペア残差のヒストリ １本ずつ原点に近いものから求めていっている。 線形ソルバーの停止条件は残差が1/10 を切ったとき。 線形ソルバーの停止条件をきつくすると、REAAMの反復回数は減っていくが、線形ソルバーの時間の増加を補うほどではないので帰って全体が遅くなる。 Best Time は tol_solver = 1/10の時

48 6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
原点付近の固有値を２０個求める 比較テスト　小問題　[4^4 格子 beta=5.0 クエンチ Matrix market より] 加速率の比較(Best同士の比較) IRAM,KSRAM: (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=1/10 高精度(誤差10^-4以下)の固有ペアを求めるときRESRAM は　IRAM/KSRAM より速い。 倍精度ぎりぎりのときは特に倍速い 本格的に大きな格子ではどうか？（並列化必要）

49 6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
原点付近の固有値を２０個求める 比較テスト　小問題　[4^4 格子 beta=5.0 クエンチ Matrix market より] 加速率の比較(Best同士の比較) IRAM,KSRAM: (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=1/10 高精度(誤差10^-4以下)の固有ペアを求めるときRESRAM は　IRAM/KSRAM より速い。 倍精度ぎりぎりのときは特に倍速い 本格的に大きな格子ではどうか？（並列化必要）

50 6. 固有値問題 比較テスト 中問題 [16^3x32 格子] 原点付近の固有値を２０個求める 線形ソルバー
KSRAM: Nested BiCGStab(単精度部分: BiCGStab/GCRO-DR) RESRAM: 単精度 BiCGStab/GCRO-DR 計算時間の比較(Best同士の比較) (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=0.5 すべての必要精度で、RESRAMが速いようだ。

51 6. 固有値問題 比較テスト 中問題 [16^3x32 格子] 原点付近の固有値を２０個求める 加速率の比較(Best同士の比較)
KSRAM: (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=0.5 すべての必要精度で、RESRAMが速いようだ。大体倍ぐらい速い。 QCD(Wilson-Dirac型)では確かにRESRAMは速そうだ。 Hermit 版はどうか？ Jacobi-Davidson法も同様に速いのかも？ 線形ソルバーを低精度に出来る ＝＞アクセラレータでさらに加速？

52 7．まとめ 格子QCD: クォーク行列の逆、行列式などの評価の高速化が重要 Wilson-Dirac型クォーク： さまざまな改良：
係数行列の形は簡単＝プログラム上の最適化では限界がすぐ来る＝＞アルゴリズムの改良は重要 さまざまな改良： 前処理：領域分割、Deflation、MultiGrid、 精度を落として加速？　並列度を上げる？　アクセラレータ？ 固有ペアの計算 これらの改良の組み合わせ 紹介し切れなかったこと 行列関数(Overlap型クォークなど)の評価法 分子動力学の改良 マルコフチェーンモンテカルロ法の改良？