---
name: "minres"
section: "4rheolef"
platform: "linux"
package: "librheolef-dev"
version: "7.0-2+b1"
release: "Debian 10 (Buster)"
platforms: ["linux"]
url: "https://new.carta.tech/minres/4rheolef"
---

# minres(4rheolef)

## NAME

**minres** -- conjugate gradient algorithm.

## SYNOPSIS

```
  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);

```

## EXAMPLE

The simplest call to **minres** has the folling form:

```
  solver_option sopt;
  sopt.max_iter = 100;
  sopt.tol = 1e-7;
  int status = minres (a, x, b, eye(), sopt);

```

## DESCRIPTION

**minres** solves the symmetric positive definite linear system Ax=b using the Conjugate Gradient method.

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:

- **x**

  approximate solution to Ax = b

- **sopt.n\_iter**

  the number of iterations performed before the tolerance was reached

- **sopt.residue**

  the residual after the final iteration

## NOTE

**minres** 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 **IML++ 1.2** iterative method library, **http://math.nist.gov/iml++**.

## IMPLEMENTATION

```
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;
}

```

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