格子量子色力学におけるソルバーについて
Yusuke Namekawa (KEK)

Contents
1Introduction

22Solvers in lattice QCD

53Benchmark results

94

Additional hot topics with multiple right hand side

125

Summary

20

Yusuke Namekawa(KEK)– 1 / 21 –HPC-Phys meeting1Introduction

All materials are made from quarks and leptons
cf. Kanamori-san’s talk at the 2nd HPC-Phys meeting• Theory of the strong interaction among quarks is called”Quantum ChromoDynamics(QCD)”

http://higgstan.com/ ← the designer got PhD on particle physics experiment
Yusuke Namekawa(KEK)– 2 / 21 –HPC-Phys meeting[Quantum ChromoDynamics(QCD)]• Theory(Lagrangian) is known, but is difficult to be solved analytically
1 2
LQCD = q̄(iD/ − m)q − G4
♦ One of Millennium Problems http://www.claymath.org/millennium-problems→ You will win one million USD, if you solve this problem♦ (cf. one of Millennium Problems on Poincare conjecture has
already been solved)• Numerical simulation of QCD on discretized spacetime (lattice QCD)
is possible♦Ax = b plays the central role
→ Solver is important
cf. Kanamori-san’s and Ishikawa-san’s talks at the 2nd, 3rdHPC-Phys meetings

Yusuke Namekawa(KEK)– 3 / 21 –http://www-het.ph.tsukuba.ac.jp
HPC-Phys meeting[Concrete form of A for Ax = b in lattice QCD]
• Concrete form of A depends on the fermion formulation♦ One choice is Wilson-type fermion(9-point stencil in 4-dimension,
complex non-symmetric large sparse matrix)
cf. Kanamori-san’s and Ishikawa-san’s talks at the 2nd, 3rd HPC-Phys meetings• Condition number K(A) becomes larger for smaller quark mass mquark
cf. Ishikawa-san’s talk at the 3rd HPC-Phys meeting♦ K(A(mud )) = O(2700), K(A(ms )) = O(100), ms /mud ∼ 27
A(x, y)=:mquark=K(A)∝4 noX†δx,y − κ(1 − γµ )Uµ (x)δx+µ,y + (1 + γµ )Uµ (x − µ)δx−µ,yµ=1
complex n × n non-symmetric matrix,1 1− (const)2
κ1n ∼ 1010

for a typical lattice QCD

mquarkYusuke Namekawa(KEK)– 4 / 21 –HPC-Phys meeting2Solvers in lattice QCD

Major solvers in lattice QCD are tabulated• There are many solver algorithms for lattice QCD→ Only solvers in the table are explained• There are many open sources for lattice QCD→ Only open sources in the table are explained• (Preconditioners are not covered in this talk)SolverCG Hestenes,Stiefel(1952)
BiCGStab van der Vorst(1992)
BiCGStab(L) Sleijpen,Fokkema(1993)
BiCGStab(DS-L) Miyauchi et al.(2001)
BiCGStab(IDS-L) Itoh,Namekawa(2003)
GMRES(m) Saad,Schultz(1986)
MultiGrid A.Brandt(1977)
Yusuke Namekawa(KEK)Open source
Bridge++Bridge++, CCSQCDSolverBenchBridge++
Bridge++Bridge++Bridge++DDalphaAMG– 5 / 21 –HPC-Phys meeting

[Lattice QCD code Bridge++ (our open source code)]• Bridge++ is a code set for numerical simulations of lattice gauge
theories including QCD→ Ver.1.5.1 has been released in Aug 2019
• Major solvers(BiCGStab series,CG,GMRES(m)) are covered• Project members:
Y.Akahoshi (YITP), S.Aoki (YITP), T.Aoyama (KEK), I.Kanamori (R-CCS), K.Kanaya (Tsukuba),
H.Matsufuru (KEK), Y.Namekawa (KEK), H.Nemura (RCNP), Y.Taniguchi (Tsukuba)♦ I have been the chairperson since 2016

Yusuke Namekawa(KEK)– 6 / 21 –HPC-Phys meeting

[CCS QCD SolverBench]• CCS QCD SolverBench is a benchmark BiCGStab program of QCD
developed by another CCS(Univ of Tsukuba)→ Ver.0.999(rev.248) has been released in Sep 2017
• BiCGStab with even-odd preconditioning is employed• Project members:
K-I.Ishikawa (Hiroshima), Y.Kuramashi (Tsukuba), A.Ukawa (Tsukuba), T.Boku (Tsukuba)

https://www.ccs.tsukuba.ac.jp/qcd/

Yusuke Namekawa(KEK)– 7 / 21 –HPC-Phys meeting[DDalphaAMG]• DDalphaAMG is a multigrid solver program in lattice QCD→ Ver.1701 has been released in Jan 2017
→ Ported to K-computer in Apr 2018 Ishikawa,Kanamori(2018)• Adaptive Algebraic MultiGrid(αAMG) algorithm with Domain
Decomposed(DD) smoother is employed• Project members:
M.Rottmann, A.Strebel, S.Heybrock, S.Bacchio, B.Leder, I.Kanamori

https://github.com/DDalphaAMG
https://github.com/i-kanamori/DDalphaAMG/tree/K

Yusuke Namekawa(KEK)– 8 / 21 –HPC-Phys meeting3Benchmark results[CG vs BiCGStab series, GMRES(m) by Bridge++]• For mud (up-down quark mass), which requires a huge Krylov space,
BiCGStab series gain 30-40%, while GMRES(m=2–16) shows no gain
• For ms (strange quark mass), which requires not so large Krylov space,
BiCGStab series and GMRES(m) gain a factor of 3
♦ Prescription is added to BiCGStab for better stability
Sleijpen and van der Vorst(1995)

0.8

1.0
163 × 32Wilson, mud

0.6
0.4
0.2
0.0

BiCGStab
BiCGStab(L=2)
BiCGStab(DS-L)
BiCGStab(IDS-L)
GMRES(m=2)

Yusuke Namekawa(KEK)

Nmult(solver) / Nmult(CGNR)

Nmult(solver) / Nmult(CGNR)

1.0

0.8
0.6

163 × 32
Wilson, mstrange

BiCGStab
BiCGStab(L=2)
BiCGStab(DS-L)
BiCGStab(IDS-L)
GMRES(m=2)

0.4
0.2
0.0

– 9 / 21 –HPC-Phys meeting[CG vs MG(MultiGrid)] Babich et al.(2010)• For mud (up-down quark mass), which requires a huge Krylov space,multigrid gains a factor of 3
• For ms (strange quark mass), which requires not so large Krylov space,multigrid has no gain due to its overhead♦ Memory cost of multigrid is larger than that of CG by a factor of 4–5
physcriticalcritical♦ NB. m∝ (mbarequark − mquark ) with mquark = −0.4175
quarkYusuke Namekawa(KEK)– 10 / 21 –HPC-Phys meeting[Nested BiCGStab with precond(SAP + SSOR) vs multigrid] Ishikawa,Kanamori(2018)
Similar results are obtained on K-computer• For mud (up-down quark mass), which requires a huge Krylov space,multigrid gains a factor of 2 over the baseline BiCGStab• For ms (strange quark mass), which requires not so large Krylov space,multigrid has no gain due to its overhead♦ The best solver depends on the target systemelapsed time [sec.]

elapsed time [sec.]

1200

250
baselinesetup
12 solvessetup
12 solves

1000

baseline
new kappa
12 solves
new kappa
12 solves

200

800
150
600
100
400
50

200
00
baselsine

AMG

AMG:tunedbaselsine

Up-down quark case

Yusuke Namekawa(KEK)AMG

AMG:tuned

Strange quark case– 11 / 21 –HPC-Phys meeting4

Additional hot topics with multiple right hand sideAxnrhs=bnrhs:=n × n matrixwhereAn ∼ 1010 for a typical lattice QCD∀nrhs = 1, 2, ...
• Block solver(multiple right hand side solver) O’Leary(1980)• Truncated solver Collins,Bali,Schäfer(2007)• Deflation de Forcrand(1996),Lüscher(2007)

Yusuke Namekawa(KEK)– 12 / 21 –HPC-Phys meeting[Block solver(multiple right hand side solver)] O’Leary(1980)

AX = B

instead of Ax = b

where
A:=

n × n matrix,

X, B

:=

n × nrhs matrixn ∼ 1010 for a typical lattice QCD∀nrhs = 1, 2, ...

• The philosophy is sharing Krylov space for multiple right hand sides♦ Practical advantage is better use of cache, which increase thesustained speed by a factor of 2-5• Two problems are known→ Next pageYusuke Namekawa(KEK)– 13 / 21 –HPC-Phys meeting[Block solver(continued)]• There are some attempts in lattice QCD
de Forcrand(1996), Sakurai et al.(2010), Tadano et al.(2010), Nakamura et al.(2011), Birk andFrommer(2012,2014), Clark et al.(2018), de Forcrand and Keegan(2018)♦ Problem 1 : naive Block solver has a gap between true andrecursion residuals→ Improved versions are proposed Dubrulle(2001), Tadano et al.(2009), ...
♦ Problem 2 : Block solver often fails to converge (breakdown andstagnation), though it can be tamed in part by QR decomposition
Dubrulle(2001), Nakamura et al.(2011), ...

→ We do not employ the block solver in a large scale simulation

Yusuke Namekawa(KEK)– 14 / 21 –HPC-Phys meeting[Block solver(continued)]• Block solver(blockCGrQ) gains a factor of 2-5, if it converged

DP := Double Precision
MP := Mixed Precision• Mixed precision is usually faster, but
it is not for a larger number of rhs,
probably due to less stability

Clark et al.(2018)

Yusuke Namekawa(KEK)– 15 / 21 –HPC-Phys meeting[Truncated solver] Collins,Bali,Schäfer(2007)• Truncated solver := many approximate solver results corrected byexact solver result• (cf. all mode averaging := truncated solver + low-mode averaging)
Blum et al.(2012)♦ Truncated solver leads to a factor of 10 speed up for an expectation value constructed from the solution xDEexactO[x]=DO

improvedE[x] ,1

hOi :=

NsampleNsampleXOi

i=1

whereO

improved

[x]

exact
Axn
rhs
approxAx ′nrhs=approxexact(O[x1] − O[x1]) +1

NapproxX

Napproxn′rhsapproxO[x ′]nrhs=2=bn, strict stopping condition (ex. 10−16 )rhs=b ′, loose stopping condition (ex. truncated at Niter = 50)
nrhs′∀nrhs , ∀nrhs = 1, 2, ...

Yusuke Namekawa(KEK)′larger gain for nrhs < nrhs

– 16 / 21 –HPC-Phys meeting[Truncated solver(continued)]• Truncated solver (+ low mode averaging) leads to O(10) speed up♦ NB. care is needed for the choice of the truncation(ex. Niter = 50). Too aggressive choice gives a wrong result.

LMArCostw/o defl.

0,5AMA0mN
mN
m=0.005 m=0.01
24c
24c

GA
GA
m=0.005 m=0.01
24c
24cmN
GA
m=0.001 m=0.001
32cID
32cID

Shintani et al.(2014)

Yusuke Namekawa(KEK)– 17 / 21 –HPC-Phys meeting[Deflation] de Forcrand(1996),Lüscher(2007),...
• Deflation := eigenvectors + solver for the remaining part
♦ Deflation is independent of nrhs i.e. larger nrhs gives larger gain
♦ The gain is a factor of 2-8, though deflation needs overhead andlarge memory consumption of eigenvector estimation
Axn

rhs
Aφi=bnrhs=

λ i φi ,=solverxn+
rhs

i = 1, ..., Ndeflation

Thenxnrhs

Ndeflation
Xi,j=1−1
φi A(φj , bn)ij
rhs

wheresolver
Pdeflation Axn
rhs
Pdeflation xn
rhs==

Yusuke Namekawa(KEK)

Pdeflation bnrhs
xn

rhs−NdeflationXi,j=1– 18 / 21 –−1
Aφi A)(φj , bnijrhs

HPC-Phys meeting[Deflation(continued)]• The gain is a factor of 2-8, though deflation gives overhead♦ NB. the best choice of Ndeflation depends on the system

Lüscher(2007)

Yusuke Namekawa(KEK)– 19 / 21 –HPC-Phys meeting5Summary

Overview of solvers in lattice QCD was presented• Major solvers are covered by open sources(Bridge++, CCSQCDSolverBench, DDalphaAMG, ...)• Benchmark results show the best solver depends on the physics♦ multigrid is best for mud (requiring a huge Krylov space)♦ BiCGStab series and GMRES(m) is faster for ms (requiring notso large Krylov space)• Additional hot topics with multiple right hand side are explained♦ Block solver(multiple right hand side solver) gains a factor of 2-5,
though it often fails to converge♦ Truncated solver leads to O(10) speed up, though too aggressive
truncation gives a wrong result♦ Deflation gains a factor of 2-8, though it needs overhead and large
memory consumption of eigenvector estimationYusuke Namekawa(KEK)– 20 / 21 –HPC-Phys meeting[Not covered in this talk]• Preconditioner♦ Even-odd(red/black), SAP(Schwarz Alternating Procedure), ILU,
SSOR, ...
[Advertise new supercomputer at KEK(SX-AURORA,156.8 TFlop)]• Unfortunately KEK supercomputer had been terminated since 2017,
but is renewal in 2019 http://scwww.kek.jp/• Tuning for discrete vector accelerator leads to O(100) speed up
60.0
50.0

Etime[s]40.0

OFP
SX-AURORA(default)
Bridge++3
8 × 16 SX-AURORA(optimized)
BiCGStab with Wilson, mud

30.0
20.0
10.0
0.0

Yusuke Namekawa(KEK)– 21 / 21 –HPC-Phys meetingAppendix

Yusuke Namekawa(KEK)– 22 / 21 –HPC-Phys meeting[Table of elementary particles and interactions]

http://higgstan.com/ ← the designer got PhD on particle physics experiment

Yusuke Namekawa(KEK)– 23 / 21 –HPC-Phys meeting