I IntroductionII Adaptive MultigridIII Formulating an algorithmIV Numerical ResultsV Concluding remarks ReferencesarXiv:1005.3043v2 [hep-lat] 22 Jun 2010Adaptive multigrid algorithm for the lattice Wilson-DiracoperatorR. Babich,1, 2J. Brannick,3R. C. Brower,1, 2M. A. Clark,4T. A. Manteuffel,5S. F. McCormick,5J. C. Osborn,6and C. Rebbi1, 21Center for Co mputational Science, Boston University,3 Cummington Street, Boston, MA 02215, USA2Department of Physics, Boston University,590 Commonwealth Avenue, B oston, MA 02215, USA3Department of Mathematics, The Pennsylvania State University,230 McAllister Building, University Park, PA 16802, USA4Harvard- Smithsonian Center for Astrophysic s,60 Garden Stree t, Cambridge, MA 02138, USA5Department of Applied Mathematics, Campus Box 526,University of Colorado at Boulder, Boulder, CO 80309, USA6Argonne Lea dership Computing Facility,Argonne National Laboratory, Argonne, IL 60439, USA(Dated: May 14, 2010)AbstractWe present an adaptive multigrid solver for application to the non-Hermitian Wilson-Diracsystem of QCD. The key components leading to the success of our proposed algorithm are the useof an adaptive projection onto coarse grids that preserves the near null space of the system matrixtogether with a simplified form of the correction based on the so-called γ5-Hermitian symmetry ofthe Dirac operator. We demonstrate that the algorithm nearly eliminates critical slowing dow n inthe chiral limit and th at it has weak dependence on the lattice volume.PACS numbers: 11.15.Ha, 12.38.Gc1I. INTRODUCTIONPerhaps the most severe computationa l challenge facing the lattice a pproach to quantumchromodynamics is the divergent increase in cost as o ne approaches the chiral limit requiredfor the experimental values of the up and down quark ma sses. (Similar difficulties confrontfield theories conjectured for physics beyond the standard model as well.) The cause iswell known: as the fermion mass approaches zero, the Dirac operator becomes singular(Re(λmin)→ 0), causing “critical slowing down” of the standard Krylov solvers typicallyused to find the propagators. This is unavo idable for all single-grid solvers. Improvingconvergence with a suitable preconditioning has been a main topic of research in latt iceQCD for many years but has, until recently, met very limited success in practice.Eigenvector deflation [1, 2] is a popular technique for accelerating solver convergenceand is generally successful provided sufficiently many eigenvectors are used in the deflationprocess; exact deflation approaches are, however, expected to scale as the square of thelattice volume O(V2) and, thus, become ineffective for large volumes. An alternative is thelocal deflation approach o f [3].Here approximate eigenvectors are used in the deflation process, and due to the localcoherence (see below) of the low modes of the Dirac operator, only a volume-independentnumber of low-mode prototypes are required. As a result, an effective deflation of theoperator is achieved with a computational effort growing approximately like V rather thanV2.Here we present an adaptive multigrid (MG) solver f or t he Dirac equationD(U)ψ = b , (1)whereDx,y(U) = (4 + m)δx,y−4Xµ=1[1 − γµ2Uµxδx+ˆµ,y+1 + γµ2Uµ†x−ˆµδx−ˆµ,y] (2)is the Wilson lattice discretization of the Dirac operator. This is expressed (implicitly) asthe tensor product of 4 × 4 Dirac gamma matrices γµand 3 × 3 SU(3) gauge matrices Uµ(x)on the nearest neighbor links (x, y) of a hypercubic spacetime lattice. While this matrixis not Hermitian, it satisfies γ5-Hermiticity (D†= γ5Dγ5); the corresponding Hermitian2matrix, H = γ5D, is maximally indefinite. The eigenvalues of D are complex and satisfyRe(λmin) > 0 for physical values of the simulation parameters.In a previous work [4], we presented an algorithm for solving the no r ma l equations ob-tained from the Wilson-Dirac system in the context of 2 dimensions, with a U(1) gaugefield. Here, we extend this approa ch to directly solve the Wilson-Dirac system and applythe resulting algorithm to the full 4-dimensional SU( 3) problem.II. ADAPTIVE MULTIGRIDThe “low” modes, eigenmodes with small-in-magnitude eigenvalues of the system matrix,are typically those r esponsible for t he poor convergence suffered by standard iterative solvers(relaxation or Krylov methods). As the o perator becomes singular, the error in the iterativelycomputed solution quickly becomes dominated by these modes. In the free field theory, theseslow-to-converge modes are geometrically smoo th and, hence, can be well represented on acoarse grid using fewer degrees of freedom. Moreover, these smooth modes on the finegrid now a gain become rough (high frequency) modes on the coarse grid. This observationmotivated the classical geometric MG approach, in which simple local averag ing is used torestrict residuals to the coarse grid and linear interpolation is used to t r ansfer corrections(obtained from solving the coarse-grid error equation) to the fine grid. We hereafter denotethe interpolation op era tor by P and restriction operator by R.Given a Hermitian positive definite (HPD) operat or A, taking the restriction operatoras R = P†and the coarse-grid o perator as Ac= P†AP gives the optimal (in an energy-norm sense) two-grid correction. It is natural to extend this recursively by defining theproblem on coarser and coarser grids until the degrees of freedom have been reduced enoughto permit an exact solve. When combined with m pre-relaxations (before restriction) andn post-relaxations (after prolongation) on each level, we arrive at the usual V (m, n)-cycle.Such an MG process is known to eliminate critical slowing down f or discretized elliptic PDEproblems, scaling as O(V ) [5].Explicitly, the error propagation operator for the two-grid solver with a single post-relaxation smoother S is given byET G= S(I − P (P†AP )−1P†A). (3)3The performance of the MG algorithm is r elat ed to range(P ) and how well this approxi-mates the slow-to-converge modes of the chosen relaxation procedure. Given a convergentsmoother, the two-grid algorithm can be shown to converge (i.e., kET GkA< 1) providedthat range(P ) approximates eigenvectors with error proportional to the size of their corre-sponding eigenvalues.For the Wilson-Dir ac system in the interacting theory, the low modes are not geomet-rically smooth, and so classical MG approaches, which assume the slow-to-converge erroris