---
name: "pminres"
section: "4rheolef"
platform: "linux"
package: "librheolef-dev"
version: "6.7-1+b4"
release: "Debian 9 (Stretch)"
platforms: ["linux"]
url: "https://new.carta.tech/pminres"
---

# pminres(4rheolef)

## NAME

**pminres** -- conjugate gradient algorithm.

## SYNOPSIS

```
    template <class Matrix, class Vector, class Preconditioner, class Real>
    int pminres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
      int &max_iter, Real &tol, odiststream *p_derr=0);

```

## EXAMPLE

The simplest call to 'pminres' has the folling form:

```
    size_t max_iter = 100;
    double tol = 1e-7;
    int status = pminres(a, x, b, EYE, max_iter, tol, &derr);

```

## DESCRIPTION

**pminres** 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

- **max\_iter**

  the number of iterations performed before the tolerance was reached

- **tol**

  the residual after the final iteration

## NOTE

**pminres** 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, class Real, class Size>
int pminres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
  Size &max_iter, Real &tol, odiststream *p_derr = 0, std::string label = "minres")
{
  Vector b = M.solve(Mb);
  Real norm_b = sqrt(fabs(dot(Mb,b)));
  if (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));
  if (p_derr) (*p_derr) << "[" << label << "] # norm_b=" <<norm_b<< std::endl;
  if (p_derr) (*p_derr) << "[" << label << "] #iteration residue/norm_b" << std::endl;
  if (p_derr) (*p_derr) << "[" << label << "] 0 " << norm_r/norm_b << std::endl;
  if (beta2 < 0 || norm_r <= tol*norm_b) {
    tol = norm_r/norm_b;
    max_iter = 0;
    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 (Size n = 1; n <= max_iter; n++) {
    // 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 ("pminres: machine precision problem");
        tol = norm_r/norm_b;
        max_iter = n;
        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;
    if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << norm_r/norm_b << std::endl;
    if (norm_r <= tol*max(1.0,norm_b)) {
      tol = norm_r/norm_b;
      max_iter = n;
      return 0;
    }
  }
  tol = norm_r/norm_b;
  return 1;
}

```