Reduced Systems for Three-Dimensional Elliptic Equations with Variable Coefficients.
We consider large sparse nonsymmetric linear  systems arising from finite difference discretization of three-dimensional (3D) convection-diffusion equations with variable coefficients.  We show that performing one step of cyclic reduction yields a system of equations which is well conditioned and for which  fast convergence can be obtained.  A certain block ordering strategy is applied,  and analytical results concerning symmetrizability  conditions and bounds on convergence rates are given.  The analysis is accompanied by  numerical examples.
Introduction
. Consider the following three-dimensional (3D) convection-
di#usion equation
on a
domain# R 3 , subject to Dirichlet, Neumann, or mixed boundary condi-
tions, where all the functions in (1.1) are trivariate, and p, q, r > 0 on # Several
discretization schemes are possible. See Morton [11] for a comprehensive survey on
numerical solution of the convection-di#usion problem. In this work we use a seven-point
discretization technique as a starting point and extend the analysis of Elman
and Golub [4], done for the two-dimensional (2D) variable coe#cient case, and the
analysis of Greif and Varah [8], [9] for the 3D problem with constant coe#cients to
the 3D problem with variable coe#cients.
Let h denote the width of a uniform mesh. In the description that follows we use
the notation G i,j,k # G(ih, jh, kh), where G is a trivariate function. The seven-point
discretization is done as follows (see, e.g., [4] for the analogous 2D case). For the first
term in (1.1) we have
and an analogous discretization is performed for (qu y ) y and (ru z ) z . For the convective
terms su x , tu y , and vu z we use either upwind or centered di#erence schemes.
Let F denote the corresponding di#erence operator, scaled by h 2 , and denote the
values of the associated computational molecule by a i,j,k , b i,j,k , c i,j,k , d i,j,k , e i,j,k ,
f i,j,k , and g i,j,k , in the following manner: if (i, j, k) is a gridpoint not next to the
boundary, then
F
# Received by the editors November 10, 1997; accepted for publication (in revised form) by Z.
August 20, 1998; published electronically August 3, 1999.
Department of Computer Science, The University of British Columbia, Vancouver, BC, V6T 1Z4
Canada (greif@cs.ubc.ca).
x
y
z
e
f
c a d

(a) seven-point operator (b) reduced operator
Fig. 1.1. Computational molecules of the unreduced and the reduced operators.
The computational molecule is graphically illustrated in Figure 1.1(a) (in the figure
the subscripts are dropped).
If centered di#erences are used to discretize the convective terms, the values of
the computational molecule are given by
a
uses upwind schemes, then the type of scheme depends on the sign of the
convective terms. Assuming that s, t, and v do not change sign in the domain, if they
are positive one can use the backward scheme, and if they are negative one can use
the forward scheme. Discretizing using backward di#erences yields
a
and for forward di#erences (1.4) needs to be modified in an obvious manner.
The sparsity structure of the matrix representing the system of equations depends
on the ordering of the unknowns. A common strategy is the red/black ordering, which
is depicted in Figure 1.2: the gridpoints are colored using two colors in a checkerboard
fashion, and the points that correspond to one of the colors (say, red) are numbered
first. In this case the corresponding matrix can be written as
# u (r)
REDUCED SYSTEMS FOR 3D ELLIPTIC EQUATIONS 31
x
y
z
43 44
53 54
26
28
29 30Fig. 1.2. Red/black ordering of the 3D grid.
where both B and E are diagonal. In (1.5) superscripts (r) and (b) are attached to
denote the associated colors. A simple process of block Gaussian elimination leads to
a smaller system, for the black points only, which is called a reduced system [2]:
Since B is diagonal, the matrix of (1.6) is sparse. In the 3D case the corresponding
di#erence operator has a computational molecule which consists of 19 points, as illustrated
in Figure 1.1(b). Once the solution for the black points is computed, the
solution for the red points corresponds to solving a diagonal system and thus is readily
obtained. Moving from system (1.5) to system (1.6) amounts to performing one
step of cyclic reduction [2]. This procedure can be repeated until a small system of
equations is obtained, which can be solved directly. An overview of the idea of cyclic
reduction and several references are given in [5, pp. 177-180].
The elimination of half of the unknowns is accompanied by permutation of the
matrix (equivalently, reordering of the unknowns). Once the permuted reduced system
is formed, an iterative method can be used to find the solution. The procedure
of performing one step of cyclic reduction for a non-self-adjoint problem and solving
the resulting system using an iterative solver was extensively investigated by Elman
and Golub for 2D problems [2], [3], [4]. They showed that one step of cyclic reduction
leads to systems with valuable properties, such as symmetrizability by a real diagonal
nonsingular matrix for a large set of the underlying PDE coe#cients (which is e#ec-
tively used to derive bounds on the convergence rates of iterative solvers), and fast
convergence. In [2] the univariate case, which is naturally more transparent, is used
to illustrate the advantages of this technique. Many of the highly e#ective techniques
presented and used by Elman and Golub can be generalized to the 3D case and will
be mentioned throughout this paper.
An outline of the rest of this paper follows. In section 2 we introduce the cyclically
reduced operator. In section 3 we discuss block orderings and present a family of
orderings for the reduced grid. In section 4 we present symmetrization results. In
section 5 we use the results of section 4 to derive bounds on convergence rates for
block stationary methods. In section 6 we present numerical examples, which include
solving the systems using Krylov subspace solvers. Finally, in section 7 we draw some
conclusions.
2. The reduced operator. One step of cyclic reduction for a 3D model problem
with constant coe#cients has been described in [9], where full details on the
construction of the matrix are given. Convergence analysis and techniques for solving
the resulting system of equations using block stationary methods are described in [8],
[9]. In the case of constant coe#cients, the values of the computational molecule associated
with a given gridpoint do not depend on the point's coordinates, as opposed
to the variable coe#cient case. Nevertheless, the sparsity structure of the reduced
matrix is identical in both cases.
The explicit di#erence equation associated with the reduced operator for the 3D
variable coe#cient case can be found in [7] and should be used for constructing the
reduced matrix. The alternative of performing the matrix products in (1.6) might be
significantly more costly, especially in the 3D case, and in particular, in programming
environments where vectorization is crucial (e.g., Matlab).
Consider the following constant coe#cient model problem: p(x, y, y, z)
y, z) # 1, s(x, y, y, y, scaling by ah 2 ,
the di#erence equation has the form
R u
-2cf
Denote the continuous operator corresponding to this model problem
# .
The reduced operator can be derived directly as a discretization scheme of the original
PDE, with O(h 2 ) correction terms in the case of centered di#erence discretization and
O(h) correction terms if upwind discretization is used. This can be done by means
analogous to the techniques used for the 2D case (see Elman and Golub [3] and Golub
and Tuminaro [6]). Consider the centered di#erence discretization. Expanding (2.2)
in a multivariate Taylor expansion about the gridpoint (ih, jh, kh) yields, after scaling
by 2ah 2 ,
R
The above computation was carried out using Maple V. The O(h 2 ) terms contain,
among other terms, the expression - h 2
which can be thought
of as addition of artificial viscosity to the original equation. The reduced right-hand
side is equal to w i,j,k with an O(h 2 ) error. Gaussian elimination yields the following
right-hand side,
a
a
a
a
a
a
REDUCED SYSTEMS FOR 3D ELLIPTIC EQUATIONS 33
(a) Natural (b) Red/black
Fig. 3.1. Two possible orderings of the block grid. Each point in these grids corresponds to a
one-dimensional (1D) block of gridpoints in the underlying 3D grid.
whose Taylor expansion about the gridpoint (ih, jh, kh), after scaling by 2ah 2 , is
evaluated at the point (ih, jh, kh). This is
another similarity with the 2D case [3].
3. Block ordering strategies for the reduced grid. The question of ordering
is of major importance, as a good ordering strategy can lead to fast convergence. An
excellent overview of the literature that deals with ordering strategies is found in a
recent report by Benzi, Szyld, and van Duin [1]. For 3D problems it seems useful
to consider the ordering of blocks of unknowns, rather than "pointwise" ordering.
Such a strategy could be particularly useful for the cyclically reduced problem, as
the reduced grid is somewhat irregular. Instead of ordering the unknowns directly in
the 3D reduced grid, the problem of ordering is divided into two parts. First, define
blocks of gridpoints and order them in a tensor-product 2D "block grid." Once the
ordering of the blocks is determined, the task of the "inner" ordering in each of the
blocks is relatively simple.
We can define an x-oriented "1D block of gridpoints" by referring to a set of
gridpoints whose collection of all x-coordinate values include all the possible values
{jh} on the grid. A simple example is a single horizontal line of gridpoints in a tensor-product
3D grid. Similarly, y-oriented and z-oriented 1D blocks can be defined. Once
the 1D blocks of gridpoints are defined, a block computational molecule can be defined
as follows.
Definition 3.1. For a certain given 3D grid and a 1D block of gridpoints in it,
the associated block computational molecule is defined as the computational molecule
in the corresponding block grid. That is, its components are the 1D blocks in the
block grid, each of which contains at least one gridpoint which belongs to the (point)
computational molecule associated with the 3D problem.
Using the above, we can now easily define di#erent families and types of orderings.
For example, a certain ordering strategy is a natural block ordering strategy relative
to the 1D blocks of gridpoints if these blocks are ordered in the block grid using
natural lexicographic ordering. Similarly, one can define a red/black block ordering
strategy, and so on (see Figure 3.1).
Below we focus on a particular family of orderings for the reduced grid, which
we call the two-plane ordering. This ordering corresponds to defining each of the 1D
blocks of gridpoints as a collection of 2n gridpoints from two horizontal lines and two
adjacent planes (here n is the number of gridpoints in a single line in the original
34 CHEN GREIF
x
y
(a) 2PNxy (b) 2PRBxy
Fig. 3.2. Two types of two-plane ordering.
x
y
Fig. 3.3. Possible orientations of the 1D blocks of gridpoints in the set of natural two-plane
orderings.
unreduced grid). A single member of this family was introduced in [9]. In Figure 3.2
two members of the family are depicted: natural two-plane ordering and red/black
two-plane ordering. For notational convenience, we label them "2PN" and "2PRB,"
respectively. Two additional letters are added in order to distinguish between di#erent
orientations of the 2n-item 1D blocks of gridpoints. Let us illustrate this for the
specific case depicted in Figure 3.2(a). Here Indices 1-8, 9-16, 17-24, and 25-
are each an x-oriented 1D block. The block grid is of size 2-2 and its components
are ordered in natural lexicographic fashion. Each of the sets of indices 1-16 and
17-32 forms an x-y-oriented pair of planes. Hence the name 2PNxy. In Figure 3.3
other orientations of the blocks in the natural two-plane ordering are depicted.

Figures

3.4(a) and (b) illustrate what blocks are associated with a single gridpoint.
In these figures, each "X" corresponds to a 1D block which contains at least one
gridpoint in the "point" computational molecule associated with the selected gridpoint
(the one that is at the center of the computational molecule). As is evident, the
REDUCED SYSTEMS FOR 3D ELLIPTIC EQUATIONS 35
x
x
x
x
x
x
x
x
x
x
(a) even
h (c) computational molecule
Fig. 3.4. Block computational molecule associated with the two-plane ordering.
(a) natural (b) red/black based on 2D blocks
Fig. 3.5. Sparsity structures of two matrices which belong to the family of two-plane orderings.
structure depends on the parity of this gridpoint's index. The block computational
molecule (Figure 3.4(c)) is obtained by taking the union of all the 1D blocks associated
with each of the gridpoints in the block, and thus it is identical to the computational
molecule of the classical nine-point operator. This allows one to conclude, e.g., that
the reduced matrix does not have block property A [14] relative to partitioning into
blocks. On the other hand, applying a 4-color scheme to the blocks of
gridpoints can be e#ective for parallelization.
The same ideas as above can be applied to 2D blocks of gridpoints. The block
grid in this case is univariate. The reduced matrix is consistently ordered relative
to partitioning associated with 2D blocks. In Figure 3.5 the sparsity structures of
matrices corresponding to natural two-plane ordering and red/black block
ordering relative to 2D blocks are depicted.
In order to illustrate the e#ectiveness of the two-plane ordering, we present in

Figure

3.6 a single 2D block of a natural two-plane matrix vs. one that corresponds
to ("point") natural lexicographic ordering. The matrices in the figure are associated
with a 12 - 12 - 12 grid. As is evident, the main diagonal block of the two-plane
matrix is more dense. Compared to the natural lexicographic ordering, there are
significantly more nonzero entries in the block diagonal submatrix whose bandwidth
does not depend on n. As a result, direct preconditioner solves will be more e#cient
due to less fill-in. If a stationary scheme such as block Jacobi is used, by [13, Thm.
3.15] it is guaranteed that if the reduced matrix is an M-matrix, then the convergence
of a scheme associated with the two-plane matrix is faster than that of a scheme with
36 CHEN GREIF
(a) lexicographic (b) natural two-plane
Fig. 3.6. A zoom on 2D blocks of the matrices corresponding to two ordering strategies of the
reduced grid.
the lexicographic ordering. (The circumstances in which the reduced matrix is an
M-matrix in the constant coe#cient case are discussed in [9].)
4. Symmetrization of the reduced matrix. In order to obtain bounds on
convergence rates for block stationary methods, we consider the following technique,
suggested and used in [2], [3], [4] and also e#ectively applied in [8], [9]: if there exists
a real diagonal matrix Q, such that -
then for a splitting
and an analogous splitting of -
namely -
we have
C)
D)
A bound for the spectral radius of the iteration matrix can thus be obtained by
evaluating the expression on the right-hand side of (4.1). Denote the entries of S
by {s i,j } n 3 /2
. Since S is sparse and has a block structure, a small amount of work
is needed in order to find Q-by requiring for the entries of -
S, which we denote by
{-s
s j,i . Since matrices that correspond to di#erent orderings are
merely symmetric permutations of one another, we can pick a matrix that corresponds
to a specific ordering strategy and do all the work for it. This will result in obtaining
general symmetrization conditions for the reduced matrices (regardless of the ordering
used). Thus we pick the specific ordering strategy 2PNxz.
Let q # denote the #th diagonal entry in Q. Then
and the symmetry conditions can be expressed as
s j,#
It is su#cient to look at 2n - 2n blocks. We start with the main diagonal blocks.
We have to examine all the entries of the main block that appear in the #th row of
REDUCED SYSTEMS FOR 3D ELLIPTIC EQUATIONS 37
the matrix, namely s #-4 , s #-3 , . , s # , . , s #+4 . For s #-4 , if # mod 2n # 5 or
is equal to 0, then # - 4 corresponds to the (i - 2, j, point. Thus
#, -
a -
and from this it follows that
a i-1,j,k
a
c i-1,j,k c i,j,k
d i-2,j,k d i-1,j,k
In this case the values associated with the center of the computational molecule
(namely, a i,j,k ) are canceled, but this happens only for rows that involve the
gridpoints. Applying the same procedure
to the rest of the entries of the main diagonal block, we obtain the following:
a
a i-1,j,k
a
a i-1,j,k
c i,j+1,k e i,j,k
a i,j+1,k
a i-1,j,k
b i,j+1,k d i-1,j+1,k
a i,j+1,k
a i-1,j,k
a i-1,j,k
a i,j-1,k
d i-1,j,k e i-1,j-1,k
a i-1,j,k
a i,j-1,k
a i-1,j,k
a i,j,k+1
a i-1,j,k
a i,j,k+1
a
a i,j+1,k
a
a i,j+1,k
a
a i,j-1,k
a
a i,j-1,k
As is evident, (4.5)-(4.11) overdetermine the nonzero values of the matrix Q.
Indeed, (4.9)-(4.11) are su#cient to determine all the diagonal entries, except the first
entry in each 2n-2n block, which at this stage can be arbitrarily chosen. We have to
make sure, therefore, that (4.5)-(4.8) are consistent with these three equations, and
this requirement imposes some additional conditions. In the constant coe#cient case
there is unconditional consistency. The problematic nature of the variable coe#cient
case can be demonstrated simply by looking at one of the consistency conditions.
Consider a gridpoint (i, j, whose associated index, #, satisfies # mod
Applying (4.9) to # - 1 means looking at the row corresponding to the (i,
gridpoint, and multiplying (4.9), applied to # - 1, by (4.11) results in an equation
for q #
, which should be consistent with (4.8). There are three additional
consistency conditions for the main block and then eight additional conditions for
the rest of the blocks of the reduced matrix. In the consistency condition, if we
equate variables that belong to the same location in the computational molecule,
we find that su#cient conditions for the above-mentioned consistency conditions to
hold are b
similarly c , the consistency
condition becomes
d i-1 e j-1
which is obviously satisfied. The actual meaning of these conditions is that the continuous
problem is separable.
The analysis for o#-diagonal blocks is identical, and the following additional conditions
are obtained:
The two equations above determine the rest of the entries of the matrix, and only
the first entry in the symmetrizer can be determined arbitrarily.
Last, in order for the symmetrizer to be real, we must require that the products
have the same sign.
All that has been said can be summarized in the following theorem, which demonstrates
another point of similarity between the 2D [4] and the 3D problems with
variable coe#cients.
Theorem 4.1. Suppose the operator of (1.1) is separable. If c i d i-1 , b j e j-1 , and
are all nonzero and have the same sign for all i, j, and k, then there exists a
real nonsingular diagonal matrix Q such that Q -1 SQ is symmetric.
The symmetrized computational molecule can be derived without actually performing
the similarity transformation. For example, the symmetrized value corresponding
to - c i,j,k c i-1,j,k
a i-1,j,k
is
a i-1,j,k
, and so on. The symmetrization
operation should not actually be performed in order to solve the linear system,
as the symmetrizing matrix has entries that are unbounded as h goes to zero. The
symmetrization should be done for the mere purpose of convergence analysis.
5. Bounds on convergence rates for block stationary methods. The reduced
matrix which corresponds to the family of natural two-plane orderings is of
size (n 3 /2) - (n 3 /2) and can be thought of as a block tridiagonal matrix, relative to
S i,j are block tridiagonal matrices with respect to 2n - 2n blocks.
REDUCED SYSTEMS FOR 3D ELLIPTIC EQUATIONS 39
In [8] two partitionings are considered: a partitioning into 1D blocks (2n - 2n
blocks) and a partitioning into 2D blocks (of size n 2
In order to find a bound,
if symmetrization is possible, then the strategy in [4] can be applied. In Theorem
5.1, we refer to the ordering strategy 2PNxz. Since other orderings are symmetric
permutations of 2PNxz, finding the bounds for other ordering strategies discussed in
this paper is straightforward.
Theorem 5.1. Suppose the continuous problem is separable and c i+1 d i , b j+1 e j ,
and f k+1 g k are all positive and bounded by # x , # y , and # z , respectively. Suppose also
that a i,j,k # for all i, j, and k. Denote -
. Then the spectral radii of the
iteration matrices associated with the block Jacobi scheme which correspond to 1D
splitting and 2D splitting are bounded by #
# and #
respectively, where #, and
# are defined as follows:
(5.2a)
(5.2c)
Proof. The proof follows by using the technique of [4, pp. 346-347]. The conditions
stated in the theorem guarantee that the matrix is symmetrizable. Denote the reduced
matrix by S and the symmetrized matrix by -
S. Suppose S # is obtained by modifying
S in the following manner: replace each occurrence of c i and d i by # x , replace each
occurrence of b j and e j by # y , replace each occurrence of f k and g k by # z , and
replace each occurrence of a i,j,k by #. Denote by S the splitting, which
is analogous, as far as sparsity structure is concerned, to the splitting
the 1D splitting, the matrix D # is block diagonal with semibandwidth 4, its sparsity
structure is identical to that of -
D, and moreover, it is componentwise smaller than or
equal to the entries of -
D. By [9, Lem. 3.3], D # is an irreducible diagonally dominant
M-matrix.
The matrix C # is nonnegative and satisfies C # -
C. Thus the Perron-Frobenius
theorem [13, p. 30] can be used to obtain an upper bound on the convergence rate
for this splitting. Since the matrix S # can now be referred to as a symmetrized
version of a matrix that is associated with a constant coe#cient case, the bound
on the convergence rate is readily obtained from [9, Thm. 3.15]. For the 2D splitting
the procedure is completely analogous and the bound is obtained from
[8, Thm. 3.6].
We remark that estimates for the convergence rates of block Gauss-Seidel and
block SOR schemes can be obtained by using the "near property A" analysis presented
in [7], [8].
6. Numerical experiments. In the examples that follow, we begin with some
results which validate the convergence analysis of section 5 for stationary methods.
We then compare the performance of Krylov subspace solvers for the reduced and the
unreduced systems. The experiments were performed on an SGI Origin 2000 machine.
The program is written in Matlab 5.
6.1. Test problem 1. Consider the separable problem
y, z)

Table
Comparison between the computed spectral radii of the block Jacobi iteration matrices and the
bounds, using centered di#erences for the two splittings, with
Splitting 1D 2D
y, z) is
constructed so that the solution is
u(x, y,
For notational convenience, let
2 . Suppose h is
su#ciently small and centered di#erence discretization is performed. Then
hence
If
k the bounds are obtained
in an identical manner. The center of the computational molecule is a = 6. In terms
of the PDE coe#cients, the condition on # means that the convergence analysis of
section 5 is applicable if the PDE coe#cients are O(n).
If the above conditions hold, the matrix is symmetrizable. Using the notation of
the previous sections, let -
S be the symmetrized matrix, let #
be a modified version of -
S, such that
each occurrence of c i+1 d i , b j+1 e j , and f k+1 g k in -
S is replaced by the upper bounds,
namely # x , # y , and # z , respectively. Since S # -
is a symmetrized version of
a matrix corresponding to the constant coe#cient case, and by [9, Lem. 3.3], it is a
diagonally dominant M-matrix. Using Theorem 5.1, the bounds on the convergence
rate of the block Jacobi method are given in Table 6.1. As is evident, the bounds are
tight even for small n. It should be noted, however, that as the PDE coe#cients grow
larger, the bounds are not expected to be as tight, as the inequalities in (6.4) become
less e#ective.
Next, the spectral radii of the block Jacobi and block Gauss-Seidel iteration
matrices and an approximation to the optimal SOR parameter have been computed
for both the upwind scheme and the centered scheme. The relaxation parameter was
computed according to the formula 2
. It is only an approximation to the
optimal SOR parameter since the matrix is only "close" to being consistently ordered
[8]. The experiments were done on an 8 - 8 - 8 grid (512 gridpoints). In Tables 6.2
and 6.3 the superscripts c and u stand for centered or upwind, respectively, R and U
stand for reduced system and unreduced system, respectively, and the subscripts J
and GS stand for Jacobi and Gauss-Seidel, respectively. We present the results for
REDUCED SYSTEMS FOR 3D ELLIPTIC EQUATIONS 41

Table
Spectral radius of the block Jacobi, block Gauss-Seidel, and approximate optimal relaxation
parameter for the reduced system, using both upwind and centered di#erences and 1D splitting.
# u (R) # c

Table
Spectral radii of the block Jacobi, block Gauss-Seidel, and optimal relaxation parameter for the
unreduced system, using both upwind and centered di#erences with 1D splitting.
two di#erent cases, one of which has convection of moderate size, and the other with
large convection, for which the upwind scheme is more e#ective than the centered
scheme.
Note that for the unreduced system, in most cases the matrix satisfies all the
conditions required for Young's SOR analysis [14]; thus, the spectral radius of the
Gauss-Seidel matrix and the optimal relaxation parameter can be computed from
the spectral radius of the block Jacobi matrix. By comparing Tables 6.2 and 6.3
it is evident that stationary solvers for the reduced system converge faster than for
the unreduced system. In one case there is convergence for the reduced system and
divergence for the unreduced system.
Moving to consider Krylov subspace solvers, in Table 6.4 we make a comparison
between the performance of solvers for the two systems. The stopping criterion was
relative residual smaller than 10 -10 . The method that is used is nonpreconditioned
Bi-CGSTAB. The table presents information on the complete process, namely, construction
of the systems and the iterative solves. The increase in iteration counts
as the grid is refined agrees with theory, at least if one assumes that for this well
conditioned and mildly nonsymmetric system, the condition number is of magnitude
O(h -2 ) and the convergence rate is similar to that of the conjugate gradient method
for symmetric positive definite systems. When one step of cyclic reduction is applied,
the savings become more dramatic as the systems grow larger. An explanation for
this is that the construction of the reduced system, which requires significantly more
floating point operations compared with the construction of the unreduced system,
becomes a less significant factor in the overall computation as the grid becomes finer.
In general, since the iterative solve is the costly component of the computation, it is
significant that the number of iterations until convergence of the unreduced solver is
larger by a factor of approximately 2 compared with the reduced solver. Figure 6.1
illustrates the saving and the convergence behavior for this problem.
In

Table

6.5 we provide some numerical evidence which suggests that the good
performance of reduced solvers is due to e#ective preconditioning of the original ma-
trix. In the table, estimates of the condition numbers for
with upwind di#erence discretization are presented. The estimates were obtained using
Matlab's command ``condest.'' The factor of approximately 2 has been obtained
for several additional cases that have been tested. More observations on the condition

Table
Comparison between the performance of the unreduced and the reduced solvers for increasing
mesh size, using nonpreconditioned Bi-CGSTAB, for
Iterations Mflops Time (sec.)
Unreduced Reduced Unreduced Reduced Unreduced Reduced
iterations
relative
residual
Fig. 6.1. Relative residuals for nonpreconditioned Bi-CGSTAB applied to linear systems arising
from discretization of the test problem, with grid. The
residual associated with the reduced solver is the lower curve.
number of the reduced matrix can be found in [7].
6.2. Test problem 2. Consider the nonseparable problem
y, z)
y, z) is
constructed so that the solution is (6.2). For this problem the convergence analysis of
section 5 does not apply. Results are given in Table 6.6. The experiments were done
so that the tensor-product grid has 13, 824 gridpoints. GMRES(5) [12]
was used, preconditioned by ILU with drop tolerance of 10 -3 . The stopping criterion
was ||r i ||/||r 0 In all cases that have been tested, setting up and solving
the reduced system is faster compared with setting up and solving the unreduced
system. It should be noted that the CPU times are a#ected by a long preconditioner
setup time. When convection dominates, the centered scheme performs poorly. (In
general, this scheme su#ers numerical instability when the Reynolds numbers are large
[10].) Additional numerical experiments indicate that both solvers have the property
that for centered di#erence discretization the solver of mildly nonsymmetric systems
converges faster than the solver of a close-to-symmetric system. This phenomenon
was proved analytically for stationary methods for the constant coe#cient case in [2]
(2D) and in [9] (3D).
REDUCED SYSTEMS FOR 3D ELLIPTIC EQUATIONS 43

Table
Comparison between estimates of condition numbers of the unreduced matrix (denoted by U)
vs. the reduced matrix (R), for

Table
Comparison of CPU times (seconds) for setting up and solving the unreduced system and the
reduced system using ILU+GMRES. R and U stand for reduced and unreduced, respectively, and the
subscripts c and u stand for centered and upwind, respectively.
50 91.2 43.5 48.7 32.7
100 148.2 37.3 62.1 23.8
6.3. Test problem 3. Consider the nonseparable problem
with Neumann boundary conditions u conditions
on the unit cube. w was constructed so that the exact solution
is u(x, y, cos(#z). Here we compare the performance of a few
Krylov subspace solvers; thus, the focus is on the actual iterative solve time, once the
systems and the preconditioners were set up.
The results in Table 6.7 are for a 20 - 20 - 20 grid. ILU(0) was used as a
preconditioner. The stopping criterion was ||r i ||/||r 0 In all cases the reduced
solver converges faster than the unreduced solver. The factor of approximately 2 can
be a good indication for the gain in performing one step of cyclic reduction in much
finer grids. Bi-CGSTAB is slightly faster than CGS. These two schemes are faster than
BiCG. The di#erences in performance between the solvers are qualitatively similar for
the reduced and the unreduced systems.

Table
Iteration counts and solving times for three Krylov solvers.
Method Iterations Time (sec.)
Unreduced Reduced Unreduced Reduced
CGS
7. Concluding remarks. A cyclically reduced operator for a 3D convection-
di#usion equation with variable coe#cients has been derived. Block orderings have
been discussed, some solving techniques for the reduced system have been examined,
44 CHEN GREIF
and numerical experiments illustrate the fact that the reduced system is easier to
solve than the unreduced system.
The results presented in this work show that one step of cyclic reduction can
be e#ectively used as a preconditioning technique for solving the convection-di#usion
equation with variable coe#cients. The questions of parallelism and applications are
topics for further investigation.

Acknowledgments

. I am very grateful to Jim Varah for carefully reading this
manuscript and for his valuable advice. I would also like to thank Gene Golub, who
has held helpful discussions with me and has pointed out some useful references. Many
thanks to the two referees for their helpful comments and suggestions.



--R

Orderings for incomplete factorization preconditioning of nonsymmetric problems.
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
Matrix Computations
Cyclic Reduction/Multigrid
Analysis of Cyclic Reduction for the Numerical Solution of Three-Dimensional Convection-Di#usion Equations
Block stationary methods for nonsymmetric cyclically reduced systems arising from three-dimensional elliptic equations SIAM J
Iterative solution of cyclically reduced systems arising from discretization of the three-dimensional convection di#usion equation
Don't suppress the wiggles
Numerical Solution of Convection-Di#usion Problems
Iterative Methods for Sparse Linear Systems
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
cyclic reduction;variable coefficients;3D elliptic problems