bicg.h File Reference

Go to the source code of this file.

Functions
template<class Matrix , class Vector , class Preconditioner , class Real >
int	BiCG (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)

Function Documentation

template<class Matrix , class Vector , class Preconditioner , class Real >

int BiCG	(	const Matrix &	A,
		Vector &	x,
		const Vector &	b,
		const Preconditioner &	M,
		int &	max_iter,
		Real &	tol
	)

Definition at line 25 of file bicg.h.

{
  Real resid;
  Vector rho_1(1), rho_2(1), alpha(1), beta(1);
  Vector z, ztilde, p, ptilde, q, qtilde;

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

  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);
    ztilde = M.trans_solve(rtilde);
    rho_1(0) = dot(z, rtilde);
    if (rho_1(0) == 0) { 
      tol = norm(r) / normb;
      max_iter = i;
      return 2;
    }
    if (i == 1) {
      p = z;
      ptilde = ztilde;
    } else {
      beta(0) = rho_1(0) / rho_2(0);
      p = z + beta(0) * p;
      ptilde = ztilde + beta(0) * ptilde;
    }
    q = A * p;
    qtilde = A.trans_mult(ptilde);
    alpha(0) = rho_1(0) / dot(ptilde, q);
    x += alpha(0) * p;
    r -= alpha(0) * q;
    rtilde -= alpha(0) * qtilde;

    rho_2(0) = rho_1(0);
    if ((resid = norm(r) / normb) < tol) {
      tol = resid;
      max_iter = i;
      return 0;
    }
  }

  tol = resid;
  return 1;
}

Here is the call graph for this function: