Block Stationary Methods for Nonsymmetric Cyclically Reduced Systems Arising from Three-Dimensional Elliptic Equations.
We consider a three-dimensional convection-diffusion model problem and examine systems of equations arising from performing one step of cyclic reduction on an equally spaced mesh, discretized using the seven-point operator. We present two ordering strategies and analyze block splittings of the resulting matrices. If the matrices are consistently ordered relative to a given partitioning, Young's analysis for the block Gauss--Seidel and block SOR methods can be applied. We compare partitionings for which this property holds with ones where the matrices do not have Property A yet still give rise to an efficient solution process. Bounds on convergence rates are derived and the work involved in solving the systems is estimated.
Introduction
. Consider the three-dimensional (3D) convection-di#usion equation
with constant coe#cients
y, z)
on the unit
subject to Dirichlet-type boundary con-
ditions. We focus on applying seven-point finite di#erence discretizations, for example
centered di#erences to the di#usive terms, and centered di#erences or first-order
upwind approximations to the convective terms. Let us define n and h so that n 3 is
the number of unknowns and is the mesh size, and let F denote the
corresponding di#erence operator, after scaling by h 2 , so that for a gridpoint u i,j,k
not next to the boundary we have
F
If we denote the mesh Reynolds numbers by
-h,
then the values of the components of the computational molecule are given by
# Received by the editors March 3, 1997; accepted for publication (in revised form) by Z. Strako-s
February 27, 1998; published electronically July 9, 1999.
http://www.siam.org/journals/simax/20-4/31771.html
Gates Building, Stanford University, Stanford, CA 94305 (greif@
sccm.stanford.edu).
# Department of Computer Science, University of British Columbia, Vancouver, BC, V6T 1Z4,
Canada (varah@cs.ubc.ca).
(a) lexicographic ordering (b) red/black ordering
Fig. 1.1. Sparsity patterns of the matrices corresponding to two possible orderings of the unknown

if centered di#erence approximations of the first derivatives are used and by
if backward first-order accurate schemes are used.
The sparsity pattern of the underlying matrix depends on the ordering of the
unknowns. In Fig. 1.1 the sparsity patterns associated with two possible ordering
strategies are illustrated. The natural lexicographic ordering in (a) is one where the
unknowns are numbered rowwise and then planewise. The red/black ordering in (b)
means we color the gridpoints using two colors, in a checkerboard fashion, and then
number all the points that correspond to one of the colors first.
As is evident from Fig. 1.1(b), if we split the matrix into four blocks of the same
size, we can see that the two diagonal blocks are diagonal matrices. This means that
the matrix has Property A [24]. A cheap and relatively simple process of elimination
of all the points that correspond to one color (say, red) leads to a smaller system of
equations, whose associated matrix is the Schur complement of the original matrix,
and is still fairly sparse. This procedure amounts to performing one step of cyclic
reduction. Notice that in general both the original and the reduced matrices are
nonsymmetric.
The cyclic reduction step can be repeated until a small system of equations is
obtained, which can then be solved directly. This procedure is called complete cyclic
reduction. It has been studied in several papers, mainly for symmetric systems arising
from two-dimensional (2D) self-adjoint elliptic problems. A general overview of the
algorithm and a list of references can be found in [12]. Early papers that present and
analyze the algorithm are those of Hockney [17], Buneman [2], and Buzbee, Golub,
and Nielson [4]. Buzbee et al. [3] use cyclic reduction for solving the Poisson equation
on irregular regions; Concus and Golub [6] discuss 2D nonseparable cases. Application
of cyclic reduction to matrices with arbitrary dimensions is done by Sweet [20], [21].
presents a fast O(n 2 ) algorithm and discusses its stability and e#ciency.
One step of cyclic reduction for symmetric positive definite systems is analyzed
by Hageman and Varga in [16] and later by Hageman, Luk, and Young [15], where it
is shown that the reduced solver generally converges faster than the unreduced solver.
In [1], Axelsson and Gustafsson use cyclic reduction in conjunction with the conjugate
gradient method. Elman and Golub have conducted an extensive investigation for 2D
elliptic non-self-adjoint problems [8], [9], [10] and have shown that one step of cyclic
reduction leads to systems with several valuable properties, such as symmetrizability
for a large set of the underlying PDE coe#cients, which is e#ectively used to derive
bounds on the convergence rates of iterative solvers, and fast convergence.
Preliminary analysis for the non-self-adjoint 3D model problem (1.1) has been
done by the authors in [14], where one step of 3D cyclic reduction has been described
in detail, and a block Jacobi solver has been analyzed, which is based on a certain
block splitting (referred to as 1D splitting throughout this paper), in conjunction with
what we called a two-plane ordering strategy.
The computational molecule of the reduced operator consists of 19 points, located
on 5 parallel planes. Let R denote the reduced di#erence operator, after scaling by
ah 2 . Then for an interior gridpoint, u i,j,k , we have
R u
-2cf
The following results hold for any ordering strategy. See [14] for the proofs, which
have been obtained by using the techniques of Elman and Golub [8], [9].
Theorem 1.1. The reduced matrix can be symmetrized by a real diagonal similarity
transformation if and only if the products bcde, befg, and cdfg are positive.
Theorem 1.2. If be, cd, fg > 0, then both the reduced matrix and the symmetrized
reduced matrix are diagonally dominant M-matrices.
In this paper our purpose is to extend the analysis initiated in [14] and examine
block stationary methods as solvers for the reduced system. In section 2 we present the
ordering strategies that are examined. In section 3 two block splittings are presented,
and bounds on convergence rates are derived. In section 4 we analyze the reduced
system in the context of consistently ordered matrices. In section 5 the amount of
computational work involved in solving the linear systems is estimated, a comparison
of the reduced system with the unreduced system is conducted, and some numerical
results which validate our analysis and illustrate the fast convergence of the reduced
system are given. Finally, in section 6 we conclude.
2. Orderings for the reduced system. We consider two ordering strategies
for the reduced grid. The two-plane ordering has been described in detail in [14]. It
corresponds to ordering the unknowns by gathering blocks of 2n gridpoints from two
horizontal lines and two adjacent planes. This ordering strategy is depicted in Fig.
2.1(a). In the figure, the numbers are the indices of the gridpoints, which are to be
expressed below by # in (2.1) and (2.15).
The connection between the index of a gridpoint, #, and its coordinate values
(i, j,
h ) is given below. The term fix is borrowed from MATLAB and
means rounding to the nearest integer toward zero.
x
y
17 192123252729316(a) two-plane (b) two-line
Fig. 2.1. Two suggested ordering strategies for the reduced grid (in the figure the unreduced grid
is of size 4 - 4 - 4).
(2.1c)
See [14] for specification of the matrix entries.
An alternative to the two-plane ordering is a straightforward generalization to
three dimensions of the two-line ordering used by Elman and Golub in [9]. It is illustrated
in Fig. 2.1(b). The reduced matrix for this ordering strategy is block pentadiagonal

Each S i,j is (n 2 /2) - (n 2 /2) and is a combination of n
uncoupled matrices, each of
size n - n.
The diagonal matrices {S j,j } are themselves block tridiagonal. Each submatrix is
of size n - n and its diagonal block is
stands for the value along the
main diagonal of the matrix, which is given by a 2
for gridpoints not
next to the boundary. See [14] for specification of the main diagonal's values associated
with gridpoints next to the boundary.
For the superdiagonal and the subdiagonal blocks of the matrices S j,j we have
the following irregular tridiagonal structure, which depends on whether j is even or
odd. The superdiagonal matrices are given by
-2de
-2ce -e 2
-2de .
-2ce -e 2
-2de
1042 CHEN GREIF AND JAMES VARAH
if j is odd or
-2ce -e 2
-2de
-2ce -e 2
if j is even.
The subdiagonal matrices are
-2bd
if j is odd or
-2bd
-2bd .
-2bd
if j is even.
The superdiagonal and the subdiagonal blocks of S, S j,j-1 are block tridiagonal:
(a) two-plane (b) two-line
Fig. 2.2. Sparsity patterns of the reduced matrices associated with the two ordering strategies
(the matrices correspond to 6 - 6 - 6 grids). Each square corresponds to an n 2
gridpoints form two coupled planes in the reduced grid).
Finally, the matrices S j,j-2 and S j,j+2 are diagonal:
2.
The connection between the gridpoint's index and its coordinate values is given
by
(2.15c)
The sparsity patterns of the matrices corresponding to two-plane ordering and
two-line ordering are depicted in Fig. 2.2.
(a) 1D splitting (b) 2D splitting
Fig. 3.1. Sparsity patterns of the block diagonal matrices associated with the block Jacobi split-
ting, for the two suggested block splittings, using the two-plane ordering strategy.
3. Block splittings and bounds on convergence rate. For the two ordering
strategies presented in section 2 the matrices can be expressed as block tridiagonal,
of the form
S is an (n 3 /2) - (n 3 /2) matrix. In the case of two-plane ordering, each block S i,j is
of size n 2
and is block tridiagonal with respect to 2n - 2n blocks. In the case of
two-line ordering, each block S i,j is of size (n 2 /2) - (n 2 /2) and is block tridiagonal
with respect to n - n blocks.
In solving the reduced system using a stationary method, various splittings are
possible. We consider two obvious ones, based on dimension. We use the term 1D
splitting for a splitting which is based on partitioning the matrix into O(n) blocks
(2n-2n blocks for the two-plane ordering and n-n blocks for the two-line ordering).
A 2D splitting is one which is based on partitioning the matrix into O(n 2 ) blocks
blocks for the two-plane ordering and (n 2 /2) - (n 2 /2) blocks for the two-line
ordering). Notice that the 1D splitting for both ordering strategies is essentially
associated with blocks of gridpoints that are x-oriented. However, the 2D splitting
for the two-line ordering corresponds to x-y oriented planes of gridpoints, whereas
for the two-plane ordering it corresponds to x-z oriented planes of gridpoints. (These
observations can be deduced by referring to Fig. 2.1.) Di#erent orientations can be
obtained by simply reordering the unknowns so that the roles of x, y, and z are
interchanged.
The sparsity patterns of the block diagonal parts of the splittings associated with
the block Jacobi scheme are depicted in Fig. 3.1.
We now compare the orderings. We have the following useful result.
Theorem 3.1. If be, cd, fg > 0, then for the 1D splitting the spectral radius of the
Jacobi iteration matrix associated with two-plane ordering is smaller than the spectral
radius of the iteration matrix associated with the two-line ordering.
Proof. By Theorem 1.2 the matrices are M-matrices. Each ordering strategy
produces a matrix which is merely a symmetric permutation of a matrix associated
with the other ordering. Suppose splitting of the two-plane
ordering matrix and S is a 1D splitting for the two-line order-
ing. There exists a permutation matrix P such that P T Consider the
splitting P T . It is straightforward to show by examining
the matrix entries that P T . The latter are both nonnegative matrices;
therefore by [23, Thm. 3.15] it follows that 0 < #(M
The same result applies to 2D splitting, provided that the orientation of the planes
of gridpoints is identical for both ordering strategies. The proof for this is identical
to the proof of Theorem 3.1.
The results indicated in Theorem 3.1 can be observed in Fig. 3.2. It is interesting
to observe that the superiority of the two-plane ordering carries over to the case
be, cd, fg < 0, which corresponds to the region of mesh Reynolds numbers larger
than 1 (for which the PDE is considered convection-dominated). We remark, however,
that the amount of computational work per each iteration is somewhat higher for the
system which corresponds to two-plane ordering. In Fig. 3.2 a few cross sections of
mesh Reynolds numbers are examined. For example, graph (a) corresponds to flow
with the same velocity in x, y, and z directions. Graph (b) corresponds to flow only
in x and y directions, and no convection in z direction, and so on. (See (1.3) for
definitions of #, and #.)
We now derive bounds on convergence rates. Below we shall attach the subscripts
1 and 2 to matrices associated with the 1D splitting and 2D splitting, respectively.
Since two-plane ordering gives rise to a more-e#cient solution procedure than two-line
ordering, we focus on it.
Denote the two splittings for the block Jacobi scheme by
In [14] we have shown that if be, cd, and fg have the same sign then a real diagonal
nonsingular symmetrizer can be found, and thus (since the symmetrizer is diagonal)
the sparsity patterns of the original nonsymmetric matrix and the symmetrized matrix
are identical. Let us attach the hat sign to a matrix to denote application of the
similarity transformation that symmetrizes it. That is, for a given matrix X and a
diagonal symmetrizer Q, Q -1 XQ is to be denoted by -
X.
The matrices -
are similar to the original iteration matrices
respectively, and thus have the same spectral radii. Following
Elman and Golub's strategy [8], [9], the symmetric matrix can be handled more easily
as far as computing the spectral radius is concerned, since we can use the following:
2.
The results presented below are for the case be, cd, fg > 0, using two-plane
ordering. These conditions are equivalent to |#| < 1 if centered di#erences
are used to discretize the convective terms. No restriction on the magnitude of the
mesh Reynolds numbers is imposed if upwind di#erences are used. For these values
tight bounds for the spectral radius of the iteration matrix can be obtained.
For -
D 1 the minimal eigenvalue has been found in [14, Thm. 3.8], and the relevant
part of this theorem is quoted below.
Proposition 3.2. The minimal eigenvalue of -
D 1 is
A lower bound for -
D 1 is given by the following proposition.
Fig. 3.2. Spectral radii of iteration matrices versus mesh Reynolds numbers for the block Jacobi
scheme, using 1D splitting and centered di#erence discretization. The broken lines correspond to
two-plane ordering. The solid lines correspond to two-line ordering.
Proposition 3.3. The minimal eigenvalue of -
D 1 is bounded from below by
-#, where

Table
Comparison between the computed spectral radius and the bound, for the 2D splitting, with
Scheme Upwind Centered
n # bound ratio # bound ratio
The proof for this part follows from [14, Lem. 3.10], where it is shown that the
spectral radius of -
D 1 is bounded by #.
Combining Propositions 3.2 and 3.3, and applying Rayleigh quotients to the matrices
we obtain the following lemma.
Lemma 3.4. The minimal eigenvalue of -
D 2 is bounded from below by #, where
# and # are the expressions given in (3.3) and (3.4).
The bound for -
C 2 can be obtained by combining [14, Lems. 3.11-3.13], as follows.
Lemma 3.5. The spectral radius of the matrix -
C 2 is bounded by
# .
Finally, Lemmas 3.4 and 3.5 lead to the following theorem.
Theorem 3.6. The spectral radii of the iteration matrices D
are bounded by #
# and #
respectively, where #, and # are defined in (3.3),
(3.4), and (3.5), respectively.
Corollary 3.7. If be, cd, fg > 0 then the block Jacobi iteration converges for
both the 1D and 2D splittings.
Proof. For this we can use Varga's result on M-matrices [23, Thm. 3.13]. Alter-
natively, Taylor expansions of the bounds given in Theorem 3.6 are given by
and
and thus are smaller than 1.
In

Table

3.1 we give some indication on the quality of the bound for the 2D
splitting. Results with a similar level of accuracy have been obtained and presented
in [14] for the 1D splitting. As can be observed, the bounds are tight and become
tighter as n increases, which suggests that they are asymptotic to the spectral radii.
We now discuss other stationary methods, namely, Gauss-Seidel and SOR. Relative
to a given partitioning, if the reduced matrix is consistently ordered, then it
is straightforward to apply Young's analysis, and the bounds in Theorem 3.6 can be
used for estimating the rate of convergence of the Gauss-Seidel and SOR schemes.
The reader is referred to [19, Defs. 4.3 and 4.4] for definitions of Property A and consistent
ordering. As stated in [19], a matrix that is consistently ordered has Property
1048 CHEN GREIF AND JAMES VARAH
conversely, a matrix with Property A can be permuted so that it is consistently
ordered. We mentioned in the introduction that the matrix of the unreduced system
has Property A. For the reduced system, we have the following observations.
Proposition 3.8. The reduced matrix associated with two-line ordering, SL , does
not have Property A relative to 1D or 2D partitionings.
Proof. Let S i,j denote the (i, j)th n-n block of SL , and let Q be an (n 2 /2)-(n 2 /2)
matrix, whose entries satisfy q be an
the (i, j)th (n 2 /2) - (n 2 /2) block submatrix of
S is nonzero, and t is a pentadiagonal matrix and thus
does not have Property A. Since Q can be referred to as a partitioning of T into
also does not have Property A.
Proposition 3.9. The reduced matrix associated with two-plane ordering, SP ,
does not have Property A relative to 1D partitioning.
Proof. Let S i,j denote the (i, j)th 2n - 2n block of SP , and let Q be an (n 2
whose entries satisfy q otherwise. It is
straightforward to see that the nonzero pattern of Q is identical to that of the matrix
associated with using a nine-point operator for a 2D grid. Since the latter does not
have Property A relative to partitioning into 1 - 1 matrices, the result follows.
On the other hand, we have the following proposition.
Proposition 3.10. The reduced matrix associated with two-plane ordering, SP ,
has Property A and, moreover, is consistently ordered relative to 2D partitioning.
Proof. The matrix is block tridiagonal relative to this partitioning [24].
For the SOR scheme we have the following result, which is completely analogous
to Elman and Golub's result for the 2D problem [9, Thm. 4].
Theorem 3.11. Let L# denote the block SOR operator associated with 2D splitting
and using two-plane ordering. If either be, cd, fg > 0 or < 0, then the choice
minimizes #(L# ) with respect to #, and #(L# - 1.
The proof of this theorem is essentially identical to the proof of Elman and Golub
in [9, Thm. 4] and follows from Young [24, Chap. 14, Sects. 5.2 and 14.3]. The algebraic
details on how to pick the signs of the diagonal symmetrizer so that the symmetrized
block diagonal part of the splitting is a diagonally dominant M-matrix are omitted.
That #(D -1
known by Corollary 3.7. The reduced matrix is consistently
ordered by Proposition 3.10.
A way to approximately determine an optimal relaxation parameter for the case
be, cd, fg > 0 is to replace #(D by the bound for it (given in Theorem 3.6) in
the expression for # in Theorem 3.11. If the bound for the block Jacobi scheme is
tight, then the estimate of # is fairly accurate.
Proposition 3.12. Suppose be, cd, fg > 0. For the system associated with 2D
splitting and for h su#ciently small, the choice
approximately minimizes #(L# ). The spectral radius of the iteration matrix is approximately
- 1.
The Taylor expansion of the estimate for the optimal relaxation parameter is
given by
Fig. 4.1. The sparsity pattern of the matrix C d .
From (3.9) it follows that the estimated asymptotic rate of convergence of the block
SOR scheme is approximately the second term in (3.9) (with the negative sign re-
moved) and is thus O(h).
4. Near-Property A for 1D splitting of the two-plane matrix. Although
the matrix associated with two-plane ordering does not have Property A relative to
the 1D partitioning, some interesting observations can be made: As before, let {S i,j }
denote the n 2
blocks of the reduced matrix. Each block S i,j is a block tridiagonal
matrix relative to 2n - 2n blocks. We attach superscripts to mark how far a block
diagonal is from the main block diagonal, and we define
See [14] for specification of the entries of these matrices. As in section 3 (with a slight
change in notation), let C be the 1D splitting of the matrix, and define -
so that
C +C d ). The matrix D-
C has Property A, but SP does not. Let
us examine the matrix that prevents SP from having Property A, namely, C d . It is
an extremely sparse matrix, and the magnitude of the nonzero values in this matrix
is bounded by 2 if be, cd, fg > 0. The nonzero pattern of C d is depicted in Fig. 4.1.
We wish to estimate how far the reduced matrix SP is from having block Property
A, relative to the 1D partitioning. Let us denote the upper part and the lower part of
C d by U d and L d , respectively, and let -
U and -
L be the upper part and lower part of
C, respectively. Then the spectral radius of the block Gauss-Seidel matrix satisfies
significantly larger than the other norms in the above inequality,
which means that the spectral radius of the Gauss-Seidel iteration matrix associated
with the two-plane ordering can be estimated by replacing the two-plane matrix by
C, which does have Property A and thus is easier to analyze. Alternatively, the
following observation has been obtained by numerical experiments:
(a) #GS vs # 2
centered (b) #GS vs # 2
a
Fig. 4.2. "Near Property A" for the 1D splitting.
Young's analysis can be applied directly to both D -
C and D - C d (both have
Property A), and thus an approximate relationship between the eigenvalues of the
block Jacobi iteration matrix and the eigenvalues of the block Gauss-Seidel iteration
matrix can be obtained.
For be, cd, fg > 0 we have observed that the spectral radius of the block Jacobi
iteration matrix satisfies
J #GS .
The first two graphs in Fig. 4.2 illustrate this phenomenon numerically. The broken
lines in graphs (a) and (b) correspond to the square of the spectral radius of the
iteration matrix associated with block Jacobi, for a 256 - 256 matrix. The solid lines
correspond to the spectral radius of the block Gauss-Seidel iteration matrix. As can
be seen, the curves are almost indistinguishable. This phenomenon becomes more
dramatic as the systems become larger.
Some analysis can be done using Varga's work on extensions of the theory of p-
cyclic matrices [22], [23, Sect. 4.4]. (In this paper we are concerned only with
Recall [23, Def. 4.2], which defines a set S of matrices as follows. The square matrix
satisfies the following properties:
diagonal entries.
2. B is irreducible and convergent, i.e., 0 < #(B) < 1.
3. B is symmetric.
If be, cd, fg > 0, the reduced matrix SP is a diagonally dominant M-matrix
which can be symmetrized, and -
D 1/2 is well defined. Define -
I -
D -1/2 .
Applying block Jacobi to the original reduced system is analogous to applying
point Jacobi to -
S, in the sense that the spectra of the iteration matrices associated
with both systems are identical. The iteration matrix associated with -
S is
D -1/2 . Showing that the matrix B belongs to the set S defined above is
easy and is omitted. Let L be the lower part of B. Define MB (#L
and mB (#(MB (#)). Let
#(B)# 1/2
Then we have [23, Thm. 4.7] (with a slight modification so as to match the
terminology used in this paper), as follows.
Theorem 4.1. Let B # S. Then, hB (# 1 if and only if B is consistently
ordered.
In some sense hB (ln #) measures the departure of the matrix B from having block
Property A. For matrices that are not consistently ordered, the following result applies
[23, Thm. 4.6].
Theorem 4.2. If B # S, then either hB (# 1 for all real #, or hB (#) is strictly
increasing for # 0. Moreover, for any #= 0,

Figure

4.2(c) demonstrates how close the function hB is to 1 for the reduced
matrix when 1D partitioning is used and provides another way to illustrate the near-
Property A of the matrix. In the figure, the function hB is computed for a symmetrized
block Jacobi 256 - 256 matrix, where
We can now analyze the Gauss-Seidel and SOR schemes. Recall [23, Thm. 4.8]
(slightly modified), as follows.
Theorem 4.3. Let LB,# denote the SOR iteration matrix. If B # S then the
Gauss-Seidel iteration matrix, which corresponds to the case
with equality possible only if B is consistently ordered.
This is a sharpened form of the Stein-Rosenberg theorem [23]. Applying this
theorem to our reduced matrix, we have the following theorem.
Theorem 4.4. If the bound for the block Jacobi iteration matrix tends to the
actual spectral radius as h # 0, then the spectral radius of the block Gauss-Seidel
iteration matrix coincides with the square of the bound for the spectral radius of the
block Jacobi iteration matrix up to O(h 2 ) terms.
Proof. Since the iteration matrix B has the same spectral radius as D -1 C, where
D-C, we can use the bound for the 1D iteration matrix, which was presented
in Theorem 3.6. For simplicity of notation, denote this bound by #. Clearly, since
1052 CHEN GREIF AND JAMES VARAH
20.10.30.50.70.9Fig. 4.3. Spectral radius of the SOR iteration matrix versus the relaxation parameter. The
uppermost curve corresponds to 1D splitting for the unreduced system, and then we have, in order,
2D splitting for the unreduced system, 1D splitting for the reduced system, and 2D splitting for the
reduced system.
Since # has a Taylor expansion of the form 1 - ch 2
that #
2-# and # 2 have the same Taylor expansion up to O(h 2 ) terms, of the form
in terms of the PDE coe#cients,
and the same for # 2 . It has been shown that the bound for #(B) is extremely tight as
so we can replace the spectral radii by the bounds for the spectral radii
in Theorem 4.3 to obtain the desired result.
The actual meaning of this result is that for systems of equations that are large
enough, the matrix nearly has Property A relative to 1D partitioning, at least as far
as the convergence properties of the block Gauss-Seidel scheme are concerned. Since
the solution process for small mesh Reynolds numbers is more e#cient for the 1D
splitting, compared to the 2D splitting, as we shall see in section 5, it was our aim to
overcome the di#culty of not being able to apply Young's analysis directly.
For the block SOR scheme, the upper bound for the spectral radius is given in
[23, Thm. 4.9] as # - 1 and is not tight. However, it is numerically evident that
the bound for the Jacobi iteration matrix can be e#ectively used to estimate the
optimal SOR parameter. In Fig. 4.3 we can observe that the behavior for the 1D
splitting is qualitatively identical to the behavior of two-cyclic consistently ordered
matrices. Here we present results for centered di#erence discretization of the problem
with 0.5. The reduced matrix is 256 - 256. In the figure we also present
the behavior of the SOR iteration matrix of the unreduced system.
5. Computational work and numerical experiments. Having done some
analysis, in this section we examine which of the 1D and 2D solvers is more e#cient
overall and show that the reduced system is superior to the unreduced system.
5.1. Aspects of computational work. If be, cd, fg > 0, then by [23, Thm.
3.15] or by (3.6) and (3.7), it is evident that the spectral radius of the iteration matrix
associated with the 2D splitting is smaller than that of the 1D iteration matrix.
However, inverting D 1 involves less computational work than inverting D 2 . We now
compare these two solution procedures.
We begin with the block Jacobi scheme. Asymptotically, there is a fixed ratio
of 1.8 between the rate of convergence of the two splittings (see (3.6) and (3.7)). In
rough terms, this number characterizes the ratio between number of iterations until
convergence for the two solvers.
As far as the computational work per iteration is concerned, if D
are the LU decompositions of the matrices of the systems that are to be
solved in each iteration, we can assume that the number of operations per iteration is
approximately the number of nonzeros in L i +U i plus the number of nonzeros in the
other part of the splitting. In order to avoid costly fill-in using Gaussian elimination
for (whose band is sparse), we use instead a technique of inner-outer iterations.
denote the number of iterations for the schemes associated with
the 1D splitting and the 2D splitting, respectively. Let us also define cost functions
as follows: c 1 (n) and c 2 (n) represent the overall number of floating point operations
for each of the solvers, and c in (n) represents the cost of the inner solve. Then
(5.1a)
The term nz(X) stands for the number of nonzeros of a matrix X, and S stands for
the reduced matrix.
Proposition 5.1. For n large enough, the scheme associated with the 2D splitting
is cheaper than the one associated with the 1D splitting only if c in (n) < 15n 3 .
Proof. If n is large enough we can use the relation k1
refer only to the
leading power of n in the expressions for c 1 (n) and c 2 (n). So doing, it follows that
c in (n)
and the result stated in the proposition readily follows.
What is left now is to examine the amount of work involved in solving the inner
system of equations. A natural choice of a splitting for this system is D
It is straightforward to show the following by Propositions 3.2 and 3.3.
Proposition 5.2. If block Jacobi based on the splitting D
is used, then the spectral radius of the inner iteration matrix, namely, I -D
bounded by #
are defined in (3.3) and (3.4).
For considering methods that are faster than block Jacobi for the inner system,
we have the following useful result.
Proposition 5.3. The inner matrix is block consistently ordered relative to 1D
partitioning.
Proof. The inner matrix is block tridiagonal relative to this partitioning.
We are now ready to prove the main result of this subsection.
Proposition 5.4. If be, cd, fg > 0, then if 1D splitting is used in solving the
inner system, the cost of solving it is higher than 15n 3 floating point operations, for
block Jacobi as well as block Gauss-Seidel and block SOR, and thus, for n large enough
and the methods considered in this paper, the 1D solver is faster than the 2D solver.
1054 CHEN GREIF AND JAMES VARAH
Proof. The Taylor expansion of the bound in Proposition 5.2 is
For h small enough, we can simply examine the leading term: The bound is approximately9 if block Jacobi is used, and since by Proposition 5.3 the matrix is consistently
ordered, Young's analysis shows that the spectral radius is approximately
81 if block
Gauss-Seidel is used and approximately 0.055 if block SOR with the optimal relaxation
parameter is used. For both of these schemes each iteration costs about 7n 3
floating point operations. Since reducing the inital error by a factor of 10 m takes
roughly is the spectral radius of the associated iteration
matrix, it follows that even for the block SOR scheme with the optimal relaxation
parameter, which is the fastest scheme considered here, after two iterations the error
is reduced only by a factor of approximately 10 2.5 , which is obviously far from satis-
factory. Thus the iteration count is larger than 2, and the cost of inner solve is larger
than 15n 3 floating point operations.
We remark that an inexact inner solve can also be considered (see, for example,
Elman and Golub's paper on inexact Uzawa algorithms [11]), but this is beyond the
scope of this work.
It is our conclusion that the solver associated with 1D splitting is more e#cient
than the one associated with the 2D splitting if upwind di#erences are used or if
centered di#erences with mesh Reynolds numbers smaller than 1 in magnitude are
used.
5.2. Comparison with the unreduced system. One step of cyclic reduction
results in a more complicated di#erence operator compared to the original, unreduced
system, and a grid which is more di#cult to handle as far as ordering of the unknowns
is concerned. Moreover, the unreduced matrix is block consistently ordered relative
to both 1D and 2D splittings (we refer to the straightforward one-line and one-plane
partitionings as the basis for 1D and 2D splittings in case of the unreduced system)
and thus Young's analysis can be easily applied. One could ask, therefore, what the
advantages of using cyclic reduction are. In this subsection we illustrate the superiority
of the reduced system over the unreduced system.
We start with the block Jacobi scheme. For the unreduced system we shall refer
to natural lexicographic ordering of the unknowns, so that the lines of gridpoints are
x-oriented and the planes are x-y oriented. We start with quoting the following result,
given in [14, Sects. 2 and 4].
Lemma 5.5. The spectral radius of the block Jacobi scheme associated with the
1D splitting for the unreduced system is
The Taylor expansion of (5.4) about
In [14] we have shown that the spectrum of the iteration matrix of the unreduced
system can be found by a sequence of diagonalizations and permutations that form
a similarity transformation of the matrix into a matrix whose associated iteration
matrix is easy to analyze, as far as its spectrum is concerned. The reader is referred
to the proof of [14, Thm. 2.1] for full details. For the 2D splitting a similar procedure
can be applied. The technique we have used is similar to the one presented in [14] and
the algebraic details are omitted.
Lemma 5.6. The spectral radius of the block Jacobi iteration matrix associated
with 2D splitting is given by
and its Taylor expansion about
The same type of analysis that has been done in the previous section, comparing the
1D splitting to the 2D splitting for the reduced system, is possible for the unreduced
system. Below we sketch the main details: Suppose inner-outer iterations are used in
solving the scheme associated with the 2D splitting. Denote, again, this splitting for
the inner system as D are now di#erent than the ones
defined in section 3). Then we have the following proposition.
Proposition 5.7. Consider the unreduced system. Suppose be, cd, fg > 0, n is
su#ciently large, and 1D splitting is used in solving the inner system. Then for the
stationary methods considered in this paper, the 1D solver is faster than the 2D solver.
Proof. The ratio between the asymptotic rate of convergence between the 1D solver
and the 2D solver is 2. The number of nonzeros of the whole matrix is approximately
7n 3 , the number of nonzeros of D 1 is approximately 3n 3 , and the number of nonzeros
of D 2 is approximately 5n 3 . Since the spectral radii for the two splittings are available,
we can find the spectral radius for the iteration matrix of the inner system. Its Taylor
expansion is given by 1
cost functions
analogous to the ones defined in section 3 for the reduced system, and using the same
line of argument, we have
c in (n)
and from this it follows that only if c in (n) < 12n 3 the 2D solver is more e#cient.
However, as in Proposition 5.4, this means at most two iterations of the inner solve
can be performed, which is not enough for the required accuracy.
Since the 1D splitting for both the reduced and the unreduced systems gives rise
to a more e#cient solve, we compare these two systems, focusing on this splitting.
See also [14, Sect. 4]. The LU decomposition for the solution of the system in each
iteration is done once and for all (see [12] for operation count) and its cost is negligible
in comparison with the amount of work done in the iterative process.
Each iteration in the reduced system costs about 10n 3 floating point operations,
whereas each iterate for the unreduced system costs approximately 7n 3 floating point
operations per iteration. Hence, the amount of computational work per iteration is
cheaper for the unreduced system by a factor of about 10/7. However, using the
asymptotic formulas (3.6) and (5.5), it is evident that the number of iterations required
for the unreduced system is larger than that required for the reduced system, and in
(a) # gs - centered (b) # gs - upwind
Fig. 5.1. Comparison between the spectral radii of the Gauss-Seidel iteration matrices of the
reduced and unreduced systems. The uppermost curve corresponds to 1D splitting for the unreduced
system, and then we have, in order, 2D splitting for the unreduced system, 1D splitting for the
reduced system, and 2D splitting for the reduced system.
the worst case, the ratio between the work required for solving the reduced system
versus the unreduced system is roughly (10/7) - (27/40), which is 27/28 and is still
smaller than 1, thus the reduced solver is more e#cient. If the convective terms are
nonzero, then this ratio becomes smaller, and in practice we have observed substantial
savings, as is illustrated in the test problem discussed in section 5.3.
Moving from comparing the block Jacobi scheme for both the reduced and the
unreduced systems to comparing Gauss-Seidel and SOR is straightforward if Young's
analysis can be used. In section 4 we showed that even though the reduced matrix is
not consistently ordered relative to 1D partitioning, it is nearly consistently ordered.
In general, convergence analysis for the Jacobi scheme does not always indicate the behavior
of the Gauss-Seidel and SOR schemes. Nevertheless, for two-cyclic consistently
ordered matrices (or matrices that are nearly so) the strong connections between the
spectra of the Jacobi iteration matrix and the Gauss-Seidel and SOR iteration matrices
[24] allow us to conclude that once the superiority of the reduced system over the
unreduced system has been shown for Jacobi, this superiority is carried over to the
other stationary schemes. Indeed, our numerical experiments verify this observation,
as is illustrated in section 5.3.
In Fig. 5.1 the superiority of the reduced system over the unreduced system for the
Gauss-Seidel scheme is illustrated numerically. The graphs were created for a small
512-point grid. It is interesting to notice that the reduced 1D Gauss-Seidel iteration
matrix is well behaved (i.e., its spectral radius is significantly smaller than 1), even
for the convection-dominated case, when centered di#erences are used. Convergence
does not occur when the block Jacobi scheme with the same values of mesh Reynolds
numbers is used. We have no bounds on convergence rates for this range of mesh
Reynolds numbers and thus cannot explain this phenomenon analytically.
The superiority of the reduced system is evident also for the SOR scheme (see Fig.
4.3). Notice that for the SOR scheme it is di#cult to determine the optimal relaxation
parameter when be, cd, and fg are negative.
We end this subsection with a remark regarding the case of convection-dominated
equations. Our convergence analysis does not cover the case of mesh Reynolds num-

Table
Comparison between iteration counts for the reduced and unreduced system, for di#erent values
of mesh Reynolds numbers. N/C marks no convergence after 2,000 iterations.
System Reduced Unreduced
centered 393 173 53 N/C 1030 444 N/C N/C
GS centered 188 77 14 322 492 198 N/C N/C
centered 36
GS upwind 219
bers that are greater than 1 in magnitude in conjunction with centered di#erence
discretization. Since the numerical solution might be oscillatory when a centered difference
scheme is used [18], analysis for this case is of less interest. Nevertheless,
Fourier analysis based on Chan and Elman's technique [5], which shows that when
one of the mesh Reynolds numbers tends to # the scheme still converges, is presented
in [13].
5.3. Test problem. Consider (1.1), where the right-hand side is such that the
solution for the continuous problem is u(x, y, sin(#z) and the
domain is the unit cube. The Dirichlet boundary conditions in this case are zero. The
performance of the solvers for this specific problem well represents the performance
for other test problems that we have examined.
We have taken the zero vector as our initial guess and have used ||r i
as a stopping criterion (here r i denotes the residual at the ith iterate). The
program stopped if the stopping criterion was not satisfied after 2,000 iterations. Our
numerical experiments were executed on an SGI Origin 2000, which has four parallel
195 MHZ processors, 512 MB RAM, and 4MB cache. The program was written in
MATLAB 5.
In the experiments that are presented, the 1D solver is used. In Table 5.1, the
grid is of size 32. The matrix of the underlying system of equations is of
size 32, 768 - 32, 768. In the table, iteration counts for the Jacobi scheme and the
Gauss-Seidel scheme are presented for four values of the PDE coe#cients and for two
discretization schemes.
The PDE coe#cients referred to in Table 5.1 are specified in (1.1). For the values
of these coe#cients in the table, namely 10, 20, 100, and 1,000, the corresponding
values of the mesh Reynolds numbers are 0.1515, 0.3030, 1.515, and 15.15. Notice that
the last two are larger than 1, and so for these values we have no analytical way of
knowing the optimal relaxation parameter and the experiments for these values were
not performed.
The following observations can be made.
1. Overall the reduced solver is substantially faster than the unreduced solver.
There are cases where the reduced solver converges whereas the unreduced
solver does not. We remark that in all cases that were examined, the CPU
time for the reduced solver was less (much less in most cases) than the CPU
time for the unreduced system.
2. For convergence is faster than for This illustrates a phe-
nomenon, which is supported by the analysis and holds also for the two-dimensional
case [9], that for small-enough mesh Reynolds numbers, the
1058 CHEN GREIF AND JAMES VARAH
"more nonsymmetric" systems converge faster than the "close to symmet-
ric" ones (close in the sense of PDE coe#cients close to zero).
3. The upwind di#erence scheme converges more slowly than the centered di#er-
ence scheme when the mesh Reynolds numbers are small in magnitude, but
convergence is extremely fast for large mesh Reynolds numbers. This applies
to both the reduced and the unreduced systems and follows from the fact that
as the PDE coe#cients grow larger, the underlying matrix is more diagonally
dominant when upwind schemes are used.
6. Concluding remarks. We have presented ordering strategies for a cyclically
reduced matrix arising from discretizing a 3D model problem with constant coe#-
cients. We have derived bounds on convergence rates for block stationary schemes
associated with what we called 1D splitting or 2D splitting. We have compared the
amount of work involved in solving the system with the suggested splittings. In gen-
eral, the 1D splitting gives rise to more-e#cient solvers. Since the matrices associated
with this splitting are not consistently ordered, we have analyzed their departure from
block Property A and have shown that, in fact, these matrices are nearly block consistently
ordered. We have shown, both analytically and numerically, that one step of
cyclic reduction results in a system which is easier to solve, compared to the original,
unreduced system.

Acknowledgments

. We would like to thank the referees for their helpful com-
ments, which substantially improved this manuscript.



--R

On the use of preconditioned conjugate gradient methods for red-black order five-point di#erence schemes
A compact non-iterative Poisson solver
The direct solution of the discrete Poisson equation on irregular regions
On direct methods for solving Poisson's equations
Fourier analysis of iterative methods for elliptic problems
Use of fast direct methods for the e
Point cyclic reductions for elliptic boundary-value problems I: The constant coefficient case
Iterative methods for cyclically reduced non-self-adjoint linear systems
Iterative methods for cyclically reduced non-self-adjoint linear systems II
Line iterative methods for cyclically reduced discrete convection-di#usion problems
Inexact and preconditioned Uzawa algorithms for saddle point problems
Matrix Computations
Analysis of Cyclic Reduction for the Numerical Solution of Three-Dimensional Convection-Di#usion Equations
Iterative solution of cyclically reduced systems arising from discretization of the three-dimensional convection di#usion equation
On the equivalence of certain iterative acceleration methods
Block iterative methods for cyclically reduced matrix equa- tions
A fast direct solution of Poisson's equation using Fourier analysis
Numerical Solution of Convection-Di#usion Problems
Iterative Methods for Sparse Linear Systems
A generalized cyclic reduction algorithm
A cyclic reduction algorithm for solving block tridiagonal systems of arbitrary dimension
A generalization of the Young-Frankel successive over-relaxation scheme
Matrix Iterative Analysis
Iterative Solution of Large Linear Systems
--TR

--CTR
M. Cheung , Michael K. Ng, Block-circulant preconditioners for systems arising from discretization of the three-dimensional convection-diffusion equation, Journal of Computational and Applied Mathematics, v.140 n.1-2, p.143-158, 1 March 2002
Liang Li , Ting-Zhu Huang , Xing-Ping Liu, Asymmetric Hermitian and skew-Hermitian splitting methods for positive definite linear systems, Computers & Mathematics with Applications, v.54 n.1, p.147-159, July, 2007
convection-diffusion;cyclic reduction;stationary methods;three-dimensional problems