f04mbf
f04mbf
© Numerical Algorithms Group, 2002.
Purpose
F04MBF Real sparse symmetric simultaneous linear equations
Synopsis
[x,itn,anorm,acond,rnorm,xnorm,inform,ifail] = f04mbf(b,aprod,msolve,...
precon<,shift,msglvl,rtol,itnlim,ifail>)
Description
F04MBF solves the system of linear equations
(A-(lambda)I)x=b (3.1)
where A is an n by n sparse symmetric matrix and (lambda) is a
scalar, which is of course zero if the solution of the equations
Ax=b
is required. It should be noted that neither A nor (A-(lambda)I)
need be positive-definite.
(lambda) is supplied as the parameter SHIFT, and allows F04MBF to
be used for finding eigenvectors of A in methods such as Rayleigh
quotient iteration, in which case (lambda) will be an
approximation to an eigenvalue of A and b an approximation to an
eigenvector of A.
The routine also provides an option to allow pre-conditioning and
this will often reduce the number of iterations required by
F04MBF.
F04MBF solves the equations by an algorithm based upon the
Lanczos process. The routine does not require A explicitly, but A
is specified via a user-supplied routine APROD which, given an n
element vector c, must return the vector z given by
z=Ac.
The pre-conditioning option is based on the following reasoning.
If A can be expressed in the form
A=I+B
where B is of rank (rho), then the Lanczos process converges (in
exact arithmetic) in at most (rho) iterations. If more generally
A can be expressed in the form
A=M+C
where M is symmetric positive-definite and C has rank (rho), then
-(1/2) -(1/2) -(1/2) -(1/2)
M AM =I+M CM
-(1/2) -(1/2)
and M AM also has rank (rho), and the Lanczos process
-(1/2) -(1/2)
applied to M AM would again converge in at most (rho)
iterations. On a computer, the number of iterations may be
greater than (rho), but the Lanczos process may still be expected
-(1/2) -(1/2)
to converge rapidly. F04MBF does not require M AM to
be formed explicitly, but implicitly solves the equations
-(1/2) -(1/2) -(1/2) 1/2
M (A-(lambda)I)M y=M b , y=M x (3.2)
with the user being required to supply a routine MSOLVE to solve
the equations
Mz=c. (3.3)
For the pre-conditioning option to be effective, it is desirable
that equations (3.3) can be solved efficiently.
If we let r denote the residual vector
r=b-(A-(lambda)I)x
corresponding to an iterate x, then, when pre-conditioning has
not been requested, the iterative procedure is terminated if it
is estimated that
||r||<=tol.||A-(lambda)I||.||x||, (3.4)
where tol is a user-supplied tolerance, ||r|| denotes the
Euclidean length of the vector r and ||A|| denotes the Frobenius
(Euclidean) norm of the matrix A. When pre-conditioning has been
requested, the iterative procedure is terminated if it is
estimated that
-(1/2) -(1/2) -(1/2) 1/2
||M r||<=tol.||M (A-(lambda)I)M ||.||M x||. (3.5)
Note that
-(1/2) -(1/2) -(1/2) -(1/2) 1/2
M r=(M b)-M (A-(lambda)I)M (M x)
-(1/2)
so that M r is the residual vector corresponding to equation
(3.2). The routine will also terminate if it is estimated that
||A-(lambda)I||.||x||>=||b||/(epsilon), (3.6)
where (epsilon) is the machine precision, when pre-conditioning
has not been requested; or if it is estimated that
-(1/2) -(1/2) 1/2 -(1/2)
||M (A-(lambda)I)M ||.||M x||>=||M b||/(epsilon)
(3.7)
when pre-conditioning has been requested. If (3.6) is satisfied
then x is almost certainly an eigenvector of A corresponding to
the eigenvalue (lambda). If (lambda) was set to 0 (for the
solution of Ax=b), then this condition simply means that A is
effectively singular.
Parameters
f04mbf
Required Input Arguments:
b (:) real
aprod function (User-Supplied)
msolve function (User-Supplied)
precon logical
Optional Input Arguments: <Default>
shift real 0
msglvl integer 0
rtol real eps
itnlim integer length(b)
ifail integer -1
Output Arguments:
x (:) real
itn integer
anorm real
acond real
rnorm real
xnorm real
inform integer
ifail integer