minimum residual algorithm
Synopsis
template <class Matrix, class Vector, class Preconditioner>
int minres (
const Matrix &
A, Vector &x,
const Vector &Mb,
const Preconditioner &
M,
const solver_option& sopt = solver_option())
int minres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Description
This function solves the symmetric positive but possibly singular linear system A*x=b
with the minimal residual method. The minres
function follows the algorithm described in
C. C. Paige and M. A. Saunders,
Solution of sparse indefinite systems of linear equations",
SIAM J. Numer. Anal., 12(4), 1975.
For more, see http://www.stanford.edu/group/SOL/software.html and also at page 60 of the PhD report:
S.-C. T. Choi,
Iterative methods for singular linear equations and least-squares problems,
Stanford University, 2006,
http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf
Example
solver_option sopt;
sopt.max_iter = 100;
sopt.tol = 1e-7;
int status = minres (A, x, b, eye(), sopt);
The fourth argument of minres
is a preconditionner: here, the eye
one is a do-nothing preconditionner, for simplicity. Finally, the solver_option
variable sopt
transmits the stopping criterion with sopt.tol
and sopt.max_iter
.
On return, the sopt.residue
and sopt.n_iter
indicate the reached residue and the number of iterations effectively performed. The return status
is zero when the prescribed tolerance tol
has been obtained, and non-zero otherwise. Also, the x
variable contains the approximate solution. See also the solver_option
for more controls upon the stopping criterion.
Implementation
This documentation has been generated from file linalg/lib/minres.h
The present template implementation is inspired from the IML++ 1.2
iterative method library, http://math.nist.gov/iml++
template <class Matrix, class Vector, class Preconditioner>
int minres (
const Matrix &
A, Vector &x,
const Vector &Mb,
const Preconditioner &
M,
const solver_option& sopt = solver_option())
{
typedef typename Vector::float_type Real;
std::string label = (sopt.label != "" ? sopt.label : "minres");
Real norm_b = sqrt(fabs(
dot(Mb,b)));
if (sopt.absolute_stopping || norm_b == Real(0.)) norm_b = 1;
Real norm_r = sqrt(fabs(beta2));
sopt.residue = norm_r/norm_b;
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl
<< "[" << label << "] 0 " << sopt.residue << std::endl;
if (beta2 < 0 || sopt.residue <= sopt.tol) {
return 0;
}
Real beta = sqrt(beta2);
Vector Mv = Mr/beta;
Real c_old = 1.;
Real s_old = 0.;
Real c = 1.;
Real s = 0.;
Vector u_old (x.ownership(), 0.);
Vector Mv_old (x.ownership(), 0.);
Vector w (x.ownership(), 0.);
Vector w_old (x.ownership(), 0.);
Vector w_old2 (x.ownership(), 0.);
for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
Mr = Mr - alpha*Mv - beta*Mv_old;
z = z - alpha*
u - beta*u_old;
if (beta2 < 0) {
sopt.residue = norm_r/norm_b;
return 2;
}
Real c_old2 = c_old;
Real s_old2 = s_old;
s_old = s;
Real rho0 = c_old*
alpha - c_old2*s_old*beta_old;
Real rho2 = s_old*
alpha + c_old2*c_old*beta_old;
Real rho1 = sqrt(sqr(rho0) + sqr(beta));
Real rho3 = s_old2 * beta_old;
w_old2 = w_old;
w_old = w;
w = (
u - rho2*w_old - rho3*w_old2)/rho1;
Mv_old = Mv;
norm_r *= s;
sopt.residue = norm_r/norm_b;
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
if (sopt.residue <= sopt.tol) return 0;
}
return 1;
}
#define dis_warning_macro(message)
Float alpha[pmax+1][pmax+1]
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
DATA::size_type size_type