tminres
Functions
tminres.hpp File Reference
#include <limits>
#include <memory>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <algorithm>

Go to the source code of this file.

Functions

template<typename Operator , typename Vector , typename Preconditioner >
int MINRES (const Operator &A, Vector &x, const Vector &b, const Preconditioner *M, double shift, int &max_iter, double &tol, bool show)
 Preconditioner Minres algorithm.

Detailed Description

Author:
U. Villa - uvilla@emory.edu
Date:
04/2012

Definition in file tminres.hpp.


Function Documentation

template<typename Operator , typename Vector , typename Preconditioner >
int MINRES ( const Operator &  A,
Vector &  x,
const Vector &  b,
const Preconditioner *  M,
double  shift,
int &  max_iter,
double &  tol,
bool  show 
)

Preconditioner Minres algorithm.

minres solves the n x n system of linear equations Ax = b or the n x n least squares problem min ||Ax - b||_2^2, where A is a symmetric matrix (possibly indefinite or singular) and b is a given vector. The dimension n is defined by length(b).

Parameters:
Aconst OPERATOR & (INPUT) A is a linear operator with a method Apply(const Vector & X, Vector & Y) which gives the product y = A*x for any given n-vector x.
xVECTOR & (INPUT/OUTPUT) "x" is the initial guess (INPUT) and the solution of the linear system (OUTPUT).
bconst VECTOR & (INPUT) "b" is the right hand side of the linear system
Mconst PRECONDITIONER * (INPUT) "M" is a linear operator (preconditioner) with a method Apply(const Vector & X, Vector & Y) such that y = M(x) solves the system My = x given any n-vector x. If M is a null pointer, preconditioning is not used. Otherwise, "M" defines a (implicitely) a positive-definite preconditioner M = C*C'.
shiftdouble (INPUT) If shift != 0, minres really solves (A - shift*I)x = b (or the corresponding least-squares problem if shift is an eigenvalue of A).
max_iterint& (INPUT/OUTPUT) "max_iter" is the maximum number of minres iterations (INPUT) and the actual number of iteration at termination (OUTPUT)
toldouble& (INPUT) "tol" is tolerance for the stopping criterion
showbool (INPUT) If "show" is true print information on screen.
Returns:
istop int "istop" is a value from [-1:9] to indicate the reason for termination. The reason is summarized in msg[istop+2] below.

When M = C*C' exists, minres implicitly solves the system

            P(A - shift*I)P'xbar = Pb,
    i.e.               Abar xbar = bbar,
    where                      P = inv(C),
                            Abar = P(A - shift*I)P',
                            bbar = Pb,

 and returns the solution      x = P'xbar.
 The associated residual is rbar = bbar - Abar xbar
                                 = P(b - (A - shift*I)x)
                                 = Pr.

Known bugs:

  • As Jeff Kline pointed out, Arnorm = ||A r_{k-1}|| lags behind rnorm = ||r_k||. On singular systems, this means that a good least-squares solution exists before Arnorm is small enough to recognize it. The solution x_{k-1} gets updated to x_k (possibly a very large solution) before Arnorm shuts things down the next iteration. It would be better to keep x_{k-1}.

Other notes:

  • ynorm = norm(x) is computed directly at the cost of an additional inner product. This can represent a bottle-neck in the parallel case. In fact it is possible to avoid norm(x) by updating certain scalar quantities (see the QLP factors in MINRES-QLP). It would add only cheap (scalar) arithmetic, but lots of extra code.
Examples:
SerialExample/ex1.cpp, SerialExample/ex2.cpp, and TrilinosExample/ex1.cpp.

Definition at line 106 of file tminres.hpp.

 All Classes Files Functions Friends