simple_stratimikos_example.cpp

Simple example program for how the Thyra::DefaultRealLinearSolverBuilder class is used to take a parameter list (e.g. read from a file) and use it so solve a linear system read in as Epetra objects from a file.

Outline

Example program source (with line numbers)

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //         Stratimikos: Thyra-based strategies for linear solvers
00005 //                Copyright (2006) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Thyra_DefaultRealLinearSolverBuilder.hpp"
00030 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00031 #include "Thyra_EpetraThyraWrappers.hpp"
00032 #include "Thyra_EpetraLinearOp.hpp"
00033 #include "EpetraExt_readEpetraLinearSystem.h"
00034 #include "Teuchos_GlobalMPISession.hpp"
00035 #include "Teuchos_VerboseObject.hpp"
00036 #include "Teuchos_XMLParameterListHelpers.hpp"
00037 #include "Teuchos_CommandLineProcessor.hpp"
00038 #include "Teuchos_StandardCatchMacros.hpp"
00039 #ifdef HAVE_MPI
00040 #  include "Epetra_MpiComm.h"
00041 #else
00042 #  include "Epetra_SerialComm.h"
00043 #endif
00044 
00045 namespace {
00046 
00047 // Helper function to compute a single norm for a vector
00048 double epetraNorm2( const Epetra_Vector &v )
00049 {
00050   double norm[1] = { -1.0 };
00051   v.Norm2(&norm[0]);
00052   return norm[0];
00053 }
00054 
00055 } // namespace
00056 
00057 int main(int argc, char* argv[])
00058 {
00059 
00060   using Teuchos::rcp;
00061   using Teuchos::RefCountPtr;
00062   using Teuchos::CommandLineProcessor;
00063 
00064   bool success = true;
00065   bool verbose = true;
00066 
00067   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00068 
00069   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00070     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00071 
00072   try {
00073     
00074     //
00075     // A) Program setup code
00076     //
00077 
00078     //
00079     // Read options from command-line
00080     //
00081     
00082     std::string     matrixFile             = "";
00083     double          tol                    = 1e-5;
00084     bool            onlyPrintOptions       = false;
00085     bool            printXmlFormat         = false;
00086 
00087     Thyra::DefaultRealLinearSolverBuilder linearSolverBuilder;
00088 
00089     CommandLineProcessor  clp(false); // Don't throw exceptions
00090 
00091     // Set up command-line options for the linear solver that will be used!
00092     linearSolverBuilder.setupCLP(&clp);
00093 
00094     clp.setOption( "matrix-file", &matrixFile
00095                    ,"Defines the matrix and perhaps the RHS and LHS for a linear system to be solved." );
00096     clp.setOption( "tol", &tol, "Tolerance to check against the scaled residual norm." );
00097     clp.setOption( "only-print-options", "continue-after-printing-options", &onlyPrintOptions
00098                    ,"Only print options and stop or continue on" );
00099     clp.setOption( "print-xml-format", "print-readable-format", &printXmlFormat
00100                    ,"Print the valid options in XML format or in readable format." );
00101     
00102     clp.setDocString(
00103       "Simple example for the use of the Stratimikos facade Thyra::DefaultRealLinearSolverBuilder.\n"
00104       "\n"
00105       "To print out just the valid options use --matrix-file=\"\" --only-print-options with --print-xml-format"
00106       " or --print-readable-format.\n"
00107       "\n"
00108       "To solve a linear system from a system read from a file use --matrix-file=\"SomeFile.mtx\""
00109       " with options read from an XML file using --linear-solver-params-file=\"SomeFile.xml\"\n"
00110       );
00111 
00112     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00113     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00114 
00115     //
00116     // Print out the valid options and stop if asked
00117     //
00118 
00119     if(onlyPrintOptions) {
00120       if(printXmlFormat)
00121         Teuchos::writeParameterListToXmlOStream(
00122           *linearSolverBuilder.getValidParameters()
00123           ,*out
00124           );
00125       else
00126         linearSolverBuilder.getValidParameters()->print(*out,1,true,false);
00127       return 0;
00128     }
00129     
00130     //
00131     // B) Epetra-specific code that sets up the linear system to be solved
00132     //
00133     // While the below code reads in the Epetra objects from a file, you can
00134     // setup the Epetra objects any way you would like.  Note that this next
00135     // set of code as nothing to do with Thyra at all, and it should not.
00136     //
00137 
00138     *out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00139     
00140 #ifdef HAVE_MPI
00141     Epetra_MpiComm comm(MPI_COMM_WORLD);
00142 #else
00143     Epetra_SerialComm comm;
00144 #endif
00145     RefCountPtr<Epetra_CrsMatrix> epetra_A;
00146     RefCountPtr<Epetra_Vector> epetra_x, epetra_b;
00147     EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, &epetra_x, &epetra_b );
00148 
00149     if(!epetra_b.get()) {
00150       *out << "\nThe RHS b was not read in so generate a new random vector ...\n";
00151       epetra_b = rcp(new Epetra_Vector(epetra_A->OperatorRangeMap()));
00152       epetra_b->Random();
00153     }
00154 
00155     if(!epetra_x.get()) {
00156       *out << "\nThe LHS x was not read in so generate a new zero vector ...\n";
00157       epetra_x = rcp(new Epetra_Vector(epetra_A->OperatorDomainMap()));
00158       epetra_x->PutScalar(0.0); // Initial guess is critical!
00159     }
00160 
00161     *out << "\nPrinting statistics of the Epetra linear system ...\n";
00162 
00163     *out
00164       << "\n  Epetra_CrsMatrix epetra_A of dimension "
00165       << epetra_A->OperatorRangeMap().NumGlobalElements()
00166       << " x " << epetra_A->OperatorDomainMap().NumGlobalElements()
00167       << "\n  ||epetraA||inf = " << epetra_A->NormInf()
00168       << "\n  ||epetra_b||2 = " << epetraNorm2(*epetra_b)
00169       << "\n  ||epetra_x||2 = " << epetraNorm2(*epetra_x)
00170       << "\n";
00171 
00172     //
00173     // C) The "Glue" code that takes Epetra objects and wraps them as Thyra
00174     // objects
00175     //
00176     // This next set of code wraps the Epetra objects that define the linear
00177     // system to be solved as Thyra objects so that they can be passed to the
00178     // linear solver.
00179     //
00180 
00181     // Create RCPs that will be used to hold the Thyra wrappers
00182     RefCountPtr<const Thyra::LinearOpBase<double> >    A;
00183     RefCountPtr<Thyra::VectorBase<double> >            x;
00184     RefCountPtr<const Thyra::VectorBase<double> >      b;
00185 
00186     // Create the Thyra wrappers
00187     if(1) {
00188       // Create an RCP directly to the EpetraLinearOp so that we can access the
00189       // right range and domains spaces to use to create the wrappers for the
00190       // vector objects.
00191       RefCountPtr<const Thyra::EpetraLinearOp>
00192         _A = rcp(new Thyra::EpetraLinearOp(epetra_A));
00193       // Create Thyra wrappers for the vector objects that will automatically
00194       // update the Epetra objects.
00195       b = Thyra::create_Vector(epetra_b,_A->spmdRange());
00196       x = Thyra::create_Vector(epetra_x,_A->spmdDomain());
00197       // Set the RCP to the base linear operator interface
00198       A = _A;
00199     }
00200 
00201     // Note that above Thyra is only interacted with in the most trival of
00202     // ways.  For most users, Thyra will only be seen as a thin wrapper that
00203     // they need know little about in order to wrap their objects in order to
00204     // pass them to Thyra-enabled solvers.
00205       
00206     //
00207     // D) Thyra-specific code for solving the linear system
00208     //
00209     // Note that this code has no mention of any concrete implementation and
00210     // therefore can be used in any use case.
00211     //
00212 
00213     // Reading in the solver parameters from the parameters file and/or from
00214     // the command line.  This was setup by the command-line options
00215     // set by the setupCLP(...) function above.
00216     linearSolverBuilder.readParameters(out.get());
00217     // Create a linear solver factory given information read from the
00218     // parameter list.
00219     RefCountPtr<Thyra::LinearOpWithSolveFactoryBase<double> >
00220       lowsFactory = linearSolverBuilder.createLinearSolveStrategy("");
00221     // Setup output stream and the verbosity level
00222     lowsFactory->setOStream(out);
00223     lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
00224     // Create a linear solver based on the forward operator A
00225     RefCountPtr<Thyra::LinearOpWithSolveBase<double> >
00226       lows = Thyra::linearOpWithSolve(*lowsFactory,A);
00227     // Solve the linear system (note: the initial guess in 'x' is critical)
00228     Thyra::SolveStatus<double>
00229       status = Thyra::solve(*lows,Thyra::NOTRANS,*b,&*x);
00230     *out << "\nSolve status:\n" << status;
00231     // Write the linear solver parameters after they were read
00232     linearSolverBuilder.writeParamsFile(*lowsFactory);
00233 
00234     //
00235     // E) Post process the solution and check the error
00236     //
00237     // Note that the below code is based only on the Epetra objects themselves
00238     // and does not in any way depend or interact with any Thyra-based
00239     // objects.  The point is that most users of Thyra can largely gloss over
00240     // the fact that Thyra is really being used for anything.
00241     //
00242 
00243     // Wipe out the Thyra wrapper for x to guarantee that the solution will be
00244     // written back to epetra_x!
00245     x = Teuchos::null;
00246 
00247     *out
00248       << "\nSolution ||epetra_x||2 = " << epetraNorm2(*epetra_x) << "\n";
00249 
00250     *out << "\nTesting the solution error ||b-A*x||/||b|| computed through the Epetra objects ...\n";
00251 
00252     // r = b - A*x
00253     Epetra_Vector epetra_r(*epetra_b);
00254     if(1) {
00255       Epetra_Vector epetra_A_x(epetra_A->OperatorRangeMap());
00256       epetra_A->Apply(*epetra_x,epetra_A_x);
00257       epetra_r.Update(-1.0,epetra_A_x,1.0);
00258     }
00259       
00260     const double
00261       nrm_r = epetraNorm2(epetra_r),
00262       nrm_b = epetraNorm2(*epetra_b),
00263       rel_err = ( nrm_r / nrm_b );
00264     const bool
00265       passed = (rel_err <= tol);
00266 
00267     *out
00268       << "||b-A*x||/||b|| = " << nrm_r << "/" << nrm_b << " = " << rel_err
00269       << " < tol = " << tol << " ? " << ( passed ? "passed" : "failed" ) << "\n";
00270 
00271     if(!passed) success = false;
00272     
00273   }
00274   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success)
00275   
00276   if (verbose) {
00277     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00278     else         *out << "\nOh no! At least one of the tests failed!\n";
00279   }
00280 
00281   return ( success ? 0 : 1 );
00282 }

Command-line options for the driver program

Echoing the command-line:

../example/simple_stratimikos_example.exe --echo-command-line --help 

Usage: ../example/simple_stratimikos_example.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
  --linear-solver-params-file        string  Name of an XML file containing parameters for linear solver options to be appended first.
                                             (default: --linear-solver-params-file="")
  --extra-linear-solver-params       string  An XML string containing linear solver parameters to be appended second.
                                             (default: --extra-linear-solver-params="")
  --linear-solver-params-used-file   string  Name of an XML file that can be written with the parameter list after it has been used on completion of this program.
                                             (default: --linear-solver-params-used-file="")
  --matrix-file                      string  Defines the matrix and perhaps the RHS and LHS for a linear system to be solved.
                                             (default: --matrix-file="")
  --tol                              double  Tolerance to check against the scaled residual norm.
                                             (default: --tol=1e-05)
  --only-print-options               bool    Only print options and stop or continue on
  --continue-after-printing-options          (default: --continue-after-printing-options)
  --print-xml-format                 bool    Print the valid options in XML format or in readable format.
  --print-readable-format                    (default: --print-readable-format)

DETAILED DOCUMENTATION:

Simple example for the use of the Stratimikos facade Thyra::DefaultRealLinearSolverBuilder.

To print out just the valid options use --matrix-file="" --only-print-options with --print-xml-format or --print-readable-format.

To solve a linear system from a system read from a file use --matrix-file="SomeFile.mtx" with options read from an XML file using --linear-solver-params-file="SomeFile.xml"


Teuchos::GlobalMPISession::GlobalMPISession(): started serial run

Example input XML file

Here is an example input XML file, called amesos.klu.xml, that specifies that the KLU solver from Amesos be used to solve the linear system:

<ParameterList>
  <Parameter name="Linear Solver Type" type="string" value="Amesos"/>
  <ParameterList name="Linear Solver Types">
    <ParameterList name="Amesos">
      <ParameterList name="Amesos Settings"/>
      <Parameter name="Solver Type" type="string" value="Klu"/>
    </ParameterList>
  </ParameterList>
  <Parameter name="Preconditioner Type" type="string" value="None"/>
</ParameterList>

Example program output

Sample output looks like when using the above amesos.klu.xml file as input with a test system:

Echoing the command-line:

../example/simple_stratimikos_example.exe --echo-command-line --matrix-file=../example/FourByFour.mtx --linear-solver-params-file=../example/amesos.klu.xml 

Teuchos::GlobalMPISession::GlobalMPISession(): started serial run

Reading linear system in Epetra format from the file '../example/FourByFour.mtx' ...

Printing statistics of the Epetra linear system ...

  Epetra_CrsMatrix epetra_A of dimension 4 x 4
  ||epetraA||inf = 1.643
  ||epetra_b||2 = 1.52593
  ||epetra_x||2 = 0

Reading parameters from XML file "../example/amesos.klu.xml" ...
  
  Solving block system using Amesos solver 10Amesos_Klu ...
  
  
  Total solve time = 0 sec

Solve status:
  solveStatus = SOLVE_STATUS_CONVERGED
  achievedTol = unknownTolerance()
  message: "Solver 10Amesos_Klu converged!"
  extraParameters: NONE

Solution ||epetra_x||2 = 1.24984

Testing the solution error ||b-A*x||/||b|| computed through the Epetra objects ...
||b-A*x||/||b|| = 0/1.52593 = 0 < tol = 1e-05 ? passed

Congratulations! All of the tests checked out!

Example program source (without line numbers)

Here is the main driver program itself without line numbers:

// @HEADER
// ***********************************************************************
// 
//         Stratimikos: Thyra-based strategies for linear solvers
//                Copyright (2006) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#include "Thyra_DefaultRealLinearSolverBuilder.hpp"
#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
#include "Thyra_EpetraThyraWrappers.hpp"
#include "Thyra_EpetraLinearOp.hpp"
#include "EpetraExt_readEpetraLinearSystem.h"
#include "Teuchos_GlobalMPISession.hpp"
#include "Teuchos_VerboseObject.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_StandardCatchMacros.hpp"
#ifdef HAVE_MPI
#  include "Epetra_MpiComm.h"
#else
#  include "Epetra_SerialComm.h"
#endif

namespace {

// Helper function to compute a single norm for a vector
double epetraNorm2( const Epetra_Vector &v )
{
  double norm[1] = { -1.0 };
  v.Norm2(&norm[0]);
  return norm[0];
}

} // namespace

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

  using Teuchos::rcp;
  using Teuchos::RefCountPtr;
  using Teuchos::CommandLineProcessor;

  bool success = true;
  bool verbose = true;

  Teuchos::GlobalMPISession mpiSession(&argc,&argv);

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

  try {
    
    //
    // A) Program setup code
    //

    //
    // Read options from command-line
    //
    
    std::string     matrixFile             = "";
    double          tol                    = 1e-5;
    bool            onlyPrintOptions       = false;
    bool            printXmlFormat         = false;

    Thyra::DefaultRealLinearSolverBuilder linearSolverBuilder;

    CommandLineProcessor  clp(false); // Don't throw exceptions

    // Set up command-line options for the linear solver that will be used!
    linearSolverBuilder.setupCLP(&clp);

    clp.setOption( "matrix-file", &matrixFile
                   ,"Defines the matrix and perhaps the RHS and LHS for a linear system to be solved." );
    clp.setOption( "tol", &tol, "Tolerance to check against the scaled residual norm." );
    clp.setOption( "only-print-options", "continue-after-printing-options", &onlyPrintOptions
                   ,"Only print options and stop or continue on" );
    clp.setOption( "print-xml-format", "print-readable-format", &printXmlFormat
                   ,"Print the valid options in XML format or in readable format." );
    
    clp.setDocString(
      "Simple example for the use of the Stratimikos facade Thyra::DefaultRealLinearSolverBuilder.\n"
      "\n"
      "To print out just the valid options use --matrix-file=\"\" --only-print-options with --print-xml-format"
      " or --print-readable-format.\n"
      "\n"
      "To solve a linear system from a system read from a file use --matrix-file=\"SomeFile.mtx\""
      " with options read from an XML file using --linear-solver-params-file=\"SomeFile.xml\"\n"
      );

    CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;

    //
    // Print out the valid options and stop if asked
    //

    if(onlyPrintOptions) {
      if(printXmlFormat)
        Teuchos::writeParameterListToXmlOStream(
          *linearSolverBuilder.getValidParameters()
          ,*out
          );
      else
        linearSolverBuilder.getValidParameters()->print(*out,1,true,false);
      return 0;
    }
    
    //
    // B) Epetra-specific code that sets up the linear system to be solved
    //
    // While the below code reads in the Epetra objects from a file, you can
    // setup the Epetra objects any way you would like.  Note that this next
    // set of code as nothing to do with Thyra at all, and it should not.
    //

    *out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
    
#ifdef HAVE_MPI
    Epetra_MpiComm comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm comm;
#endif
    RefCountPtr<Epetra_CrsMatrix> epetra_A;
    RefCountPtr<Epetra_Vector> epetra_x, epetra_b;
    EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, &epetra_x, &epetra_b );

    if(!epetra_b.get()) {
      *out << "\nThe RHS b was not read in so generate a new random vector ...\n";
      epetra_b = rcp(new Epetra_Vector(epetra_A->OperatorRangeMap()));
      epetra_b->Random();
    }

    if(!epetra_x.get()) {
      *out << "\nThe LHS x was not read in so generate a new zero vector ...\n";
      epetra_x = rcp(new Epetra_Vector(epetra_A->OperatorDomainMap()));
      epetra_x->PutScalar(0.0); // Initial guess is critical!
    }

    *out << "\nPrinting statistics of the Epetra linear system ...\n";

    *out
      << "\n  Epetra_CrsMatrix epetra_A of dimension "
      << epetra_A->OperatorRangeMap().NumGlobalElements()
      << " x " << epetra_A->OperatorDomainMap().NumGlobalElements()
      << "\n  ||epetraA||inf = " << epetra_A->NormInf()
      << "\n  ||epetra_b||2 = " << epetraNorm2(*epetra_b)
      << "\n  ||epetra_x||2 = " << epetraNorm2(*epetra_x)
      << "\n";

    //
    // C) The "Glue" code that takes Epetra objects and wraps them as Thyra
    // objects
    //
    // This next set of code wraps the Epetra objects that define the linear
    // system to be solved as Thyra objects so that they can be passed to the
    // linear solver.
    //

    // Create RCPs that will be used to hold the Thyra wrappers
    RefCountPtr<const Thyra::LinearOpBase<double> >    A;
    RefCountPtr<Thyra::VectorBase<double> >            x;
    RefCountPtr<const Thyra::VectorBase<double> >      b;

    // Create the Thyra wrappers
    if(1) {
      // Create an RCP directly to the EpetraLinearOp so that we can access the
      // right range and domains spaces to use to create the wrappers for the
      // vector objects.
      RefCountPtr<const Thyra::EpetraLinearOp>
        _A = rcp(new Thyra::EpetraLinearOp(epetra_A));
      // Create Thyra wrappers for the vector objects that will automatically
      // update the Epetra objects.
      b = Thyra::create_Vector(epetra_b,_A->spmdRange());
      x = Thyra::create_Vector(epetra_x,_A->spmdDomain());
      // Set the RCP to the base linear operator interface
      A = _A;
    }

    // Note that above Thyra is only interacted with in the most trival of
    // ways.  For most users, Thyra will only be seen as a thin wrapper that
    // they need know little about in order to wrap their objects in order to
    // pass them to Thyra-enabled solvers.
      
    //
    // D) Thyra-specific code for solving the linear system
    //
    // Note that this code has no mention of any concrete implementation and
    // therefore can be used in any use case.
    //

    // Reading in the solver parameters from the parameters file and/or from
    // the command line.  This was setup by the command-line options
    // set by the setupCLP(...) function above.
    linearSolverBuilder.readParameters(out.get());
    // Create a linear solver factory given information read from the
    // parameter list.
    RefCountPtr<Thyra::LinearOpWithSolveFactoryBase<double> >
      lowsFactory = linearSolverBuilder.createLinearSolveStrategy("");
    // Setup output stream and the verbosity level
    lowsFactory->setOStream(out);
    lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
    // Create a linear solver based on the forward operator A
    RefCountPtr<Thyra::LinearOpWithSolveBase<double> >
      lows = Thyra::linearOpWithSolve(*lowsFactory,A);
    // Solve the linear system (note: the initial guess in 'x' is critical)
    Thyra::SolveStatus<double>
      status = Thyra::solve(*lows,Thyra::NOTRANS,*b,&*x);
    *out << "\nSolve status:\n" << status;
    // Write the linear solver parameters after they were read
    linearSolverBuilder.writeParamsFile(*lowsFactory);

    //
    // E) Post process the solution and check the error
    //
    // Note that the below code is based only on the Epetra objects themselves
    // and does not in any way depend or interact with any Thyra-based
    // objects.  The point is that most users of Thyra can largely gloss over
    // the fact that Thyra is really being used for anything.
    //

    // Wipe out the Thyra wrapper for x to guarantee that the solution will be
    // written back to epetra_x!
    x = Teuchos::null;

    *out
      << "\nSolution ||epetra_x||2 = " << epetraNorm2(*epetra_x) << "\n";

    *out << "\nTesting the solution error ||b-A*x||/||b|| computed through the Epetra objects ...\n";

    // r = b - A*x
    Epetra_Vector epetra_r(*epetra_b);
    if(1) {
      Epetra_Vector epetra_A_x(epetra_A->OperatorRangeMap());
      epetra_A->Apply(*epetra_x,epetra_A_x);
      epetra_r.Update(-1.0,epetra_A_x,1.0);
    }
      
    const double
      nrm_r = epetraNorm2(epetra_r),
      nrm_b = epetraNorm2(*epetra_b),
      rel_err = ( nrm_r / nrm_b );
    const bool
      passed = (rel_err <= tol);

    *out
      << "||b-A*x||/||b|| = " << nrm_r << "/" << nrm_b << " = " << rel_err
      << " < tol = " << tol << " ? " << ( passed ? "passed" : "failed" ) << "\n";

    if(!passed) success = false;
    
  }
  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,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 );
}

Generated on Thu Sep 18 12:37:39 2008 for Stratimikos by doxygen 1.3.9.1