Power Method Examples
[Assorted Thyra Operator/Vector Example Code]

Here we show some simple examples of using code in the Thyra package with an example linear ANA algorithm for the power method for estimating the dominate eigen value of a matrix. More...

Modules
group	Templated Serial Implementation of the Power Method
	Here is an example program that shows the use of the example serial templated matrix class SerialTridiagLinearOp with the example linear ANA implementation sillyPowerMethod().
Functions
template<class Scalar>
bool	sillyPowerMethod (const Thyra::LinearOpBase< Scalar > &A, const int maxNumIters, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType tolerance, Scalar *lambda, std::ostream *out=NULL)
	Silly little example power method abstract numerical algorithm (ANA).

---

Detailed Description

Here we show some simple examples of using code in the Thyra package with an example linear ANA algorithm for the power method for estimating the dominate eigen value of a matrix.

These example programs are meant to mimic the power method example program shown in the Epetra documentation.

The power method ANA is implemented in the function sillyPowerMethod() and its implementation is shown below:

template<class Scalar>
bool sillyPowerMethod(
  const Thyra::LinearOpBase<Scalar>                              &A
  ,const int                                                     maxNumIters
  ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType   tolerance
  ,Scalar                                                        *lambda
  ,std::ostream                                                  *out          = NULL
  )
{
  typedef Teuchos::ScalarTraits<Scalar>   ST;         // We need to use ScalarTraits to support arbitrary types!         
  typedef typename ST::magnitudeType      ScalarMag;  // This is the type for a norm
  using Teuchos::arrayArg;
  if(out) *out << "\nStarting power method ...\n\n";
  // Create workspace vectors
  Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> >
    q = createMember(A.domain()),
    z = createMember(A.range()),
    r = createMember(A.range());
  // Randomize initial z
  Thyra::seed_randomize<Scalar>(0); // Make repeated runs with same data unique
  Thyra::randomize( Scalar(-ST::one()), Scalar(+ST::one()), &*z );
  // Perform iterations
  for( int iter = 0; iter < maxNumIters; ++iter ) {
    const ScalarMag z_nrm = Thyra::norm(*z);          // Compute natural norm of z
    Thyra::V_StV( &*q, Scalar(ST::one()/z_nrm), *z ); // q = (1/||z}*z 
    Thyra::apply( A, Thyra::NOTRANS , *q, &*z );      // z = A*q
    *lambda = A.range()->scalarProd(*q,*z);           // lambda = <q,z> : Approximate maximum absolute eigenvalue
    if( iter%(maxNumIters/10) == 0 || iter+1 == maxNumIters ) {
      Thyra::linear_combination(                      // r = z - lambda*q : Compute residual of eigenvalue equation
        2,arrayArg<Scalar>(ST::one(),-*lambda)()
        ,arrayArg<const Thyra::VectorBase<Scalar>*>(&*z,&*q)()
        ,ST::zero(),&*r
        );
      const ScalarMag r_nrm = Thyra::norm(*r);        // Compute natural norm of r
      if(out) *out << "Iter = " << iter << ", lambda = " << (*lambda) << ", ||A*q-lambda*q|| = " << r_nrm << std::endl;
      if( r_nrm < tolerance ) return true;            // Success!
    }
  }
  return false; // Failure
} // end sillyPowerMethod

The above templated function sillyPowerMethod() is used in the following various example implementations which use several different scalar types:

1. Templated Serial Implementation of the Power Method

---

Function Documentation

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

Silly little example power method abstract numerical algorithm (ANA). This little function is just a silly little ANA that implements the power method for estimating the dominate eigenvalue for a matrix 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: sillyPowerMethod.hpp, and sillyPowerMethod_serial.cpp. Definition at line 47 of file sillyPowerMethod.hpp.

---