ARGONNE NATIONAL LABORATORY 9700 South Cass Avenue Argonne, Illinois 60439

THE EFFECT OF MANY-COLOR ORDERINGS ON THE CONVERGENCE OF ITERATIVE METHODS ^{∗}

Mark T. Jones and Paul E. Plassmann

Mathematics and Computer Science Division Preprint MCS-P292-0292 May 1992

ABSTRACT

It is well known that orderings based on graph colorings can slow the convergence of many iterative methods. However, orderings based on graph colorings can be used to build scalable, parallel iterative methods. For example, SOR with consistent orderings, SSOR, and incomplete Cholesky can all be implemented in a scalable manner using graph colorings. It has been observed experimentally that if, for a fixed problem, the number of colors used is increased, then the SSOR PCG method will converge faster. In this paper we prove that for a model problem the multilevel SOR algorithm with consistent orderings will converge more quickly as the number of colors used is increased. We also give a formula for computing the optimal relaxation parameter and rate of convergence as a function of the number of colors. We prove a similar, but limited, result for ICCG.

These results are useful because they allow the algorithm designer to trade off parallel efficiency for fewer iterations. For example, the nine-point Laplacian can be colored using a minimum of four colors. However, a matrix ordering based on four-colorings may result in poor convergence. If eight rather than four colors are used, the convergence rate will be increased without, perhaps, unduly affecting the parallel efficiency.

Key words: Fourier analysis, ICCG algorithm, multicoloring, parallel algorithms, SOR algorithm, sparse matrices

AMS(MOS) subject classifications: 65F10, 65N22, 65Y05

∗ This work was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under Contract W-31-109-Eng-38.

1. Introduction. Many iterative methods exhibit excellent parallel efficiency if the matrices have been reordered according to a graph coloring [13]. However, for most of these iterative methods, the convergence rate after reordering is slower than that of the “natural” ordering. For example, for the standard five-point difference operator the incomplete Cholesky preconditioner has a condition number of \frac{2}{(2+\sqrt{2})π^{2}h^{2}} with the natural ordering [2] and a condition number of \frac{1}{π^{2}h^{2}} with a red/black ordering [8].

Harrar [5] has observed experimentally for the SSOR PCG method that the convergence rate improves if more colors are used to color the graph of the matrix. Using the natural ordering and the colorings described in §2, we show in Figure 1 the reduction in the residual for the ICCG method. These results are for the standard nine-point difference operator on a 128 × 128 grid. We observe that as the number of colors is increased, the number of iterations necessary to achieve the same accuracy decreases. For large problems, one can use more colors with some

Fig. 1. The convergence rate and the cost of ICCG with varying numbers of colors

sacrifice of parallel efficiency and reduce solution time by reducing the number of iterations. For example, in the right plot in Figure 1 we show the predicted times associated with these multicolorings when we use a variant of the parallel execution model given by Elman and Agrón [3] for solving the same 128 × 128 grid problem with 16 processors. For the model, we chose the parameters τ = 10, c = 10, and s = 1; the same parameters suggested by Elman and Agrón in §6 of their paper.

For two different iterative methods, a stationary method and a preconditioned method, we use Fourier analysis techniques to prove that increasing the number of colors will result in faster convergence. This analytic information can be used by practitioners to make algorithmic design trade-offs on high-performance architectures.

The first method that we examine is a multilevel version of the red/black SOR method [15]: a generalization suggested by Kuo and Levy [9]. By using a multilevel approach, a consistent ordering can be obtained for any number of colors. The optimal relaxation factor can be cheaply computed for any consistently ordered, symmetric positive definite matrix [4]. Consider the four-color example given in Figure 2. At the outermost level of the multilevel algorithm, the four colors are split into two groups to form a 2 × 2 block matrix. This 2 × 2 matrix is trivially consistently ordered, and an optimal relaxation parameter can be computed. At the second level, each block is split into two matrices to form 2 × 2 block matrices for which optimal relaxation parameters can be computed. At the innermost level, the matrices are diagonal. Given that each processor will have approximately

Fig. 2. A multilevel SOR algorithm for four colors

the same number of grid points of a given color, the load balancing throughout the algorithm will be quite good. Further discussion of this algorithm is given in §3. We note that for the specific case of four-coloring the nine-point Laplacian, an alternative to the multilevel SOR algorithm is a nonconsistently ordered SOR algorithm for which the optimal relaxation parameter can be determined by solving a quartic equation [1].

The multilevel SOR algorithm has advantages over the traditional approach of obtaining a consistent ordering by breadth-first search (BFS) and then using red/black SOR. In the BFS approach, subproblems associated with narrow strips of grid points must be solved. For the large matrices that must be solved on parallel computers, these subproblems cannot fit onto a single processor. For example, if we wish to solve a symmetric matrix generated from a finite difference stencil on a cube, the subproblems are planes of the cube. The subproblem size increases as the cube size is increased, but the memory on a single processor is fixed; for a scalable parallel implementation the subproblems must be split across processors. Therefore, one is faced once again with the problem of solving sparse matrices in parallel. In the multilevel algorithm, the parallel solution of the subproblems is inherent in the algorithm. The only operations required are sparse matrix by vector multiplication and the solution of the diagonal submatrices; both operations are parallel.

The BFS algorithm has an additional drawback when solving unstructured problems in parallel. The parallel determination of the subproblem strips is difficult. In general, no good parallel method for BFS exists; therefore, the ordering of the nodes can take an inordinate amount of time on a parallel computer. In contrast, good parallel graph coloring heuristics exist. Recently Jones and Plassmann have given a scalable, parallel heuristic for finding colorings for unstructured, sparse matrices [7].

In addition to the multilevel SOR algorithm, we analyze the conjugate gradient algorithm preconditioned by an incomplete Cholesky factorization [11, 12]. Jones and Plassmann have given scalable, parallel methods for the ICCG algorithm for solving unstructured, sparse, symmetric positive definite matrices [6].

The remainder of the paper is organized as follows. In §2 we analyze the convergence properties of the multilevel SOR algorithm with varying numbers of colors. We further discuss the multilevel algorithm and give experimental results in §3. In §4 we give a similar but more limited analysis for ICCG.

2. Analysis of Multilevel SOR. In this section, we use Fourier analysis to find the spectral radius of the multilevel SOR algorithm for the five-point Laplacian operator for any number of colors. This general analysis is restricted to a particular coloring pattern. We also give a shorter, restricted proof for a different coloring pattern for which the outer iteration uses blocks rather than strips. We believe that the extension of these results to other operators and coloring patterns is straightforward in many cases. More important, we have found these results to hold true in practical, complicated problems.

Consider an N × N grid, whose nodes are colored by using the following pattern.

\begin{matrix}n & 1 & 3 & 5 & & n − 1 & 2 & 4 \\ n − 1 & 2 & 4 & 6 & ⋯ & n & 1 & 3 \\ n & 1 & 3 & 5 & ⋯ & n − 1 & 2 & 4 \\ n − 1 & 2 & 4 & 6 & & n & 1 & 3 \end{matrix} (2.1)

We index the x-direction by j and the y-direction by k. For convenience we use the notation

(2.2)  S_{j,η} = sin(jηπh)
C_{j,η} = cos(jηπh),

where h = 1/(N + 1). Let x_{j,k} be the discretization of the scalar function x at the point (j, k) on the grid. If the color of the grid point (j, k) is c, we denote this point by x_{j,k}^{(c)} and define the Fourier transform of this point by the equation

x_{j,k}^{(c)} = ∑_{(η, ξ)∈K^{(c)}} x̂_{η,ξ}^{(c)}S_{j,η}S_{k,ξ}. (2.3)

We do not explicitly compute the set K^{(c)} here; we only note that the lowest frequency pair included is (η, ξ) = (1, 1).

Consider the action of the standard five-point Laplacian operator Δ on the scalar function x at grid point (j, k):

− \frac{h^{2}}{4} Δ(j, k)x = x_{j,k} − \frac{1}{4}(x_{j−1,k} + x_{j+1,k} + x_{j,k−1} + x_{j,k+1}). (2.4)

In the frequency domain, we note that the transformed Laplacian operator is approximately block diagonal, with a block Δ̂(η, ξ) for each frequency pair (η, ξ). By “approximately” we mean that for more than four colors, there is a term involving S_{1,η}C_{j,η}S_{k,ξ} that does not cancel in the Fourier representation of Equation 2.4. However, we are interested only in the low-frequency limit. For fixed η, we have that S_{1,η} → 0 as N → ∞. In this limit, which we assume for the rest of this section, these blocks can be expressed as

Δ̂(η, ξ) = [ \begin{matrix}M̂_{n}(η, ξ) & N̂_{n}(η, ξ) \\ N̂_{n}^{T} (η, ξ) & M̂_{n}(η, ξ) \end{matrix} ] . (2.5)

This block operator acts on the vector x̂(η, ξ) given by

x̂(η, ξ) = [ x̂_{η,ξ}^{(1)} ; x̂_{η,ξ}^{(2)} ; ⋮ ; x̂_{η,ξ}^{(n)} ] . (2.6)

Substituting Equation 2.3 into Equation 2.4, we are able to compute the two n/2 × n/2 matrices M̂_{n}(η, ξ) and N̂_{n}(η, ξ)

M̂_{n}(η, ξ) = [ \begin{matrix}1 & α_{η} & β_{ξ} \\ α_{η} & 1 & 0 & β_{ξ} \\ β_{ξ} & 0 & 1 & α_{η} & β_{ξ} \\ 0 & β_{ξ} & α_{η} & 1 & 0 & ... \\ & & β_{ξ} & 0 & ... & & β_{ξ} \\ & & & ... & & 1 & α_{η} \\ & & & & β_{ξ} & α_{η} & 1 \end{matrix} ] , (2.7)

N̂_{n}(η, ξ) = [ \begin{matrix}& & & & & 0 & β_{ξ} \\ & & & & & β_{ξ} & 0 \\ & & & & 0 \\ & & & 0 \\ & & 0 \\ β_{ξ} & 0 \\ 0 & β_{ξ} \end{matrix} ] , (2.8)

where α_{η} = −C_{1,η}/2 and β_{ξ} = −C_{1,ξ}/4.

We consider only the outer level of the multilevel SOR iteration by grouping all grid points colored by the first n/2 colors into a first group, and the remainder of the grid points into a second group. For example, with n = 8, we have the following partition of the grid.

8 1 3 5 7 2 4 6 7 2 4 6 8 1 3 5 8 1 3 5 7 2 4 6 7 2 4 6 8 1 3 5 (2.9)

Note that the grid is partitioned into strips; for n = 4 we obtain the equivalent of the standard line red/black ordering.

This multilevel SOR iteration is consistently ordered. Thus, we can compute [15] the optimal relaxation parameter, ω^{∗}, for SOR from the spectral radius, ρ_{J}, of the Jacobi iteration matrix from the equation

ω^{∗} = \frac{2}{1 + \sqrt{1 − ρ_{J}^{2}}} . (2.10)

The computation of ρ_{J} is more straightforward in the frequency domain. Because of the special structure of M̂_{n} and N̂_{n} and the approximate block diagonal structure of the Jacobi iteration matrix, we find

ρ_{J} ≈ max_{(η,ξ)} ρ(M̂_{n}^{−1}(η, ξ)N̂_{n}(η, ξ)). (2.11)

It is clear that the maximum in Equation 2.11 is obtained for the lowest frequency mode, (η, ξ) = (1, 1), where the block diagonal approximation is most accurate.

To simplify the notation in the following discussion, let M_{n} = M_{n}(1, 1), N_{n} = N_{n}(1, 1), β = β_{1}, and α = α_{1}. For n = 2^{i+1}, we define ρ_{i} = ρ(M̂_{n}^{−1}N̂_{n}). It will be useful to define the following recursion relation:

(2.12)  s_{1} = \frac{cos(πh)/2}{2 − cos(πh)}
r_{1} = s_{1}
s_{i+1} = \frac{s_{i}^{2}}{1 − r_{i}^{2}}
r_{i+1} = r_{i}(1 + s_{i+1}).

Theorem 2.1. Let r_{i} and s_{i} be defined by the recursion relations shown in Equation 2.12. The spectral radius of the Jacobi iteration satisfies ρ_{i} = r_{i} + s_{i}.

Proof: Suppose v is the eigenvector corresponding to the largest eigenvalue of M̂_{n}^{−1}N̂_{n}. By inspection, we note that the product N̂_{n}v has the special form [ x x 0 ⋯ 0 x x]^{T}. Thus, the computation of the spectral radius is equivalent to computing the sum of four components of M̂_{n}^{−1} multiplied by β,

ρ_{i} = |β [(M̂_{n}^{−1})_{1,1} + (M̂_{n}^{−1})_{1,2} + (M̂_{n}^{−1})_{1, \frac{n}{2} −1} + (M̂_{n}^{−1})_{1, \frac{n}{2}} ]| . (2.13)

The value of these components can be computed by an induction argument and use of the Sherman-Morrison-Woodbury formula (A + UV^{T})^{−1} = A^{−1} − A^{−1}U(I + V^{T} A^{−1}U)^{−1}V^{T} A^{−1} . (2.14) Let M̂_{n} = A + UV^{T} , (2.15) where the matrix A is

A = [ \begin{matrix}M̂_{n/2} & 0 \\ 0 & M̂_{n/2} \end{matrix} ] , (2.16)

and the matrices U and V are given by

(2.17)  U = [ \begin{matrix}0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{matrix} ] , V = β [ \begin{matrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix} ] ,

where U and V are \frac{n}{2} × 4 matrices. Here we show the middle 4 rows; the rest of these matrices are zero. Let S be the last two columns of M̂_{n/2}^{−1}, and let S̄ be the first two columns of M̂_{n/2}^{−1}. We note that the matrix M̂_{n/2}^{−1} is persymmetric (i.e., symmetric about the “other” diagonal), so that S̄ is the persymmetric version of S. Also, let K be the leading principal 2 × 2 matrix of βM̂_{n/2}^{−1}. We compute the following matrices,

(2.18)  A^{−1}U = [ \begin{matrix}0 & S \\ S̄ & 0 \end{matrix} ]
A^{−1}V = β [ \begin{matrix}S & 0 \\ 0 & S̄ \end{matrix} ]
V^{T} A^{−1}U = [ \begin{matrix}0 & K \\ K & 0 \end{matrix} ] .

Continuing with our computation, we have

[I + ( \begin{matrix}0 & K \\ K & 0 \end{matrix} )]^{−1} = [ \begin{matrix}K_{e} & −KK_{e} \\ −KK_{e} & K_{e} \end{matrix} ] , (2.19)

where K_{e} = (I − K^{2})^{−1}. (2.20) Combining these facts with Equation 2.14, we have

(A + UV^{T})^{−1} = A^{−1} + β [ \begin{matrix}SKK_{e}S^{T} & −SK_{e}S̄^{T} \\ −S̄K_{e}S^{T} & S̄KK_{e}S̄^{T} \end{matrix} ] . (2.21)

We observe that A is structured so that if we look at A as a block matrix with 2 × 2 blocks, each of these 2 × 2 blocks is symmetric and has equal diagonal elements. Therefore, each block, and in particular K and K_{e}, can be diagonalized when conjugated by the orthogonal matrix

Q = \frac{1}{\sqrt{2}} [ \begin{matrix}1 & 1 \\ −1 & 1 \end{matrix} ] . (2.22)

Thus, to compute ρ_{i} as given in Equation 2.13, we need only compute the largest diagonal elements of the upper left 2 × 2 block and upper right 2 × 2 block of (A + UV^{T})^{−1} when conjugated by Q. We denote these values by r_{i} and s_{i}, respectively. Direct calculation establishes the recursion relation given in Equation 2.12. □

Table 1 A comparison of the predicted and calculated values of ρ_{J} for 8 colors

Mesh Size Prediction Calculated
4 × 4  0.51433  0.50965
8 × 8  0.79573  0.79456
16 × 16  0.93520  0.93511
24 × 24  0.96919  0.96917
32 × 32  0.98213  0.98213

In Table 1, by direct calculation of ρ_{J} we show that our low frequency Fourier representation of the Laplacian is quite accurate. The values for ρ_{J} predicted by Theorem 2.1 converge rapidly to the computed values. We note that Parter [14] obtains an approximation to the dominant eigenvalue ρ_{R} of the k-line Richardson method (equivalent to the Jacobi iteration) for elliptic partial differential equations with constant coefficients. For the Laplace equation with Dirichlet boundary conditions, this approximation is given by

ρ_{R} ≈ 1 − \frac{k}{2}Λh^{2} , (2.23)

where Λ is the minimum eigenvalue of the Laplacian on the given domain. In our case, the domain is the unit square; thus we have Λ = 2π^{2}.

In Figure 3 we compare the estimate obtained by Parter in Equation 2.23 with the Fourier analysis estimate given by Theorem 2.1. The exact values for the spectral radius are indistinguishable from the Fourier values in this figure.

We note that a multilevel algorithm based on these colorings has the disadvantage that the inner iterations will require information exchange along the entire strip. For large problems it is unrealistic to expect that these strips will fit onto one processor. One can avoid interprocessor communication at all levels of the multilevel iteration, except the outer level, by using a “block” coloring. One possible

Fig. 3. Comparison of Parter’s estimate for ρ_{J} with the Fourier estimate for a 64 × 64 grid

block eight-coloring is the following.

1 4 5 8 3 2 7 6 5 8 1 4 7 6 3 2 (2.24)

Let σ_{i} be the spectral radius of the block Jacobi iteration corresponding to the outer level of block coloring with 2^{i+2} colors. We obtain the following interesting result.

Theorem 2.2. We have that asymptotically σ_{1} = ρ_{1}. That is, the spectral radius for the outer level of a block 8-coloring is asymptotically equivalent to line R/B.

Proof: Consider the Fourier representation of the 5-point Laplacian for the coloring given in Equation 2.24:

Δ̂^{B}(η, ξ) = [ \begin{matrix}M̂_{n}^{B}(η, ξ) & N̂_{n}^{B}(η, ξ) \\ N̂_{n}^{B T} (η, ξ) & M̂_{n}^{B}(η, ξ) \end{matrix} ] . (2.25)

We compute that the matrices M̂_{n}^{B} and N̂_{n}^{B} have the form

M̂_{n}^{B}(η, ξ) = [ \begin{matrix}1 & 0 & β_{η} & β_{ξ} \\ 0 & 1 & β_{ξ} & β_{η} \\ β_{η} & β_{ξ} & 1 & 0 \\ β_{ξ} & β_{η} & 0 & 1 \end{matrix} ] , (2.26)

N̂_{n}^{B}(η, ξ) = [ \begin{matrix}0 & 0 & β_{η} & β_{ξ} \\ 0 & 0 & β_{ξ} & β_{η} \\ β_{η} & β_{ξ} & 0 & 0 \\ β_{ξ} & β_{η} & 0 & 0 \end{matrix} ] , (2.27)

where β_{η} = −C_{1,η}/4 and β_{ξ} = −C_{1,ξ}/4.

Consider the computation of the spectral radius of the Jacobi iteration matrix for the lowest frequency mode. We note that conjugation by the block diagonal orthogonal matrix, with the diagonal blocks given by Equation 2.22, yields a system with spectral radius equivalent to that obtained from a line red/black analysis. □

Table 2 A comparison of the calculated values of ρ_{J} for line R/B and the block 8-coloring

Mesh Size Line R/B Block 8-Coloring
4 × 4  0.67928  0.66667
8 × 8  0.88625  0.88515
16 × 16  0.96652  0.96646
24 × 24  0.98435  0.98434

This result is verified by the computational results shown in Table 2. Finally, we note that experimentally it can also be demonstrated that σ_{i} = ρ_{i} asymptotically for i > 1.

Fig. 4. Convergence of different numbers of colors

The improvement in convergence rates as the number of colors increases is not limited to the strip and block colorings. Experimentally, we have found the improvement with every coloring pattern that we have used. For example, if we use the following coloring pattern

\begin{matrix}1 & 2 & ⋯ & \frac{n}{2} & \frac{n}{2} + 1 & ⋯ & n − 1 & n \\ \frac{n}{2} + 1 & \frac{n}{2} + 2 & ⋯ & n & 1 & ⋯ & \frac{n}{2} − 1 & \frac{n}{2} \\ 1 & 2 & ⋯ & \frac{n}{2} & \frac{n}{2} + 1 & ⋯ & n − 1 & n \\ \frac{n}{2} + 1 & \frac{n}{2} + 2 & ⋯ & n & 1 & ⋯ & \frac{n}{2} − 1 & \frac{n}{2} \end{matrix} (2.28)

on the five-point difference operator on a 128 × 128 grid, we observe in Figure 4 that as the number of colors is increased, the convergence of the multilevel SOR algorithm improves.

3. Results for Multilevel Algorithm. In this section we briefly discuss the multilevel SOR algorithm when used with the block colorings described in the second part of §2. Following Kuo and Levy [9], we consider solving the inner subproblems iteratively; we may perform only a few iterations for each subproblem. We are now concerned with the convergence of the multilevel SOR algorithm when the inner subproblems are not solved to full accuracy.

From Lanzkron et al. [10] we know that the multilevel SOR algorithm can be characterized as an inner/outer iteration method. To guarantee convergence for a two-level inner/outer iteration method, the outer iteration must be a regular splitting and the inner iteration must be weak regular. For the model problem the splitting

A = D − L − U = M − N
M = ω^{−1}(D − ωL)^{−1}
(3.29)  N = ω^{−1}(ωU + (1 − ω)D)

induced by the outer iteration is not a regular splitting because N is not positive when ω > 1. However, we expect that if the inner problems are solved to reasonable accuracy, then the method will converge.

Fig. 5. Comparison of R/B SOR with 8-color multilevel SOR with varying numbers of inner iterations

In Figure 5 we show the convergence of the multilevel SOR algorithm with different numbers of inner iterations and compare this convergence with that of the standard red/black SOR algorithm for the 5-pt stencil. We also show that, given a parallel execution time model^{1} similar to that used in Section 1, the multilevel SOR algorithm can have performance advantages over algorithms using fewer numbers

1 The results are from a 16 × 16 grid using 16 processors with model parameters s = 1, τ = 1000, and c = 0.1.

of colors. This advantage is magnified as the ratio of the message startup cost relative to the cost of computation increases.

4. Analysis of Incomplete Factorization. In this section we examine the effect of many color orderings on the condition number of the model system preconditioned by incomplete Cholesky factors. We prove the limited result that a particular four-color ordering results in a smaller condition number than a twocolor, R/B ordering. We have observed experimentally (see Figure 1) that the convergence continues to improve as the number of colors increases beyond four.

We consider a slightly different four-coloring from that used in the preceding section.

2 4 2 4 1 3 1 3 2 4 2 4 1 3 1 3 (4.30)

Consider the discrete Laplacian matrix A ordered according to this coloring. The incomplete factorization A ≈ LL^{T} results in the local operator, L_{j,k}, acting at the grid point (j, k) given by

(4.31) L_{j,k} = \{ \begin{matrix}1, & (j, k) color 1 \\ \sqrt{\frac{7}{8}} − \frac{1}{4}(E_{y} + E_{y}^{−1}), & (j, k) color 2 \\ \sqrt{\frac{7}{8}} − \frac{1}{4}(E_{x} + E_{x}^{−1}), & (j, k) color 3 \\ \sqrt{\frac{5}{7}} − \sqrt{\frac{2}{7}}(E_{x} + E_{x}^{−1} + E_{y} + E_{y}^{−1}), & (j, k) color 4 \end{matrix} ,

where E_{x} and E_{y} are the shift operators in the x and y directions. As observed by Kuo and Chan [8], this representation is accurate only far from the boundary. Therefore, the results derived from equation 4.31 are valid asymptotically.

We note that the Fourier representation of the preconditioned system, L̂^{−1}ÂL̂^{−T}, is block diagonal, using the basis given in equation 2.3. The block corresponding to the frequency pair (η, ξ) is given by

L̂_{ξ,η}^{−1}Â_{ξ,η}L̂_{ξ,η}^{−T} = ( \begin{matrix}1 & 0 & 0 & 0 \\ 0 & \frac{−2(−2+a)(2+a)}{7} & \frac{−2ab}{7} & \frac{2b−8a^{2}b}{7\sqrt{2}\sqrt{5}} \\ 0 & \frac{−2ab}{7} & \frac{−2(−2+b)(2+b)}{7} & \frac{2a−8b^{2}a}{7\sqrt{2}\sqrt{5}} \\ 0 & \frac{2b−8a^{2}b}{7\sqrt{2}\sqrt{5}} & \frac{2a−8b^{2}a}{7\sqrt{2}\sqrt{5}} & \frac{49−12a^{2}−12b^{2}}{35} \end{matrix} ) , (4.32)

where a = cos ξπh and b = cos ηπh.

If λ_{ξ,η} are the eigenvalues of block matrix in equation 4.32, the condition number of the preconditioned system is given by

κ(L^{−1}AL^{−T}) = \begin{matrix}max & ξ,η & | λ_{ξ,η} | \\ min & ξ,η & | λ_{ξ,η} | \end{matrix}. (4.33)

It is evident that the maximum and minimum eigenvalues occur for a = b = 0 and (η, ξ) = (1, 1), respectively. Thus, we find that

max_{ξ,η} | λ_{ξ,η} | = \frac{7}{5}, (4.34)

and

(4.35)  min_{ξ,η} | λ_{ξ,η} | = \frac{1}{70} [89 − 44 cos^{2}(πh) − 16 cos^{4}(πh)]

−\frac{1}{70} [81 + 8 cos^{2}(πh) − 912 cos^{4}(πh) + 1408 cos^{6}(πh) + 256 cos^{8}(πh)]^{\frac{1}{2}} .

We compute that min_{ξ,η} | λ_{ξ,η} | = \frac{56}{29}(πh)^{2} + o(h^{2}). Thus, we obtain the approximation

κ(L^{−1}AL^{−T}) ≈ \frac{29}{40π^{2}h^{2}} . (4.36)

If two colors are used, Kuo and Chan [8] prove that the condition number is asymptotically equal to \frac{1}{π^{2}h^{2}} .

In Table 3 our approximations to the eigenvalues are compared to the calculated eigenvalues for several mesh sizes. We observe that the computed values converge to predicted values with increasing mesh size.

Table 3 Comparison of calculated eigenvalues and predicted eigenvalues for ICCG

Approximate Calculated Approximate Calculated
Mesh Size Maximum λ Maximum λ Minimum λ Minimum λ
4 × 4  1.4000  1.1197  0.7623  0.4435
8 × 8  1.4000  1.3219  0.2347  0.2154
16 × 16  1.4000  1.3782  0.0659  0.0644
32 × 32  1.4000  1.3943  0.0175  0.0174

Experimental evidence similar to the graph in §1 supports the conjecture that increasing the number of colors beyond four will continue to improve the convergence rate of ICCG. However, proving this conjecture is more difficult than for the multilevel SOR algorithm. The reason for this difficulty is that the block diagonal approximation is asymptotically accurate only for the low-frequency analysis, not for the high-frequency analysis. Because we can accurately determine only the smallest eigenvalue of the preconditioner, an exact bound for the condition number when more colors are used is not yet possible.

5. Conclusion. The convergence of iterative methods when used with multicoloring can be improved if the number of colors is increased. For two specific iterative methods, multilevel SOR and ICCG, we have proven that if the number of colors is increased, then the convergence will be faster. This general result is not limited to specific colorings schemes; it appears to apply to many different coloring schemes.

Using experimental models, we have also shown that the overall execution time can be decreased when more colors are used. A trade-off exists between accelerating convergence by using more colors and improving parallel efficiency by using fewer colors. This trade off should be examined when implementing iterative methods using multicoloring.

REFERENCES

[1] L. M. Adams, R. J. LeVeque, and D. M. Young, Analysis of the SOR iteration for the 9-point Laplacian, SIAM Journal on Numerical Analysis, 25 (1988), pp. 1156–1180.

[2] T. F. Chan and H. C. Elman, Fourier analysis of iterative methods for elliptic boundary value problems, SIAM Review, 31 (1989), pp. 20–49.

[3] H. C. Elman and E. Agrón, Ordering techniques for the preconditioned conjugate gradient method on parallel computers, Computer Physics Communications, 53 (1989), pp. 253–269.

[4] L. A. Hageman and D. M. Young, Applied Iterative Methods, Academic Press, New York, 1981.

[5] D. L. Harrar II, Conjugate Gradient Methods for Red/Black Systems on Vector Computers, Ph.D. thesis, Department of Applied Mathematics, University of Virginia, 1989.

[6] M. T. Jones and P. E. Plassmann, Scalable iterative solution of sparse linear systems, Preprint MCS-P277-1191, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., 1991.

[7] M. T. Jones and P. E. Plassmann, A parallel graph coloring heuristic, SIAM Journal on Scientific and Statistical Computing (to appear).

[8] C.-C. J. Kuo and T. F. Chan, Two-color Fourier analysis of iterative algorithms for elliptic problems with red/black ordering, SIAM Journal on Scientific and Statistical Computing, 11 (1990), pp. 767–793.

[9] C.-C. J. Kuo and B. C. Levy, A two-level four-color SOR method, SIAM Journal on Numerical Analysis, 26 (1989), pp. 129–151.

[10] P. J. Lanzkron, D. J. Rose, and D. B. Szyld, Convergence of nested classical iterative methods for linear systems, Numerische Mathematik, 58 (1991), pp. 685–702.

[11] T. A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Mathematics of Computation, 34 (1980), pp. 473–497.

[12] J. Meijerink and H. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Mathematics of Computation, 31 (1977), pp. 148–162.

[13] J. M. Ortega, Introduction to Parallel and Vector Solution of Linear Systems, Plenum, New York, 1988.

[14] S. V. Parter, “Multi-line” iterative methods for elliptic difference equations and fundamental frequencies, Numerische Mathematik, 3 (1961), pp. 305–319.

[15] D. M. Young, Iterative Solution of Large Linear Systems, Academic Press, New York, 1971.