Millepede-II  V04-08-02
Functions/Subroutines
minresmodule Module Reference

MINRES solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite and/or singular. More...

Functions/Subroutines

subroutine, public minres (n, Aprod, Msolve, b, shift, checkA, precon, x, itnlim, nout, rtol, istop, itn, Anorm, Acond, rnorm, Arnorm, ynorm)
 Solution of linear equation system. More...
 

Detailed Description

MINRES solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite and/or singular.

 The software for MINRES (f90 version) is provided by SOL, Stanford University
 under the terms of the OSI Common Public License (CPL):
 http://www.opensource.org/licenses/cpl1.0.php

 Contributors:
     Chris Paige <chris@cs.mcgill.ca>
     Sou-Cheng Choi <scchoi@stanford.edu>

     Michael Saunders <saunders@stanford.edu>
     Systems Optimization Laboratory (SOL)
     Stanford University
     Stanford, CA 94305-4026, USA
     (650)723-1875

 09 Oct 2007: F90 version constructed from the F77 version.
              Initially used compiler option -r8, but this is nonstandard.
 15 Oct 2007: Test on Arnorm = ||Ar|| added to recognize singular systems.
 15 Oct 2007: Temporarily used real(8) everywhere.
 16 Oct 2007: Use minresDataModule to define dp = selected_real_kind(15).
              We need "use minresDataModule"
              at the beginning of modules AND inside interfaces.

              g95 compiles successfully with the following options:
   g95 -c -g -O0 -pedantic -Wall -Wextra -fbounds-check -ftrace=full minresModule.f90

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Function/Subroutine Documentation

◆ minres()

subroutine, public minresmodule::minres ( integer, intent(in)  n,
  Aprod,
  Msolve,
real(dp), dimension(n), intent(in)  b,
real(dp), intent(in)  shift,
logical, intent(in)  checkA,
logical, intent(in)  precon,
real(dp), dimension(n), intent(out)  x,
integer, intent(in)  itnlim,
integer, intent(in)  nout,
real(dp), intent(in)  rtol,
integer, intent(out)  istop,
integer, intent(out)  itn,
real(dp), intent(out)  Anorm,
real(dp), intent(out)  Acond,
real(dp), intent(out)  rnorm,
real(dp), intent(out)  Arnorm,
real(dp), intent(out)  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.

 Parameters
 ----------

 n       input      The dimension of the matrix A.
 b(n)    input      The rhs vector b.
 x(n)    output     Returns the computed solution x.

 Aprod   external   A subroutine defining the matrix A.
                       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.

                       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 ).

 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 > 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 >= 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/.

 16 Oct 2007: Added a stopping rule for singular systems,
              as derived in Sou-Cheng Choi's PhD thesis.
              Note that ||Ar|| small => r is a null vector for A.
              Subroutine minrestest2 in minresTestModule.f90
              tests this option.  (NB: Not yet working.)
-------------------------------------------------------------------

Definition at line 294 of file minresModule.f90.

References precon().

Referenced by mminrs(), and solglo().