Extracting algorithm of an Iterative-method module

[Adam2018]

Dear all,

I'm not familiar with vector and matrices as classes. I appreciate your help with extracting the *correct* underlying algorithm of the below
module. Thanks in Advance, Adam

//*****************************************************************
// Iterative template routine -- Preconditioned Richardson
//
// IR solves the unsymmetric linear system Ax = b using 
// Iterative Refinement (preconditioned Richardson iteration).
//
// 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
//  
//*****************************************************************

template < class Matrix, class Vector, class Preconditioner, class Real >
int 
IR(const Matrix &A, Vector &x, const Vector &b,
   const Preconditioner &M, int &max_iter, Real &tol)
{
  Real resid;
  Vector z;

  Real normb = norm(b);
  Vector r = b - A*x;

  if (normb == 0.0) 
    normb = 1;
  
  if ((resid = norm(r) / normb) <= tol) {
    tol = resid;
    max_iter = 0;
    return 0;
  }
  
  for (int i = 1; i <= max_iter; i++) {
    z = M.solve(r);
    x += z;
    r = b - A * x;
    
    if ((resid = norm(r) / normb) <= tol) {
      tol = resid;
      max_iter = i;
      return 0;
    }
  }

  tol = resid;
  return 1;
}

[jonnin]

matrix is not a c++ thing. Someone created a matrix class.
vector (note the spelling/caps!) is NOT a math vector or physics vector. It is just a container that is a wrapper for array type storage. Vector (again, spelling/caps!) as seen here however is also some sort of non-c++ thing that was created by this author.

It is difficult to get the algorithm from this code. I don't know what solve() does, and I am not familiar with this approach generally (ironically, and strangely, while I did 2 decades of linear algebra coding, we did not need ax = b at all, it was all weird stuff like AX +XB = C (all matrices) ... I never had to code a solver for it...

I suggest you hit google on this approach (look at iterative refinement for ax = b) and get the algorithm THAT way, and understand it at the math level. THEN review this code to see how it accomplishes that.

At a guess, the entire program you posted IS the relevant algorithm, apart from the function header and return codes, which if you pulled it apart would need to be adjusted.