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 "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 }

Generated on Thu Sep 18 12:37:35 2008 for Stratimikos Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1