A Parallel SSOR Preconditioner for Lattice QCD

TL;DR

This paper introduces a parallel SSOR preconditioner for lattice QCD that significantly reduces iteration counts and CPU time in simulations, especially on parallel computers like Quadrics QH4.

Contribution

It proposes a new parallelizable SSOR preconditioning scheme based on local lexicographic ordering, improving efficiency in lattice QCD computations.

Findings

Achieves a factor of 2 reduction in iterations with BiCGstab.

Realizes a 1.7 times CPU-time reduction on 512-processor Quadrics QH4.

Effective implementation using the Eisenstat-trick and machine-specific optimizations.

Abstract

We present a parallelizable SSOR preconditioning scheme for Krylov subspace iterative solvers which proves to be efficient in lattice QCD applications involving Wilson fermions. Our preconditioner is based on a locally lexicographic ordering of the lattice points. In actual hybrid Monte Carlo applications with the bi-conjugate gradient stabilized method BiCGstab, we achieve a gain factor of about 2 in the number of iterations compared to conventional odd-even preconditioning. Whether this translates into similar reductions in run time will depend on the parallel computer in use. We discuss implementation issues using the `Eisenstat-trick' and machine specific advantages of the method for the APE100/Quadrics parallel computer. In a full QCD simulation with Wilson fermions on a 512-processor Quadrics QH4 we find a gain in cpu-time of a factor of 1.7 over odd-even preconditioning for a 24^3 x 40 lattice.

Equations

$$Mx=\varphi ,$$

$$V_1^{-1}M{V}_2^{-1}\tilde{x}=\tilde{\varphi },\text{mbox}where\tilde{\varphi }=V_1^{-1}\varphi ,\tilde{x}=V_{2}x.$$

$$M=I-L-U$$

$$V_1=I-L,V_2=I-U.$$

$$v=V_2^{-1}r,u=V_1^{-1}(r-v),w=v+u.$$

$$r_i=(I-L){\tilde{r}}_i.$$

$$(I-L)x^{k+1/2}=U{x}^k+\varphi ,(I-U)x^{k+1}=L{x}^{k+1/2}+\varphi ,k=0,1,\dots$$

$$M_{x,y}$$

$$M=\left(\begin{array}[]{cc}\mbox{{\bf 1}}&-\kappa D_{oe}\\ -\kappa D_{eo}&\mbox{{\bf 1}}\end{array}\right)$$

$$I-L=\left(\begin{array}[]{cc}\mbox{{\bf 1}}&0\\ -\kappa D_{eo}&\mbox{{\bf 1}}\end{array}\right)\mbox{ with }(I-L)^{-1}=\left(\begin{array}[]{cc}\mbox{{\bf 1}}&0\\ \kappa D_{eo}&\mbox{{\bf 1}}\end{array}\right)$$

$$I-U=\left(\begin{array}[]{cc}\mbox{{\bf 1}}&-\kappa D_{oe}\\ 0&\mbox{{\bf 1}}\end{array}\right)\mbox{ with }(I-U)^{-1}=\left(\begin{array}[]{cc}\mbox{{\bf 1}}&\kappa D_{oe}\\ 0&\mbox{{\bf 1}}\end{array}\right).$$

$$(I-L)^{-1}M(I-U)^{-1}=\left(\begin{array}[]{cc}\mbox{{\bf 1}}&0\\ 0&\mbox{{\bf 1}}-\kappa^{2}D_{eo}D_{oe}\end{array}\right),$$

$$y_x=p_x+\kappa \left(\sum _{\mu =1}^4{m}_{x,x-\mu }^{+}{y}_{x-\mu }\right),$$

$$\text{mbox}stopiterationif\frac{||Mx-\varphi |{|}_2}{||\varphi |{|}_2}<ϵ,$$

Full text

1 Introduction

The computation of the effects of fermionic forces (from virtual quark-antiquark creation and annihilation processes) onto the vacuum structure of quantum chromodynamics (QCD) presents a severe bottleneck to the numerical evaluation of this fundamental theory of strong interactions.

The stochastic sampling of the vacuum state by the standard Hybrid Monte Carlo algorithm (HMC) e. g. implies the continual computation of fermionic Greens functions on a stochastic source $\phi$. In terms of the discretized Dirac operator $M$, the fermionic Greens function (or ‘quark propagator’) is nothing but the solution of the linear equation

[TABLE]

where $M$ is a huge sparse matrix, of rank $r=3\times 4\times V$ and $V$ is the volume of the underlying 4-dimensional space-time lattice. There are different discretizations of the Dirac operator in use. For the study of weak interaction processes, the Wilson fermion scheme [1] appears to be the most attractive.

Today, parallel computers have matured and provide sufficient compute power to drive simulations of full QCD with Wilson fermions to the level of 100 and more Teraflops-hours. This performance allows for sufficient sampling of QCD vacuum configurations on reasonably sized lattices. As a consequence there is topical interest to (i) pioneer new algorithms for generating QCD configurations on one hand [2, 3, 4] and (ii) achieve progress in improving HMC (as the method of choice) on the other hand111It appears promising to find that this complementary strategy of tackling the fermion problem in numerical field theory joins both physicists and applied mathematicians..

Meanwhile, the stabilized bi-conjugate gradient algorithm (BiCGstab) [5] has been established as an efficient inverter in lattice QCD applications [6, 7, 8, 9] since it requires less iterations and less computing time as compared to the minimal residual algorithm or CG on the normal equations: improvements of about 50 % (with respect to the conjugate gradient) can be achieved in the regime of small quark masses. Numerical studies in Refs. [6, 7, 10] indicate that further progress is now to be expected through preconditioning.

In the past, various approaches for preconditioning the Dirac matrix $M$ have been taken:

The first attempt started from the perception that in the weak coupling limit the quark propagator is diagonal in momentum space and hence Fourier transformation would act as a useful preconditioner. Such ‘Fourier acceleration’ [11], however, needs smoothness of fields and hence gauge fixing. This handicaps the computation with a substantial overhead which at the end of the day eats up most of the previous gain. 2. 2.

The state-of-the-art222Non-Krylov-subspace based methods like multigrid as of today did not achieve a level of maturity to be competitive with inverters like CG in lattice gauge theory. This is due to the disordered nature of the gluonic field configurations that act as fluctuating coefficients within the discretized differential operator stencils. preconditioning approach in lattice QCD333The transformation of $M$ from natural into odd-even order might be looked upon as a particular SSOR preconditioning, see Section 3. rests upon the odd-even decomposition of the matrix $M$ [12] that yields an efficiency gain by a factor of 2 - 3 when inverting $M$ in the context of actual fermionic force and quenched propagator computations [13, 6]. It is crucial that this second preconditioning approach lends itself easily to parallelization. 3. 3.

A third and promising method has been evaluated by Oyanagi in the 80s [14]. He used the standard incomplete LU (ILU) factorization of the matrix $M$ based on the natural, globally lexicographic ordering of the lattice points. For Wilson fermions with Wilson parameter $r=1$, the ILU preconditioner is identical to the SSOR preconditioner with respect to that ordering.

Ref. [14] found ILU preconditioning (for the conjugate residual method) to outperform the odd-even decomposition on vector machines, achieving a large gain factor (in terms of iteration numbers) over conventional unpreconditioned methods. As it stands, Oyanagi’s method works satisfactorily on vector machines. However, on local memory or grid-oriented parallel computers, this preconditioner can hardly be implemented efficiently, since parallelism can only be achieved by working concurrently on lattice points lying on the same diagonal hyper-plane. Hockney [15, 16] reports that the run times on several parallel and vector computers actually degrade as compared to the odd-even preconditioning.

Using the so-called ‘Eisenstat-trick’ (see Section 2), the ILU preconditioning can be implemented in a more efficient manner than in [14, 15]444Incorporating the Eisenstat trick into the Oyanagi ILU preconditioner could thus lead to a revision of some of the conclusions from [15], so that a more detailed numerical study is advisable..

In this paper, which is based on the results in [17], we intend to open the stage for use of general parallel SSOR preconditioning techniques in lattice QCD based on appropriate orderings of the lattice points. Our approach may be regarded as a generalization of the odd-even (or red-black) ordering to a multiple color layout [18] or, alternatively, as a localization of the globally lexicographic ordering. The parallel implementation will be carried out in combining two steps:

Partition the whole lattice into equally shaped sublattices and introduce a lexicographic ordering on each sublattice. 2. 2.

During the SSOR preconditioning step, in parallel sweep through all sublattices in lexicographic order.

In Section 2 we discuss the SSOR preconditioner in detail and we give efficient implementations for the SSOR preconditioned BiCGstab and MR methods involving the Eisenstat-trick. Section 3 explains how different orderings of the lattice points lead to variants of the SSOR preconditioning procedure. In particular, this section introduces the locally lexicographic ordering through which we obtain the improvements over odd-even preconditioning.

In Section 4, the parallel implementation in conjunction with sub-blocking is discussed.

Section 5 presents the behavior of our new preconditioning method for Hybrid Monte Carlo applications and its performance on APE100/Quadrics parallel computers. On lattices of size $8^{3}\times 16$ we test the efficiency of the method numerically where we find improvement factors of about 2 in terms of iteration numbers compared to the standard odd-even method. We discuss the effects of granularity and further report about first experiences with the preconditioner in a large scale simulation using a 512-node APE100/Quadrics on a $24^{3}\times 40$ lattice, in the chirally sensible region.

2 SSOR Preconditioning

When preconditioning (1), we take two non-singular matrices $V_{1}$ and $V_{2}$ which act as a left and a right preconditioner, respectively, i.e. we consider the new system

[TABLE]

(In the following, preconditioned quantities will be denoted by tildes.) We could now apply BiCGstab, as given in Algorithm 2 for the original system, (or any other Krylov subspace method) directly to (2), replacing each occurrence of $M$ and $\phi$ by $V_{1}^{-1}MV_{2}^{-1}$ and $\tilde{\phi}$, resp.

However, this would yield the preconditioned iterates $\tilde{x}^{k}$ together with preconditioned residuals. Therefore, one usually reformulates the algorithm incorporating an implicit back-transformation to the unpreconditioned quantities. The resulting algorithm is similar to Algorithm 2, requiring two additional systems with matrix $V=V_{1}V_{2}$ and two systems with matrix $V_{1}$ to be solved in each iterative step (see [5]).

The purpose of preconditioning is to reduce the number of iterations and the computing time necessary to achieve a given accuracy. This means that $V$ has to be a sufficiently good approximation to the inverse of $M$ (thus decreasing the number of iterations) while solving systems with $V,V_{1}$ should be sufficiently cheap (since these solves represent the overhead of the preconditioned upon the basic method).

In the present work we are only interested in the SSOR (or, more precisely, symmetric Gauß-Seidel) preconditioner. Consider the decomposition of $M$ into its diagonal, strictly lower and strictly upper triangular parts, i.e.

[TABLE]

with $L,U$ strictly lower and upper triangular matrices, respectively. Then the SSOR preconditioner is given by

[TABLE]

For the SSOR preconditioner we have $V_{1}+V_{2}-M=I$. This important relation can be exploited through the ‘Eisenstat-trick’ [19], because we now have $V^{-1}MV_{2}^{-1}=V_{2}^{-1}+V_{1}^{-1}(I-V_{2}^{-1})$, so that the matrix vector product $w=V_{1}^{-1}MV_{2}^{-1}r$ can economically be computed via

[TABLE]

In this manner we get the standard algorithmic formulation of the SSOR preconditioned BiCGstab method, stated as Algorithm 2.

The Eisenstat trick is not restricted to the BiCGstab method but can be applied in any Krylov subspace method. As another example, we state the SSOR preconditioned minimal residual (MR) method as Algorithm 2.

Note that these algorithms use the preconditioned residuals $\tilde{r}_{i}$ which are related to the unpreconditioned residuals $r_{i}=\phi-Mx_{i}$ via

[TABLE]

To save computational costs, a stopping criterion for Algorithm 2 will usually be based on $\tilde{r}_{i}$. Upon successful stopping, one can then compute $r_{i}$ and test for convergence using $r_{i}$. If the solution is not yet accurate enough, one continues the iteration with a stronger stopping criterion. In our numerical experiments to be reported in Section 5 it turned out that such additional steps were never necessary.

Also note that multiplications with $M$ are completely absent in this formulation, the only matrix operations being solves with each of the matrices $I-L$ and $I-U$. Since these matrices are triangular, the solves can be done directly via forward or backward substitution, respectively. For example, denoting the components of a vector $y$ by $(y)_{i}$ and the entries of $L$ by $(L)_{ij}$ and similarly for $U$, the forward solve $(I-L)y=p$ and backward solve $(I-U)z=p$ become simply

[TABLE]

forward solve

[TABLE]

backward solve

Due to the sparsity pattern of $M$, most of the entries in $L$ and $U$ are zero so that only a few $j$ actually contribute to the sums over $j$. We will come back to this point in more detail in the next section.

It is important to note that for $i$ fixed, the number of non-zero entries of $L$ and $U$ involved when updating $y_{i}$ in the forward together with $z_{i}$ in the backward solve is just the number of non-zero entries of the $i$-th row of the matrix $M$. Therefore, in terms of computational cost, a forward followed by a backward solve is quite precisely as expensive as a multiplication with $M$ (which is required in the unpreconditioned method). Hence, due to the Eisenstat-trick, in terms of matrix vector operations, each step of the SSOR preconditioned BiCGstab (or MR) method requires the same work as in the unpreconditioned method.

We finally note that the ‘classical’ symmetric Gauß-Seidel iteration

[TABLE]

represents a sequence of alternating forward and backward solves. This explains the terminology SSOR preconditioner where SSOR stands for symmetric successive over-relaxation, the over-relaxed variant of the symmetric Gauß-Seidel method. See [20], e.g.

3 Orderings

When writing down equation (1) with Wilson fermion matrix $M$555The 4 matrices $\gamma_{\mu}$ are $4\times 4$ Dirac matrices, and the $U_{\mu}$ are $3\times 3$ SU(3) matrices that represent the gluonic degrees of freedom on the lattice.

[TABLE]

we have the freedom to choose any ordering scheme for the lattice points $x$. Different orderings yield different matrices $M$, however, which are permutationally similar to each other, i.e. one matrix can be retrieved from the other via the transformation $M\rightarrow P^{\dagger}MP$ with a permutation matrix $P$. In general, $P^{\dagger}LP$ and $P^{\dagger}UP$ will not represent the strictly lower or strictly upper triangular part of $P^{\dagger}MP$. Consequently, the SSOR preconditioned matrix for $P^{\dagger}MP$ will not be permutationally similar to that of $M$, so that the quality of the SSOR preconditioner will usually depend on the ordering scheme chosen. One purpose of this paper is to explain that the odd-even preconditioning and the Oyanagi preconditioning can both be regarded as the SSOR preconditioning belonging to particular orderings of the grid points666It is shown in [17] that the SSOR preconditioner for the Wilson fermion matrix (with Wilson parameter $r=1$) is identical to the ILU preconditioner for any ordering of the lattice points. We will then introduce a new ‘locally lexicographic’ ordering which is adapted for parallel computation and for which the SSOR preconditioned system can be solved more rapidly than with odd-even preconditioning with respect to both, the number of iterations as well as actual run time on particular parallel computers.

Consider an arbitrary numbering (ordering) of the lattice points. For a given grid point $x$, the corresponding row in the matrix $L$ or $U$ contains exactly the coupling coefficients of those nearest neighbors of $x$ which have been numbered before or after $x$, resp. Therefore, a generic formulation of the forward solve for this ordering is given by Algorithm 3. The backward solves are done similarly, now running through the grid points in reverse order and taking those grid points $x\pm\mu$ which were numbered after (instead of before) $x$. Due to this analogy, we will restrict our discussion to the forward solves in the sequel.

Odd-even ordering.

As a first specific example, let us consider the familiar odd-even ordering where all odd lattice points are numbered before the even ones. From our generic formulation in Algorithm 3 we then get Algorithm 3 for the forward solve.

In traditional QCD computations, the odd-even preconditioning is not implemented by using the above formulation of the forward (and backward) solve in Algorithm 2. This is due to the fact that—very exceptionally—for this particular ordering the inverses of $I-L$ and $I-U$ can be determined directly: With the odd-even ordering, the matrix $M$ has the form

[TABLE]

so that

[TABLE]

and

[TABLE]

Hence

[TABLE]

where $\mbox{{\bf 1}}-\kappa^{2}D_{eo}D_{oe}$ is called the matrix of the odd-even reduced system. Very exceptionally again, in this particular case, the preconditioned matrix $(I-L)^{-1}M(I-U)^{-1}$ thus has such a simple structure that it is possible to apply Algorithm 1 directly, using the preconditioned matrix explicitly. Since the matrix $(I-L)^{-1}M(I-U)^{-1}$ is $2\times 2$ block diagonal with the first diagonal block being the identity, the computations for the second diagonal block are completely decoupled from those of the first block and they are identical to those which are performed when applying the algorithm for the odd-even reduced system directly. This is precisely what is done in traditional odd-even preconditioning. So traditional odd-even preconditioning is equivalent to SSOR preconditioning for the odd-even ordered system.

So far, odd-even preconditioning is generally considered as the only successful preconditioner in a parallel computing environment. For typical, realistic configurations, it gains a factor 2-3 in the numbers of iterations and also in computing time as compared to solving the unpreconditioned system.

Globally lexicographic ordering.

Assume now that $M$ is given with respect to the natural (lexicographic) ordering of the lattice points. This means that grid point $x=(i_{1},i_{2},i_{3},i_{4})$ is numbered before $x^{\prime}=(i_{1}^{\prime},i_{2}^{\prime},i_{3}^{\prime},i_{4}^{\prime})$ if and only if ($i_{4}<i_{4}^{\prime}$) or ($i_{4}=i_{4}^{\prime}$ and $i_{3}<i_{3}^{\prime}$) or ($i_{4}=i_{4}^{\prime},i_{3}=i_{3}^{\prime}$ and $i_{2}<i_{2}^{\prime}$) or ($i_{4}=i_{4}^{\prime},i_{3}=i_{3}^{\prime},i_{2}=i_{2}^{\prime}$ and $i_{1}<i_{1}^{\prime}$). The corresponding forward solve is given as Algorithm 3.

The SSOR preconditioning for the lexicographic ordering, applied to the MR and other conjugate residual methods, was considered by Oyanagi [14]. He showed that it yields a further improvement over odd-even preconditioning as far as the number of iterations is concerned. However, its parallel implementation is more difficult and less efficient, since only grid points lying on the same diagonal hyper-plane can be worked upon in parallel in the forward and backward solves. Hockney [16] reports that on the 592 processor ACPMAPS at Fermi-Lab, the run time (without the Eisenstat trick) actually degrades as compared to odd-even preconditioning. See also [15].

Locally lexicographic ordering.

Unlike the lexicographical and the odd-even ordering, the ordering we propose now is adapted to the parallel computer used to solve equation (1). We assume that the processors of the parallel computer are connected as a $p_{1}\times p_{2}\times p_{3}\times p_{4}$ 4-dimensional grid. Note that this includes lower dimensional grids by setting some of the $p_{i}$ to 1. The total number of processors is $p=p_{1}p_{2}p_{3}p_{4}$. The space-time lattice can be matched to the processor grid in an obvious natural manner, producing a local lattice of size $n^{loc}_{1}\times n^{loc}_{2}\times n^{loc}_{3}\times n^{loc}_{4}$ with $n^{loc}_{i}=n_{i}/p_{i}$ on each processor. Here, for simplicity, we assume that each $p_{i}$ divides $n_{i}$ and that we have $n^{loc}_{i}\geq 2$ for $i=1,\ldots,4$.

Let us partition the whole lattice into $n^{loc}=n^{loc}_{1}n^{loc}_{2}n^{loc}_{3}n^{loc}_{4}$ groups. Each group corresponds to a fixed position of the local grid and contains all grid points appearing at this position within their respective local grid. Associating a color with each of the groups, we can interpret this process as a coloring of the lattice points. For the 2-dimensional case, such a coloring is depicted in Figure 1 (with $n^{loc}_{1}=n^{loc}_{2}=4$), where the 16 colors are denoted by the letters $a$ – $q$.

We now consider an ordering where all points of color $b$ are ordered after all those of color $a$ if on the local grids the position belonging to color $a$ is lexicographically less than that of color $b$. In Figure 1, this corresponds to the alphabetic ordering of the colors $a$ – $q$. Such an ordering is termed locally lexicographic. From now on we will use the prefix ‘$ll$-’ in expressions like ‘$ll$-first’, ‘$ll$-ordering’ to make clear that they refer to the locally lexicographic ordering.

Since we assumed $n^{loc}_{i}\geq 2$ for all $i$, all nearest neighbors of a given grid point have colors different from that point. This implies that when performing the forward and backward solves in Algorithm 2, grid points having the same color can be worked upon in parallel, thus yielding an optimal parallelism of $p$, the number of processors. A formulation of the $ll$-forward solve is given as Algorithm 3. Here, we use ‘$\leq_{ll}$’ as a symbol for ‘$ll$-less than’.

For grid points lying in the ‘interior’ of each local grid, we have $x-\mu\leq_{ll}x\leq_{ll}x+\mu$ for $\mu=1,\ldots,4$. The update thus becomes

[TABLE]

whereas on the ‘local boundaries’ we will have between 0 (for the $ll$-first point) and 8 (for the $ll$-last point) summands to add to $p_{x}$. The second sum in the update formula then contains those grid points which belong to neighboring processors. The arrows in Figure 1 illustrate this situation for the 2-dimensional analogue, where 2 grid points contribute to the update at locally interior points (colors $f,g,k,l$) and 0 to 4 on the local boundaries. (An arrow from point $x\pm\mu$ to $x$ means that $x\pm\mu$ appears in the sum for updating point $x$.)

With the lexicographic ordering on the whole lattice, a forward solve followed by a backward solve has the effect that the information at any grid point will have spread out to all other grid points. This yields a heuristic justification of why this ordering is superior (as far as convergence properties are concerned) to the odd-even ordering, where the information is passed to the (second-) nearest neighbors only. The $ll$-ordering can be regarded as a compromise between these two extreme cases, spreading information to all points within a local grid. The $ll$-forward and $ll$-backward solves parallelize efficiently and better than with the lexicographic ordering on the whole grid. The parallelism achieved is $p$ and thus less than with the odd-even ordering, but it is optimal since we have $p$ processors. If we change the number of processors, the $ll$-ordering, and consequently the properties of the corresponding SSOR preconditioner will change, too. Heuristically, we expect the convergence properties to degrade as the size of the local grid becomes smaller but to always remain better than with the odd-even preconditioner (as long as $n_{i}^{loc}\geq 2$ for $i=1,\ldots,4$).

4 Parallel Implementation

In this section we give a detailed discussion on the realization of the $ll$-forward solve on a parallel computer. In order to render all necessary communication transparent we assume a model parallel computer where communication is done via explicit message passing between processors.

As in the previous section, we assume that the processors are connected as a 4-dimensional (or lower dimensional) grid and that the space-time lattice is mapped onto that grid in the obvious manner, yielding a local subgrid on each processor. The local coordinates of a grid point $x$ within a subgrid are denoted by the 4-tuple $x=(i_{1},i_{2},i_{3},i_{4})$. For processor $\pi$ of the processor grid we denote its neighbors in dimension $\mu$ by $\pi\pm\mu$. Incorporating all necessary data transmission into Algorithm 3 via message passing statements—sending messages as soon as possible while receiving them as late as possible—we end up with Algorithm 4 which formulates the $ll$-forward solve for processor $\pi$.

Note that if processor $\pi$ issues a ‘send’ to transfer $y_{x}$ to processor $\pi-\mu$ because of $i_{\mu}=1$, the corresponding ‘receive’ (for $i_{\mu}=n_{\mu}^{loc}$) in processor $\pi-\mu$ reads this data as $y_{x^{\prime}}$ with $x^{\prime}=x+\mu$. If $\mu=1$, exactly $n_{1}^{loc}-2$ updates are done between a matching send/receive pair. If $\mu=2$ this number increases to $n^{loc}_{1}(n_{2}^{loc}-2)$, for $\mu=3$ it becomes $n^{loc}_{1}n^{loc}_{2}(n^{loc}_{3}-2)$ and for $\mu=4$ it is $n^{loc}_{1}n^{loc}_{2}n^{loc}_{3}(n^{loc}_{4}-2)$. Therefore, if messages can be handled in parallel to the computation, the message passing can be hidden behind the computation so that the communication overhead will be small. Also note that if $p_{\mu}=1$ the send/receive pairs for direction $\mu$ have to be dropped in Algorithm 4.

Algorithm 4 has to communicate the local boundary values all one by one. This is inherent in the $ll$-ordering and can become a disadvantage on machines which have a large start-up time for communication, irrespective of the message length. In these cases, it is preferable to issue only a small number of (large) messages. Then the following outer product’ of the $ll$-ordering with an odd-even ordering can be used: Divide each local grid into two halves $G_{1},G_{2}$ by bisecting perpendicularly to direction 4 (say), where $n^{loc}_{4}$ is supposed to be a multiple of 2. Identify odd and even processors on the processor grid, neglecting the 4-th coordinate, i.e., processor $\pi=(j_{1},j_{2},j_{3},j_{4})$ is odd if and only if $j_{1}+j_{2}+j_{3}$ is odd. On each even processor, we first take a locally lexicographic ordering on $G_{1}$ and then on $G_{2}$, on odd processors we reverse the role of $G_{1}$ and $G_{2}$. The forward solve with respect to this ordering consists of two half steps. First, the odd processors update $G_{1}$ while the even processors work on $G_{2}$. No communication is necessary there. Then the updated values on the local boundaries have to be communicated (requiring a total of 7 large messages, one for direction 4 and two for each other direction). Now the second half step can be performed, where odd processors update $G_{2}$ and even processors update $G_{1}$. An illustration of this $oe\times ll$–ordering for the 2-dimensional case is given in Figure 2.

Implementing the $ll$-SSOR preconditioning requires to distinguish between grid sites along a given direction. Assume in a forward solve that a grid point is located at the ‘left’ boundary of the local lattice. This site is not updated with data from the left neighbor, and this fact has to be inquired by a conditional statement or a case statement. The second case is that a grid site is located at a ‘right’ boundary where the right neighbor has to be communicated from the neighbor sublattice (in this situation the neighbor processor). The third possibility is that the site is not a boundary grid point. This leads to 81 different cases a site can be associated with in four dimensions.

5 Results

Our numerical tests of the locally lexicographic SSOR preconditioner were performed on APE100/Quadrics machines, a SIMD parallel architecture optimized for fast floating point arithmetic on block data-structures like $3\times 3$ SU(3) matrices or any other data structure that can be blocked similarly. Arithmetic on the Quadrics is compliant with the IEEE 754 single precision standard. We had access to a 32-node Quadrics Q4 at Wuppertal University and a 512-node Quadrics QH4 at the computer center of ENEA/Casaccia Rome.

In order to speed up the performance on the Quadrics it is favorable (as it is the case on most high speed parallel computers) to partition the code into code blocks i. e. sections that can be computed entirely from the registers without recourse to the local or remote memory, thus using the pipelining features of the processors efficiently. This requires to modify Algorithm 4 by moving all if-statements out of the four-fold nested loop over $i_{1},\ldots,i_{4}$ resulting in a sequence of 81 nested loops, one for each possible combination of the three cases $i_{\mu}=1,1<i_{\mu}<n_{\mu}^{loc},i_{\mu}=n_{\mu}^{loc}$ for $\mu=1,\ldots,4$.

Evidence of improvement.

In Fig. 3 first evidence of improvement by $ll$-SSOR preconditioning of BiCGstab is presented.

Shown are three different iteration numbers of a valence mass “trajectory” for a given thermalized full-QCD configuration at $\beta=5.6$ and $\kappa_{\mbox{\tiny sea}}=0.156$ corresponding to an intermediate pion mass. The standard BiCGstab solution of (1) is compared with the odd-even preconditioned method and the new locally lexicographic ordering scheme. The measurements are taken on a lattice of size $8^{3}\times 16$, hence the computational problem being of granularity $g=N/p=256$ on the Quadrics Q4 with 32 nodes. In this exploratory implementation, as discussed in the previous sections, the sublattice administrated by a given processor is used as the local lattice in the $ll$-SSOR scheme. At this place, the convergence stopping criteria,

[TABLE]

are applied on the three different residuals used in the three variant algorithms, respectively, with equal numerical value $\epsilon=10^{-8}$. We verified that the ‘true’ residuals $Mx-\phi$, calculated explicitly from the solution $x$, eventually were smaller for the preconditioned variants of the BiCGstab solver. The gain in iterations is more than a factor of 2 for $\kappa_{\mbox{\tiny valence}}=\kappa_{\mbox{\tiny sea}}$ as is the case within the entire valence $\kappa$ range considered. Compared to the unpreconditioned BiCGstab the gain is even a factor of 4.

The improvement in terms of cpu-time, on the Q4, is smaller than the gain in iteration numbers since the $ll$-SSOR preconditioning leads to certain breaks in code blocks due to the 81 different cases mentioned above. Thus, pipelined arithmetic is not exploited as well as in the other cases. For the given choice of parameters we find an improvement of about a factor of about 1.5 compared to standard $oe$-preconditioning throughout the valence mass range considered.

The convergence behavior of the three methods is demonstrated in Fig. 4 on a typical configuration at the same parameter values for $\kappa_{\mbox{\tiny sea}}=\kappa_{\mbox{\tiny val}}=0.156$.

In addition to the acceleration, the convergence behavior turns out to be much smoother for the $ll$-preconditioned BiCGstab than in the unpreconditioned case. We have plotted both the accumulated and the ‘true’ residual. The latter approaches a horizontal line for $\|r\|_{2}\approx 10^{-6}$ which indicates the limits of single precision accuracy. For $ll$-preconditioning, however, the final residual that can be achieved appears to be somewhat more accurate than for pure BiCGstab.

Local lattice size dependency.

It is important to assess the dependency of the improvement due to $ll$-SSOR preconditioning on the size of the local lattice. As our Q4 is equipped with a fixed number of 32 processors, we emulated larger numbers of local lattices by sub-blocking the system on each processor. The smallest local lattice size we can investigate is a $2^{4}$ lattice, apart from the unpreconditioned case (which we treat as having local lattice size 1) and the $oe$-realization of the solution of (1) (local lattice size 2).

The dependency of the number of iterations on the local lattice size is depicted in Fig. 5 on an equally sized lattice as used before and at equal physical parameters but with $\kappa_{\mbox{\tiny sea}}=0.1575$.

As expected, smaller local lattices lead to a less efficient preconditioning. However, the number of iterations at first increases only slightly if we move from a local lattice size of $256$ to $16$. This behavior might of course be dependent on $\kappa_{\mbox{\tiny sea}}$, hence the crossover to the $oe$-result (the diamond in Fig. 5) might occur at a larger local lattice size for $\kappa_{\mbox{\tiny sea}}$ closer to $\kappa_{c}$.

Speeding up QCD computations.

Next we discuss improvements due to the $ll$-SSOR preconditioning implemented in an actual Hybrid Monte Carlo simulation of full QCD with Wilson fermions. We depict in Fig. 6a a time series of iteration numbers that are each averages over a molecular dynamics trajectory of a length of 125 steps. The measurements are taken on a system of size $24^{3}\times 40$ at $\beta=5.6$ and $\kappa_{\mbox{\tiny sea}}=0.1575$.

The first part of the curve has been generated using standard $oe$-preconditioning for the BiCGstab inverter. The sudden step in iterations is due to the fact that we switched the preconditioner to $ll$-SSOR at this location. The iteration numbers show a quite stable value in both branches of the curve and hence nicely demonstrate the stability of the preconditioning on different background field configurations within the canonical ensemble. In the actual large scale simulation, the $ll$-SSOR preconditioning scheme indeed gives an improvement of nearly a factor of 2 as to the iteration numbers for a quite light pion mass [21].

We remark that we have verified in our simulation that the educated guess procedure as advocated by Brower et al. [22] leads to an improvement that is additive to the gain as achieved by the use of $ll$-SSOR preconditioning.

The time behavior of the system on the 512-node APE100/Quadrics QH4 is depicted in Fig. 6b. We remark that on this machine topology with $8\times 8\times 8$ nodes on a 3-dimensional grid, we arrive at a local grid size of $3\times 3\times 24\times 5$ for the $24^{3}\times 40$ lattice. On a strict SIMD architecture without local addressing, as is the case on the Quadrics, $oe$-schemes are quite tricky to implement and lead to a substantial loss in performance in the range of 30 %. Therefore $ll$-SSOR can help a lot to overcome the SIMD odd-even limitations. The total improvement as to cpu-time is about a factor of 1.7.

Details on our implementation and the experience with $ll$-SSOR preconditioning in the context of large scale HMC simulations in full QCD will be published elsewhere [23].

6 Conclusion and Outlook

We have presented a new local grid point ordering scheme that allows to carry through efficient preconditioning of Krylov subspace solvers like BiCGstab in the context of computations of Greens functions for lattice quantum chromodynamics with Wilson fermions on parallel computers.

Using the Eisenstat-trick we can avoid the explicit multiplication by the fermionic matrix $M$ and we remain with a forward and a backward solve to invert the upper and lower triangular parts of $M$.

In actual Hybrid Monte Carlo simulations with Wilson fermions at small pion masses and on large lattices, we verify an improvement of a factor of 2 in BiCGstab concerning iteration numbers as compared to the state-of-the-art odd-even preconditioned scheme. In this manner, our code is thus speeded up on a 512-node APE100/Quadrics by a factor of about 1.7 at $\beta=5.6$ and $\kappa_{\mbox{\tiny sea}}=0.1575$.

We are confident that the $ll$-SSOR scheme will help to drive future simulations with dynamical Wilson fermions deeper into the chiral region dominated by a small pion mass.

Acknowledgments

Part of the numerical test in this paper have been carried out on the Quadrics QH2 at ENEA/Casaccia/Roma. Th. L. and G. R. thank Profs. N. Cabibbo and A. Mathis for their kind hospitality and Dotts. R. Cannata and R. Guadagni for very valuable support. We are grateful to our friends from the SESAM group for important discussions and their patience while we were blocking the Quadrics Q4 at Wuppertal University and HLRZ/KFA-Jülich. We acknowledge helpful discussions with P. Fiebach from the department of Mathematics, University of Wuppertal and G. Hockney from Fermi-Lab, Batavia, IL.

The work of Thomas Lippert is supported by the Deutsche Forschungsgemeinschaft DFG under grant No. Schi 257/5-1.