.TH "minres" 5rheolef "Version 7.2" "rheolef" \" -*- nroff -*-
.ad l
.nh
.SH NAME
minres \- minimum residual algorithm (rheolef-7\&.2)
.PP
 
.SH "SYNOPSIS"
.PP
.PP
.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())
.fi
.PP
.SH "DESCRIPTION"
.PP
This function solves the symmetric positive but possibly \fIsingular\fP linear system \fCA*x=b\fP with the minimal residual method\&. The \fCminres\fP function follows the algorithm described in 
.PP
.nf
    C\&. C\&. Paige and M\&. A\&. Saunders,
    Solution of sparse indefinite systems of linear equations',
    SIAM J\&. Numer\&. Anal\&., 12(4), 1975\&.

.fi
.PP
 For more, see http://www.stanford.edu/group/SOL/software.html and also at page 60 of the PhD report: 
.PP
.nf
   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

.fi
.PP
 
.SH "EXAMPLE"
.PP
.PP
.nf
    solver_option sopt;
    sopt\&.max_iter = 100;
    sopt\&.tol = 1e-7;
    int status = minres (A, x, b, eye(), sopt);
.fi
.PP
 The fourth argument of \fCminres\fP is a preconditionner: here, the \fBeye(5)\fP one is a do-nothing preconditionner, for simplicity\&. Finally, the \fBsolver_option(4)\fP variable \fCsopt\fP transmits the stopping criterion with \fCsopt\&.tol\fP and \fCsopt\&.max_iter\fP\&.
.PP
On return, the \fCsopt\&.residue\fP and \fCsopt\&.n_iter\fP indicate the reached residue and the number of iterations effectively performed\&. The return \fCstatus\fP is zero when the prescribed tolerance \fCtol\fP has been obtained, and non-zero otherwise\&. Also, the \fCx\fP variable contains the approximate solution\&. See also the \fBsolver_option(4)\fP for more controls upon the stopping criterion\&.
.SH "IMPLEMENTATION"
.PP
This documentation has been generated from file linalg/lib/minres\&.h
.PP
The present template implementation is inspired from the \fCIML++ 1\&.2\fP iterative method library, http://math.nist.gov/iml++
.PP
.PP
.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())
.fi
.PP
 
.PP
.nf
{
  // 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) {
    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;
}

.fi
.PP
 
.SH AUTHOR
Pierre  Saramito  <Pierre.Saramito@imag.fr>
.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.