Epetra Implemenation of the CG 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 sillyCgSolve(). More...
This example program is contained in the source file

 ./thyra/example/sillyCgSolve_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 a & -1 \\ -1 & 2 a & -1 \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 a & -1 \\ & & & -1 & 2 a \end{array}\right] \]

where $a$ is an adjustable diagonal scale factor that makes the matrix $A$ more or less well conditioned.

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

The above matrix-building function Thyra::createTridiagEpetraLinearOp() is then called along with the sillyCgSolve() function 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;
    double diagScale                  = 1.001;
    bool   useWithNonEpetraVectors    = false;
    double tolerance                  = 1e-4;
    int    maxNumIters                = 300;

    CommandLineProcessor  clp(false); // Don't throw exceptions
    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( "diag-scale", &diagScale, "Scaling of the diagonal to improve conditioning." );
    clp.setOption( "use-with-non-epetra-vectors", "use-with-epetra-vectors", &useWithNonEpetraVectors
                   , "Use non-epetra vectors with Epetra operator or not." );
    clp.setOption( "tol", &tolerance, "Relative tolerance for linear system solve." );
    clp.setOption( "max-num-iters", &maxNumIters, "Maximum of CG iterations." );

    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 CG example using Epetra implementation\n***\n" << std::scientific;

    //
    // (D) Setup a simple linear system with tridagonal operator:
    //
    //       [  a*2   -1                ]
    //       [ -1    a*2  -1            ]
    //  A =  [         .   .    .       ]
    //       [            -1  a*2    -1 ]
    //       [                 -1   a*2 ]
    //
    // (D.1) Create the tridagonal Epetra matrix operator
    if(verbose) *out << "\n(1) Constructing tridagonal Epetra matrix A_epetra of global dimension = " << globalDim << " ...\n";
    RefCountPtr<Epetra_Operator>
      A_epetra = createTridiagEpetraLinearOp(
        globalDim
#ifdef HAVE_MPI
        ,mpiComm
#endif
        ,diagScale,verbose,*out
        );
    // (D.2) Wrap the Epetra_Opertor in a Thyra::EpetraLinearOp object
    RefCountPtr<Thyra::LinearOpBase<double> >
      A = rcp(new Thyra::EpetraLinearOp(A_epetra));
    // (D.3) Create RHS vector b and set to a random value
    RefCountPtr<const Thyra::VectorSpaceBase<double> >
      b_space = A->range();
    // (D.4) Create the RHS vector b and initialize it to a random vector
    RefCountPtr<Thyra::VectorBase<double> > b = createMember(b_space);
    Thyra::seed_randomize<double>(0);
    Thyra::randomize( -1.0, +1.0, &*b );
    // (D.5) Create LHS vector x and set to zero
    RefCountPtr<Thyra::VectorBase<double> > x = createMember(A->domain());
    Thyra::assign( &*x, 0.0 );
    //
    // (E) Solve the linear system with the silly CG solver
    //
    result = sillyCgSolve(*A,*b,maxNumIters,tolerance,&*x,verbose?&*out:0);
    if(!result) success = false;
    //
    // (F) Check that the linear system was solved to the specified tolerance
    //
    RefCountPtr<Thyra::VectorBase<double> > r = createMember(A->range());                     
    Thyra::assign(&*r,*b);                                        // r = b
    Thyra::apply(*A,Thyra::NOTRANS,*x,&*r,-1.0,+1.0);             // r = -A*x + r
    const double r_nrm = Thyra::norm(*r), b_nrm =Thyra:: norm(*b);
    const double rel_err = r_nrm/b_nrm, relaxTol = 2.0*tolerance;
    result = rel_err <= relaxTol;
    if(!result) success = false;
    if(verbose)
      *out
        << "\n||b-A*x||/||b|| = "<<r_nrm<<"/"<<b_nrm<<" = "<<rel_err<<(result?" <= ":" > ")
        <<"2.0*tolerance = "<<relaxTol<<": "<<(result?"passed":"failed")<<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.

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

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


Generated on Thu Sep 18 12:37:49 2008 for Epetra to Thyra Adapter Software by doxygen 1.3.9.1