simple_stratimikos_example.cpp

Simple example program for how the Stratimikos::DefaultLinearSolverBuilder 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 "Stratimikos_DefaultLinearSolverBuilder.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::RCP;
00062   using Teuchos::CommandLineProcessor;
00063   typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
00064 
00065   bool success = true;
00066   bool verbose = true;
00067 
00068   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00069 
00070   Teuchos::RCP<Teuchos::FancyOStream>
00071     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00072 
00073   try {
00074 
00075     
00076     //
00077     // A) Program setup code
00078     //
00079 
00080     //
00081     // Read options from command-line
00082     //
00083     
00084     std::string     matrixFile             = "";
00085     double          tol                    = 1e-5;
00086     bool            onlyPrintOptions       = false;
00087     bool            printXmlFormat         = false;
00088     bool            showDoc                = true;
00089 
00090     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00091 
00092     CommandLineProcessor  clp(false); // Don't throw exceptions
00093 
00094     // Set up command-line options for the linear solver that will be used!
00095     linearSolverBuilder.setupCLP(&clp);
00096 
00097     clp.setOption( "matrix-file", &matrixFile
00098                    ,"Defines the matrix and perhaps the RHS and LHS for a linear system to be solved." );
00099     clp.setOption( "tol", &tol, "Tolerance to check against the scaled residual norm." );
00100     clp.setOption( "only-print-options", "continue-after-printing-options", &onlyPrintOptions
00101                    ,"Only print options and stop or continue on" );
00102     clp.setOption( "print-xml-format", "print-readable-format", &printXmlFormat
00103                    ,"Print the valid options in XML format or in readable format." );
00104     clp.setOption( "show-doc", "hide-doc", &showDoc
00105                    ,"Print the valid options with or without documentation." );
00106     
00107     clp.setDocString(
00108       "Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.\n"
00109       "\n"
00110       "To print out just the valid options use --matrix-file=\"\" --only-print-options with --print-xml-format"
00111       " or --print-readable-format.\n"
00112       "\n"
00113       "To solve a linear system from a system read from a file use --matrix-file=\"SomeFile.mtx\""
00114       " with options read from an XML file using --linear-solver-params-file=\"SomeFile.xml\"\n"
00115       );
00116 
00117     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00118     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00119 
00120     //
00121     // Print out the valid options and stop if asked
00122     //
00123 
00124     if(onlyPrintOptions) {
00125       if(printXmlFormat)
00126         Teuchos::writeParameterListToXmlOStream(
00127           *linearSolverBuilder.getValidParameters()
00128           ,*out
00129           );
00130       else
00131         linearSolverBuilder.getValidParameters()->print(
00132           *out,PLPrintOptions().indent(2).showTypes(true).showDoc(showDoc)
00133           );
00134       return 0;
00135     }
00136 
00137     
00138     //
00139     // B) Epetra-specific code that sets up the linear system to be solved
00140     //
00141     // While the below code reads in the Epetra objects from a file, you can
00142     // setup the Epetra objects any way you would like.  Note that this next
00143     // set of code as nothing to do with Thyra at all, and it should not.
00144     //
00145 
00146     *out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00147     
00148 #ifdef HAVE_MPI
00149     Epetra_MpiComm comm(MPI_COMM_WORLD);
00150 #else
00151     Epetra_SerialComm comm;
00152 #endif
00153     RCP<Epetra_CrsMatrix> epetra_A;
00154     RCP<Epetra_Vector> epetra_x, epetra_b;
00155     EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, &epetra_x, &epetra_b );
00156 
00157     if(!epetra_b.get()) {
00158       *out << "\nThe RHS b was not read in so generate a new random vector ...\n";
00159       epetra_b = rcp(new Epetra_Vector(epetra_A->OperatorRangeMap()));
00160       epetra_b->Random();
00161     }
00162 
00163     if(!epetra_x.get()) {
00164       *out << "\nThe LHS x was not read in so generate a new zero vector ...\n";
00165       epetra_x = rcp(new Epetra_Vector(epetra_A->OperatorDomainMap()));
00166       epetra_x->PutScalar(0.0); // Initial guess is critical!
00167     }
00168 
00169     *out << "\nPrinting statistics of the Epetra linear system ...\n";
00170 
00171     *out
00172       << "\n  Epetra_CrsMatrix epetra_A of dimension "
00173       << epetra_A->OperatorRangeMap().NumGlobalElements()
00174       << " x " << epetra_A->OperatorDomainMap().NumGlobalElements()
00175       << "\n  ||epetraA||inf = " << epetra_A->NormInf()
00176       << "\n  ||epetra_b||2 = " << epetraNorm2(*epetra_b)
00177       << "\n  ||epetra_x||2 = " << epetraNorm2(*epetra_x)
00178       << "\n";
00179 
00180 
00181     //
00182     // C) The "Glue" code that takes Epetra objects and wraps them as Thyra
00183     // objects
00184     //
00185     // This next set of code wraps the Epetra objects that define the linear
00186     // system to be solved as Thyra objects so that they can be passed to the
00187     // linear solver.
00188     //
00189 
00190     RCP<const Thyra::LinearOpBase<double> >
00191       A = Thyra::epetraLinearOp( epetra_A );
00192     RCP<Thyra::VectorBase<double> >
00193       x = Thyra::create_Vector( epetra_x, A->domain() );
00194     RCP<const Thyra::VectorBase<double> >
00195       b = Thyra::create_Vector( epetra_b, A->range() );
00196 
00197     // Note that above Thyra is only interacted with in the most trival of
00198     // ways.  For most users, Thyra will only be seen as a thin wrapper that
00199     // they need know little about in order to wrap their objects in order to
00200     // pass them to Thyra-enabled solvers.
00201 
00202       
00203     //
00204     // D) Thyra-specific code for solving the linear system
00205     //
00206     // Note that this code has no mention of any concrete implementation and
00207     // therefore can be used in any use case.
00208     //
00209 
00210     // Reading in the solver parameters from the parameters file and/or from
00211     // the command line.  This was setup by the command-line options
00212     // set by the setupCLP(...) function above.
00213     linearSolverBuilder.readParameters(out.get());
00214     // Create a linear solver factory given information read from the
00215     // parameter list.
00216     RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
00217       lowsFactory = linearSolverBuilder.createLinearSolveStrategy("");
00218     // Setup output stream and the verbosity level
00219     lowsFactory->setOStream(out);
00220     lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
00221     // Create a linear solver based on the forward operator A
00222     RCP<Thyra::LinearOpWithSolveBase<double> >
00223       lows = Thyra::linearOpWithSolve(*lowsFactory,A);
00224     // Solve the linear system (note: the initial guess in 'x' is critical)
00225     Thyra::SolveStatus<double>
00226       status = Thyra::solve(*lows,Thyra::NOTRANS,*b,&*x);
00227     *out << "\nSolve status:\n" << status;
00228     // Write the linear solver parameters after they were read
00229     linearSolverBuilder.writeParamsFile(*lowsFactory);
00230 
00231 
00232     //
00233     // E) Post process the solution and check the error
00234     //
00235     // Note that the below code is based only on the Epetra objects themselves
00236     // and does not in any way depend or interact with any Thyra-based
00237     // objects.  The point is that most users of Thyra can largely gloss over
00238     // the fact that Thyra is really being used for anything.
00239     //
00240 
00241     // Wipe out the Thyra wrapper for x to guarantee that the solution will be
00242     // written back to epetra_x!  At the time of this writting this is not
00243     // really needed but the behavior may change at some point so this is a
00244     // good idea.
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     {
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

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

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        std::stringName of an XML file containing parameters for linear solver options to be appended first.
                                             (default: --linear-solver-params-file="")
  --extra-linear-solver-params       std::stringAn XML string containing linear solver parameters to be appended second.
                                             (default: --extra-linear-solver-params="")
  --linear-solver-params-used-file   std::stringName 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                      std::stringDefines 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)
  --show-doc                         bool    Print the valid options with or without documentation.
  --hide-doc                                 (default: --show-doc)

DETAILED DOCUMENTATION:

Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.

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"


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 isUsed="true" 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:

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

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 


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 Amesos_Klu ...
  
  
  Total solve time = 0 sec

Solve status:
  solveStatus = SOLVE_STATUS_CONVERGED
  achievedTol = unknownTolerance()
  message:
    Solver Amesos_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|| = 1.11022e-16/1.52593 = 7.27571e-17 < 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 "Stratimikos_DefaultLinearSolverBuilder.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::RCP;
  using Teuchos::CommandLineProcessor;
  typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;

  bool success = true;
  bool verbose = true;

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

  Teuchos::RCP<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;
    bool            showDoc                = true;

    Stratimikos::DefaultLinearSolverBuilder 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.setOption( "show-doc", "hide-doc", &showDoc
                   ,"Print the valid options with or without documentation." );
    
    clp.setDocString(
      "Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.\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,PLPrintOptions().indent(2).showTypes(true).showDoc(showDoc)
          );
      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
    RCP<Epetra_CrsMatrix> epetra_A;
    RCP<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.
    //

    RCP<const Thyra::LinearOpBase<double> >
      A = Thyra::epetraLinearOp( epetra_A );
    RCP<Thyra::VectorBase<double> >
      x = Thyra::create_Vector( epetra_x, A->domain() );
    RCP<const Thyra::VectorBase<double> >
      b = Thyra::create_Vector( epetra_b, A->range() );

    // 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.
    RCP<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
    RCP<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!  At the time of this writting this is not
    // really needed but the behavior may change at some point so this is a
    // good idea.
    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);
    {
      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 Sun Nov 23 12:18:26 2008 for Stratimikos by  doxygen 1.3.9.1