cg(5rheolef)

✨ cg — bestie edition ✨

okay so this is the cg page from rheolef-7.2, section 5 of the manual (yes, cg(5rheolef), that's the whole vibe). sit down, grab a boba, we're solving linear systems 💅

---

🎀 NAME 🎀

cg — the conjugate gradient algorithm (from rheolef-7.2). that's it, that's the name, she's simple but she's powerful.

---

✨ SYNOPSIS ✨

here's the signature, no cap, copy it exactly bestie:

template <class Matrix, class Vector, class Vector2, class Preconditioner>
int cg (const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M,
        const solver_option& sopt = solver_option())

---

🦋 EXAMPLE 🦋

okay here's a lil taste, a lil sample, a lil amuse-bouche:

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

see how cute that is?? 100 iterations, tolerance 1e-7, and eye() doing absolutely nothing as a preconditioner (we love a low-effort king 👑).

---

💖 DESCRIPTION 💖

so. what does she DO. this function solves the symmetric positive definite linear system

A*x=b

with the conjugate gradient method. not me being obsessed with an iterative solver 😭. the cg function follows the algorithm described on p. 15 in:

R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst, Templates for the solution of linear systems: building blocks for iterative methods, SIAM, 1994.

(the citation is literally right there, i'm not making it up, i would never 💀)

the fourth argument of cg is a preconditionner — and here, the eye(5) one is a do-nothing preconditionner, for simplicity. it's giving "i'm here but i'm not doing anything" energy and honestly? valid.

finally, the solver_option(4) variable sopt transmits the stopping criterion with sopt.tol and sopt.max_iter. those two are the main characters here, fr fr.

on return:

• sopt.residue and sopt.n_iter indicate the reached residue and the number of iterations effectively performed. receipts!! 🧾
• 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. that's the whole point bestie, x is the glow-up.

see also the solver_option(4) for more controls upon the stopping criterion. iykyk.

---

🔥 IMPLEMENTATION 🔥

this documentation has been generated from file linalg/lib/cg.h

the present template implementation is inspired from the IML++ 1.2 iterative method library, http://math.nist.gov/iml++

and now, the sacred text, verbatim, do not touch, do not edit, do not even breathe on it:

template <class Matrix, class Vector, class Vector2, class Preconditioner>
int cg (const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M,
        const solver_option& sopt = solver_option())

{
  typedef typename Vector::size_type  Size;
  typedef typename Vector::float_type Real;
  std::string label = (sopt.label != "" ? sopt.label : "cg");
  Vector b = M.solve(Mb);
  Real norm2_b = dot(Mb,b);
  if (sopt.absolute_stopping || norm2_b == Real(0)) norm2_b = 1;
  Vector Mr = Mb - A*x;
  Real  norm2_r = 0;
  if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl;
  Vector p;
  for (sopt.n_iter = 0; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
    Vector r = M.solve(Mr);
    Real prev_norm2_r = norm2_r;
    norm2_r = dot(Mr, r);
    sopt.residue = sqrt(norm2_r/norm2_b);
    if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
    if (sopt.residue <= sopt.tol) return 0;
    if (sopt.n_iter == 0) {
      p = r;
    } else {
      Real beta = norm2_r/prev_norm2_r;
      p = r + beta*p;
    }
    Vector Mq = A*p;
    Real alpha = norm2_r/dot(Mq, p);
    x  += alpha*p;
    Mr -= alpha*Mq;
  }
  return 1;
}

look at that loop, look at that beta, look at that alpha — she's doing math and she's doing it beautifully 💜

---

👑 AUTHOR 👑

Pierre Saramito <Pierre.Saramito@imag.fr>

(we stan, thank you for your service 🫡)

---

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

ok so real talk for a sec: that last bit means no warranty, like at all, to the extent permitted by law. she's free software bestie, you can change her, you can share her, but if she misbehaves that's on you 🙈

---

and that's the page!! this came from rheolef, Version 7.2, and the original footer said cg(5rheolef) — so yeah, we're the section 5 girly, version 7.2, and we ate. periodt 💅✨🎀🦋