CG Examples
[Assorted Thyra Operator/Vector Example Code]

Here we show some simple examples of using code in the Thyra package with an example ANA algorithm for the iterative solution of symmetric positive-definite linear systems using the conjugate gradient (CG) method. More...

Modules
group	Templated Serial Implementation of the CG Method
	Here is an example program that shows the use of the example serial templated matrix class SerialTridiagLinearOp with the example linear ANA implementation sillyCgSolve().
group	Templated MPI Implementation of the CG Method
	Here is an example program that shows the use of the example MPI templated matrix class MPITridiagLinearOp with the example linear ANA implementation sillyCgSolve().
Functions
template<class Scalar>
bool	sillyCgSolve (const Thyra::LinearOpBase< Scalar > &A, const Thyra::VectorBase< Scalar > &b, const int maxNumIters, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType tolerance, Thyra::VectorBase< Scalar > *x, std::ostream *out=NULL)
	Silly little example unpreconditioned CG solver.

---

Detailed Description

Here we show some simple examples of using code in the Thyra package with an example ANA algorithm for the iterative solution of symmetric positive-definite linear systems using the conjugate gradient (CG) method.

The CG ANA is implemented in the function sillyCgSolve() and its implementation is shown below:

template<class Scalar>
bool sillyCgSolve(
  const Thyra::LinearOpBase<Scalar>                              &A
  ,const Thyra::VectorBase<Scalar>                               &b
  ,const int                                                     maxNumIters
  ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType   tolerance
  ,Thyra::VectorBase<Scalar>                                     *x
  ,std::ostream                                                  *out          = NULL
  )
{
  using Teuchos::RefCountPtr;
  typedef Teuchos::ScalarTraits<Scalar>   ST;         // We need to use ScalarTraits to support arbitrary types.         
  typedef typename ST::magnitudeType      ScalarMag;  // This is the type returned from a vector norm.
  // Validate input
  TEST_FOR_EXCEPT(x==NULL);
  THYRA_ASSERT_LINEAR_OP_VEC_APPLY_SPACES("sillyCgSolve()",A,Thyra::NOTRANS,*x,&b); // A*x - b agree?
  Teuchos::EVerbosityLevel vl = Teuchos::VERB_MEDIUM; // Set the verbosity level
  if(out) *out << "\nStarting CG solver ...\n" << std::scientific << "\ndescribe A:\n"<<describe(A,vl)
               << "\ndescribe b:\n"<<describe(b,vl)<<"\ndescribe x:\n"<<describe(*x,vl)<<"\n";
  // Get the vector space (domain and range spaces should be the same)
  RefCountPtr<const Thyra::VectorSpaceBase<Scalar> > space = A.domain();
  // Compute initial residual : r = b - A*x
  RefCountPtr<Thyra::VectorBase<Scalar> > r = createMember(space);                     
  Thyra::assign(&*r,b);                                              // r = b
  Thyra::apply(A,Thyra::NOTRANS,*x,&*r,Scalar(-ST::one()),ST::one());// r = -A*x + r
  const ScalarMag r0_nrm = Thyra::norm(*r);      // Compute ||r0|| = sqrt(<r0,r0>) for convergence test
  if(r0_nrm == ST::zero()) return true;          // Trivial RHS and initial LHS guess?
  // Create workspace vectors and scalars
  RefCountPtr<Thyra::VectorBase<Scalar> > p = createMember(space), q = createMember(space);
  Scalar rho_old;
  // Perform the iterations
  for( int iter = 0; iter <= maxNumIters; ++iter ) {
    // Check convergence and output iteration
    const ScalarMag r_nrm = Thyra::norm(*r);          // Compute ||r|| = sqrt(<r,r>)
    const bool isConverged = r_nrm/r0_nrm <= tolerance;
    if( iter%(maxNumIters/10+1) == 0 || iter == maxNumIters || isConverged ) {
      if(out) *out << "Iter = " << iter << ", ||b-A*x||/||b-A*x0|| = " << (r_nrm/r0_nrm) << std::endl;
      if( r_nrm/r0_nrm < tolerance ) return true;     // Converged to tolerance, Success!
    }
    // Compute iteration
    const Scalar rho = space->scalarProd(*r,*r);      // <r,r>              -> rho
    if(iter==0) Thyra::assign(&*p,*r);                // r                  -> p   (iter == 0)
    else Thyra::Vp_V( &*p, *r, Scalar(rho/rho_old) ); // r+(rho/rho_old)*p  -> p   (iter  > 0)
    Thyra::apply(A,Thyra::NOTRANS,*p,&*q);            // A*p                -> q
    const Scalar alpha = rho/space->scalarProd(*p,*q);// rho/<p,q>          -> alpha
    Thyra::Vp_StV( x,   Scalar(+alpha), *p );         // +alpha*p + x       -> x
    Thyra::Vp_StV( &*r, Scalar(-alpha), *q );         // -alpha*q + r       -> r
    rho_old = rho;                                    // rho                -> rho_old (remember rho for next iter)
  }
  return false; // Failure
} // end sillyCgSolve

This above templated function is used in the following various example implementations which use several different scalar types:

1. Templated Serial Implementation of the CG Method

2. Templated MPI Implementation of the CG Method

---

Function Documentation

template<class Scalar> bool sillyCgSolve (  const Thyra::LinearOpBase< Scalar > &  A, const Thyra::VectorBase< Scalar > &  b, const int  maxNumIters, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType  tolerance, Thyra::VectorBase< Scalar > *  x, std::ostream *  out = NULL )

Silly little example unpreconditioned CG solver. This little function is just a silly little ANA that implements the CG (conjugate gradient) method for solving symmetric positive definite systems using the foundational Thyra operator/vector interfaces. This function is small and is meant to be looked at so study its implementation by clicking on the below link to its definition. Examples: sillyCgSolve.hpp, sillyCgSolve_mpi.cpp, and sillyCgSolve_serial.cpp. Definition at line 47 of file sillyCgSolve.hpp.

---