CG Class Reference

#include <CG.h>

Inheritance diagram for CG:

Collaboration diagram for CG:

Public Member Functions
	CG ()
virtual	~CG ()
virtual MFRPtr	solve (LinearSystem &ls)
virtual void	cleanup ()
virtual void	smooth (LinearSystem &ls)
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
	CG (const CG &)

Detailed Description

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

Definition at line 20 of file CG.h.

Constructor & Destructor Documentation

CG::CG

(

)

Definition at line 11 of file CG.cpp.

       :
  preconditioner(0)
{}

CG::~CG ( )

virtual

Definition at line 15 of file CG.cpp.

{
}

CG::CG ( const CG &  )

private

Member Function Documentation

void CG::cleanup ( )

virtual

Implements LinearSolver.

Definition at line 20 of file CG.cpp.

References LinearSolver::cleanup(), and preconditioner.

{
  preconditioner->cleanup();
}

void CG::smooth ( LinearSystem &  ls )

virtual

Implements LinearSolver.

Definition at line 142 of file CG.cpp.

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

MFRPtr CG::solve ( LinearSystem &  ls )

virtual

Implements LinearSolver.

Definition at line 26 of file CG.cpp.

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

{
  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> p;
  shared_ptr<MultiField> z(dynamic_pointer_cast<MultiField>(x->newClone()));
  shared_ptr<MultiField> q(dynamic_pointer_cast<MultiField>(x->newClone()));

  MFRPtr rho;
  MFRPtr rhoPrev;

  MFRPtr alpha;
  MFRPtr beta;

#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++)
  {
      z->zero();
      ls.replaceDelta(z);
      ls.replaceB(r);

      preconditioner->smooth(ls);


      rhoPrev = rho;
      rho = r->dotWith(*z);

      rho->reduceSum();

      if (!p)
      {
          p = dynamic_pointer_cast<MultiField>(z->newCopy());
      }
      else
      {
          MFRPtr beta = (*rho) / (*rhoPrev);
          *p *= *beta;
          *p += *z;
      }

      matrix.multiply(*q,*p);

      MFRPtr ptq = p->dotWith(*q);
      ptq->reduceSum();

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

      x->msaxpy(*alpha,*p);
      r->msaxpy(*alpha,*q);

      MFRPtr 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

LinearSolver* CG::preconditioner

Definition at line 31 of file CG.h.

Referenced by cleanup(), and solve().

---

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

• CG.h
• CG.cpp