[HELP] User-defined header in my CPP program

[jemeripol]

Hi guys,

I need a bit of your help again.

I have found this template, named CG.h (which obviously is a header file)

in my main.cpp,

how am i going to use this code to function as it is designed?

any idea will be appreciated. Thank you. :)

Here is the code:

//*****************************************************************
// Iterative template routine -- CG
//
// CG solves the symmetric positive definite linear
// system Ax=b using the Conjugate Gradient method.
//
// CG follows the algorithm described on p. 15 in the 
// SIAM Templates book.
//
// 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 CG(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)
{
  Real resid;
  Vector p, z, q;
  Vector alpha(1), beta(1), rho(1), rho_1(1);

  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);
    rho(0) = dot(r, z);
    
    if (i == 1)
      p = z;
    else {
      beta(0) = rho(0) / rho_1(0);
      p = z + beta(0) * p;
    }
    
    q = A*p;
    alpha(0) = rho(0) / dot(p, q);
    
    x += alpha(0) * p;
    r -= alpha(0) * q;
    
    if ((resid = norm(r) / normb) <= tol) {
      tol = resid;
      max_iter = i;
      return 0;     
    }

    rho_1(0) = rho(0);
  }
  
  tol = resid;
  return 1;
}

[TwilightSpectre]

The same as normal.

Assuming you have some set of matrix / vector classes (note - physics vector, not std::vector) you just include the file with #include "CG.h" and call the function with appropriate parameters.

[jemeripol]

Can you give me an example?

something like

int main()
{
      //What can i write here? :)
}

Can i call it like a normal function?

like ?

CG(a[2][2]... )

[TwilightSpectre]

Sure. Here is an example (with some made-up terms, because I'm not sure exactly what maths library it is using):

#include "someMathLibrary.h"
#include "CG.h"

int main() {
    // Note I'm making up all these types. 
    Mat3f A;
    Vec2f x;
    Vec2f b;
    Preconditionor M;
    float result;

    CG(A, x, b, M, 100, result);
    // the result is now stored in 'x' and 'result'.
    // Test the return value of 'CG' to determine if it failed or not.

    return 0;
}

[jemeripol]

From what I am seeing,
A, x, and b are now declared as a variable of certain type, which is recognized only by the header CG.h or someMathLibrary.h..

Which means, you can now use it like normal variables?

Therefore, I need to study the construction of the user-defined headers to use it. :)

Wow. Thank you the help sir. :)
I owe you much..
I'll try to make some codes from it.