.\"
.de Id
..
.de Sp
.if n .sp
.if t .sp 0.4
..
.TH minres 4rheolef "rheolef-7.0" "rheolef-7.0" "rheolef-7.0"
.\" label: /*D:minres
.SH NAME
\fBminres\fP -- conjugate gradient algorithm.
.\" skip: @findex minres
.\" skip: @cindex conjugate gradient algorithm
.\" skip: @cindex iterative solver
.\" skip: @cindex preconditioner
.SH SYNOPSIS

.\" begin_example
.Sp
.nf
  template <class Matrix, class Vector, class Preconditioner, class Real>
  int minres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
    const solver_option& sopt = solver_option);
.Sp
.fi
.\" end_example
.SH EXAMPLE

The simplest call to \fBminres\fP has the folling form:
.\" begin_example
.Sp
.nf
  solver_option sopt;
  sopt.max_iter = 100;
  sopt.tol = 1e-7;
  int status = minres (a, x, b, eye(), sopt);
.Sp
.fi
.\" end_example
.SH DESCRIPTION

\fBminres\fP solves the symmetric positive definite linear
system Ax=b using the Conjugate Gradient method.
.PP
The return value indicates convergence within max_iter (input)
iterations (0), or no convergence within max_iter iterations (1).
Upon successful return, output arguments have the following values:
.\" begin table
.\" start item
.TP
.B x
approximate solution to Ax = b

.\" start item
.TP
.B sopt.n_iter
the number of iterations performed before the tolerance was reached

.\" start item
.TP
.B sopt.residue
the residual after the final iteration
.\" end table
.SH NOTE

\fBminres\fP follows the algorithm described in
"Solution of sparse indefinite systems of linear equations", C. C. Paige
and M. A. Saunders, SIAM J. Numer. Anal., 12(4), 1975.
For more, see http://www.stanford.edu/group/SOL/software.html and also the
PhD "Iterative methods for singular linear equations and least-squares problems",
S.-C. T. Choi, Stanford University, 2006,
http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf at page 60.
The present implementation style is inspired from \fBIML++ 1.2\fP iterative method library,
\fBhttp://math.nist.gov/iml++\fP.
.\" skip start:AUTHOR:

.\" skip start:DATE:

.\" skip start:METHODS:
.\" END
.SH IMPLEMENTATION
.\" begin_example
.Sp
.nf
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())
{
  // Size &max_iter, Real &tol, odiststream *p_derr = 0
  typedef typename Vector::size_type  Size;
  typedef typename Vector::float_type Real;
  std::string label = (sopt.label != "" ? sopt.label : "minres");
  Vector b = M.solve(Mb);
  Real norm_b = sqrt(fabs(dot(Mb,b)));
  if (sopt.absolute_stopping || norm_b == Real(0.)) norm_b = 1;
  Vector Mr = Mb - A*x;
  Vector z = M.solve(Mr);
  Real beta2 = dot(Mr, z);
  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) {
    dis_warning_macro ("beta2 = " << beta2 << " < 0: stop");
    return 0;
  }
  Real beta = sqrt(beta2);
  Real eta = beta;
  Vector Mv = Mr/beta;
  Vector  u =  z/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++) {
    // Lanczos
    Mr = A*u;
    z = M.solve(Mr);
    Real alpha = dot(Mr, u);
    Mr = Mr - alpha*Mv - beta*Mv_old;
    z  =  z - alpha*u  - beta*u_old;
    beta2 = dot(Mr, z);
    if (beta2 < 0) {
        dis_warning_macro ("minres: machine precision problem");
        sopt.residue = norm_r/norm_b;
        return 2;
    }
    Real beta_old = beta;
    beta = sqrt(beta2);
    // QR factorisation
    Real c_old2 = c_old;
    Real s_old2 = s_old;
    c_old = c;
    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;
    // Givens rotation
    c = rho0 / rho1;
    s = beta / rho1;
    // update
    w_old2 = w_old;
    w_old  = w;
    w = (u - rho2*w_old - rho3*w_old2)/rho1;
    x += c*eta*w;
    eta = -s*eta;
    Mv_old = Mv;
     u_old = u;
    Mv = Mr/beta;
     u =  z/beta;
    // check residue
    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;
}
.Sp
.fi
.\" end_example
.SH COPYRIGHT
Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito@imag.fr> GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>.
This is free software: you are free to change and redistribute it.  There is NO WARRANTY, to the extent permitted by law.