cg.h

Go to the documentation of this file.

/* $Id: cg.h,v 1.10.2.3 2019/05/28 11:13:02 garloff Exp $ */

#ifndef TBCI_SOLVER_CG_H
#define TBCI_SOLVER_CG_H

#include "tbci/basics.h"
#include "tbci/solver/precond.h"

NAMESPACE_TBCI


INST3(template <typename T, Matrix<T>, Vector<T> > class NN friend \
int CG (const Matrix<T>&, Vector<T>&, const Vector<T>&,\
        const Preconditioner_Sig<T, Matrix<T> >&, unsigned int&, double&);)
INST3(template <typename T, BdMatrix<T>, Vector<T> > class NN friend \
int CG (const BdMatrix<T>&, Vector<T>&, const Vector<T>&,\
        const Preconditioner_Sig<T, BdMatrix<T> >&, unsigned int&, double&);)


template < typename T, typename SysMatrix, typename SysVector >
int CG (const SysMatrix &A, SysVector &x, const SysVector &b,
        const Preconditioner_Sig<T, SysMatrix> &M, unsigned int &max_iter, double &tol)
{
  double resid;
  unsigned int dim = A.rows();
  SysVector p(dim), z(dim), q(dim);
  typename SysVector::value_type alpha, beta, rho, rho_1(0);

  double normb = b.fabs();
  if (normb == 0.0)   normb = 1.0; else normb = 1.0 / normb;

  SysVector r(b - A*x);

  
  if ((resid = r.fabs() * normb) <= tol) 
  { 
    tol = resid;
    max_iter = 0;
    return 0;
  }

//Aenderung Marco:
  for (unsigned int i = 1; i <= max_iter; i++) 
  {
    z = M.solve(r);
    rho = dot(r,z);
    
    if (i == 1)
      p = z;
    else {
      beta = rho/ rho_1;
      p = z + beta * p;
    }
    
    q = A*p;
    alpha = rho/ dot(p,q);
    
    x += alpha * p;
    r -= alpha * q;
    
    if ((resid = r.fabs() * normb) <= tol) 
    {
      tol = resid;
      max_iter = i;
      return 0;     
    }

    rho_1 = rho;
  }
  
  tol = resid;
  return 1;
}

  
INST3(template <typename T, Matrix<T>, Vector<T> > class NN friend \
int CG2 (const Matrix<T>&, Vector<T>&, const Vector<T>&,\
        const Preconditioner_Sig<T, Matrix<T> >&, unsigned int&, double&);)
INST3(template <typename T, BdMatrix<T>, Vector<T> > class NN friend \
int CG2 (const BdMatrix<T>&, Vector<T>&, const Vector<T>&,\
        const Preconditioner_Sig<T, BdMatrix<T> >&, unsigned int&, double&);)

template < typename T, typename SysMatrix, typename SysVector>
int CG2(const SysMatrix &A, SysVector &x, const SysVector &b,
        const Preconditioner_Sig<T, SysMatrix> &M, unsigned int &max_iter, double &tol)
{
  double resid;
  unsigned int dim = A.rows();
  SysVector p(dim), z(dim), q(dim);
  typename SysVector::value_type alpha, beta, rho, rho_1(0);

  double normb = b.fabs();
  if (normb == 0.0)   normb = 1; else normb = 1 / normb;

  SysVector r(b - A*x);

  
  if ((resid = r.fabs() * normb) <= tol) 
  { 
    tol = resid;
    max_iter = 0;
    return 0;
  }

//Aenderung Marco:
  for (unsigned int i = 1; i <= max_iter; i++) 
  {
    z = M.solve(r);
    rho = (r*z);
    
    if (i == 1)
      p = z;
    else {
      beta = rho/ rho_1;
      p = z + beta * p;
    }
    
    q = A*p;
    alpha = rho/ (p*q);
    
    x += alpha * p;
    r -= alpha * q;
    
    if ((resid = r.fabs() * normb) <= tol) 
    {
      tol = resid;
      max_iter = i;
      return 0;     
    }

    rho_1 = rho;
  }
  
  tol = resid;
  return 1;
}

NAMESPACE_END

#endif  /* TBCI_SOLVER_CG_H */