Preserving Symmetry in Preconditioned Krylov Subspace Methods.
We consider the problem of solving a linear system A when A is nearly symmetric and when the system is preconditioned by a symmetric positive definite matrix M. In the symmetric case, we can recover symmetry by using M-inner products in the conjugate gradient (CG) algorithm. This idea can also be used in the nonsymmetric case, and near symmetry can be preserved similarly. Like CG, the new algorithms are mathematically equivalent to split preconditioning but do not require M to be factored. Better robustness in a specific sense can also be observed. When combined with truncated versions of iterative methods, tests show that this is more effective than the common practice of forfeiting near-symmetry altogether.
Introduction
Consider the solution of the linear system
by a preconditioned Krylov subspace method. Assume at first that A is symmetric positive
definite (SPD) and let M be an SPD matrix that is a preconditioner for the matrix A. Then
one possibility is to solve either the left-preconditioned system
This work was supported by the NSF under contract NSF/CCR-9214116, by the ONR under contract
ONR-N00014-92-J-1890, and by the ARO under contract DAAL-03-9-C-0047 (Univ. of Tenn. subcontract
ORA 4466.04, Amendment 1 and 2).
y Department of Mathematics, University of California Los Angeles, 405 Hilgard Avenue, Los Angeles,
California, 90095-1555 (chan@math.ucla.edu and myeung@math.ucla.edu).
z Department of Computer Science and Minnesota Supercomputer Institute, University of Minnesota,
4-192 EE/CSci Bldg., 200 Union St., S.E., Minneapolis, Minnesota, 55455-0154 (chow@cs.umn.edu and
saad@cs.umn.edu).
or the right-preconditioned system
Both of these systems have lost their symmetry. A remedy exists when the preconditioner
M is available in factored form, e.g., as an incomplete Choleski factorization,
in which case a simple way to preserve symmetry is to 'split' the preconditioner between left
and right, i.e., we could solve
which involves a symmetric positive definite matrix. This can also be done when M can
be factored as . Unfortunately, the requirement that M be available in
factored forms is often too stringent.
However, this remedy is not required. As is well-known, we can preserve symmetry by
using a different inner product. Specifically, we observe that M \Gamma1 A is self-adjoint for the
M-inner product,
since we have
If we rewrite the CG-algorithm for this new inner product, denoting by r the
original residual and by z the residual for the preconditioned system, we would
obtain the following algorithm, see, e.g., [4].
Algorithm 1.1 Preconditioned Conjugate Gradient
1. Start: Compute
2. Iterate: For convergence do,
(a)
(b) x
(c) r
(d) z
Note that even though M-inner products are used, we do not need to multiply by M which
may not be available explicitly; we only need to solve linear systems with the matrix co-efficient
M . It can also be seen that with a change of variables, the iterates produced by
this algorithm are identical to both those of CG with split preconditioning, and CG with
right-preconditioning using the M \Gamma1 -inner product. All three of these algorithms are thus
equivalent.
The question we now raise is the following. When A is only nearly symmetric, it is
often the case that there exists a preconditioner M which is SPD. In this situation, the use
of either of the forms (2) or (3) is just as unsatisfactory as in the fully symmetric case.
Indeed, whatever degree of symmetry was available in A is now entirely lost. Although
the above remedy based on M-inner products is always used in the symmetric case, it is
rather surprising that this problem is seldom ever mentioned in the literature for the nearly
symmetric case. In the nonsymmetric case, when M exists in factored form, some form of
balancing can also be achieved by splitting the preconditioner from left and right. However,
there does not seem to have been much work done in exploiting M-inner products even when
M is not available in factored form. This dichotomy between the treatment of the symmetric
and the nonsymmetric cases is what motivated this study.
Ashby et. al. have fully considered the case of using alternate inner products when the
matrix A is symmetric. Work of which we are aware that consider the use of alternate inner
products when A is near-symmetric are Young and Jea [9] and Meyer [5]. In the latter, the
product (with M SPD) is used with ORTHOMIN and ORTHODIR.
This paper is organized as follows. In Section 2, it is shown how alternative inner products
may be used to preserve symmetry in GMRES. Section 3 considers the use of truncated
iterative methods when the preconditioned system is close to being symmetric. This has
been hypothesized by many authors, for example, Axelsson [2] and Meyer [5]. In Section 4 we
consider the symmetrically preconditioned Bi-CG algorithm. Section 5 tests these algorithms
numerically using a Navier-Stokes problem parameterized by the Reynolds number, and thus
nearness to symmetry. We conclude this paper in Section 6.
preconditioning in GMRES
When A is nearly symmetric, split preconditioning may be used to preserve the original
degree of symmetry. Alternatively, left-preconditioning with the M-inner product, or right-
preconditioning with the M \Gamma1 -inner product may be used. These latter two alternatives
will be developed for the Arnoldi process, and used as the basis of 'symmetric' preconditioned
versions of GMRES. Like CG, it will be shown that these symmetric versions are
mathematically equivalent to split preconditioning, but do not require the preconditioner to
be symmetrically factored. We begin by exploring the options and implementation issues
associated with left symmetric preconditioning.
2.1 Left symmetric preconditioning
The GMRES algorithm is based on the Arnoldi process. Without preconditioning, the
Arnoldi algorithm based on the classical Gram-Schmidt process is as follows.
Algorithm 2.1 Arnoldi-Classical Gram-Schmidt
1. Choose a vector v 1 of norm 1.
2. For
3.
4. Compute
5.
7. If h
8.
9. EndDo
Consider here the case when A is left-preconditioned, i.e., the matrix A involved in the
algorithm is to be replaced by some preconditioner, which we
assume to be SPD. We wish to define a procedure to implement the above algorithm for the
using M-inner products, and if possible, avoid additional matrix-vector
products, e.g., with M . Once this is accomplished, we will define the corresponding GMRES
procedure.
In the preconditioned case, it is customary to define the intermediate vectors in the
product
which is then preconditioned to get
We now reformulate the operations of the above algorithm for with the M-inner
product. Only the computation of the inner products changes. In the classical Gram-Schmidt
formulation, we would first compute the scalars h ij in Line 4 of Algorithm 2.1,
Then we would modify the vector z j to obtain the next Arnoldi vector (before normalization),
To complete the orthonormalization step we must normalize the final "
z j . Because of the
M-orthogonality of " z j versus all previous v i 's we observe that
Thus, the desired M-norm can be computed according to (9) and then computing
One potentially serious difficulty with the above procedure is that the inner product
as computed by (9) may be negative in the presence of round-off. There are two
remedies. First, we can compute this M-norm explicitly at the expense of an additional
matrix-vector multiplication with M , i.e., from
As was pointed out earlier, this is undesirable, since the operator M is often not available
explicitly. Indeed in many cases, only the preconditioning operation M \Gamma1 is available from a
sequence of operations, as is the case for multigrid preconditioning. Another difficulty with
computing h ij with (7) is that it is not immediately amenable to a modified Gram-Schmidt
implementation. Indeed, consider the first step of a hypothetical modified Gram-Schmidt
step, which consists of M-orthonormalizing z against v 1 ,
As was observed, the inner product (z; is equal to (w; v 1 ) which is computable. Now
we need M "
z to compute in the modified Gram-Schmidt process. However,
no such vector is available, and we can only compute the (z; v i ) M of classical Gram-Schmidt.
An alternative is to save the set of vectors Mv i (again, not computed by multiplying by
M) which would allow us to accumulate inexpensively both the vector "
z j and the vector "
via the relation
which is obtained from (8). Now the inner product ("z is given by
In this form, this inner product is guaranteed to be nonnegative as desired. This leads to
the following algorithm.
Algorithm 2.2 Arnoldi-Classical Gram-Schmidt and M-inner products
1. Choose a vector w 1 such that v has M-norm 1.
2. For
3.
4. Compute
5. "
w .
As is noted, the above algorithm requires that we save two sets of vectors: the v j 's and
the w i 's. The v i 's form the needed Arnoldi basis, and the w i 's are required when computing
the vector "
w j in Line 5. If we do save these two sets of vectors we can also now easily
formulate the algorithm with the modified Gram-Schmidt version of the Arnoldi procedure.
Algorithm 2.3 Arnoldi-Modified Gram-Schmidt and M-inner products
1. Choose a vector w 1 such that v has M-norm 1.
2. For
3.
4. For
5.
7. EndDo
8.
9. h
11.
12. EndDo
2.2 Right symmetric preconditioning
The matrix AM \Gamma1 is self-adjoint with the M \Gamma1 -inner product. The situation for right-
preconditioning with this inner product is much simpler, mainly because M \Gamma1 z is available
when z needs to be normalized in the M \Gamma1 norm. However, M \Gamma1 z is normally computed at
the next iteration in the standard Arnoldi algorithm; a slight reorganization of the Arnoldi-
Modified Gram-Schmidt algorithm yields the following.
Algorithm 2.4 Arnoldi-Modified Gram-Schmidt and M \Gamma1 -inner products
1. Choose a vector v 1 of M \Gamma1 -norm 1, compute w
2. For
3.
4. For
5.
7. EndDo
8.
Note that the preconditioned vector is computed in Line 8, while in the standard algorithm
it is computed before Line 3. Again, both the v's and the w's need to be saved, where
in this case.
The additional storage of the w's, however, makes this algorithm naturally `flexible,' i.e.,
it accommodates the situation where M varies at each step as when M \Gamma1 v is the result
of some unspecified computation. If M \Gamma1 is not a constant operator, then a basis for the
right-preconditioned Krylov subspace cannot be constructed from the v's alone. However,
the vectors w
do form a basis for this subspace, where M \Gamma1
denotes the preconditioning
operation at the j-th step. The use of this extra set of vectors is exactly how the
standard flexible variant of GMRES is implemented [6].
2.3 Using M-inner products in GMRES
The vectors v i form an orthonormal basis of the Krylov subspace. In the following we
denote by the matrix whose column vectors are the vectors v i produced by
the Arnoldi-Modified Gram-Schmidt algorithm with M-inner products (Algorithm 2.3). A
similar notation will be used for the matrix Wm . We also denote by
Hm the (m+1) \Theta m upper
Hessenberg matrix whose nonzero entries h ij are defined by the algorithm. Hm denotes the
portion of
Hm . These matrices satisfy a number of relations similar to the ones
obtained from using the standard Euclidean inner product.
Proposition 2.1 The following properties are satisfied by the vectors v i and w i in Algorithm
2.3:
Hm ,
2.
Hm ,
3. V T
4. W T
5. If A is Hermitian then Hm is Hermitian and tridiagonal.
Consider now the implementation of a GMRES procedure based on the orthogonalization
process of Algorithm 2.3. Since we are using M-inner products we should be able to minimize
the M-norm of the residual vectors M all vectors of the affine subspace
where z as the first coordinate
vector. Then we have,
Hm y
Therefore we have the equality,
Hm
A result of the equality (12) is that we can minimize the M-norm of the (preconditioned)
residual vector M \Gamma1 (b\GammaAx) by simply minimizing the 2-norm of fie 1 \gamma
Hm y as in the standard
GMRES algorithm.
Algorithm 2.5 Left-preconditioned GMRES with M-inner products
1. Compute
2. For
3.
4. For
5.
7. EndDo
8.
9. h
11.
12. EndDo
13. Compute the minimizer y m of kfie 1 \gamma
14. Compute the approximate solution
15. If satisfied Stop; else set x
An equality similar to (12) can be shown for the right-preconditioned case with M
inner products. We can summarize with the following theorem, which we state without
proof.
Theorem 2.1 The approximate solution xm obtained from the left-preconditioned GMRES
algorithm with M-inner products minimizes the residual M-norm kM \Gamma1 (b \gamma Ax)kM over all
vectors of the affine subspace x 0 +Km in which
where z . Also, the approximate solution xm obtained from the right-preconditioned
GMRES algorithm with M \Gamma1 -inner products minimizes the residual M \Gamma1 -norm kb
over the same affine subspace.
2.4 Equivalence of the algorithms
We can show that both left and right symmetric preconditioning are mathematically equivalent
to split preconditioning. In the latter case, M must be factored into and we
solve
Denoting by B the preconditioned matrix the GMRES procedure applied to
the above system for the u variable, minimizes the residual norm,
over all vectors u in the space u 0 +K (u)
m with
in which
~
Note that the variables u and x are related by As a result, this procedure
minimizes
over all x in the space x 0 +K (x)
m with
We now make the following observation. For any k  0 we have,
where z . Indeed, this can be easily proved by induction. Hence, the space K (x)
m is
identical with the space
which is nothing but (13). In noting that
we have proved the following result.
Theorem 2.2 Let . Then the approximate solution obtained by GMRES applied
to the split preconditioned system (14) is identical with that obtained from the GMRES
algorithm for the left preconditioned system (2) using the M-inner product.
Again, the same statement can be made about right-preconditioning. All that must be
noticed is that the same minimization (16) is taking place, and that the minimization is over
the same subspace in each of the left, right, and split preconditioning options [7, Sec. 9.3.4].
We emphasize in particular that it is the split preconditioned residual that is minimized in
all three algorithms.
3 Truncated iterative methods
Truncated iterative methods are an alternative to restarting, when the number of steps required
for convergence is large and the computation and storage of the Krylov basis becomes
excessive. When A is exactly symmetric, a three-term recurrence governs the vectors in the
Arnoldi process, and it is only necessary to orthogonalize the current Arnoldi vector against
the previous two vectors. If A is nearly symmetric, an incomplete orthogonalization against
a small number of previous vectors may be advantageous over restarted methods. The advantage
here may offset the cost of maintaining the extra set of vectors to maintain the
initial degree of symmetry. The incomplete Arnoldi procedure outlined below stores only
the previous k Arnoldi vectors, and orthogonalizes the new vectors against them. It differs
from the full Arnoldi procedure only in Line 4, which would normally be a loop from 1 to j.
It can be considered to be the full Arnoldi procedure when k is set to infinity.
Algorithm 3.1 Incomplete Arnoldi Procedure
1. Choose a vector v 1 of norm 1.
2. For do
3.
4. For do
5.
7. EndDo
8. h
9. If h
11. EndDo
The truncated version of GMRES uses this incomplete Arnoldi procedure and is called
Quasi-GMRES [3]. The practical implementation of this algorithm allows the solution to be
updated at each iteration, and is thus called a 'direct' version, or DQGMRES [8].
To suggest that truncated iterative methods may be effective in cases of near symmetry,
we study the asymptotic behavior of the iterates of DQGMRES as the coefficient matrix A
varies from nonsymmetry to (skew) symmetry. We first decompose A as
in which S is symmetric or skew symmetric, and set
We will first establish asymptotic relations among the variables in the incomplete and full
Arnoldi procedures. Then we will apply the incomplete procedure to A, and the full procedure
to S, using the superscripts I and F to distinguish between the variables appearing in
the two procedures. (Note that since S is (skew) symmetric, the full procedure on S is the
same as the incomplete procedure with k  2.)
Moreover, if we denote the degree of the minimal polynomial of v F
1 with respect to S by
, then h F
. In the proof of the following lemma, we
also use " v I
j to denote the vectors w I and w F obtained at the end of Line 7 in the
incomplete and complete Arnoldi procedures.
Lemma 3.1 Assume the truncation parameter k  2. If v I
h I
Proof. The proof is by induction on the index j. By Lines 5 and 6 of the Arnoldi procedure,
h I
we have
h I
and hence the lemma holds for 1. Assume that the lemma has been proved for
On that hypothesis, we prove it for of the Arnoldi procedure,
I
=h I
which yields that
I
by the induction hypothesis. Therefore
for the w I and w F in Line 3 of the Arnoldi procedures. Using another induction on the index
i in Lines 5 and 6, and the induction hypothesis on j and, in the mean time, noting that
h I
for
From the last equation,
h I
and then the induction step is complete. QED.
We now turn to the DQGMRES algorithm. Consider the linear system
and denote by x G
m and x Q
m the approximate solutions by the GMRES and DQGMRES algo-
rithms, respectively. Let  be the degree of the minimal polynomial of the vector b \gamma Sx 0
with respect to S. A result of the lemma can be stated as follows.
Theorem 3.1 Given the same initial guess x 0 to GMRES and DQGMRES with k  2, then
at any given step m with
Proof. By the definitions of DQGMRES and GMRES, we have
H I
H I
H I
and
where
(b
(b
we have by the lemma,
H I
and therefore the desired equation holds. QED.
If we let xA be the exact solution to be the exact solution to
it is obvious that O("). Since, on the other hand, x G
immediately have
the following corollary.
Corollary 3.1 For any initial guess x 0 and any k  2,
The corollary suggests that we may use DQGMRES with small k when A is nearly
symmetric or nearly skew symmetric.
4 Symmetric preconditioning in Bi-CG
The Bi-CG algorithm is based on Lanczos biorthogonalization. Both left-symmetric and
right-symmetric preconditioning are relatively straightforward, and no extra vectors are re-
quired. For reference, Algorithm 4.1 gives the right-preconditioned Bi-CG algorithm with
preconditioner M . The symmetric right-preconditioned Bi-CG algorithm (right-preconditioned
using M \Gamma1 -inner products) is developed immediately afterward.
Algorithm 4.1 Right-preconditioned Bi-CG
1. Compute
0 such that (r
2. 3. For convergence Do:
4. ff
5. x
6. r
7. r
8.
9.
11. EndDo
Note in Line 7 of the above algorithm that the preconditioned coefficient matrix of the dual
system is (AM the dual residual r   is the residual of
which is a left-preconditioned version of some linear system with A T .
To develop the symmetric right-preconditioned Bi-CG, M \Gamma1 -inner products are used in
Algorithm 4.1 above. However, the preconditioned coefficient matrix of the dual system
must be the adjoint of AM \Gamma1 in the M \Gamma1 inner product. This is A T M \Gamma1 as shown by
The dual system thus involves the coefficient matrix A T M \Gamma1 . Algorithm 4.2 gives the symmetric
right-preconditioned Bi-CG algorithm with preconditioner M .
Algorithm 4.2 Right-preconditioned Bi-CG with M \Gamma1 -inner products
1. Compute
0 such that (M \Gamma1 r
2. 3. For convergence Do:
4. ff
5. x
6. r
7. r
8.
9.
11. EndDo
Like GMRES, both left and right symmetric preconditioned versions of Bi-CG are equivalent
to the split preconditioned version, and this can be shown by a change of variables.
However, in both left and right symmetric preconditioned versions, the exact, rather than
the split preconditioned residual is available.
The unpreconditioned Bi-CG algorithm cannot have a serious breakdown if A is SPD and
r
0 is chosen to be r 0 . This is because r
and the vectors Ap
never become orthogonal. In fact, the cosine
can be bounded below by the reciprocal of the condition number of A.
Similarly, in the symmetric right-preconditioned version of Bi-CG, if both A and M are
SPD, and r
We measure the cosines rather than the quantities (Ap
because the
and r vectors have magnitudes going to 0 as the algorithms progress. Recall that in the
case when (M
we have a lucky breakdown.
For the case of regular right- or left-preconditioning, or if r
in the symmetrically
preconditioned cases, then no such lower bounds as the above exist, and the algorithms are
liable to break down.
When A is near-symmetric, it is our hypothesis that the probability of breakdown is
lower in the symmetrically preconditioned cases, and this will be shown by experiment in
the next section.
5 Numerical Experiments
Section 5.1 tests the idea of using symmetric preconditionings with truncated iterative meth-
ods. Section 5.2 tests the breakdown behavior of symmetrically preconditioned Bi-CG.
5.1 Truncated iterative methods
To test the idea of using symmetric preconditionings with truncated iterative methods for
nearly symmetric systems, we selected a standard fluid flow problem where the degree of
symmetry in the matrices is parameterized by the Reynolds number. The flow problem
is the two-dimensional square-lid driven cavity, and was discretized by the Galerkin Finite
Element method. Rectangular elements were used, with biquadratic basis functions for
velocities, and linear discontinuous basis functions for pressure. We considered a segregated
solution method for the Navier-Stokes equations, where the velocity and pressure variables
are solved separately; the matrices arising from a fully-coupled solution method are otherwise
indefinite. In particular, we considered the expression of the conservation of momentum,
where u denotes the vector of velocity variables, p denotes the pressure variable, and Re is
the Reynolds number. The boundary conditions for the driven cavity problem over the unit
square are on the top edge of the square, and on the other three sides
and the corners. The reference pressure specified at the bottom-left corner is 0.
The matrices are the initial Jacobians at each Newton iteration, assuming a zero pressure
distribution. For convenience, however, we chose the right-hand sides of the linear systems
to be the vector of all ones. A mesh of 20 by 20 elements was used, leading to momentum
equation matrices of order 3042 and having 91204 nonzero entries. The nodes corresponding
to the boundaries were not assembled into the matrix, and the degrees of freedom were
numbered element by element. For Reynolds number 0., the matrix is SPD, and is equal to
the symmetric part of the matrices with nonzero Reynolds number.
We generated matrices with Reynolds number less than 10, which gives rise to the nearly
symmetric case. For Reynolds number 1., the degree of symmetry measured by
has value 7:5102 \Theta 10 \Gamma4 and this measure increases linearly with the Reynolds number (at
least up to Re=10).
In the numerical experiments below, we show the number of matrix-vector products
consumed by GMRES(k) and DQGMRES(k) to reduce the actual residual norm to less than
\Gamma6 of the original residual norm, with a zero initial guess. Several values of k are used. A
dagger (y) in the tables indicates that there was no convergence within 500 matrix-vector
products. The incomplete Choleski factorization IC(0) of the Re=0 problem was used as the
preconditioner in all the problems.
For comparison, we first show in Table 1, the results using standard right-preconditioning.

Table

2 shows the results using right-preconditioning with M \Gamma1 inner products, or equiva-
lently, split preconditioning. In the latter case, care was taken to assure that GMRES did
not stop before the actual residual norm was within twice the tolerance. For DQGMRES,
since an accurate residual norm estimate is not available within the algorithm, the exact
residual norm was computed and used for the stopping criterion for the purpose of this com-
parison. The right-preconditioned methods have a slight advantage in this comparison (by
as many as 20 Mat-Vec's), since they directly minimize the actual residual norm, whereas
the symmetrically preconditioned methods minimize a preconditioned residual norm.
6 214 128 73 446 361 y 94 97 258 197 100 94

Table

1: Mat-Vec's for convergence for right-preconditioned methods.

Table

2: Mat-Vec's for convergence for symmetric right-preconditioned methods.
The results in Table 1 show the irregular performance of DQGMRES(k) for these small
values of k when the preconditioned system is not symmetric. The performance is entirely
regular in Table 2, where the preconditioned system is near symmetric. For Reynolds numbers
up to 3, the systems are sufficiently symmetric so that DQGMRES(2) behaves the same
as DQGMRES with much larger k. The performance remains regular until beyond Reynolds
number 7, when the number of steps to convergence begins to become irregular, like in the
right-preconditioned case.
GMRES with either right or symmetric preconditioning does not show any marked difference
in performance; apparently the symmetry of the preconditioned system is not as
essential here for this problem. However, the results do show that DQGMRES(k) with small
values of k may perform as well, in terms of number of steps, as GMRES(k) with large values
of k, particularly if near-symmetry is preserved. Since the former is much more efficient, the
combination of preserving symmetry and truncated iterative methods may result in a much
more economical method, as well as the more regular behavior shown above.
We also performed the same experiments with orthogonal projection methods, namely
the Full Orthogonalization Method (FOM) and its truncated variant, the Direct Incomplete
Orthogonalization Method (DIOM) [7]. The results were very similar to the results above,
and are not shown here. Indeed, the development of the algorithms and the theory above is
identical for these methods.
For interest, we also performed tests where an ILU(0) preconditioner was constructed for
each matrix and compared right and split preconditioning. For the near-symmetric systems
here, there was very little difference in these results compared to using IC(0) constructed
from the Re=0 case for all the matrices. Thus the deterioration in performance as the
Reynolds number increases is not entirely due to a relatively less accurate preconditioner,
but is more due to the increased nonsymmetry and non-normality of the matrices. Although
the eigenvalues of the preconditioned matrices are identical, their eigenvectors and
hence their degree of non-normality may change completely. Unfortunately, it is difficult to
quantitatively relate non-normality and convergence.
5.2 Breakdown behavior of Bi-CG
To test the breakdown behavior of Bi-CG, MATLAB was used to generate random matrices of
order 300 with approximately 50 percent normally distributed nonzero entries. The matrices
were adjusted so that
A /
A
i.e., the symmetric part was shifted so that the lowest eigenvalue was 10 \Gamma5 and then " times
the skew-symmetric part was added back. The parameter " was altered to get varying degrees
of nonsymmetry.
For each " that we tested, 100 matrices were generated, and the smallest value of the
cosines corresponding to the denominators in the algorithms were recorded. In the right-
preconditioned case, we recorded the minimum of
for all j, and for the symmetric right-preconditioned case, we recorded the minimum of
for all j. The relative residual norm reduction was 10 \Gamma9 when the iterations were stopped.
The initial guesses were 0, and r
was set to r 0 . IC(0) of the symmetric part was used as
the preconditioner.

Table

3 shows the frequencies of the size of minimum cosines for the right-preconditioned
(first row of each pair of rows) and the symmetrically-preconditioned cases (second row of
each pair of rows). For example, all 100 minimum cosines were between
the symmetrically-preconditioned case. The average number of Bi-CG steps and the average
minimum cosine is also shown. The last column, labeled 'better', shows the number of times
that the minimum cosine was higher in the improved algorithm.
The Table shows that the right-preconditioned algorithm can produce much smaller
cosines, indicating a greater probability for breakdown. The difference between the algorithms
is less as the degree of nonsymmetry is increased. For there is almost no
difference in the breakdown behavior of the algorithms. The Table shows that the number of
Bi-CG steps is not significantly reduced in the new algorithm, nor is the average minimum
cosine of the modified algorithm significantly increased. It is the probability that a small
cosine is not encountered that is better.
It is important to note that this behavior only applies when r
0 is set to r 0 . When r
0 is
chosen randomly, there is no gain in the symmetrically-preconditioned algorithm, as shown
in

Table

4.

Table

5 shows the number of steps and the minimum cosines for the two algorithms
applied to the driven cavity problem described in Section 5.1 above. Figure 1 shows a plot of
the minimum cosines as the two algorithms progress for the Note that the
minimum cosines are higher and much smoother in the symmetrically-preconditioned case.
In the Re = 7 problem, the cosines are still higher, but the smoothness is lost.
6 Conclusions
When solving linear systems with matrices that are very close to being symmetric, this paper
has shown that it is possible to improve upon the standard practice of using a (nonsymmetric)
preconditioner for that matrix along with a left- or right-preconditioned iterative method.
The original degree of symmetry may be maintained by using a symmetric preconditioner
and an alternate inner product (or split preconditioning, if appropriate). By combining this
idea with truncated iterative methods, solution procedures that converge more quickly and
require less storage are developed. The truncated methods also seem to become more robust
with the truncation parameter k when near-symmetry is maintained. The Bi-CG algorithm
also seems to be more robust with respect to serious breakdown when near-symmetry is
maintained.
" steps 3e-6 1e-5 3e-5 1e-4 3e-4 1e-3 3e-3 1e-2 3e-2 average better
(0.) 31.77 100 1.87 74

Table

3: Frequencies of minimum cosines for right-preconditioned (first row of each pair of
rows) and symmetrically-preconditioned (second row of each pair of rows) Bi-CG.
" steps 3e-6 1e-5 3e-5 1e-4 3e-4 1e-3 3e-3 1e-2 3e-2 average better

Table

4: Frequencies of minimum cosines when r  is chosen randomly.
Bi-CG steps min cosines
Re. right symm right symm
3 73 72 2.02e-4 9.07e-3

Table

5: Steps and minimum cosines for the driven cavity problem.

Figure

1: Minimum cosines in right-preconditioned Bi-CG (solid line) and symmetrically-
preconditioned Bi-CG (dashed line) for the

Acknowledgments

The first and third authors wish to acknowledge the support of RIACS, NASA Ames under
Contract NAS 2-13721, where this collaboration originated. The second and third authors
also wish to acknowledge the support of the Minnesota Supercomputer Institute.



--R

A taxonomy for conjugate gradient methods.
Conjugate gradient type methods for unsymmetric and inconsistent systems of linear equations.

An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix
The concept of special inner products for deriving new conjugate gradient-like solvers for non-symmetric sparse linear systems
A flexible inner-outer preconditioned GMRES algorithm
Iterative Methods for Sparse Linear Systems.
DQGMRES: a direct quasi-minimal residual algorithm based on incomplete orthogonalization
Generalized conjugate-gradient acceleration of nonsym- metrizable iterative methods
--TR
preconditioned iterative methods;inner-products;incomplete orthogonalization