Epetra Implemenation of the Power Method
[Assorted Epetra to Thyra Operator/Vector Adapter Example Code]

Here is an example program that shows the use of the Epetra adapter subclasses with the example linear ANA implementation sillyPowerMethod(). More...
This example program is contained in the source file

 ./thyra/example/sillyPowerMethod_epetra.cpp 

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

This example program is broken down into a few pieces in an attempt to be more understandable.

The first piece of example code we show is the function Thyra::createTridiagEpetraLinearOp(). This function creates an Epetra_CrsMatrix object that represents 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]. \]

After the Epetra_CrsMatrix object is created, it is used to create a Thyra::EpetraLinearOp object that adapts it into a Thyra::LinearOp object.

The implementation of the function Thyra::createTridiagEpetraLinearOp() is shown below:

The above matrix-building function Thyra::createTridiagEpetraLinearOp() is then called in the following MPI-enabled main() driver program:

int main(int argc, char *argv[])
{
  using Teuchos::CommandLineProcessor;
  using Teuchos::RefCountPtr;
 
  bool success = true;
  bool verbose = true;
  bool result;
  
  //
  // (A) Setup and get basic MPI info
  //
  
  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
  const int procRank = Teuchos::GlobalMPISession::getRank();
  const int numProc  = Teuchos::GlobalMPISession::getNProc();
#ifdef HAVE_MPI
  MPI_Comm mpiComm = MPI_COMM_WORLD;
#endif
  
  //
  // (B) Setup the output stream (do output only on root process!)
  //

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

    //
    // (C) Read in commandline options
    //

    int    globalDim                  = 500;
    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( "global-dim", &globalDim, "Global 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( globalDim < 2, std::logic_error, "Error, globalDim=" << globalDim << " < 2 is not allowed!" );

    if(verbose) *out << "\n***\n*** Running power method example using Epetra implementation\n***\n" << std::scientific;

    //
    // (D) Setup the operator and run the power method!
    //

    //
    // (1) Setup the initial tridiagonal operator
    //
    //       [  2  -1             ]
    //       [ -1   2  -1         ]
    //  A =  [      .   .   .     ]
    //       [          -1  2  -1 ]
    //       [             -1   2 ]
    //
    if(verbose) *out << "\n(1) Constructing tridagonal Epetra matrix A of global dimension = " << globalDim << " ...\n";
    RefCountPtr<Epetra_Operator>
      A_epetra = createTridiagEpetraLinearOp(
        globalDim
#ifdef HAVE_MPI
        ,mpiComm
#endif
        ,1.0,verbose,*out
        );
    // Wrap in an Thyra::EpetraLinearOp object
    RefCountPtr<Thyra::LinearOpBase<double> >
      A = rcp(new Thyra::EpetraLinearOp(A_epetra));
    //
    if( verbose && dumpAll ) *out << "\nA =\n" << *A; // This works even in parallel!
    
    //
    // (2) Run the power method ANA
    //
    if(verbose) *out << "\n(2) Running the power method on matrix A ...\n";
    double  lambda      = 0.0;
    double  tolerance   = 1e-3;
    int     maxNumIters = 10*globalDim;
    result = sillyPowerMethod<double>(*A,maxNumIters,tolerance,&lambda,(verbose?&*out:NULL));
    if(!result) success = false;
    if(verbose) *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
    
    //
    // (3) Increase dominance of first eigenvalue
    //
    if(verbose) *out << "\n(3) Scale the diagonal of A by a factor of 10 ...\n";
    scaleFirstDiagElement( 10.0, &*A );
    
    //
    // (4) Run the power method ANA again
    //
    if(verbose) *out << "\n(4) Running the power method again on matrix A ...\n";
    result = sillyPowerMethod<double>(*A,maxNumIters,tolerance,&lambda,(verbose?&*out:NULL));
    if(!result) success = false;
    if(verbose) *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
    
  }
  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 driver program should be very strightforward to follow and generates the exact same numerical results as in the double case in the tempalted serial version.

Notice however that the above driver program calls a helper function called scaleFirstDiagElement() which hides the details of scaling the first diagonal entry of the tridiagonal matrix operator. This helper function is not too tricking and takes the form:

void scaleFirstDiagElement( const double diagScale, Thyra::LinearOpBase<double> *A )
{
  using Teuchos::RefCountPtr;
  TEST_FOR_EXCEPT(A==NULL);
  // (A) Get at the underlying Epetra_Operator object that the EpetraLinearOp
  // object directly maintains.
  const RefCountPtr<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*A);
  // (B) Perform a dynamic cast to Epetra_CrsMatrix.
  // Note, the dyn_cast<>() template function will throw std::bad_cast
  // with a nice error message if the cast fails! 
  Epetra_CrsMatrix &crsMatrix = Teuchos::dyn_cast<Epetra_CrsMatrix>(*epetra_op);
  // (C) Query if the first row of the matrix is on this processor and if
  // it is get it and scale it.
  if(crsMatrix.MyGlobalRow(0)) {
    const int numRowNz = crsMatrix.NumGlobalEntries(0);
    TEST_FOR_EXCEPT( numRowNz != 2 );
    int returndNumRowNz; double rowValues[2]; int rowIndexes[2];
    crsMatrix.ExtractGlobalRowCopy( 0, numRowNz, returndNumRowNz, &rowValues[0], &rowIndexes[0] );
    TEST_FOR_EXCEPT( returndNumRowNz != 2 );
    for( int k = 0; k < numRowNz; ++k) if (rowIndexes[k] == 0) rowValues[k] *= diagScale;
    crsMatrix.ReplaceGlobalValues( 0, numRowNz, rowValues, rowIndexes );
  }
} // end scaleFirstDiagElement()

One important thing to notice about the above helper function is how the underlying Epetra_CrsMatrix object is accessed given the Thyra::LinearOp interface to the Thyra::EpetraLinearOp object that was created by Thyra::createTridiagEpetraLinearOp(). Actually, most of the C++ magic is performed in the helper function Thyra::get_Epetra_Operator() which is then followed by a simple dynamic cast.

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

 ./thyra/example/sillyPowerMethod_epetra.exe 

where ./ is the base build directory for epetra (e.g. ???/Trilinos/$BUILD_DIR/packages/epetra).

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

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_epetra.exe --echo-command-line --help are:

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


Generated on Sun May 20 13:06:11 2007 for Epetra to Thyra Adapter Software by doxygen 1.3.9.1