Main Page | Modules | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Class Members | Examples | Search for 

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 ExampleTridiagSerialLinearOp with the example linear ANA implementation sillyCgSolve().
group  Templated SPMD Implementation of the CG Method
 Here is an example program that shows the use of the example SPMD templated matrix class ExampleTridiagSpmdLinearOp with the example linear ANA implementation sillyCgSolve() or silliestCgSolve().

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.
template<class Scalar>
bool sillierCgSolve (const Thyra::LinearOpBase< Scalar > &A_in, const Thyra::VectorBase< Scalar > &b_in, const int maxNumIters, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType tolerance, Thyra::VectorBase< Scalar > *x_inout, std::ostream *out=NULL)
 Silly little example unpreconditioned CG solver (calls templated code).
template<class Scalar>
bool silliestCgSolve (Thyra::ConstLinearOperator< Scalar > const &A, Thyra::ConstVector< Scalar > const &b, const int maxNumIters, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType tolerance, Thyra::Vector< Scalar > &x, std::ostream *out=NULL)
 Silly little example unpreconditioned CG solver that uses handles.

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
  )
{
  // Create some typedefs and some other stuff to make the code cleaner
  typedef Teuchos::ScalarTraits<Scalar> ST; typedef typename ST::magnitudeType ScalarMag;
  const Scalar one = ST::one(), zero = ST::zero(); using Thyra::NOTRANS;
  typedef Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<Scalar> > VectorSpacePtr;
  typedef Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > VectorPtr;
  // Validate input
  TEST_FOR_EXCEPT(x==NULL);
  THYRA_ASSERT_LINEAR_OP_VEC_APPLY_SPACES("sillyCgSolve()",A,Thyra::NOTRANS,*x,&b); // Does A*x - b agree?
  Teuchos::EVerbosityLevel vl = Teuchos::VERB_MEDIUM;
  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";
  // Initialization
  VectorSpacePtr space = A.domain();
  VectorPtr r = createMember(space);
  V_V(&*r,b); apply(A,NOTRANS,*x,&*r,Scalar(-one),one); // r = -A*x + b
  const ScalarMag r0_nrm = norm(*r);
  if(r0_nrm==zero) return true;
  VectorPtr 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 = norm(*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; // Success!
    }
    // Compute iteration
    const Scalar rho = scalarProd(*r,*r);         // <r,r>              -> rho
    if(iter==0) V_V(&*p,*r);                      // r                  -> p   (iter == 0)
    else Vp_V( &*p, *r, Scalar(rho/rho_old) );    // r+(rho/rho_old)*p  -> p   (iter  > 0)
    apply(A,NOTRANS,*p,&*q);                      // A*p                -> q
    const Scalar alpha = rho/scalarProd(*p,*q);   // rho/<p,q>          -> alpha
    Vp_StV( x,   Scalar(+alpha), *p );            // +alpha*p + x       -> x
    Vp_StV( &*r, Scalar(-alpha), *q );            // -alpha*q + r       -> r
    rho_old = rho;                                // rho                -> rho_old (remember for next iter)
  }
  return false; // Failure
} // end sillyCgSolve

Another version of this CG algorithm is demonstrated in the below function:

template<class Scalar>
bool silliestCgSolve(
  Thyra::ConstLinearOperator<Scalar>                           const& A
  ,Thyra::ConstVector<Scalar>                                  const& b
  ,const int                                                          maxNumIters
  ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType        tolerance
  ,Thyra::Vector<Scalar>                                            & x
  ,std::ostream                                                     * out = NULL
  )
{
  // Create some typedefs and inject some names into local namespace
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType ScalarMag;
  const Scalar one = ST::one(), zero = ST::zero();
  using Thyra::VectorSpace;
  using Thyra::Vector;
  // Initialization of the algorithm
  const VectorSpace<Scalar>     space = A.domain();
  Vector<Scalar>                r = b - A*x;
  ScalarMag                     r0_nrm = norm(r);
  if(r0_nrm==zero) return true;
  Vector<Scalar>                p(space), q(space);
  Scalar                        rho_old = -one;
  // Perform the iterations
  for( int iter = 0; iter <= maxNumIters; ++iter ) {
    // Check convergence and output iteration
    const ScalarMag    r_nrm = norm(r);
    const bool         isConverged = ( (r_nrm/r0_nrm) <= tolerance );
    if( ( iter%(maxNumIters/10+1) == 0 || iter == maxNumIters || isConverged ) && out )
      *out<<"Iter = "<<iter<<", ||b-A*x||/||b-A*x0|| = "<<(r_nrm/r0_nrm)<<std::endl;
    if( r_nrm/r0_nrm < tolerance ) return true; // Success!
    // Compute the iteration
    const Scalar      rho = inner(r,r);
    if(iter==0)       copyInto(r,p);
    else              p = Scalar(rho/rho_old)*p + r;
    q = A*p;
    const Scalar      alpha = rho/inner(p,q);
    x += Scalar(+alpha)*p;
    r += Scalar(-alpha)*q;
    rho_old = rho;
  }
  return false; // Failure
} // end silliestCgSolve

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

  1. Templated Serial Implementation of the CG Method

  2. Templated SPMD 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 48 of file sillyCgSolve.hpp.

template<class Scalar>
bool sillierCgSolve (  const Thyra::LinearOpBase< Scalar > &  A_in,
const Thyra::VectorBase< Scalar > &  b_in,
const int  maxNumIters,
const typename Teuchos::ScalarTraits< Scalar >::magnitudeType  tolerance,
Thyra::VectorBase< Scalar > *  x_inout,
std::ostream *  out = NULL
) 
 

Silly little example unpreconditioned CG solver (calls templated code).

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:
sillierCgSolve.hpp, and sillyCgSolve_serial.cpp.

Definition at line 50 of file sillierCgSolve.hpp.

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

Silly little example unpreconditioned CG solver that uses handles.

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:
sillierCgSolve.hpp, and silliestCgSolve.hpp.

Definition at line 49 of file silliestCgSolve.hpp.

Generated on Thu Sep 18 12:32:32 2008 for Thyra Operator/Vector Support by doxygen 1.3.9.1