Source-linked AI summary
Adaptive Aggregation Based Domain Decomposition Multigrid for the Lattice Wilson Dirac Operator
Andreas Frommer, Karsten Kahl, Stefan Krieg, Björn Leder, Matthias Rottmann
TL;DR
Lattice QCD requires repeated solution of very large Dirac linear systems, but conventional Krylov methods exhibit critical slowing down near physically relevant parameters. The paper combines domain decomposition smoothing with adaptive aggregation-based multigrid, improving setup and enabling recursive multigrid; experiments report substantial speed-ups over established approaches.
Problem
Standard Krylov methods and common preconditioners retain unfavorable scaling and critical slowing down for large lattice-QCD Dirac systems near physical parameters.
Method
The method combines a domain decomposition smoother with non-smoothed adaptive aggregation and improves the setup for constructing a recursive multigrid hierarchy.
Results
22.2 times faster than BiCGStab for the smallest mass parameter, with a factor of 3.9 remaining for the largest value tested.
Takeaways & Limitations
The approach provides an efficient parallel multigrid preconditioner for realistic lattice-QCD systems and can be applied recursively rather than only as a two-grid method.
Takeaways & Limitations
SAP alone cannot remedy unfavorable Krylov-method scaling with respect to system size and quark mass.
Abstract
from arXiv · showhide
In lattice QCD computations a substantial amount of work is spent in solving discretized versions of the Dirac equation. Conventional Krylov solvers show critical slowing down for large system sizes and physically interesting parameter regions. We present a domain decomposition adaptive algebraic multigrid method used as a precondtioner to solve the "clover improved" Wilson discretization of the Dirac equation. This approach combines and improves two approaches, namely domain decomposition and adaptive algebraic multigrid, that have been used seperately in lattice QCD before. We show in extensive numerical test conducted with a parallel production code implementation that considerable speed-up over conventional Krylov subspace methods, domain decomposition methods and other hierarchical approaches for realistic system sizes can be achieved.
1. Introduction.
The paper develops a domain decomposition adaptive algebraic multigrid method for large lattice-QCD Dirac systems, combining complementary hierarchical approaches to address critical slowing down and enable efficient parallel computation.
- Lattice QCD repeatedly requires solving very large sparse linear systems, consuming substantial supercomputer resources.
- For Wilson-discretized Dirac systems, standard Krylov methods suffer critical slowing down near physical mass parameters and lattice spacing a → 0.Existing odd-even preconditioning, deflation, and domain decomposition methods provide speed-ups but do not change this scaling behavior.
- Prior multigrid methods faced difficulties from the random gauge fields, while adaptive algebraic multigrid and inexact deflation provided distinct recent approaches.The paper positions its method at the intersection of these developments.
- The proposed method combines a domain decomposition smoother with non-smoothed aggregation and improves the adaptive setup for constructing the multigrid hierarchy.It is designed for natural parallelization on hybrid architectures and reduces the time required for hierarchy construction.
- Unlike inexact deflation, the method solves the little Dirac system only to low accuracy, allowing recursive application beyond a two-grid method.The paper presents this as turning domain decomposition into a true multigrid method.
2. Lattice Quantum Chromodynamics.
Lattice QCD discretizes four-dimensional quantum chromodynamics on a hyper-cubic lattice, producing very large sparse Dirac systems whose solution is computationally expensive. The paper introduces the Wilson discretization, its clover improvement, and domain-decomposition behavior relevant to multigrid construction.
- QCD predictions requiring non-perturbative computation are obtained by discretizing the theory and solving the resulting lattice systems numerically.
- Wilson discretization uses gauge links and covariant finite differences to represent the continuum Dirac operator on the lattice.
- SAP reduces errors associated with larger eigenvalues well but leaves small-eigenvalue components relatively unaffected, motivating its use as a multigrid smoother.
3. Algebraic Multigrid Methods.
Algebraic multigrid separates error reduction between a smoother and a coarse-grid correction, with aggregation-based operators designed to represent near-kernel components. The paper develops adaptive, domain-aware constructions for lattice QCD, including Γ5-compatible transfer operators and an eigenvector-based setup.
- Multigrid alternates smoothing with coarse-grid correction, recursively applying this two-level procedure until the remaining operator is small enough for direct solution.The coarse correction targets error components that the smoother does not sufficiently reduce.
- The coarse space should approximate the near kernel of D through sparse prolongation and restriction operators, keeping the coarse problem computationally cheaper than the original system.The restriction and interpolation ranges target complementary left and right near-kernel components.
- Local coherence motivates aggregation: low-eigenvalue eigenvectors approximately coincide across many lattice blocks and can therefore be represented from their blockwise components.Aggregation-based transfer operators exploit this structure to construct a smaller coarse representation.
- 3.1. Aggregation-based Intergrid Transfer Operators.: Aggregation partitions fine variables into lattice-block-based aggregates, each inducing N coarse variables while preserving a nearest-neighbor-like coarse connection structure.Interpolation is built by decomposing test vectors over aggregates and locally orthonormalizing their restricted components.
- 3.1. Aggregation-based Intergrid Transfer Operators.: The restriction can be tied to interpolation through R = (Γ5P)^H, while Γ5-compatible aggregates enable the special choice R = P^H.Γ5-compatibility requires each aggregate to contain only spin 0 and 1 variables or only spin 2 and 3 variables.
- 3.2. Petrov-Galerkin Approach in Lattice QCD.: Although singular-vector-based interpolation is theoretically possible, the paper uses eigenvector-based adaptivity because singular-vector approximation increased setup time without improving solver performance numerically.Adaptive setup generates test vectors by smoothing homogeneous equations to expose components that current smoothers or multigrid methods fail to reduce.
4. Multigrid and Inexact Deflation.
The paper recasts inexact deflation in algebraic multigrid terms and contrasts its composition with a multigrid V-cycle. The key distinction is that multigrid can solve the coarse system inaccurately, reducing per-iteration cost and enabling recursive extension.
- Inexact deflation: Inexact deflation constructs a subspace from locally supported test vectors and forms the coarse operator Dc = P^HDP, paralleling aggregation-based interpolation.Its aggregates are not Γ5-compatible, so Γ5 symmetry is not preserved on the coarse grid.
- Inexact deflation: The associated projectors provide a coarse-grid correction and split the solution into coarse and inexactly deflated components.The remaining component is solved with a right-preconditioned Krylov method such as GCR.
- Composition: Multigrid combines SAP smoothing with coarse-grid correction in a V-cycle, whereas inexact deflation applies SAP to the deflated system.Because SAP is non-stationary, flexible Krylov methods are required when these preconditioners are used.
- Comparison: Multigrid permits low-accuracy coarse solves without noticeably degrading the preconditioner, while inexact deflation requires comparable accuracy for the coarse system in each iteration.Recursive application of the two-grid approach can further reduce the cost and produce a true multigrid method.
- Adaptive setup: The adaptive setup uses random test vectors and approximate inverse applications to amplify components associated with the near kernel.The update can be interpreted as a V-cycle iteration or a bootstrap AMG-style test of the coarse correction.
5. DD-αAMG.
DD-αAMG combines domain-decomposition smoothing with aggregation-based adaptive algebraic multigrid. It uses a V-cycle as a flexible GMRES preconditioner and modifies the adaptive setup to improve efficiency.
- Method: DD-αAMG combines a domain decomposition smoother with non-smoothed aggregation for the Wilson Dirac operator.The method uses SAP, implemented as ν red-black Schwarz iterations.
- Coarse operator: The adaptive coarse system is Dc = P^HDP, with P built from test vectors over lattice-block-based, Γ5-compatible aggregates.The aggregation blocks need not match those used by the SAP smoother.
- Krylov preconditioning: A standard V-cycle with ν post-smoothing steps is used as a right preconditioner for flexible GMRES.The method is not run as a stand-alone V-cycle solver.
- Adaptive setup: The adaptive setup modifies inexact deflation by using the current test vector as the initial guess for the second update stage.Extensive testing identified this modification as the most efficient among the tested setup choices.
6. Numerical Results.
The numerical study implements DD-αAMG in parallel C/MPI code and evaluates its setup, solver behavior, configuration dependence, and scaling. Results show substantial speedups and motivate balancing setup cost against the number of right-hand sides.
- Implementation: The implementation uses C and MPI, focuses mainly on a two-grid method, and includes preliminary results for a recursively extended multigrid version.The tests were run on the Juropa cluster using production-oriented parallel code.
- Preconditioning: Odd-even preconditioning improves the Schur-complement condition and typically reduces iteration counts and execution time by a factor of 2–3.The multigrid coarse solves use odd-even preconditioned restarted GMRES with restart length 30.
- Parameters: The default parameter set uses four local iterative steps, a coarse residual reduction factor of 20, and an overall residual reduction factor of 10^10.These stopping criteria are reported for the coarse solve and full FGMRES iteration, respectively.
- Setup evaluation: The setup iteration count should be tuned to the number of right-hand sides because increasing setup cost can reduce solver cost while worsening coarse-grid conditioning.Two setup steps may minimize total time for one solve, whereas five may be preferable for many right-hand sides.
- Scaling studies: The study varies mass, system size, and processor count to assess solver and parallel scaling under changing spectral conditions.Increasing volume raises the density of small eigenvalues, while the mass parameter changes the operator conditioning.
- Comparison with BiCGStab: 3.5 is the speed-up factor over BiCGStab in total timing on the ill-conditioned test case.Excluding setup, the speed-up is 29; approximately 77% of DD-αAMG solve time is spent on the coarse grid.
- Configuration dependence: For ninv ≥4, DD-αAMG iteration counts vary only marginally across six independent configurations.The reported averages use 20 runs, and fluctuations nearly disappear once at least four setup iterations are used.
Mass Scaling.
DD-αAMG is substantially less sensitive to the mass parameter than BiCGStab in the reported scaling study. Its coarse-grid iteration count nevertheless increases strongly toward the smallest mass.
- Timing: 22.2 is the DD-αAMG speed-up over BiCGStab at the smallest mass parameter, while 3.9 remains at the largest.These figures compare timing for one right-hand side.
- Mass dependence: BiCGStab becomes 18.5 times more expensive from the largest to the smallest mass, whereas DD-αAMG timing increases by a factor of 3.2.The DD-αAMG iteration count increases by only 1.4, but its coarse-grid iteration count increases by 8.0.
System Size Scaling.
As lattice size increases at fixed mass parameter and lattice spacing, BiCGStab iteration counts grow with spatial extent, while DD-αAMG remains nearly constant in both iterations and time.
- BiCGStab iteration counts appear to scale with N_s, nearly doubling from N_s = 16 to N_s = 32.
- DD-αAMG maintains an almost constant iteration count and time across the tested lattice sizes.
Weak Scaling.
The weak-scaling study increases lattice size with process count while keeping the local lattice size fixed; DD-αAMG scales well to thousands of processes, with a machine-specific communication penalty at 1,024 processes. Comparisons also show benefits over inexact deflation and AMG, especially for multilevel configurations and high core counts.
- Weak Scaling: DD-αAMG weak scaling increases lattices from 16^4 on one node to 128^2 · 64^2 on 1,024 nodes while keeping 16 · 8^3 local lattice size per core.The test uses 100 DD-αAMG iterations with ninv = 5 in setup.
- Weak Scaling: The scaling study fixes 50 odd-even preconditioned GMRES coarse-grid iterations, corresponding to the same number of 100 MPI Allreduce operations.
- Weak Scaling: DD-αAMG exhibits the usual log(p) communication dependence and an exceptional increase from 512 to 1,024 processes caused by slower MPI Allreduce operations on Juropa.
- Weak Scaling: Apart from the Juropa-specific increase, DD-αAMG scales well up to 8,192 processes.
- Comparison with Inexact Deflation: For an ill-conditioned 48^4 lattice on 2,592 cores, ninv = 5 makes DD-αAMG twice as fast as the fastest inexact-deflation solver, while ninv = 2 minimizes one-right-hand-side setup-plus-solve time.The reported setup-plus-solve times are 6.48s for DD-αAMG and 10.89s for inexact deflation.
- Multilevel Comparison: Three-level DD-αAMG improves solver time over two levels by a factor of 2.5 to 3 for small m0, while four levels are also evaluated.The small-m0 regime corresponds to the physically interesting parameter region described in the passage.
- Comparison with AMG: Compared with AMG, DD-αAMG has slightly fewer solve iterations, lower setup cost, and solve times 2 to 3 times smaller on 8,192 cores.The passage attributes the high-core advantage to lower global communication from the domain decomposition smoother.