 |
Millepede-II
V04-00-00_preview
|
|
Millepede-II
V04-00-00_preview
|
MINRES algorithm. More...
Go to the source code of this file.
Functions/Subroutines | |
| subroutine | minres (n, b, r1, r2, v, w, w1, w2, x, y, aprod, msolve, checka, precon, shift, nout, itnlim, rtol, istop, itn, anorm, acond, rnorm, ynorm) |
| Solution of linear equation system. | |
MINRES algorithm.
Definition in file minres.f90.
| subroutine minres | ( | integer | n, |
| double precision, dimension(n) | b, | ||
| double precision, dimension(n) | r1, | ||
| double precision, dimension(n) | r2, | ||
| double precision, dimension(n) | v, | ||
| double precision, dimension(n) | w, | ||
| double precision, dimension(n) | w1, | ||
| double precision, dimension(n) | w2, | ||
| double precision, dimension(n) | x, | ||
| double precision, dimension(n) | y, | ||
| , external | aprod, | ||
| , external | msolve, | ||
| logical | checka, | ||
| logical | precon, | ||
| double precision | shift, | ||
| integer | nout, | ||
| integer | itnlim, | ||
| double precision | rtol, | ||
| integer | istop, | ||
| integer | itn, | ||
| double precision | anorm, | ||
| double precision | acond, | ||
| double precision | rnorm, | ||
| double precision | ynorm | ||
| ) |
Solution of linear equation system.
------------------------------------------------------------------
MINRES is designed to solve the system of linear equations
Ax = b
or the least-squares problem
min || Ax - b ||_2,
where A is an n by n symmetric matrix and b is a given vector.
The matrix A may be indefinite and/or singular.
1. If A is known to be positive definite, the Conjugate Gradient
Method might be preferred, since it requires the same number
of iterations as MINRES but less work per iteration.
2. If A is indefinite but Ax = b is known to have a solution
(e.g. if A is nonsingular), SYMMLQ might be preferred,
since it requires the same number of iterations as MINRES
but slightly less work per iteration.
The matrix A is intended to be large and sparse. It is accessed
by means of a subroutine call of the form
SYMMLQ development:
call Aprod ( n, x, y )
which must return the product y = Ax for any given vector x.
More generally, MINRES is designed to solve the system
(A - shift*I) x = b
or
min || (A - shift*I) x - b ||_2,
where shift is a specified scalar value. Again, the matrix
(A - shift*I) may be indefinite and/or singular.
The work per iteration is very slightly less if shift = 0.
Note: If shift is an approximate eigenvalue of A
and b is an approximate eigenvector, x might prove to be
a better approximate eigenvector, as in the methods of
inverse iteration and/or Rayleigh-quotient iteration.
However, we're not yet sure on that -- it may be better
to use SYMMLQ.
A further option is that of preconditioning, which may reduce
the number of iterations required. If M = C C' is a positive
definite matrix that is known to approximate (A - shift*I)
in some sense, and if systems of the form My = x can be
solved efficiently, the parameters precon and Msolve may be
used (see below). When precon = .true., MINRES will
implicitly solve the system of equations
P (A - shift*I) P' xbar = P b,
i.e. Abar xbar = bbar
where P = C**(-1),
Abar = P (A - shift*I) P',
bbar = P b,
and return the solution x = P' xbar.
The associated residual is rbar = bbar - Abar xbar
= P (b - (A - shift*I)x)
= P r.
In the discussion below, eps refers to the machine precision.
eps is computed by MINRES. A typical value is eps = 2.22d-16
for IEEE double-precision arithmetic.
Parameters
----------
n input The dimension of the matrix A.
b(n) input The rhs vector b.
r1(n) workspace
r2(n) workspace
v(n) workspace
w(n) workspace
w1(n) workspace
w2(n) workspace
x(n) output Returns the computed solution x.
y(n) workspace
Aprod external A subroutine defining the matrix A.
For a given vector x, the statement
call Aprod ( n, x, y )
must return the product y = Ax
without altering the vector x.
Msolve external An optional subroutine defining a
preconditioning matrix M, which should
approximate (A - shift*I) in some sense.
M must be positive definite.
For a given vector x, the statement
call Msolve( n, x, y )
must solve the linear system My = x
without altering the vector x.
In general, M should be chosen so that Abar has
clustered eigenvalues. For example,
if A is positive definite, Abar would ideally
be close to a multiple of I.
If A or A - shift*I is indefinite, Abar might
be close to a multiple of diag( I -I ).
NOTE. The program calling MINRES must declare
Aprod and Msolve to be external.
checkA input If checkA = .true., an extra call of Aprod will
be used to check if A is symmetric. Also,
if precon = .true., an extra call of Msolve
will be used to check if M is symmetric.
precon input If precon = .true., preconditioning will
be invoked. Otherwise, subroutine Msolve
will not be referenced; in this case the
actual parameter corresponding to Msolve may
be the same as that corresponding to Aprod.
shift input Should be zero if the system Ax = b is to be
solved. Otherwise, it could be an
approximation to an eigenvalue of A, such as
the Rayleigh quotient b'Ab / (b'b)
corresponding to the vector b.
If b is sufficiently like an eigenvector
corresponding to an eigenvalue near shift,
then the computed x may have very large
components. When normalized, x may be
closer to an eigenvector than b.
nout input A file number.
If nout .gt. 0, a summary of the iterations
will be printed on unit nout.
itnlim input An upper limit on the number of iterations.
rtol input A user-specified tolerance. MINRES terminates
if it appears that norm(rbar) is smaller than
rtol * norm(Abar) * norm(xbar),
where rbar is the transformed residual vector,
rbar = bbar - Abar xbar.
If shift = 0 and precon = .false., MINRES
terminates if norm(b - A*x) is smaller than
rtol * norm(A) * norm(x).
istop output An integer giving the reason for termination...
-1 beta2 = 0 in the Lanczos iteration; i.e. the
second Lanczos vector is zero. This means the
rhs is very special.
If there is no preconditioner, b is an
eigenvector of A.
Otherwise (if precon is true), let My = b.
If shift is zero, y is a solution of the
generalized eigenvalue problem Ay = lambda My,
with lambda = alpha1 from the Lanczos vectors.
In general, (A - shift*I)x = b
has the solution x = (1/alpha1) y
where My = b.
0 b = 0, so the exact solution is x = 0.
No iterations were performed.
1 Norm(rbar) appears to be less than
the value rtol * norm(Abar) * norm(xbar).
The solution in x should be acceptable.
2 Norm(rbar) appears to be less than
the value eps * norm(Abar) * norm(xbar).
This means that the residual is as small as
seems reasonable on this machine.
3 Norm(Abar) * norm(xbar) exceeds norm(b)/eps,
which should indicate that x has essentially
converged to an eigenvector of A
corresponding to the eigenvalue shift.
4 Acond (see below) has exceeded 0.1/eps, so
the matrix Abar must be very ill-conditioned.
x may not contain an acceptable solution.
5 The iteration limit was reached before any of
the previous criteria were satisfied.
6 The matrix defined by Aprod does not appear
to be symmetric.
For certain vectors y = Av and r = Ay, the
products y'y and r'v differ significantly.
7 The matrix defined by Msolve does not appear
to be symmetric.
For vectors satisfying My = v and Mr = y, the
products y'y and r'v differ significantly.
8 An inner product of the form x' M**(-1) x
was not positive, so the preconditioning matrix
M does not appear to be positive definite.
If istop .ge. 5, the final x may not be an
acceptable solution.
itn output The number of iterations performed.
Anorm output An estimate of the norm of the matrix operator
Abar = P (A - shift*I) P', where P = C**(-1).
Acond output An estimate of the condition of Abar above.
This will usually be a substantial
under-estimate of the true condition.
rnorm output An estimate of the norm of the final
transformed residual vector,
P (b - (A - shift*I) x).
ynorm output An estimate of the norm of xbar.
This is sqrt( x'Mx ). If precon is false,
ynorm is an estimate of norm(x).
------------------------------------------------------------------
MINRES is an implementation of the algorithm described in
the following reference:
C. C. Paige and M. A. Saunders (1975),
Solution of sparse indefinite systems of linear equations,
SIAM J. Numer. Anal. 12(4), pp. 617-629.
------------------------------------------------------------------
MINRES development:
1972: First version, similar to original SYMMLQ.
Later lost @#%*!
Oct 1995: Tried to reconstruct MINRES from
1995 version of SYMMLQ.
30 May 1999: Need to make it more like LSQR.
In middle of major overhaul.
19 Jul 2003: Next attempt to reconstruct MINRES.
Seems to need two vectors more than SYMMLQ. (w1, w2)
Lanczos is now at the top of the loop,
so the operator Aprod is called in just one place
(not counting the initial check for symmetry).
22 Jul 2003: Success at last. Preconditioning also works.
minres.f added to http://www.stanford.edu/group/SOL/.
FUTURE WORK: A stopping rule is needed for singular systems.
We need to estimate ||Ar|| as in LSQR. This will be
joint work with Sou Cheng Choi, SCCM, Stanford.
Note that ||Ar|| small => r is a null vector for A.
Michael A. Saunders na.msaunders@na-net.ornl.gov
Department of MS&E saunders@stanford.edu
Stanford University
Stanford, CA 94305-4026 (650) 723-1875
------------------------------------------------------------------
Subroutines and functions
USER Aprod , Msolve
BLAS1 daxpy , dcopy , ddot , dnrm2 } These are all in
Utilities daxpy2, dload2, dscal2 } the file minresblas.f
Functions
Definition at line 289 of file minres.f90.
References daxpy(), dcopy(), ddot(), dload2(), dnrm2(), dscal2(), and precon().
1.8.1