minres(5rheolef)

📦
version 7.2 · 4 versions
, in debian trixie, ubuntu noble, ubuntu questing

💖 minres — the minimum residual algorithm 💖

okay so bestie, welcome to minres(5rheolef), your new fave page in the rheolef manual. this one's all about a solver that just keeps going until the residual is basically gone. main character energy fr. ✨

🌸 SYNOPSIS 🌸

like, here's the whole vibe in one glance:

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())

🦋 DESCRIPTION 🦋

so this function solves the symmetric positive but possibly singular linear system A*x=b using the minimal residual method. yes, "possibly singular" — she can handle it, no cap. 💅

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.

and for more, see http://www.stanford.edu/group/SOL/software.html — plus, if you want the deep lore, page 60 of this 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 🎀

here's her in action, living her best life:

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 — and here the eye(5) one is a do-nothing preconditionner, for simplicity. honestly? iconic. 🫶

finally, the solver_option(4) variable sopt transmits the stopping criterion with sopt.tol and sopt.max_iter.

on return, 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(4) 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++

okay here's the whole thing, verbatim, because we do not gatekeep the source 🧚:

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) {
    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;
}
lowkey important: if beta2 < 0 mid-loop she throws a dis_warning_macro saying "minres: machine precision problem", sets sopt.residue = norm_r/norm_b and returns 2. that's the "bestie the machine said no" moment. 😭

👑 AUTHOR 👑

Pierre Saramito <Pierre.Saramito@imag.fr>

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.

and that's the whole page, bestie!! rheolef, Version 7.2, minres(5rheolef) — slay, now go solve something symmetric and possibly singular 💅✨🌸

📦 where this page comes from