Millepede-II  V04-07-02
Functions/Subroutines
minresqlpmodule Module Reference

MINRESQLP 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 minresqlp (n, Aprod, b, shift, Msolve, disable, nout, itnlim, rtol, maxxnorm, trancond, Acondlim, x, istop, itn, rnorm, Arnorm, xnorm, Anorm, Acond)
 Solution of linear equation system or least squares problem. More...
 
subroutine, public symortho (a, b, c, s, r)
 SymOrtho: Stable Householder reflection. More...
 

Detailed Description

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

"A" is really (A - shift*I), where shift is an input real scalar.

 09 Sep 2013: Version 27
-------------------------------------------------------------------

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

-------------------------------------------------------------------

 Authors:
     Sou-Cheng Choi <sctchoi@uchicago.edu>
     Computation Institute (CI)
     University of Chicago
     Chicago, IL 60637, USA

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

 Contributor:
     Christopher Paige <paige@cs.mcgill.ca>

   See also: Makefile

Function/Subroutine Documentation

◆ minresqlp()

subroutine, public minresqlpmodule::minresqlp ( integer(ip), intent(in)  n,
  Aprod,
real(dp), dimension(n), intent(in)  b,
real(dp), intent(in), optional  shift,
optional  Msolve,
logical, intent(in), optional  disable,
integer(ip), intent(in), optional  nout,
integer(ip), intent(in), optional  itnlim,
real(dp), intent(in), optional  rtol,
real(dp), intent(in), optional  maxxnorm,
real(dp), intent(in), optional  trancond,
real(dp), intent(in), optional  Acondlim,
real(dp), dimension(n), intent(out)  x,
integer(ip), intent(out), optional  istop,
integer(ip), intent(out), optional  itn,
real(dp), intent(out), optional  rnorm,
real(dp), intent(out), optional  Arnorm,
real(dp), intent(out), optional  xnorm,
real(dp), intent(out), optional  Anorm,
real(dp), intent(out), optional  Acond 
)

Solution of linear equation system or least squares problem.

------------------------------------------------------------------

 MINRESQLP  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 roughly the same
        number of iterations as MINRESQLP but less work per iteration.
        But if a low-accuracy solution is adequate, MINRESQLP will
        terminate sooner.

     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 roughly the same number of iterations as
        MINRESQLP but slightly less work per iteration.

     3. If A is indefinite and well-conditioned, and Ax = b has a
        solution, i.e., it is not a least-squares problem, MINRES might
        be preferred as it requires the same number of iterations as
        MINRESQLP 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

                call Aprod ( n, x, y )

     which must return the product y = Ax for any given vector x.


     More generally, MINRESQLP is designed to solve the system

                (A - shift*I) x = b
     or
         min || (A - shift*I) x - b ||_2,

     where shift is a specified real scalar.  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 are not yet sure on that -- it may be better
     to use SYMMLQ.

     In this documentation, ' denotes the transpose of
     a vector or a matrix.

     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 parameter Msolve may be used (see below).
     When an external procedure Msolve is supplied, MINRESQLP 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 MINRESQLP.  A typical value is eps = 2.22d-16
     for IEEE double-precision arithmetic.

     Parameters
     ----------
     Some inputs are optional, with default values described below.
     Mandatory inputs are n, Aprod, and b.
     All outputs other than x are optional.

     n       input      The dimension of the matrix or operator A.

     b(n)    input      The rhs vector b.

     x(n)    output     Returns the computed solution  x.

     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.
                        An extra call of Aprod is
                        used to check if A is symmetric.

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

     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. Default to 0.

     nout    input      A file number. The calling program must open a file
                        for output using for example:
                          open(nout, file='MINRESQLP.txt', status='unknown')
                        If nout > 0, a summary of the iterations
                        will be printed on unit nout. If nout is absent or
                        the file associated with nout is not opened properly,
                        results will be written to 'MINRESQLP_tmp.txt'.
                        (Avoid 0, 5, 6 because by convention stderr=0,
                        stdin=5, stdout=6.)

     itnlim  input      An upper limit on the number of iterations. Default to 4n.

     rtol    input      A user-specified tolerance.  MINRESQLP terminates
                        if it appears that norm(rbar) is smaller than
                              rtol*[norm(Abar)*norm(xbar) + norm(b)],
                        where rbar = bbar - Abar xbar,
                        or that norm(Abar*rbar) is smaller than
                              rtol*norm(Abar)*norm(rbar).

                        If shift = 0 and Msolve is absent, MINRESQLP
                        terminates if norm(r) is smaller than
                              rtol*[norm(A)*norm(x) + norm(b)],
                        where r = b - Ax,
                        or if norm(A*r) is smaller than
                              rtol*norm(A)*norm(r).

                        Default to machine precision.

     istop   output     An integer giving the reason for termination...
               0        Initial value of istop.

               1        beta_{k+1} < eps.
                        Iteration k is the final Lanczos step.

               2        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 Abar. Also, x = (1/alpha1) b
                        is a solution of Abar x = b.

                        Otherwise (if Msolve is present), 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.

               3        b = 0, so the exact solution is x = 0.
                        No iterations were performed.

               4        Norm(rbar) appears to be less than
                        the value  rtol * [ norm(Abar) * norm(xbar) + norm(b) ].
                        The solution in  x  should be an acceptable
                        solution of Abar x = b.

               5        Norm(rbar) appears to be less than
                        the value  eps * norm(Abar) * norm(xbar).
                        This means that the solution is as accurate as
                        seems reasonable on this machine.

               6        Norm(Abar rbar) appears to be less than
                        the value  rtol * norm(Abar) * norm(rbar).
                        The solution in x should be an acceptable
                        least-squares solution.

               7        Norm(Abar rbar) appears to be less than
                        the value  eps * norm(Abar) * norm(rbar).
                        This means that the least-squares solution is as
                        accurate as seems reasonable on this machine.

               8        The iteration limit was reached before convergence.

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

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

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

               12       xnorm has exceeded maxxnorm or will exceed it
                        next iteration.

               13       Acond (see below) has exceeded Acondlim or 0.1/eps,
                        so the matrix Abar must be very ill-conditioned.

               14       | gamma_k | < eps.
                        This is very likely a least-squares problem but
                        x may not contain an acceptable solution yet.

               15       norm(Abar x) < rtol * norm(Abar) * norm(x).
                        If disable = .true., then a null vector will be
                        obtained, given rtol.

                        If istop >= 7, 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).

     xnorm   output     An estimate of the norm of xbar.
                        This is sqrt( x'Mx ).  If Msolve is absent,
                        xnorm is an estimate of norm(x).

   maxxnorm  input      An upper bound on norm(x). Default value is 1e7.

   trancond  input      If trancond > 1, a switch is made from MINRES
                        iterations to MINRES-QLP iterations when
                        Acond > trancond.
                        If trancond = 1, all iterations are MINRES-QLP
                        iterations.
                        If trancond = Acondlim, all iterations are
                        conventional MINRES iterations (which are
                        slightly cheaper).
                        Default to 1e7.

   Acondlim  input      An upper bound on Acond. Default value is 1e15.

   disable   input      All stopping conditions are disabled except
                        norm(Ax) / norm(x) < tol. Default to .false..

------------------------------------------------------------------

     MINRESQLP is an implementation of the algorithm described in
     the following references:

     Sou-Cheng Choi,
     Iterative Methods for Singular Linear Equations and Least-
     Squares Problems, PhD dissertation, ICME, Stanford University,
     2006.

     Sou-Cheng Choi, Christopher Paige, and Michael Saunders,
     MINRES-QLP: A Krylov subspace method for indefinite or
     singular symmetric systems, SIAM Journal of Scientific
     Computing 33:4 (2011) 1810-1836.

     Sou-Cheng Choi and Michael Saunders,
     ALGORITHM & DOCUMENTATION: MINRES-QLP for singular Symmetric and Hermitian
     linear equations and least-squares problems, Technical Report,
     ANL/MCS-P3027-0812, Computation Institute,
     University of Chicago/Argonne National Laboratory, 2012.

     Sou-Cheng Choi and Michael Saunders,
     ALGORITHM xxx: MINRES-QLP for singular Symmetric and Hermitian
     linear equations and least-squares problems,
     ACM Transactions on Mathematical Software, to appear, 2013.

     FORTRAN 90 and MATLAB implementations are
     downloadable from
            http://www.stanford.edu/group/SOL/software.html
            http://home.uchicago.edu/sctchoi/

------------------------------------------------------------------

     MINRESQLP development:
     14 Dec 2006: Sou-Cheng's thesis completed.
                  MINRESQLP includes a stopping rule for singular
                  systems (using an estimate of ||Ar||) and very many
                  other things(!).
                  Note that ||Ar|| small => r is a null vector for A.
     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 minresqlpDataModule to define
                  dp = selected_real_kind(15).
                  We need "use minresqlpDataModule" at the
                  beginning of modules AND inside interfaces.
     06 Jun 2010: Added comments.
     12 Jul 2011: Created complex version zminresqlpModule.f90
                  from real version minresqlpModule.f90.
     23 Aug 2011: (1) Tim Hopkins ran version 17 on the NAG Fortran compiler
                  We removed half a dozen unused variables in MINRESQLP
                  and also local var sgn_a and sgn_b in SMMORTHO,
                  as they result in division by zero for inputs a=b=0.
                  (2) Version 18 was submitted to ACM TOMS for review.
     20 Aug 2012: Version 19:
                  (1) Added optional inputs and outputs, and
                  default values for optional inputs.
                  (2) Removed inputs 'checkA' and 'precon'.
                  (3) Changed slightly the order of parameters in the
                  MINRESQLP API.
                  (4) Updated documentation.
                  (5) Fixed a minor bug in printing x(1) in iteration
                      log during MINRES mode.
                  (6) Made sure MINRESQLP is portable in both single
                      and double precison.
                  (7) Fixed a bug to ensure the 2x2 Hermitian reflectors
                      are orthonormal. Make output c real.
     24 Apr 2013: istop = 12 now means xnorm just exceeded maxxnorm.
     28 Jun 2013: likeLS introduced to terminate with big xnorm
                  only if the problem seems to be singular and inconsistent.
     08 Jul 2013: (1) dot_product replaces ddotc.
     04 Aug 2013: If present(maxxnorm), use maxxnorm_ = min(maxxnorm, one/eps).
     09 Sep 2013: Initialize relresl and relAresl to zero.
------------------------------------------------------------------

Definition at line 411 of file minresqlpModule.f90.

References minresqlpdatamodule::debug, minresqlpblasmodule::dnrm2(), minresqlpdatamodule::eps, minresqlpdatamodule::one, minresqlpdatamodule::prcsn, symortho(), and minresqlpdatamodule::zero.

Referenced by mminrsqlp(), and solgloqlp().

◆ symortho()

subroutine, public minresqlpmodule::symortho ( real(dp), intent(in)  a,
real(dp), intent(in)  b,
real(dp), intent(out)  c,
real(dp), intent(out)  s,
real(dp), intent(out)  r 
)

SymOrtho: Stable Householder reflection.

  USAGE:
     SymOrtho(a, b, c, s, r)

  INPUTS:
    a      first element of a two-vector  [a; b]
    b      second element of a two-vector [a; b]

  OUTPUTS:
    c      cosine(theta), where theta is the implicit angle of reflection
    s      sine(theta)
    r      two-norm of [a; b]

  DESCRIPTION:
     Stable Householder reflection that gives c and s such that
        [ c  s ][a] = [r],
        [ s -c ][b]   [0]
     where r = two norm of vector [a, b],
        c = a / sqrt(a**2 + b**2) = a / r,
        s = b / sqrt(a**2 + b**2) = b / r.
     The implementation guards against overlow in computing sqrt (a**2 + b**2).


  REFERENCES:
    Algorithm 4.9, stable unsymmetric Givens rotations in
     Golub and van Loan's book Matrix Computations, 3rd edition.

  MODIFICATION HISTORY:
    20/08/2012: Fixed a bug to ensure the 2x2 Hermitian reflectors
                are orthonormal.
    05/27/2011: Created this file from Matlab SymGivens2.m

  KNOWN BUGS:
     MM/DD/2004: description

  AUTHORS: Sou-Cheng Choi, CI, University of Chicago
           Michael Saunders, MS&E, Stanford University

  CREATION DATE: 05/27/2011

Definition at line 1334 of file minresqlpModule.f90.

References minresqlpdatamodule::debug, minresqlpdatamodule::one, minresqlpdatamodule::realmin, and minresqlpdatamodule::zero.

Referenced by minresqlp().