# 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:

## 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  1.4.7