Physics Letters B 692 (2010) 157–160

Contents lists available at ScienceDirect

Physics Letters B
www.elsevier.com/locate/physletb

Distance preconditioning for lattice Dirac operators

G.M. de Divitiis a,b , R. Petronzio a,b , N. Tantalo b,c,∗
a
Università di Roma “Tor Vergata”, I-00133 Rome, Italy
b
INFN sezione di Roma “Tor Vergata”, I-00133 Rome, Italy
c
Centro Enrico Fermi, I-00184 Rome, Italy

a r t i c l e i n f o a b s t r a c t

Article history: We propose a preconditioning of the Dirac operator based on the factorisation of a predefined function
Received 23 June 2010 related to the decay of the propagator with the distance. We show that it can improve the accuracy
Accepted 14 July 2010 of correlators involving heavy quarks at large distances and accelerate the computation of light quark
Available online 24 July 2010
propagators.
Editor: A. Ringwald
© 2010 Elsevier B.V. Open access under CC BY license.
Keywords:
Lattice QCD
Preconditioned inversion of Dirac operators

1. Introduction allowed by the computer architecture. Depending on the values

of the quark mass the solution of Eq. (1) poses different numeri-
A key ingredient of lattice QCD simulations is the inversion cal problems. For light quarks the matrix ( D [U ] + M )x, y is badly
of the Dirac operator which enters the generation of unquenched conditioned and its numerical inversion requires a big number of
gauge field configurations and the computation of hadronic observ- iterations. At the other extreme, the number of iterations required
ables. One needs to solve numerically a linear system of the form for heavy quark masses is small but there may be problems with
  the numerical accuracy resulting for the time-slices far away from
D [U ] + M x, y S ( y ) = η(x) (1) the source ( y 0  0). Indeed Eq. (2) is a global condition while
y for heavy quark propagators the time-slices far away from the
where D [U ] is the chosen discretization of the massless interact- source are exponentially suppressed by a factor of the order of
ing Dirac operator, M is the quark mass in lattice units, η(x) is exp(− M y 0 ) and give a negligible contribution to the norm on the
a source vector that is different from zero on a single time-slice left side of Eq. (2). When this problem arises one cannot trust nu-
(that without any loss we shall assume to be at x0 = 0). The so- merical results at large times and it becomes impossible to extract
lution S ( y ) is obtained by iterative numerical algorithms, solvers, physical informations by fitting the leading exponentials contribut-
devised to invert so-called sparse matrices, like the matrices that ing to correlation functions.
result from the discretization of differential equations by finite dif- In order to alleviate both difficulties, we propose a precondi-
ferences methods. In this Letter we shall not discuss the details of tioning of the Dirac operator that factorises from the propagator
any particular solver (see Ref. [1] for a complete review and for an a function aiming to modify its leading decay with the distance.1
updated list of references). For any solver one checks if the condi- The simplest choice is to factorize a function α ( y 0 ), to solve nu-
tion merically the preconditioned equation, and to restore the original
 
   propagator by multiplying each time-slice for 1/α ( y 0 ). α ( y 0 ) is
 D [U ] + M x, y S ( y ) − η(x) < r
n
(2) defined such that all the different time-slices give comparable con-

y tributions to the calculation of the residue in the preconditioned
case. Our preconditioning is inspired to what is usually done in
is satisfied. Here S n ( y ) is the tentative solution at iteration num-
deriving the Eichten and Hill [2] lattice HQET action but of course
ber n, the norm is any good norm in field space and r, the residue,
does not introduce any approximation. Indeed, the choice above is
is the global numerical accuracy requested for the solution. Typi-
suited for heavy quark propagators, while for light quark masses
cally r is a small number of the order of the arithmetic precision
we will introduce a generalisation of the factorised function.

* Corresponding author.
1
E-mail address: nazario.tantalo@roma2.infn.it (N. Tantalo). See Refs. [3,4] for different approaches.

0370-2693 © 2010 Elsevier B.V. Open access under CC BY license.

doi:10.1016/j.physletb.2010.07.031
158 G.M. de Divitiis et al. / Physics Letters B 692 (2010) 157–160

Fig. 1. The red points correspond to the correlation function −C PP ( y 0 ) obtained by inverting the lattice Dirac equation in the free theory with M  0.5 and with a residue
r = 10−11 . The black points correspond to the same quantity but have been obtained with a residue r = 10−6 . The blue points correspond to the correlation function
 ( y ) obtained by solving the preconditioned lattice Dirac equation with M  0.5 and α = 0.4. The two black lines correspond to r 2 for the two values of the residue
−C PP 0 0
used in the calculations. We use logarithmic scale on the y-axis. (For interpretation of the references to color in this figure legend, the reader is referred to the web version
of this Letter.)

Fig. 2. The red and black sets of points are the effective masses of the corresponding correlators shown in Fig. 1. The blue set of points is the effective mass of the
corresponding correlator in the top panel multiplied for the restoration factor as explained in the text. (For interpretation of the references to color in this figure legend, the
reader is referred to the web version of this Letter.)

2. Preconditioning heavy quark propagators the results obtained in the preset case with r = 10−11 can be
considered as exact. If instead the time extent of the lattice is
We work with the O (a)-improved Wilson lattice Dirac opera- rather large the brute force approach cannot be considered be-
tor but the numerical problems that we address arise also with cause the required residues would be smaller than what is allowed
alternative discretisations of the continuum action and the pro- on double-precision architectures, even in the case of moderately
posed solution can as well be easily implemented in those cases. heavy quarks. In the case under consideration, by choosing a loose
We have tested our preconditioning scheme both for heavy and precision, i.e. a residue r = 10−6 , we make the numerical problem
for light quark masses and in the free and in the interacting case. evident and we show that also such an “extreme” situation can
We start with the results for the heavy quarks. We first want to be recovered by using our proposal. Moreover, we notice that a
pick up a case where the problem arises. As an example, we have residue r = 10−6 is the smallest allowed on single-precision archi-
calculated the correlation function tectures that presently are considerably much faster than double-
precision ones.
  
C PP ( y 0 ) = − tr S † ( y ) S ( y ) (3) We now come to the proposed solution. We redefine the quark
fields and the propagators as follows
y

by solving Eq. (1) for a heavy quark propagator of mass M  0.5 S ( y) = α( y0 ) S  ( y)

in the free theory for different choices of the residue. More pre-
cisely the red points in Fig. 1 have been obtained with a residue
η( y ) = α ( y 0 )η ( y ) (4)
r = 10−11 while the black points with a residue r = 10−6 and the Once the previous expressions are inserted in Eq. (1) we get the
two black lines correspond to the squares of these two values of r. preconditioned system
As is clearly visible from Fig. 1, and from Fig. 2 where we show
 
the effective masses of the correlations shown in Fig. 1, the black D  [U ] + M x, y S  ( y ) = η (x) (5)
points start to deviate from the red ones for y 0  18, i.e. when
y
the correlator, which in this case is just the square module of the
propagator, becomes smaller than the square of the “loose” residue that we solve numerically in place of Eq. (1). In order to write the
r = 10−6 . preconditioned Dirac operator it is sufficient to modify the forward
If the time extent of the lattice is not too large the prob- and backward lattice covariant derivatives in the time direction ac-
lem can be solved by brute force by lowering the residue and cording to
 G.M. de Divitiis et al. / Physics Letters B 692 (2010) 157–160 159

Fig. 3. The red points correspond to the effective mass of the correlation function −C PP ( y 0 ) obtained by inverting the lattice Dirac equation in the interacting theory with
amh  0.35 and with a residue r = 10−11 . The black points correspond to the same quantity but have been obtained with a residue r = 10−6 . The blue points correspond
pcac

to the effective mass of the restored correlation function −C PP  ( y ) obtained by solving the preconditioned lattice Dirac equation with am pcac
 0.35 and α0 = 0.4. (For
0 h
interpretation of the references to color in this figure legend, the reader is referred to the web version of this Letter.)

∇0 S ( y ) = U 0 ( y ) S ( y + 0̂) − S ( y ) Table 1
Gauge configurations have been generated with n f = 2 dynamical O (a)-improved
α ( y 0 + 1)
→ U 0 ( y ) S  ( y + 0̂) − S  ( y ) Wilson quarks with c sw = 1.90952. The figures in the last column correspond to
α( y0 ) the average of the number of iterations required to invert the Dirac equation in the
† † unitary theory by using the SAP+GCR inverter for several values of the precondi-
∇0 S ( y ) = S ( y ) − U 0 ( y − 0̂) S ( y − 0̂) tioning parameter α0 . The values corresponding to α0 = 0.0 have been obtained
without using our preconditioning technique.
α ( y 0 − 1) †
→ S ( y) − U 0 ( y − 0̂) S  ( y − 0̂) (6) β L0 L1 L2 L3 ksea r α0 Iterations
α( y0 )
D5 5 .3 48 × 243 0.13625 10−11 0.0 175
Particular care has to be used at the boundaries of the lattice in D5 5 .3 48 × 243 0.13625 10−11 0.4 141
order to respect the boundary conditions originally satisfied by
E3 5 .3 64 × 323 0.13605 10−10 0.0 99
the quark fields. If as in the case of Fig. 1 S ( y ) satisfies anti- E3 5 .3 64 × 323 0.13605 10−10 0.2 78
periodic boundary conditions along the time direction, it follows E3 5 .3 64 × 323 0.13605 10−10 0.4 69
from Eq. (4) that
E4 5 .3 64 × 323 0.13610 10−10 0.0 115
E4 5 .3 64 × 323 0.13610 10−10 0.2 91
S ( y + L 0 0̂) = − S ( y ) E4 5 .3 64 × 323 0.13610 10−10 0.4 81
α( y0 )
S  ( y + L 0 0̂) = − S  ( y) (7) E5 5 .3 64 × 323 0.13625 10−10 0.0 194
α ( y 0 + L 0 0̂) E5 5 .3 64 × 323 0.13625 10−10 0.2 153
E5 5 .3 64 × 323 0.13625 10−10 0.4 141
The blue points in Fig. 1 correspond to the correlation function
  † 
 
C PP ( y0 ) = − tr S ( y) S  ( y) (8) mass of about amhPCAC  0.35. The unpreconditioned correlators de-
y cay approximately as fast as in the free theory case and from the
difference of the black (unpreconditioned, r = 10−6 ) and red (un-
obtained by solving Eq. (5) with the loose residue r = 10−6 but
preconditioned, r = 10−11 ) sets of data we see the same distortion
after having factorized the function
of Fig. 2. The blue points have been obtained by solving Eq. (5)
  after having factorized α ( y 0 ) with α0 = 0.4 and by restoring the
α ( y 0 ) = cosh α0 ( y 0 − L 0 /2) (9)
results according to Eq. (10). Also in the interacting theory the blue
by setting α0 = 0.4. As expected, the preconditioned correlator points are identical to red points though they have been obtained
stays above the line of the loose precision residue and the “ex- with the same loose residue r = 10−6 used to obtain the black
act” result can be back recovered as follows points.
We close this section by observing that our preconditioning
 2
C PP ( y 0 ) = α ( y 0 ) C PP ( y0 ) (10) technique may be particularly useful when working with the
Schrödinger Functional [5,6] formulation of the theory. In this
In Fig. 2 the blue points correspond to the effective mass of the case, contrary to the case of periodic boundary conditions along
preconditioned correlator after the “restoration” of Eq. (10) and fall the time direction, the correlators decay exponentially over the
exactly on top of the red ones in spite of the fact that they have whole time extent of the lattice and one has to choose very small
been obtained with the same loose precision that affected the non- residues also in computing relatively light quark propagators. We
preconditioned black points. have performed several successful experiments with our precondi-
In Fig. 3 we show the same plot as in Fig. 2 but in the inter- tioning technique also in the Schrödinger Functional case by using
acting theory. The gauge ensamble used correspond to the entry α (x0 ) = exp(−α0 x0 ).
E5 in Table 1. The size of the lattice is L 0 L 1 L 2 L 3 = 64 × 323 and
the hopping parameter of the sea quarks is ksea = 0.13625 corre-
3. Preconditioning light quark propagators
sponding to a PCAC quark mass of about amPCAC sea  0.07. The data
shown in Fig. 3 correspond to a pseudoscalar–pseudoscalar correla-
tor, as in the free theory case, of two degenerate heavy quarks with In this section we shall briefly discuss how the ideas developed
hopping parameters kh = 0.125 corresponding to a PCAC quark and discussed in the previous section can be used to accelerate
160 G.M. de Divitiis et al. / Physics Letters B 692 (2010) 157–160

the numerical calculation of light quark propagators. We start our preconditioned system without compromising the numerical accu-
discussion by generalizing Eq. (4) as follows racy of the solution, an operation that should be performed on
double-precision computer architectures.
S ( y ) = β( y 0 , y 1 , y 2 , y 3 ) S  ( y ) In Table 1 we quantify the gain in computational time that can
be achieved by showing the number of iterations of the SAP+GCR
η( y ) = β( y 0 , y 1 , y 2 , y 3 )η ( y ) (11) solver required to solve the lattice Dirac equation for light quarks
with and without our preconditioning. The SAP+GCR solver has
In the following we shall consider the particular choice
been introduced and explained in details by the author in Ref. [7].
The gauge ensambles used to perform this test have been gener-

3
β( y 0 , y 1 , y 2 , y 3 ) = α( yμ ) ated within the CLS agreement [8] with the parameters given in
the table. In the case of the E-lattices the SAP+GCR solver has
μ=0
been ran on 128 processors of a cluster of PC’s by dividing the

3
1 global lattices into blocks of 44 points. In the case of the D-lattice
= (12) the SAP+GCR solver has been ran on 32 processors of a cluster
cosh[α0 ( y μ − L μ /2)]
μ=0 of PC’s by dividing the global lattices into blocks of 64 points. The
table shows that by increasing the value of the parameter α0 the
and the preconditioned lattice Dirac operator can be obtained as
number of iterations goes down with a time gain that can eas-
easily as before by changing all the covariant derivatives according
ily reach the 30%. In the case under discussion, we checked that
to
higher values of α0 would induce a “heavy quark” like behavior
and produce distorted results for the reasons discussed at length
∇μ S ( y ) = U μ ( y ) S ( y + μ̂) − S ( y ) in the previous section.
α ( y μ + 1) The method discussed in this Letter can be generalised by
→ U μ ( y ) S  ( y + μ̂) − S  ( y ) adding some Dirac structure in the factorised function, an option
α( yμ )
presently under investigation.
† †
∇μ S ( y ) = S ( y ) − U μ ( y − μ̂) S ( y − μ̂)
Acknowledgements
α ( y μ − 1) †
→ S ( y) − U μ ( y − μ̂) S  ( y − μ̂) (13)
α( yμ ) We thank our colleagues of the CLS community for sharing the
efforts required to generate the dynamical gauge ensambles used
and by changing accordingly the boundary conditions in all direc-
in this study.
tions as done in Eqs. (7) for the time direction.
An important difference of the present case with respect to the
References
one discussed in the previous section is that the restoration of the
true propagator must be performed before making the contractions [1] M. Luscher, arXiv:1002.4232 [hep-lat].
needed to build correlation functions by using the first of Eqs. (11). [2] E. Eichten, B.R. Hill, Phys. Lett. B 234 (1990) 511.
Here the preconditioning is not to gain precision, but to ac- [3] A. Juttner, M. Della Morte, PoS LAT2005 (2006) 204, arXiv:hep-lat/0508023.
celerate the convergence of the inversion. Therefore, by applying [4] B. Joo, R.G. Edwards, M.J. Peardon, arXiv:0910.0992 [hep-lat].
[5] M. Luscher, R. Narayanan, P. Weisz, U. Wolff, Nucl. Phys. B 384 (1992) 168,
Eqs. (11), (12) and (13) to the calculation of a light quark prop-
arXiv:hep-lat/9207009.
agator one aims to make the propagator to decay faster than the [6] S. Sint, Nucl. Phys. B 421 (1994) 135, arXiv:hep-lat/9312079.
original unpreconditioned operator. By judiciously choosing the pa- [7] M. Luscher, Comput. Phys. Commun. 156 (2004) 209, arXiv:hep-lat/0310048.
rameter α0 it is possible to change the condition number of the [8] https://twiki.cern.ch/twiki/bin/view/CLS/WebHome.