cgs.h

Go to the documentation of this file.

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

#ifndef TBCI_SOLVER_CGS_H
#define TBCI_SOLVER_CGS_H  

#include "tbci/basics.h"

#ifdef SOLVER_LOG
# include <fstream>
# include <iomanip>
#endif

NAMESPACE_TBCI

/* NOTE: Template parameter T and T_tol should really be double (or 
 * similar: float, long double are OK), as fabs() is used, which 
 * returns double. This could be changed to abs(), but will then cause 
 * some trouble for complex numbers.
 */


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


template < typename T, typename SysMatrix, typename SysVector >
int CGS(const SysMatrix &A, SysVector &x, const SysVector &b,
        const Preconditioner_Sig<T, SysMatrix> &M, unsigned int &max_iter, 
        double &tol, const unsigned off = 0)
{
  const unsigned int dim = A.rows();
  SysVector p(dim), phat(dim), q(dim), qhat(dim), vhat(dim), u(dim), uhat(dim);
  typename SysVector::value_type rho_1, rho_2((typename SysVector::value_type)1), alpha, beta;

  double residsqr;
  double tolsqr = TBCI__ sqr (tol);
  double normbsqr = b.fabssqr();
  if (normbsqr < 1e-32) 
        normbsqr = 1e32; 
  else 
        normbsqr = 1.0 / normbsqr;

  SysVector r(b - A*x);
  SysVector rtilde(r);
  //M.update (A);
  
  if ((residsqr = r.fabssqr() * normbsqr) <= tolsqr) 
  {
    tol = MATH__ sqrt(residsqr);
    max_iter = 0;
    return 0;
  }
#ifdef SOLVER_LOG
  STD__ ofstream cnvg;
  if (off) cnvg.open ("cgs_cnvg.gnu", STD__ ios::app);
  else cnvg.open ("cgs_cnvg.gnu");
#endif

  for (unsigned int i = 1; i <= max_iter; i++) 
  {
#ifdef SOLVER_LOG
    cnvg << i+off << "\t" << MATH__ sqrt(residsqr) << "\n";
#endif
    rho_1 = dot(rtilde, r);
    if (rho_1 == (typename SysVector::value_type)0)
    {
      tol = MATH__ sqrt (residsqr);
      return 2;
    }
    if (i == 1) 
    {
      u = r;
      p = u;
    } else 
    {
      beta = rho_1 / rho_2;
      u = r + beta * q;
      p = u + beta * (q + beta * p);
    }
    phat = M.solve(p);
    vhat = A*phat;
    alpha = rho_1 / dot(rtilde, vhat);
    q = u - alpha * vhat;
    uhat = M.solve(u + q);
    x += alpha * uhat;
    qhat = A * uhat;
    r -= alpha * qhat;
    rho_2 = rho_1;
    if ((residsqr = r.fabssqr() * normbsqr) < tolsqr) 
    {
      tol = MATH__ sqrt(residsqr);
      max_iter = i;
#ifdef SOLVER_LOG
      cnvg << i+off << "\t" << MATH__ sqrt(residsqr) << STD__ endl; cnvg.close ();
#endif
      return 0;
    }
  }
  
  tol = MATH__ sqrt(residsqr);
  return 1;
}

NAMESPACE_END

#endif  /* TBCI_SOLVER_CGS_H */