cheby.h

Go to the documentation of this file.

// $Id: cheby.h,v 1.9.2.4 2019/05/28 11:13:02 garloff Exp $

#ifndef TBCI_SOLVER_CHEBY_H
#define TBCI_SOLVER_CHEBY_H

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

NAMESPACE_TBCI

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


/*statt:
template < class Matrix, class Vector, class Preconditioner, class Real,
           class Type >
int 
CHEBY(const Matrix &A, Vector &x, const Vector &b,
      const Preconditioner &M, int &max_iter, Real &tol,
      Type eigmin, Type eigmax)*/

#ifndef HAVE_BCXX_TYPENAME_BUG
# define VALUE_TYPE typename SysVector::value_type
#else 
# define VALUE_TYPE SysVector::value_type
#endif


template < typename T, typename SysMatrix, typename SysVector >
int CHEBY(const SysMatrix &A, SysVector &x, const SysVector &b,
          const Preconditioner_Sig<T, SysMatrix> &M, unsigned int &max_iter, double &tol,
          T eigmin, T eigmax)
{
  double resid;
  unsigned int dim = A.rows();
  T alpha(0), beta, c, d;
  SysVector p(dim), q(dim), z(dim);

  double normb = b.fabs();
  SysVector r(b - A * x);

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

  c = (eigmax - eigmin) / 2.0;
  d = (eigmax + eigmin) / 2.0;

  for (unsigned int i = 1; i <= max_iter; i++) 
  {
    z = M.solve(r);

    if (i == 1) 
    {
      p = z;
      alpha = 2.0 / d;
    } else 
    {
      beta = c * alpha / 2.0;       // calculate new beta
      beta = beta * beta;
      alpha = 1.0 / (d - beta);     // calculate new alpha
      p = z + (VALUE_TYPE)beta * p;             // update search direction
    }

    q = A * p;
    x += VALUE_TYPE (alpha) * p;                 // update approximation vector
    r -= VALUE_TYPE (alpha) * q;                 // compute residual

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

  tol = resid;
  return 1;                         // no convergence
}

NAMESPACE_END

#endif  /* TBCI_SOLVER_CHEBY_H */