BCGStab Class Reference

#include <BCGStab.h>

Inheritance diagram for BCGStab:

Collaboration diagram for BCGStab:

Public Member Functions
	BCGStab ()
virtual	~BCGStab ()
virtual MFRPtr	solve (LinearSystem &ls)
virtual void	cleanup ()
virtual void	smooth (LinearSystem &ls)
int	getTotalIterations () const
	DEFINE_TYPENAME ("BCGStab")
Public Member Functions inherited from LinearSolver
	LinearSolver ()

Public Attributes
LinearSolver *	preconditioner
Public Attributes inherited from LinearSolver
int	nMaxIterations
int	verbosity
double	relativeTolerance
double	absoluteTolerance

Private Member Functions
	BCGStab (const BCGStab &)

Private Attributes
int	_totalIterations

Detailed Description

Solve a linear system using stabilized bi conjugate-gradient method.

Definition at line 20 of file BCGStab.h.

Constructor & Destructor Documentation

BCGStab::BCGStab

(

)

Definition at line 11 of file BCGStab.cpp.

                 :
  preconditioner(0),
  _totalIterations(0)
{}

BCGStab::~BCGStab ( )

virtual

Definition at line 16 of file BCGStab.cpp.

{
}

BCGStab::BCGStab ( const BCGStab &  )

private

Member Function Documentation

void BCGStab::cleanup ( )

virtual

Implements LinearSolver.

Definition at line 21 of file BCGStab.cpp.

References LinearSolver::cleanup(), and preconditioner.

{
  preconditioner->cleanup();
}

BCGStab::DEFINE_TYPENAME

(

"BCGStab"

)

int BCGStab::getTotalIterations ( ) const

inline

Definition at line 31 of file BCGStab.h.

{ return _totalIterations;}

void BCGStab::smooth ( LinearSystem &  ls )

virtual

Implements LinearSolver.

Definition at line 173 of file BCGStab.cpp.

{
  throw CException("cannot use BCGStab as preconditioner");
}

MFRPtr BCGStab::solve ( LinearSystem &  ls )

virtual

Implements LinearSolver.

Definition at line 27 of file BCGStab.cpp.

References _totalIterations, LinearSolver::absoluteTolerance, MultiFieldMatrix::computeResidual(), LinearSystem::getB(), LinearSystem::getBPtr(), LinearSystem::getDelta(), LinearSystem::getDeltaPtr(), LinearSystem::getMatrix(), MultiField::getOneNorm(), LinearSystem::getResidual(), MultiField::msaxpy(), MultiFieldMatrix::multiply(), MultiField::newCopy(), LinearSolver::nMaxIterations, preconditioner, LinearSolver::relativeTolerance, LinearSystem::replaceB(), LinearSystem::replaceDelta(), LinearSolver::smooth(), LinearSolver::verbosity, and MultiField::zero().

{
  const MultiFieldMatrix& matrix = ls.getMatrix();

  // original system is in delta form
  shared_ptr<MultiField> x(ls.getDeltaPtr());
  shared_ptr<MultiField> bOrig(ls.getBPtr());

  matrix.computeResidual(ls.getDelta(),ls.getB(),ls.getResidual());

  MFRPtr rNorm0(ls.getResidual().getOneNorm());

  shared_ptr<MultiField> r(dynamic_pointer_cast<MultiField>(ls.getResidual().newCopy()));
  shared_ptr<MultiField> rTilda(dynamic_pointer_cast<MultiField>(ls.getResidual().newCopy()));

  shared_ptr<MultiField> p;
  shared_ptr<MultiField> pHat(dynamic_pointer_cast<MultiField>(x->newClone()));
  shared_ptr<MultiField> v(dynamic_pointer_cast<MultiField>(x->newClone()));
  shared_ptr<MultiField> t(dynamic_pointer_cast<MultiField>(x->newClone()));

  MFRPtr rho;
  MFRPtr rhoPrev;

  MFRPtr alpha;
  MFRPtr omega;

#ifndef  FVM_PARALLEL
  if (verbosity >0)
    cout << 0 << ": " << *rNorm0 << endl;
#endif

#ifdef  FVM_PARALLEL
  if (verbosity >0 && MPI::COMM_WORLD.Get_rank() == 0)
    cout << 0 << ": " << *rNorm0 << endl;
#endif

  for(int i = 0; i<nMaxIterations; i++)
  { 
      _totalIterations++;
      rhoPrev = rho;
      rho = r->dotWith(*rTilda);

      rho->reduceSum();

      if (!p)
      {
          p = dynamic_pointer_cast<MultiField>(r->newCopy());
      }
      else
      {
          MFRPtr rhoRatio = (*rho) / (*rhoPrev);
          MFRPtr alphaByOmega = (*alpha) / (*omega);
          MFRPtr beta = (*rhoRatio) * (*alphaByOmega);
          p->msaxpy(*omega,*v);
          *p *= *beta;
          *p += *r;
      }

      pHat->zero();
      ls.replaceDelta(pHat);
      ls.replaceB(p);

      preconditioner->smooth(ls);

      matrix.multiply(*v,*pHat);

      MFRPtr rtv = rTilda->dotWith(*v);
      rtv->reduceSum();

      alpha = (*rho)/(*rtv);

      x->msaxpy(*alpha,*pHat);
      r->msaxpy(*alpha,*v);

      MFRPtr rNorm = r->getOneNorm();

      if (*rNorm < absoluteTolerance)
      {
          break;
      }

      shared_ptr<MultiField> sHat = pHat;
      sHat->zero();
      ls.replaceDelta(sHat);
      ls.replaceB(r);

      preconditioner->smooth(ls);

      matrix.multiply(*t,*sHat);

      MFRPtr tdotr = t->dotWith(*r);
      MFRPtr tdott = t->dotWith(*t);

      tdotr->reduceSum();
      tdott->reduceSum();

      omega = (*tdotr) / (*tdott);

      x->msaxpy(*omega,*sHat);
      r->msaxpy(*omega,*t);

      rNorm = r->getOneNorm();

      MFRPtr normRatio(rNorm->normalize(*rNorm0));

#ifndef FVM_PARALLEL
      if (verbosity >0)
        cout << i+1 << ": " << *rNorm << endl;
#endif

#ifdef  FVM_PARALLEL
  if (verbosity >0 && MPI::COMM_WORLD.Get_rank() == 0)
    cout << i+1 << ": " << *rNorm << endl;
#endif

      if (*rNorm < absoluteTolerance || *normRatio < relativeTolerance)
        break;


  }
  ls.replaceDelta(x);
  ls.replaceB(bOrig);
#ifdef FVM_PARALLEL
  x->sync();
#endif

  matrix.computeResidual(ls.getDelta(),ls.getB(),ls.getResidual());
  MFRPtr rNormn(ls.getResidual().getOneNorm());

#ifndef FVM_PARALLEL
  if (verbosity >0)
    cout << "n" << ": " << *rNormn << endl;
#endif


#ifdef  FVM_PARALLEL
  if (verbosity >0 && MPI::COMM_WORLD.Get_rank() == 0)
    cout << "n" << ": " << *rNormn << endl;
#endif



  return rNorm0;
}

Member Data Documentation

int BCGStab::_totalIterations

private

Definition at line 40 of file BCGStab.h.

Referenced by solve().

LinearSolver* BCGStab::preconditioner

Definition at line 35 of file BCGStab.h.

Referenced by cleanup(), and solve().

---

The documentation for this class was generated from the following files:

• BCGStab.h
• BCGStab.cpp