This paper was written for The European Centre for Medium Range Weather Forecasting Workshop on Developments in Numerical Methods for Very High Resolution Global Models, June 2000

Wathen, A. J.

English

Oxford University Research Archive (ORA)

Preconditioned iterative methods for implicit equations

Mathematical Institute Eprints Archive

Preonditioned iterative methods for impliit
equations
A. J. Wathen
Oxford University Computing Laboratory
These notes are about methods for nding the solution x of the linear
simultaneous equations
Ax = b (1)
where A is a known real or omplex nnmatrix (whih is assumed invertible)
and b is a known n-vetor.
The oeÆient matrix would generally be onsidered to be either full,
banded or sparse. For any of these types of matries, good numerial al-
gorithms based upon systemati (Gaussian) elimination exist. (For general
software, you need look no further than the NAG library. Hardware-spei
routines exist on many advaned arhiteture mahines.) Suh algorithms
terminate with (in general) an aurate numerial answer after a preditable
ammount of omputation. The basi diÆulty is that for large n suh pre-
ditable requirements in terms of CPU time, omputer memory and dis
storage may be ompletely unavailable whether you work on a workstation
or a superomputer. (See for example the table of work estimates on page 388
of [2℄). However, iterative solution methods are often available and pratia-
ble for suh large problems espeially when the oeÆient matrix is sparse.
In the ontext of dynami simulations suh as arise in weather predition
(as well as other areas inluding petroleum reservoir modelling, semion-
dutor proess modelling and nanial option priing) suh systems usually
derive from an impliit time disretisation method, the simplest of whih
would be the bakwards Euler method. Linearisation via Newtons method
or Piard iteration might also be involved. In suh impliit time-dependent
omputations the oeÆient matrix always has the form
1
t
I + L or I +tL
1
in the ase of nite dierene replaement of the spatial terms or
1
t
^
M + L or
^
M +tL
where
^
M is the `mass matrix' in the ase of nite element approximation
in the spatial variables - the appearane of
^
M rather than I is of no great
signiane sine it is known that simple diagonal saling of
^
M renders the
solution of systems involving
^
M trivial sine diag(
^
M)
 1
^
M is (spetrally) very
lose to the identity (see [23℄). For approximations based on expansions
in terms of globally-dened funtions (suh as arise in spetral methods),
the orthogonality properties of the basis funtions are important. Here L
represents the spatial dierential opertor terms.
This is a simple truth, but from the point of view of solving the matrix
system at eah time step it indiates that the hoie of time-step t an
have a profound eet: if t is small then the oeÆient matrix looks more
and more like the identity matrix whereas larger values of t make the
oeÆient matrix loser to a steady-state operator. For diret (elimination)
algorithms the hoie of t does not aet the work in solving suh systems
(only the auray of the omputed results) however for iterative methods,
solution with small t generally takes very few iterations; many iterations
might be required in the ase of a larger time-step. Of ourse the hoie
of timestep is usually determined by the time sale of a simulation and any
auray/stability onstraints.
The two leading lasses of iterative methods are those based on the Multi-
grid approah, and the Conjugate Gradient method and its generalisations.
The multigrid method was originally derived for grid-based problems suh
as arise from nite dierene or nite element disretisation of partial dier-
ential equations. For ertain model problems extremely rapid onvergene
rates are proved and observed with this method. For problems with no ob-
vious geometri sequene of (multi-)grids, there has been attention direted
towards nding a mimiking algebrai proedure (Algebrai Multigrid), but
suh methods are generally less rapidly onvergent. Muh researh on Multi-
grid methods ontinues - I am not an expert on it! A general referene for
Multigrid methods is [14℄.
Methods of Conjugate Gradient (CG) type (so-alled Krylov subspae
methods) are however more widely appliable, and onvergene rates in many
ases are fast. It is methods of this type that I mainly intend to review in
these notes.
2
In all the methods of this type, referene to the oeÆient matrix, A
is only required for the omputation of matrixvetor produts, Ap (and
sometimes also A
T
p) where p is a known n-vetor. For a sparse matrix these
produts an usually be enoded so that only alulation on the non-zeros is
required.
In many appliations, preonditioning of the original linear system is em-
ployed. Indeed after an explosion of new Krylov subspae methods in the
early 1990's it was realised that nding and eetive preonditioner is muh
more important in general than the spei hoie of iterative method. The
idea here is that if a preonditioning matrix M an be easily omputed suh
that onvergene of the appropriate iterative method is faster for the preon-
ditioned system
M
 1
Ax = M
 1
b (2)
than the original, AND systems of equations withM as oeÆient matrix are
very rapidly solvable (suh a solution will be required at eah iteration), then
the overall solution time will be shorter. The inverse ofM is not required; the
equation (2) is only shematially equivalent to the preonditioned methods
that are used. Indeed a good preonditioner might be a `proess' suh as a
single multigrid yle whih only for theoretial purposes an be onsidered
as a matrix.
A simple but not neessary motivation for a good preonditioner is that
it be an easily invertible approximation of the oeÆient matrix A. Alterna-
tively one may have an approximate inverse so that preonditioning involves
only a matrixvetor multipliation at eah iteration.
There are many approahes to the onstrution of preonditioners. For
partiular appliations problems, there is sometimes a spei and inherently
simpler approximation whih `aptures the main physis' and is more read-
ily invertible whih may be employed as a preonditioner. For example one
might have omputed a triangular fatorisation but need solution of many
linear systems with slightly dierent oeÆient matries - the original fa-
torisation may then be used as a good preonditioner. This type of situation
might arise for example in nonlinear/time-dependent problems.
For problems with partiular mathematial struture, there are often ways
to utilise that struture for preonditioning - for example for problems with
onstraints there have been some very suessful reent advanes whih are
based on approximation of subproblems (see [15℄,[12℄).
For self-adjoint and positive denite ellipti partial dierential equation
3
problems, suessful preonditioners based on domain deomposition are also
onveniently parallelisable. There is muh work on these methods - see
the sequene of proeedings of whih the latest is an be obtained from
http://www.gre.a.uk/ dd11/dd11/proeedings.html.
Popular and general preonditioners for sparse systems an be onstruted
using inomplete triangular fatorisation of A ([13℄,[9℄) or by (parallel) on-
strution of sparse approximate inverses ([3℄,[10℄). The quality of these pre-
onditioners varies when employed for dierent appliation problems and
it an be expensive (even prohibitive) to onstrut them for very large sale
problems. Less generally appliable preonditioners have also been suggested
whih work well on ertain lasses of problems (see [25℄,[11℄,[21℄) or the spe-
ial issue of Numerial Linear Algebra and Appliations assoiated with the
onferene on Preonditioned held last summer in Minneapolis (see [16℄).
We now turn to iterative methods.
For real symmetri (or omplex hermitian) and positive denite systems,
the (Hestenes-Steifel) Conjugate Gradient method ([7℄, page 370) is to set
initial vetors x
0
= 0; p
0
= 0, and for eah iteration
ompute residual r
k 1
= b  Ax
k 1
ompute onjugate searh diretion p
k
= r
k 1
+ 
k
p
k 1
with 
k
=
r
T
k 1
r
k 1
r
T
k 2
r
k 2
update solution x
k
= x
k 1
+ 
k
p
k
with 
k
=
r
T
k 1
r
k 1
p
T
k
Ap
k
:
The residual is usually omputed reursively as
r
k 1
= r
k 2
  
k 1
Ap
k 1
so that only one matrixvetor multipliation is required per iteration.
The most readily obtained onvergene bound for this method is that the
exat solution x and the CG iterates x
k
satisfy
kx
k
  xk
A
 2
 
1 
p

1 +
p

!
k
kx
0
  xk
A
(3)
where  is the spetral ondition number of A (the ratio of largest to smallest
eigenvalues of A), and kyk
2
A
= y
T
Ay. Thus rapid onvergene is guarenteed if
4
 is small. This is the motivation for preonditioning: to redue . If a better
guess for x
0
is available (for example the solution at a previous time-step),
then this an be used and (3) indiates that few iterations are expeted to
rene the solution. Contrast this with diret solution proedures whih take
the same time whatever is known about the solution.
For real symmetri (or omplex hermitian) indenite matries and other
matries whose eigenvalues all lie on a single line in the omplex plane, there
exist Conjugate Gradient-like methods (for example the method of Conju-
gate Residuals of whih probably the best of the various implementations is
MINRES [17℄) with similar onvergene properties but whih require a little
more work per iteration. See the `taxonomy' [1℄. Similarly to Conjugate
Gradients, rapid onvergene an be guarenteed for problems with lusted
eigenvalues - the preise theoretial onvergene bounds/estimates are more
ompliated, see [8℄,[24℄.
For non-symmetri matries, there are a number of ompetitive CG-like
methods. None has the monotoni error- or residual-redution and onstant
work per iteration like in the symmetri ase, and onvergene an be more
errati. Available onvergene bounds are generally not desriptive in on-
trast to symmetri problems. Numerial divergene is oasionally seen.
Nevertheless, a number of methods of this type are signiantly used in a
variety of situations.
The most longstanding of the non-symmetri methods is that based on
(Hestenes-Steifel) CG for the normal equations
A
T
Ax = A
T
b;
the oeÆient matrix of whih is symmetri and positive denite. This ap-
proah is not reommended unless (A) is small as (A
T
A) = (A)
2
and
slow onvergene as suggested by (3) is usually experiened. There are situ-
ations (suh as impliit omputations with a very small time-step or where
a high-quality preonditioner is available) where this method might be used.
The GMRES method (Generalised Minimal Residual method) ([19℄) the-
oretially is error-minimising at eah iteration, however a Gram-Shmidt-like
(in fat Arnoldi) orthogonalisation of the k
th
`searh diretion' p
k
against all
previous searh diretion vetors p
0
; p
1
; : : : ; p
k 1
is required at the k
th
itera-
tion - thus the work per iteration inreases as the number of iterations grows.
This is generally impratial in both storage as well as omputation. The
method may however be used in a yli manner: do l steps of GMRES
5
(storing and operating on l vetors) and then restart. This method is alled
GMRES(l). A popular value for l is about 10. The advantage of this method
is that onvergene an always be guarenteed by inreasing l. The disadvan-
tage is that too large a value for l may be required! Indeed, reent examples
provided by Embree ([4℄) indiate that it is quite possible for GMRES(l
1
)
to require fewer iterations and less work to obtain a solution to any onver-
gene tolerane than GMRES(l
2
) when l
1
< l
2
. Thus inreasing the restart
length may not be protable unless it is possible to employ the non-restarted
method GMRES(n).
Restarting is generally not favoured if it an be helped, and in appliations
suh as petroleum reservoir modelling where this method is widely employed
(usually an implementation alled Orthomin), the time-step is often kept
small enough that few enough iterations are required so that restating is not
neessary. In omputations I have been involved with, it is not lear that
the adaptive time-step seletion mehanism is not in fat responding to the
onvergene rate of the iterative linear solver rather than any disretisation
errors or non-onvergene of the (outer) non-linear iteration.
The BiCG method ([5℄) requires only 2 or 3 p-vetors, but omputation
of both of the produts Ap and A
T
p is needed. (The latter may be awkward
to ompute depending on the data struture for storing A). The method an
also stall and has been generally ast aside in favour of the modied methods
below.
The CGS (CG Squared) method ([20℄) is a variant of BiCG whih avoids
the alulation of the transpose produt A
T
p. It exhibits rapid nal on-
vergene, though large jumps (inreases) in the residual are usually seen at
earlier stages of the iteration.
A further variant of BiCG is the BiCG-STAB method ([22℄) whih is
similar to CGS, but has a `more monotoni' residual redution. A useful
hybrid ombination of a few steps of GMRES with BiCG-STAB leads to the
useful BiCG-STAB(l) method ([18℄). Both CGS and BiCG-STAB and its
variants are widely used in semiondutor devie simulations.
The QMR (Quasi Minimal Residual) method ([6℄) is related to both BiCG
and GMRES looked a very attrative approah when it was announed but
despite the availability of good software, it seems not to be ommonly used
exept for ertain model redution problems in iruit simulation.
The GMRES(l), CGS and BiCG-STAB(l) methods with appropriate pre-
onditioning are all signiantly used in a number of pratial situations not
indiated here.
6
Referenes
[1℄ S.F. Ashby, T.A. Manteuffel and P.E. Saylor, A taxonomy
for onjugate gradient methods, SIAM J. Numer. Anal. 27 (1990), pp.
1542-1568.
[2℄ O. Axelsson and V.A. Barker, Finite Element Solution of Bound-
ary Value Problems: Theory and Computation, Aademi Press, New
York, 1984.
[3℄ M. Benzi and M. Tuma, A sparse approximate inverse preonditioner
for nonsymmetri linear systems, SIAM J. Si. Comput. 19 (1998), pp.
968-994.
[4℄ M. Embree, Convergene of Krylov subspae methods for non-normal
matries D. Phil thesis, Oxford University, 1999.
[5℄ R. Flether, Conjugate Gradient Methods for Indenite Systems, in
Numerial Analysis, Dundee 1975, G.AWatson, ed. Springer, New York,
1976.
[6℄ R.W. Freund and N.M. Nahtigal, QMR: a quasi-minimal resid-
ual method for non-hermitian linear systems, Numer. Math. 60 (1991),
pp. 315-339.
[7℄ G.H. Golub and C.F. Van Loan, Matrix Computations, John Hop-
kins University Press, Baltimore, 1983.
[8℄ A. Greenbaum, Iterative methods for solving linear systems, SIAM,
Philadelphia, 1997.
[9℄ I. Gustafsson, A Class of First Order Fatorization Methods, BIT 18
(1978), pp. 142-156.
[10℄ M. J. Grote and T. Hukle, Parallel preonditioning with sparse
approximate inverses, SIAM J. Si. Comput. 18 (1997), pp. 838-853.
[11℄ T.J.R Hughes, I. Levit and J. Winget, An Element-by-Element
Solution Algorithm for Problems of Strutural and Solid Mehanis,
Comput. Meths. Appl. Meh. Engrg 36 (1983), pp. 241-254.
7
[12℄ C. Keller, N.I.M. Gould, and A.J. Wathen, Constraint preondi-
tioning for indenite linear systems, SIAM J. Matrix Anal. Appl. 21(4)
(2000), pp. 1300-1317.
[13℄ J.A. Meijerink and H.A. Van der Vorst, An Iterative Solution
Method for Linear Systems of whih the CoeÆient Matrix is a Sym-
metri M-Matrix, Math. Comput. 31 (1977), pp. 148-162.
[14℄ S.F. MCormik, Multigrid Methods, SIAM, Philadelphia, 1987.
[15℄ M.F. Murphy, G.H. Golub and A.J. Wathen, 2000, A note on
preonditioning for indenite linear systems, SIAM J. Si. Comput., to
appear.
[16℄ Numer. Linear Algebra Appl. (Autumn 2000 issue).
[17℄ C.C. Paige and M.A. Saunders, Solution of Sparse Indenite Sys-
tems of Linear Equations, SIAM J. Numer. Anal. 12 (1975), pp. 617-629.
[18℄ G.L.G. Sleijpen and D.R. Fokkema, BiCGstab(L) for linear equa-
tions involving unsymmetri matries with omplex spetrum, Eletroni
Trans. Numeri. Anal.,1 (1993), pp. 11-32.
[19℄ Y. Saad and M.H. Shultz, GMRES: A generalised minimal residual
algorithm for solving nonsymmetri linear systems, SIAM J. Si. Stat.
Comput. 7 (1986), pp. 856-869.
[20℄ P. Sonneveld, CGS: a fast Lanzos-type solver for non-symmetri
linear systems, SIAM J. Si. Stat. Comput. 10 (1989), pp 36-52.
[21℄ G. Strang, A proposal for Toeplitz matrix alulation, Studies Appl.
Math. 74 (1986), pp. 171-176.
[22℄ H.A. van der Vorst, Bi-CGSTAB: A fast and smoothly onvergent
variant of BiCG for the solution of non-symmetri linear systems, SIAM
J. Si. Stat. Comput. 13 (1992), pp. 631-644.
[23℄ A.J. Wathen, Realisti eigenvalue bounds for the Galerkin Mass Ma-
trix, IMA J. Numer. Anal. 7 (1987), pp. 449-457.
8
[24℄ A.J. Wathen, B. Fisher and D.J. Silvester, The onvergene of
iterative solution methods for symmetri and indenite linear systems,
in `Numerial Analysis 1997', Eds. D.F. GriÆths & G.A. Watson, Pit-
man Researh Notes in Mathematis Series, Addison Wesley Longman,
Harlow, England, 1998, pp. 230-243.
[25℄ H. Yserentant, On the Multi-level Splitting of Finite Element Spaes,
Numer. Math. 49 (1986), pp. 379-412.
9


Oxford University Research Archive

Preonditioned iterative methods for impliitequationsA. J. WathenOxford University Computing LaboratoryThese notes are about methods for nding the solution x of the linearsimultaneous equations Ax = b (1)where A is a known real or omplex nnmatrix (whih is assumed invertible)and b is a known n-vetor.The oeÆient matrix would generally be onsidered to be either full,banded or sparse. For any of these types of matries, good numerial al-gorithms based upon systemati (Gaussian) elimination exist. (For generalsoftware, you need look no further than the NAG library. Hardware-speiroutines exist on many advaned arhiteture mahines.) Suh algorithmsterminate with (in general) an aurate numerial answer after a preditableammount of omputation. The basi diÆulty is that for large n suh pre-ditable requirements in terms of CPU time, omputer memory and disstorage may be ompletely unavailable whether you work on a workstationor a superomputer. (See for example the table of work estimates on page 388of [2℄). However, iterative solution methods are often available and pratia-ble for suh large problems espeially when the oeÆient matrix is sparse.In the ontext of dynami simulations suh as arise in weather predition(as well as other areas inluding petroleum reservoir modelling, semion-dutor proess modelling and nanial option priing) suh systems usuallyderive from an impliit time disretisation method, the simplest of whihwould be the bakwards Euler method. Linearisation via Newtons methodor Piard iteration might also be involved. In suh impliit time-dependentomputations the oeÆient matrix always has the form1tI + L or I +tL1in the ase of nite dierene replaement of the spatial terms or1tM̂ + L or M̂ +tLwhere M̂ is the `mass matrix' in the ase of nite element approximationin the spatial variables - the appearane of M̂ rather than I is of no greatsigniane sine it is known that simple diagonal saling of M̂ renders thesolution of systems involving M̂ trivial sine diag(M̂) 1M̂ is (spetrally) verylose to the identity (see [23℄). For approximations based on expansionsin terms of globally-dened funtions (suh as arise in spetral methods),the orthogonality properties of the basis funtions are important. Here Lrepresents the spatial dierential opertor terms.This is a simple truth, but from the point of view of solving the matrixsystem at eah time step it indiates that the hoie of time-step t anhave a profound eet: if t is small then the oeÆient matrix looks moreand more like the identity matrix whereas larger values of t make theoeÆient matrix loser to a steady-state operator. For diret (elimination)algorithms the hoie of t does not aet the work in solving suh systems(only the auray of the omputed results) however for iterative methods,solution with small t generally takes very few iterations; many iterationsmight be required in the ase of a larger time-step. Of ourse the hoieof timestep is usually determined by the time sale of a simulation and anyauray/stability onstraints.The two leading lasses of iterative methods are those based on the Multi-grid approah, and the Conjugate Gradient method and its generalisations.The multigrid method was originally derived for grid-based problems suhas arise from nite dierene or nite element disretisation of partial dier-ential equations. For ertain model problems extremely rapid onvergenerates are proved and observed with this method. For problems with no ob-vious geometri sequene of (multi-)grids, there has been attention diretedtowards nding a mimiking algebrai proedure (Algebrai Multigrid), butsuh methods are generally less rapidly onvergent. Muh researh on Multi-grid methods ontinues - I am not an expert on it! A general referene forMultigrid methods is [14℄.Methods of Conjugate Gradient (CG) type (so-alled Krylov subspaemethods) are however more widely appliable, and onvergene rates in manyases are fast. It is methods of this type that I mainly intend to review inthese notes. 2In all the methods of this type, referene to the oeÆient matrix, Ais only required for the omputation of matrixvetor produts, Ap (andsometimes also ATp) where p is a known n-vetor. For a sparse matrix theseproduts an usually be enoded so that only alulation on the non-zeros isrequired.In many appliations, preonditioning of the original linear system is em-ployed. Indeed after an explosion of new Krylov subspae methods in theearly 1990's it was realised that nding and eetive preonditioner is muhmore important in general than the spei hoie of iterative method. Theidea here is that if a preonditioning matrix M an be easily omputed suhthat onvergene of the appropriate iterative method is faster for the preon-ditioned system M 1Ax = M 1b (2)than the original, AND systems of equations withM as oeÆient matrix arevery rapidly solvable (suh a solution will be required at eah iteration), thenthe overall solution time will be shorter. The inverse ofM is not required; theequation (2) is only shematially equivalent to the preonditioned methodsthat are used. Indeed a good preonditioner might be a `proess' suh as asingle multigrid yle whih only for theoretial purposes an be onsideredas a matrix.A simple but not neessary motivation for a good preonditioner is thatit be an easily invertible approximation of the oeÆient matrix A. Alterna-tively one may have an approximate inverse so that preonditioning involvesonly a matrixvetor multipliation at eah iteration.There are many approahes to the onstrution of preonditioners. Forpartiular appliations problems, there is sometimes a spei and inherentlysimpler approximation whih `aptures the main physis' and is more read-ily invertible whih may be employed as a preonditioner. For example onemight have omputed a triangular fatorisation but need solution of manylinear systems with slightly dierent oeÆient matries - the original fa-torisation may then be used as a good preonditioner. This type of situationmight arise for example in nonlinear/time-dependent problems.For problems with partiular mathematial struture, there are often waysto utilise that struture for preonditioning - for example for problems withonstraints there have been some very suessful reent advanes whih arebased on approximation of subproblems (see [15℄,[12℄).For self-adjoint and positive denite ellipti partial dierential equation3problems, suessful preonditioners based on domain deomposition are alsoonveniently parallelisable. There is muh work on these methods - seethe sequene of proeedings of whih the latest is an be obtained fromhttp://www.gre.a.uk/ dd11/dd11/proeedings.html.Popular and general preonditioners for sparse systems an be onstrutedusing inomplete triangular fatorisation of A ([13℄,[9℄) or by (parallel) on-strution of sparse approximate inverses ([3℄,[10℄). The quality of these pre-onditioners varies when employed for dierent appliation problems andit an be expensive (even prohibitive) to onstrut them for very large saleproblems. Less generally appliable preonditioners have also been suggestedwhih work well on ertain lasses of problems (see [25℄,[11℄,[21℄) or the spe-ial issue of Numerial Linear Algebra and Appliations assoiated with theonferene on Preonditioned held last summer in Minneapolis (see [16℄).We now turn to iterative methods.For real symmetri (or omplex hermitian) and positive denite systems,the (Hestenes-Steifel) Conjugate Gradient method ([7℄, page 370) is to setinitial vetors x0 = 0; p0 = 0, and for eah iterationompute residual rk 1 = b  Axk 1ompute onjugate searh diretion pk = rk 1 + kpk 1with k = rTk 1rk 1rTk 2rk 2update solution xk = xk 1 + kpkwith k = rTk 1rk 1pTkApk :The residual is usually omputed reursively asrk 1 = rk 2   k 1Apk 1so that only one matrixvetor multipliation is required per iteration.The most readily obtained onvergene bound for this method is that theexat solution x and the CG iterates xk satisfykxk   xkA  2  1 p1 +p!k kx0   xkA (3)where  is the spetral ondition number of A (the ratio of largest to smallesteigenvalues of A), and kyk2A = yTAy. Thus rapid onvergene is guarenteed if4 is small. This is the motivation for preonditioning: to redue . If a betterguess for x0 is available (for example the solution at a previous time-step),then this an be used and (3) indiates that few iterations are expeted torene the solution. Contrast this with diret solution proedures whih takethe same time whatever is known about the solution.For real symmetri (or omplex hermitian) indenite matries and othermatries whose eigenvalues all lie on a single line in the omplex plane, thereexist Conjugate Gradient-like methods (for example the method of Conju-gate Residuals of whih probably the best of the various implementations isMINRES [17℄) with similar onvergene properties but whih require a littlemore work per iteration. See the `taxonomy' [1℄. Similarly to ConjugateGradients, rapid onvergene an be guarenteed for problems with lustedeigenvalues - the preise theoretial onvergene bounds/estimates are moreompliated, see [8℄,[24℄.For non-symmetri matries, there are a number of ompetitive CG-likemethods. None has the monotoni error- or residual-redution and onstantwork per iteration like in the symmetri ase, and onvergene an be moreerrati. Available onvergene bounds are generally not desriptive in on-trast to symmetri problems. Numerial divergene is oasionally seen.Nevertheless, a number of methods of this type are signiantly used in avariety of situations.The most longstanding of the non-symmetri methods is that based on(Hestenes-Steifel) CG for the normal equationsATAx = AT b;the oeÆient matrix of whih is symmetri and positive denite. This ap-proah is not reommended unless (A) is small as (ATA) = (A)2 andslow onvergene as suggested by (3) is usually experiened. There are situ-ations (suh as impliit omputations with a very small time-step or wherea high-quality preonditioner is available) where this method might be used.The GMRES method (Generalised Minimal Residual method) ([19℄) the-oretially is error-minimising at eah iteration, however a Gram-Shmidt-like(in fat Arnoldi) orthogonalisation of the kth `searh diretion' pk against allprevious searh diretion vetors p0; p1; : : : ; pk 1 is required at the kth itera-tion - thus the work per iteration inreases as the number of iterations grows.This is generally impratial in both storage as well as omputation. Themethod may however be used in a yli manner: do l steps of GMRES5(storing and operating on l vetors) and then restart. This method is alledGMRES(l). A popular value for l is about 10. The advantage of this methodis that onvergene an always be guarenteed by inreasing l. The disadvan-tage is that too large a value for l may be required! Indeed, reent examplesprovided by Embree ([4℄) indiate that it is quite possible for GMRES(l1)to require fewer iterations and less work to obtain a solution to any onver-gene tolerane than GMRES(l2) when l1 < l2. Thus inreasing the restartlength may not be protable unless it is possible to employ the non-restartedmethod GMRES(n).Restarting is generally not favoured if it an be helped, and in appliationssuh as petroleum reservoir modelling where this method is widely employed(usually an implementation alled Orthomin), the time-step is often keptsmall enough that few enough iterations are required so that restating is notneessary. In omputations I have been involved with, it is not lear thatthe adaptive time-step seletion mehanism is not in fat responding to theonvergene rate of the iterative linear solver rather than any disretisationerrors or non-onvergene of the (outer) non-linear iteration.The BiCG method ([5℄) requires only 2 or 3 p-vetors, but omputationof both of the produts Ap and ATp is needed. (The latter may be awkwardto ompute depending on the data struture for storing A). The method analso stall and has been generally ast aside in favour of the modied methodsbelow.The CGS (CG Squared) method ([20℄) is a variant of BiCG whih avoidsthe alulation of the transpose produt ATp. It exhibits rapid nal on-vergene, though large jumps (inreases) in the residual are usually seen atearlier stages of the iteration.A further variant of BiCG is the BiCG-STAB method ([22℄) whih issimilar to CGS, but has a `more monotoni' residual redution. A usefulhybrid ombination of a few steps of GMRES with BiCG-STAB leads to theuseful BiCG-STAB(l) method ([18℄). Both CGS and BiCG-STAB and itsvariants are widely used in semiondutor devie simulations.The QMR (Quasi Minimal Residual) method ([6℄) is related to both BiCGand GMRES looked a very attrative approah when it was announed butdespite the availability of good software, it seems not to be ommonly usedexept for ertain model redution problems in iruit simulation.The GMRES(l), CGS and BiCG-STAB(l) methods with appropriate pre-onditioning are all signiantly used in a number of pratial situations notindiated here. 6Referenes[1℄ S.F. Ashby, T.A. Manteuffel and P.E. Saylor, A taxonomyfor onjugate gradient methods, SIAM J. Numer. Anal. 27 (1990), pp.1542-1568.[2℄ O. Axelsson and V.A. Barker, Finite Element Solution of Bound-ary Value Problems: Theory and Computation, Aademi Press, NewYork, 1984.[3℄ M. Benzi and M. Tuma, A sparse approximate inverse preonditionerfor nonsymmetri linear systems, SIAM J. Si. Comput. 19 (1998), pp.968-994.[4℄ M. Embree, Convergene of Krylov subspae methods for non-normalmatries D. Phil thesis, Oxford University, 1999.[5℄ R. Flether, Conjugate Gradient Methods for Indenite Systems, inNumerial Analysis, Dundee 1975, G.AWatson, ed. Springer, New York,1976.[6℄ R.W. Freund and N.M. Nahtigal, QMR: a quasi-minimal resid-ual method for non-hermitian linear systems, Numer. Math. 60 (1991),pp. 315-339.[7℄ G.H. Golub and C.F. Van Loan, Matrix Computations, John Hop-kins University Press, Baltimore, 1983.[8℄ A. Greenbaum, Iterative methods for solving linear systems, SIAM,Philadelphia, 1997.[9℄ I. Gustafsson, A Class of First Order Fatorization Methods, BIT 18(1978), pp. 142-156.[10℄ M. J. Grote and T. Hukle, Parallel preonditioning with sparseapproximate inverses, SIAM J. Si. Comput. 18 (1997), pp. 838-853.[11℄ T.J.R Hughes, I. Levit and J. Winget, An Element-by-ElementSolution Algorithm for Problems of Strutural and Solid Mehanis,Comput. Meths. Appl. Meh. Engrg 36 (1983), pp. 241-254.7[12℄ C. Keller, N.I.M. Gould, and A.J. Wathen, Constraint preondi-tioning for indenite linear systems, SIAM J. Matrix Anal. Appl. 21(4)(2000), pp. 1300-1317.[13℄ J.A. Meijerink and H.A. Van der Vorst, An Iterative SolutionMethod for Linear Systems of whih the CoeÆient Matrix is a Sym-metri M-Matrix, Math. Comput. 31 (1977), pp. 148-162.[14℄ S.F. MCormik, Multigrid Methods, SIAM, Philadelphia, 1987.[15℄ M.F. Murphy, G.H. Golub and A.J. Wathen, 2000, A note onpreonditioning for indenite linear systems, SIAM J. Si. Comput., toappear.[16℄ Numer. Linear Algebra Appl. (Autumn 2000 issue).[17℄ C.C. Paige and M.A. Saunders, Solution of Sparse Indenite Sys-tems of Linear Equations, SIAM J. Numer. Anal. 12 (1975), pp. 617-629.[18℄ G.L.G. Sleijpen and D.R. Fokkema, BiCGstab(L) for linear equa-tions involving unsymmetri matries with omplex spetrum, EletroniTrans. Numeri. Anal.,1 (1993), pp. 11-32.[19℄ Y. Saad and M.H. Shultz, GMRES: A generalised minimal residualalgorithm for solving nonsymmetri linear systems, SIAM J. Si. Stat.Comput. 7 (1986), pp. 856-869.[20℄ P. Sonneveld, CGS: a fast Lanzos-type solver for non-symmetrilinear systems, SIAM J. Si. Stat. Comput. 10 (1989), pp 36-52.[21℄ G. Strang, A proposal for Toeplitz matrix alulation, Studies Appl.Math. 74 (1986), pp. 171-176.[22℄ H.A. van der Vorst, Bi-CGSTAB: A fast and smoothly onvergentvariant of BiCG for the solution of non-symmetri linear systems, SIAMJ. Si. Stat. Comput. 13 (1992), pp. 631-644.[23℄ A.J. Wathen, Realisti eigenvalue bounds for the Galerkin Mass Ma-trix, IMA J. Numer. Anal. 7 (1987), pp. 449-457.8[24℄ A.J. Wathen, B. Fisher and D.J. Silvester, The onvergene ofiterative solution methods for symmetri and indenite linear systems,in `Numerial Analysis 1997', Eds. D.F. GriÆths & G.A. Watson, Pit-man Researh Notes in Mathematis Series, Addison Wesley Longman,Harlow, England, 1998, pp. 230-243.[25℄ H. Yserentant, On the Multi-level Splitting of Finite Element Spaes,Numer. Math. 49 (1986), pp. 379-412.9

ORA - Oxford University Research Archive

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.11.2936

Preconditioned iterative methods for implicit equations

Abstract

Similar works

Full text

Available Versions

Oxford University Research Archive (ORA)

Mathematical Institute Eprints Archive

Oxford University Research Archive

ORA - Oxford University Research Archive