CG Examples
[Assorted Thyra Operator/Vector Example Code]

Collaboration diagram for CG Examples:

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

 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().
 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, const Teuchos::Ptr< Thyra::VectorBase< Scalar > > &x, std::ostream &out)
 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, const Teuchos::Ptr< Thyra::VectorBase< Scalar > > &x_inout, std::ostream &out)
 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)
 Silly little example unpreconditioned CG solver that uses handles.
template<class Scalar>
void sillyModifiedGramSchmidt (Thyra::MultiVectorBase< Scalar > *V_inout, Teuchos::RCP< Thyra::MultiVectorBase< Scalar > > *R_out)
 Silly little implementation of the modified Gram-Schmidt algorithm to compute a QR factorization V=Q*R of a multi-vector V.

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,
  const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x,
  std::ostream &out
  )
{

  // 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 Teuchos::as;
  using Teuchos::RCP; using Thyra::VectorSpaceBase; using Thyra::VectorBase;
  using Thyra::NOTRANS; using Thyra::V_V; using Thyra::apply;
  

  // Validate input
  THYRA_ASSERT_LINEAR_OP_VEC_APPLY_SPACES("sillyCgSolve()", A, Thyra::NOTRANS, *x, &b);
  Teuchos::EVerbosityLevel vl = Teuchos::VERB_MEDIUM;

  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
  const RCP<const VectorSpaceBase<Scalar> > space = A.domain();
  const RCP<VectorBase<Scalar> > r = createMember(space);
  // r = -A*x + b
  V_V(r.ptr(), b); apply<Scalar>(A, NOTRANS, *x, r.ptr(), -one, one);
  const ScalarMag r0_nrm = norm(*r);
  if (r0_nrm==zero) return true;
  const RCP<VectorBase<Scalar> > p = createMember(space), q = createMember(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 << "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 = inner(*r, *r);        // <r,r>              -> rho
    if (iter==0) V_V(p.ptr(), *r);           // r                  -> p   (iter == 0)
    else Vp_V( p.ptr(), *r, rho/rho_old );   // r+(rho/rho_old)*p  -> p   (iter  > 0)
    apply<Scalar>(A, NOTRANS, *p, q.ptr());  // A*p                -> q
    const Scalar alpha = rho/inner(*p, *q);  // rho/<p,q>          -> alpha
    Vp_StV( x, +alpha, *p );                 // +alpha*p + x       -> x
    Vp_StV( r.ptr(), -alpha, *q );           // -alpha*q + r       -> r
    rho_old = rho;                           // rho                -> rho_old (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
  )
{
 
 // 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<<"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,
const Teuchos::Ptr< Thyra::VectorBase< Scalar > > &  x,
std::ostream &  out 
)

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 49 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,
const Teuchos::Ptr< Thyra::VectorBase< Scalar > > &  x_inout,
std::ostream &  out 
)

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 51 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 
)

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 51 of file silliestCgSolve.hpp.

template<class Scalar>
void sillyModifiedGramSchmidt ( Thyra::MultiVectorBase< Scalar > *  V_inout,
Teuchos::RCP< Thyra::MultiVectorBase< Scalar > > *  R_out 
)

Silly little implementation of the modified Gram-Schmidt algorithm to compute a QR factorization V=Q*R of a multi-vector V.

Parameters:
V [in/out] On input, contains the columns to compute the factorization for. On output, contains the columns of Q.
R [out] On output, contains the upper triangular matrix R.
ToDo: Finish documentation!
Examples:
sillyModifiedGramSchmidt.hpp.

Definition at line 49 of file sillyModifiedGramSchmidt.hpp.


Generated on Wed May 12 21:26:55 2010 for Thyra Operator/Vector Support by  doxygen 1.4.7