Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
simple_stratimikos_example.cpp
Go to the documentation of this file.
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 "EpetraExt_readEpetraLinearSystem.h"
00031 #include "Teuchos_GlobalMPISession.hpp"
00032 #include "Teuchos_VerboseObject.hpp"
00033 #include "Teuchos_XMLParameterListHelpers.hpp"
00034 #include "Teuchos_CommandLineProcessor.hpp"
00035 #include "Teuchos_StandardCatchMacros.hpp"
00036 #ifdef HAVE_MPI
00037 #  include "Epetra_MpiComm.h"
00038 #else
00039 #  include "Epetra_SerialComm.h"
00040 #endif
00041 
00042 
00043 namespace {
00044 
00045 
00046 // Helper function to compute a single norm for a vector
00047 double epetraNorm2( const Epetra_Vector &v )
00048 {
00049   double norm[1] = { -1.0 };
00050   v.Norm2(&norm[0]);
00051   return norm[0];
00052 }
00053 
00054 
00055 } // namespace
00056 
00057 
00058 int main(int argc, char* argv[])
00059 {
00060 
00061   using Teuchos::rcp;
00062   using Teuchos::RCP;
00063   using Teuchos::CommandLineProcessor;
00064   typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
00065 
00066   bool success = true;
00067   bool verbose = true;
00068 
00069   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00070 
00071   Teuchos::RCP<Teuchos::FancyOStream>
00072     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00073 
00074   try {
00075 
00076     
00077     //
00078     // A) Program setup code
00079     //
00080 
00081     //
00082     // Read options from command-line
00083     //
00084     
00085     std::string     matrixFile             = "";
00086     std::string     extraParamsFile        = "";
00087     double          tol                    = 1e-5;
00088     bool            onlyPrintOptions       = false;
00089     bool            printXmlFormat         = false;
00090     bool            showDoc                = true;
00091 
00092     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00093 
00094     CommandLineProcessor  clp(false); // Don't throw exceptions
00095 
00096     // Set up command-line options for the linear solver that will be used!
00097     linearSolverBuilder.setupCLP(&clp);
00098 
00099     clp.setOption( "matrix-file", &matrixFile
00100                    ,"Defines the matrix and perhaps the RHS and LHS for a linear system to be solved." );
00101     clp.setOption( "tol", &tol, "Tolerance to check against the scaled residual norm." );
00102     clp.setOption( "extra-params-file", &extraParamsFile, "File containing extra parameters in XML format.");
00103     clp.setOption( "only-print-options", "continue-after-printing-options", &onlyPrintOptions
00104                    ,"Only print options and stop or continue on" );
00105     clp.setOption( "print-xml-format", "print-readable-format", &printXmlFormat
00106                    ,"Print the valid options in XML format or in readable format." );
00107     clp.setOption( "show-doc", "hide-doc", &showDoc
00108                    ,"Print the valid options with or without documentation." );
00109     
00110     clp.setDocString(
00111       "Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.\n"
00112       "\n"
00113       "To print out just the valid options use --matrix-file=\"\" --only-print-options with --print-xml-format"
00114       " or --print-readable-format.\n"
00115       "\n"
00116       "To solve a linear system read from a file use --matrix-file=\"SomeFile.mtx\""
00117       " with options read from an XML file using --linear-solver-params-file=\"SomeFile.xml\"\n"
00118       );
00119 
00120     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00121     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00122 
00123     //
00124     // Print out the valid options and stop if asked
00125     //
00126 
00127     if(onlyPrintOptions) {
00128       if(printXmlFormat)
00129         Teuchos::writeParameterListToXmlOStream(
00130           *linearSolverBuilder.getValidParameters()
00131           ,*out
00132           );
00133       else
00134         linearSolverBuilder.getValidParameters()->print(
00135           *out,PLPrintOptions().indent(2).showTypes(true).showDoc(showDoc)
00136           );
00137       return 0;
00138     }
00139 
00140     
00141     //
00142     // B) Epetra-specific code that sets up the linear system to be solved
00143     //
00144     // While the below code reads in the Epetra objects from a file, you can
00145     // setup the Epetra objects any way you would like.  Note that this next
00146     // set of code as nothing to do with Thyra at all, and it should not.
00147     //
00148 
00149     *out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00150     
00151 #ifdef HAVE_MPI
00152     Epetra_MpiComm comm(MPI_COMM_WORLD);
00153 #else
00154     Epetra_SerialComm comm;
00155 #endif
00156     RCP<Epetra_CrsMatrix> epetra_A;
00157     RCP<Epetra_Vector> epetra_x, epetra_b;
00158     EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, &epetra_x, &epetra_b );
00159 
00160     if(!epetra_b.get()) {
00161       *out << "\nThe RHS b was not read in so generate a new random vector ...\n";
00162       epetra_b = rcp(new Epetra_Vector(epetra_A->OperatorRangeMap()));
00163       epetra_b->Random();
00164     }
00165 
00166     if(!epetra_x.get()) {
00167       *out << "\nThe LHS x was not read in so generate a new zero vector ...\n";
00168       epetra_x = rcp(new Epetra_Vector(epetra_A->OperatorDomainMap()));
00169       epetra_x->PutScalar(0.0); // Initial guess is critical!
00170     }
00171 
00172     *out << "\nPrinting statistics of the Epetra linear system ...\n";
00173 
00174     *out
00175       << "\n  Epetra_CrsMatrix epetra_A of dimension "
00176       << epetra_A->OperatorRangeMap().NumGlobalElements()
00177       << " x " << epetra_A->OperatorDomainMap().NumGlobalElements()
00178       << "\n  ||epetraA||inf = " << epetra_A->NormInf()
00179       << "\n  ||epetra_b||2 = " << epetraNorm2(*epetra_b)
00180       << "\n  ||epetra_x||2 = " << epetraNorm2(*epetra_x)
00181       << "\n";
00182 
00183 
00184     //
00185     // C) The "Glue" code that takes Epetra objects and wraps them as Thyra
00186     // objects
00187     //
00188     // This next set of code wraps the Epetra objects that define the linear
00189     // system to be solved as Thyra objects so that they can be passed to the
00190     // linear solver.
00191     //
00192 
00193     RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp( epetra_A );
00194     RCP<Thyra::VectorBase<double> > x = Thyra::create_Vector( epetra_x, A->domain() );
00195     RCP<const Thyra::VectorBase<double> > 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 
00215     // Augment parameters if appropriate
00216     if(extraParamsFile.length()) {
00217       Teuchos::updateParametersFromXmlFile( "./"+extraParamsFile,
00218         &*linearSolverBuilder.getNonconstParameterList() );
00219     }
00220 
00221     // Create a linear solver factory given information read from the
00222     // parameter list.
00223     RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
00224       linearSolverBuilder.createLinearSolveStrategy("");
00225 
00226     // Setup output stream and the verbosity level
00227     lowsFactory->setOStream(out);
00228     lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
00229 
00230     // Create a linear solver based on the forward operator A
00231     RCP<Thyra::LinearOpWithSolveBase<double> > lows =
00232       Thyra::linearOpWithSolve(*lowsFactory, A);
00233 
00234     // Solve the linear system (note: the initial guess in 'x' is critical)
00235     Thyra::SolveStatus<double> status =
00236       Thyra::solve<double>(*lows, Thyra::NOTRANS, *b, x.ptr());
00237     *out << "\nSolve status:\n" << status;
00238 
00239     // Write the linear solver parameters after they were read
00240     linearSolverBuilder.writeParamsFile(*lowsFactory);
00241 
00242 
00243     //
00244     // E) Post process the solution and check the error
00245     //
00246     // Note that the below code is based only on the Epetra objects themselves
00247     // and does not in any way depend or interact with any Thyra-based
00248     // objects.  The point is that most users of Thyra can largely gloss over
00249     // the fact that Thyra is really being used for anything.
00250     //
00251 
00252     // Wipe out the Thyra wrapper for x to guarantee that the solution will be
00253     // written back to epetra_x!  At the time of this writting this is not
00254     // really needed but the behavior may change at some point so this is a
00255     // good idea.
00256     x = Teuchos::null;
00257 
00258     *out
00259       << "\nSolution ||epetra_x||2 = " << epetraNorm2(*epetra_x) << "\n";
00260 
00261     *out << "\nTesting the solution error ||b-A*x||/||b|| computed through the Epetra objects ...\n";
00262 
00263     // r = b - A*x
00264     Epetra_Vector epetra_r(*epetra_b);
00265     {
00266       Epetra_Vector epetra_A_x(epetra_A->OperatorRangeMap());
00267       epetra_A->Apply(*epetra_x,epetra_A_x);
00268       epetra_r.Update(-1.0,epetra_A_x,1.0);
00269     }
00270       
00271     const double
00272       nrm_r = epetraNorm2(epetra_r),
00273       nrm_b = epetraNorm2(*epetra_b),
00274       rel_err = ( nrm_r / nrm_b );
00275     const bool
00276       passed = (rel_err <= tol);
00277 
00278     *out
00279       << "||b-A*x||/||b|| = " << nrm_r << "/" << nrm_b << " = " << rel_err
00280       << " < tol = " << tol << " ? " << ( passed ? "passed" : "failed" ) << "\n";
00281 
00282     if(!passed) success = false;
00283     
00284   }
00285   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success)
00286   
00287   if (verbose) {
00288     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00289     else         *out << "\nOh no! At least one of the tests failed!\n";
00290   }
00291 
00292   return ( success ? 0 : 1 );
00293 
00294 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines