Templated Serial Implementation of the Power Method
[Power Method Examples]

Here is an example program that shows the use of the example serial templated matrix class ExampleTridiagSerialLinearOp with the example linear ANA implementation sillyPowerMethod(). More...

Classes

class  ExampleTridiagSerialLinearOp< Scalar >
 Simple example subclass for serial tridiagonal matrices. More...

Detailed Description

Here is an example program that shows the use of the example serial templated matrix class ExampleTridiagSerialLinearOp with the example linear ANA implementation sillyPowerMethod().

This example program is contained in the source file:

 ./example/operator_vector/sillyPowerMethod_serial.cpp 

where ./ is the base source directory for Thyra (i.e. ???/Trilinos/packages/Thyra).

The class ExampleTridiagSerialLinearOp that derives from the base class Thyra::SpmdLinearOpBase is quite simple and its complete implementation looks like:

template<class Scalar>
class ExampleTridiagSerialLinearOp : public Thyra::SpmdLinearOpBase<Scalar> {
private:

  Thyra::Index         dim_;
  std::vector<Scalar>  lower_;   // size = dim - 1    
  std::vector<Scalar>  diag_;    // size = dim
  std::vector<Scalar>  upper_;   // size = dim - 1

public:

  using Thyra::SpmdLinearOpBase<Scalar>::euclideanApply;

  ExampleTridiagSerialLinearOp() : dim_(0) {}

  ExampleTridiagSerialLinearOp( const Thyra::Index dim, const Scalar lower[], const Scalar diag[], const Scalar upper[] )
    { this->initialize(dim,lower,diag,upper); }
  
  void initialize(
    const Thyra::Index              dim      // >= 2
    ,const Scalar                   lower[]  // size == dim - 1
    ,const Scalar                   diag[]   // size == dim
    ,const Scalar                   upper[]  // size == dim - 1
    )
    {
      TEST_FOR_EXCEPT( dim < 2 );
      this->setLocalDimensions(Teuchos::null,dim,dim); // Needed to setup range() and domain()
      dim_ = dim;
      lower_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) lower_[k] = lower[k];
      diag_.resize(dim);     for( int k = 0; k < dim;   ++k ) diag_[k]  = diag[k];
      upper_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) upper_[k] = upper[k];
    }

  // Overridden form Teuchos::Describable */

  std::string description() const
    {
      return (std::string("ExampleTridiagSerialLinearOp<") + Teuchos::ScalarTraits<Scalar>::name() + std::string(">"));
    }

protected:

  // Overridden from SingleScalarEuclideanLinearOpBase

  bool opSupported(Thyra::ETransp M_trans) const { return true; }  // This class supports everything!

  // Overridden from SpmdLinearOpBase

  void euclideanApply(
    const Thyra::ETransp                         M_trans
    ,const RTOpPack::ConstSubVectorView<Scalar>  &x_in
    ,const RTOpPack::SubVectorView<Scalar>       *y_out
    ,const Scalar                                alpha
    ,const Scalar                                beta
    ) const
    {
      typedef Teuchos::ScalarTraits<Scalar> ST;
      // Get raw pointers to the values
      const Scalar *x = x_in.values();
      Scalar       *y = y_out->values();
      // Perform y = beta*y (being careful to set y=0 if beta=0 in case y is uninitialized on input!)
      Thyra::Index k = 0;
      if( beta == ST::zero() ) {
        for( k = 0; k < dim_; ++k ) y[k] = ST::zero();
      }
      else if( beta != ST::one() ) {
        for( k = 0; k < dim_; ++k ) y[k] *= beta;
      }
      // Perform y = alpha*op(M)*x 
      k = 0;
      if( M_trans == Thyra::NOTRANS ) {
        y[k] += alpha * ( diag_[k]*x[k] + upper_[k]*x[k+1] );                         // First row
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );  // Middle rows
        y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] );                       // Last row
      }
      else if( M_trans == Thyra::CONJ ) {
        y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
        y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
      }
      else if( M_trans == Thyra::TRANS ) {
        y[k] += alpha * ( diag_[k]*x[k] + lower_[k]*x[k+1] );
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] + lower_[k]*x[k+1] );
        y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] );
      }
      else if( M_trans == Thyra::CONJTRANS ) {
        y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
        y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
      }
      else {
        TEST_FOR_EXCEPT(true); // Throw exception if we get here!
      }
    }

};  // end class ExampleTridiagSerialLinearOp

The above serial matrix class is used in an example program (see runPowerMethodExample() below) that calls sillyPowerMethod(). In this example program, the matrix constructed and used is the well-known tridiagonal matrix

\[ A= \left[\begin{array}{rrrrrrrrrr} 2 & -1 \\ -1 & 2 & -1 \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{array}\right]. \]

The power method is then run on the matrix $A$ run for a number of iterations (or until convergence to some tolerance).

After this, the first diagonal element $A_{(1,1)}=2$ is then scaled to $A_{(1,1)}=20$ and the power method is run again (which much faster convergence).

The following templated function implements the example described above:

template<class Scalar>
bool runPowerMethodExample(
  const int    dim
  ,const int   maxNumIters
  ,const bool  verbose
  ,const bool  dumpAll
  )
{

  using Teuchos::OSTab;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType    ScalarMag;
  
  bool success = true;
  bool result;

  Teuchos::RefCountPtr<Teuchos::FancyOStream>
    out = (verbose ? Teuchos::VerboseObjectBase::getDefaultOStream() : Teuchos::null);

  if(verbose)
    *out << "\n***\n*** Running power method example using scalar type = \'" << ST::name() << "\' ...\n***\n" << std::scientific;

  //
  // (1) Setup the initial tridiagonal operator
  //
  //       [  2  -1             ]
  //       [ -1   2  -1         ]
  //  A =  [      .   .   .     ]
  //       [          -1  2  -1 ]
  //       [             -1   2 ]
  //
  if(verbose) *out << "\n(1) Constructing tridiagonal matrix A of dimension = " << dim << " ...\n";
  std::vector<Scalar> lower(dim-1), diag(dim), upper(dim-1);
  const Scalar one = ST::one(), two = Scalar(2)*one;
  int k = 0;
  diag[k] = two; upper[k] = -one;                        //  First row
  for( k = 1; k < dim - 1; ++k ) {
    lower[k-1] = -one; diag[k] = two; upper[k] = -one;   //  Middle rows
  }
  lower[k-1] = -one; diag[k] = two;                      //  Last row
  Teuchos::RefCountPtr<ExampleTridiagSerialLinearOp<Scalar> >
    A = Teuchos::rcp( new ExampleTridiagSerialLinearOp<Scalar>(dim,&lower[0],&diag[0],&upper[0]) );
  if( verbose && dumpAll ) *out << "\nA =\n" << *A;

  //
  // (2) Run the power method ANA
  //
  if(verbose) *out << "\n(2) Running the power method on matrix A ...\n";
  Scalar     lambda      = ST::zero();
  ScalarMag  tolerance   = ScalarMag(1e-3)*Teuchos::ScalarTraits<ScalarMag>::one();
  if(1){
    OSTab tab(out);
    result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,out.get());
    if(!result) success = false;
    if(verbose) *out << "\nEstimate of dominate eigenvalue lambda = " << lambda << std::endl;
  }

  //
  // (3) Increase dominance of first eigenvalue
  //
  if(verbose) *out << "\n(3) Increasing first diagonal entry by factor of 10 ...\n";
  diag[0] *= 10;
  A->initialize(dim,&lower[0],&diag[0],&upper[0]);
  if( verbose && dumpAll ) *out << "A =\n" << *A;

  //
  // (4) Run the power method ANA
  //
  if(verbose) *out << "\n(4) Running the power method again on matrix A ...\n";
  if(1){
    OSTab tab(out);
    result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,out.get());
    if(!result) success = false;
    if(verbose) *out << "\nEstimate of dominate eigenvalue lambda = " << lambda << std::endl;
  }
  
  return success;

} // end runPowerMethodExample()

The above templated function runPowerMethodExample() is instantiated with the following scalar types:

and is called multiple times from within the following main() program function:

int main(int argc, char *argv[])
{

  using Teuchos::CommandLineProcessor;
 
  bool success = true;
  bool verbose = true;
  bool result;

  Teuchos::RefCountPtr<Teuchos::FancyOStream>
    out = Teuchos::VerboseObjectBase::getDefaultOStream();

  try {

    //
    // Read in command-line options
    //
    
    int    dim          = 4;
    bool   dumpAll      = false;

    CommandLineProcessor  clp;
    clp.throwExceptions(false);
    clp.addOutputSetupOptions(true);
    clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
    clp.setOption( "dim", &dim, "Dimension of the linear system." );
    clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
    CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;

    TEST_FOR_EXCEPTION( dim < 2, std::logic_error, "Error, dim=" << dim << " < 2 is not allowed!" );

    int    maxNumIters  = 10*dim;
    
    // Run using float
    result = runPowerMethodExample<float>(dim,maxNumIters,verbose,dumpAll);
    if(!result) success = false;

    // Run using double
    result = runPowerMethodExample<double>(dim,maxNumIters,verbose,dumpAll);
    if(!result) success = false;

#if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
    
    // Run using std::complex<float>
    result = runPowerMethodExample<std::complex<float> >(dim,maxNumIters,verbose,dumpAll);
    if(!result) success = false;
    
    // Run using std::complex<double>
    result = runPowerMethodExample<std::complex<double> >(dim,maxNumIters,verbose,dumpAll);
    if(!result) success = false;

#endif

#ifdef HAVE_TEUCHOS_GNU_MP
    
    // Run using mpf_class
    result = runPowerMethodExample<mpf_class >(dim,maxNumIters,verbose,dumpAll);
    if(!result) success = false;

#if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
    
    // Run using std::complex<mpf_class>
    //result = runPowerMethodExample<std::complex<mpf_class> >(dim,maxNumIters,verbose,dumpAll);
    //if(!result) success = false;
    //The above commented-out code throws a floating-point exception?
 
#endif


#endif
    
  }
  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
  
  if (verbose) {
    if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
    else         *out << "\nOh no! At least one of the tests failed!\n";
  }
  
  return success ? 0 : 1;

} // end main()

The above example program is built as part of the Thyra package (unless examples where disabled at configure time) and the executable can be found at:

 ./example/operator_vector/sillyPowerMethod_serial.exe 

where ./ is the base build directory for Thyra (e.g. $TRILINOS_BUILD_DIR/packages/thyra).

This example program should run successfully with no arguments and produces the following output:

***
*** Running power method example using scalar type = 'float' ...
***

(1) Constructing tridiagonal matrix A of dimension = 4 ...

(2) Running the power method on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = 5.745785e-01, ||A*q-lambda*q|| = 4.269211e-01
  Iter = 4, lambda = 3.301265e+00, ||A*q-lambda*q|| = 5.727251e-01
  Iter = 8, lambda = 3.595284e+00, ||A*q-lambda*q|| = 1.493326e-01
  Iter = 12, lambda = 3.616292e+00, ||A*q-lambda*q|| = 4.170860e-02
  Iter = 16, lambda = 3.617903e+00, ||A*q-lambda*q|| = 1.145323e-02
  Iter = 20, lambda = 3.618024e+00, ||A*q-lambda*q|| = 3.140391e-03
  Iter = 24, lambda = 3.618033e+00, ||A*q-lambda*q|| = 8.610450e-04
  
  Estimate of dominate eigenvalue lambda = 3.618033e+00

(3) Increasing first diagonal entry by factor of 10 ...

(4) Running the power method again on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = 8.132596e+00, ||A*q-lambda*q|| = 9.248439e+00
  Iter = 4, lambda = 2.005555e+01, ||A*q-lambda*q|| = 2.941549e-03
  Iter = 8, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.368338e-06
  
  Estimate of dominate eigenvalue lambda = 2.005556e+01

***
*** Running power method example using scalar type = 'double' ...
***

(1) Constructing tridiagonal matrix A of dimension = 4 ...

(2) Running the power method on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = 5.745785e-01, ||A*q-lambda*q|| = 4.269210e-01
  Iter = 4, lambda = 3.301264e+00, ||A*q-lambda*q|| = 5.727251e-01
  Iter = 8, lambda = 3.595283e+00, ||A*q-lambda*q|| = 1.493328e-01
  Iter = 12, lambda = 3.616291e+00, ||A*q-lambda*q|| = 4.170859e-02
  Iter = 16, lambda = 3.617903e+00, ||A*q-lambda*q|| = 1.145327e-02
  Iter = 20, lambda = 3.618024e+00, ||A*q-lambda*q|| = 3.140456e-03
  Iter = 24, lambda = 3.618033e+00, ||A*q-lambda*q|| = 8.610081e-04
  
  Estimate of dominate eigenvalue lambda = 3.618033e+00

(3) Increasing first diagonal entry by factor of 10 ...

(4) Running the power method again on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = 8.132596e+00, ||A*q-lambda*q|| = 9.248439e+00
  Iter = 4, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.941630e-03
  Iter = 8, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.259420e-06
  
  Estimate of dominate eigenvalue lambda = 2.005556e+01

***
*** Running power method example using scalar type = 'std::complex<float>' ...
***

(1) Constructing tridiagonal matrix A of dimension = 4 ...

(2) Running the power method on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = (1.617793e+00,1.490116e-08), ||A*q-lambda*q|| = 1.333844e+00
  Iter = 4, lambda = (3.550094e+00,-5.820766e-09), ||A*q-lambda*q|| = 2.528435e-01
  Iter = 8, lambda = (3.612626e+00,7.217750e-09), ||A*q-lambda*q|| = 7.334734e-02
  Iter = 12, lambda = (3.617626e+00,6.402843e-10), ||A*q-lambda*q|| = 2.020937e-02
  Iter = 16, lambda = (3.618003e+00,-1.698209e-08), ||A*q-lambda*q|| = 5.542850e-03
  Iter = 20, lambda = (3.618032e+00,-8.076313e-09), ||A*q-lambda*q|| = 1.519599e-03
  Iter = 24, lambda = (3.618034e+00,-1.025001e-08), ||A*q-lambda*q|| = 4.167181e-04
  
  Estimate of dominate eigenvalue lambda = (3.618034e+00,-1.025001e-08)

(3) Increasing first diagonal entry by factor of 10 ...

(4) Running the power method again on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = (6.482740e+00,-1.490116e-08), ||A*q-lambda*q|| = 7.724599e+00
  Iter = 4, lambda = (2.005554e+01,-1.837986e-08), ||A*q-lambda*q|| = 1.717740e-02
  Iter = 8, lambda = (2.005555e+01,-3.240486e-07), ||A*q-lambda*q|| = 1.432346e-05
  
  Estimate of dominate eigenvalue lambda = (2.005555e+01,-3.240486e-07)

***
*** Running power method example using scalar type = 'std::complex<double>' ...
***

(1) Constructing tridiagonal matrix A of dimension = 4 ...

(2) Running the power method on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = (1.617793e+00,2.775558e-17), ||A*q-lambda*q|| = 1.333844e+00
  Iter = 4, lambda = (3.550094e+00,1.040834e-17), ||A*q-lambda*q|| = 2.528435e-01
  Iter = 8, lambda = (3.612625e+00,1.040834e-17), ||A*q-lambda*q|| = 7.334728e-02
  Iter = 12, lambda = (3.617625e+00,-1.116728e-17), ||A*q-lambda*q|| = 2.020930e-02
  Iter = 16, lambda = (3.618003e+00,1.081492e-17), ||A*q-lambda*q|| = 5.542757e-03
  Iter = 20, lambda = (3.618032e+00,-3.577867e-18), ||A*q-lambda*q|| = 1.519668e-03
  Iter = 24, lambda = (3.618034e+00,9.903509e-18), ||A*q-lambda*q|| = 4.166393e-04
  
  Estimate of dominate eigenvalue lambda = (3.618034e+00,9.903509e-18)

(3) Increasing first diagonal entry by factor of 10 ...

(4) Running the power method again on matrix A ...
  
  Starting power method (target tolerrance = 1.000000e-03) ...
  
  Iter = 0, lambda = (6.482739e+00,1.942890e-16), ||A*q-lambda*q|| = 7.724598e+00
  Iter = 4, lambda = (2.005554e+01,6.367994e-17), ||A*q-lambda*q|| = 1.717743e-02
  Iter = 8, lambda = (2.005556e+01,3.011987e-16), ||A*q-lambda*q|| = 1.414814e-05
  
  Estimate of dominate eigenvalue lambda = (2.005556e+01,3.011987e-16)

Congratulations! All of the tests checked out!

This example program also takes a number of command-line options. To see what the command-line options are, use the --help option. The command-line options returned from ./sillyPowerMethod_serial.exe --echo-command-line --help are:

Echoing the command-line:

../example/operator_vector/sillyPowerMethod_serial.exe --echo-command-line --help 

Usage: ../example/operator_vector/sillyPowerMethod_serial.exe [options]
  options:
  --help                                Prints this help message
  --pause-for-debugging                 Pauses for user input to allow attaching a debugger
  --echo-command-line                   Echo the command-line but continue as normal
  --output-all-front-matter     bool    Set if all front matter is printed to the default FancyOStream or not
  --output-no-front-matter              (default: --output-no-front-matter)
  --output-show-line-prefix     bool    Set if the line prefix matter is printed to the default FancyOStream or not
  --output-no-show-line-prefix          (default: --output-no-show-line-prefix)
  --output-show-tab-count       bool    Set if the tab count is printed to the default FancyOStream or not
  --output-no-show-tab-count            (default: --output-no-show-tab-count)
  --output-show-proc-rank       bool    Set if the processor rank is printed to the default FancyOStream or not
  --output-no-show-proc-rank            (default: --output-no-show-proc-rank)
  --output-to-root-rank-only    int     Set which processor (the root) gets the output.  If < 0, then all processors get output.
                                        (default: --output-to-root-rank-only=0)
  --verbose                     bool    Determines if any output is printed or not.
  --quiet                               (default: --verbose)
  --dim                         int     Dimension of the linear system.
                                        (default: --dim=4)
  --dump-all                    bool    Determines if quantities are dumped or not.
  --no-dump                             (default: --no-dump)

To see the full listing of this example program click: sillyPowerMethod_serial.cpp


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