Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
MixedOrderPhysicsBasedPreconditioner.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 "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00031 #include "Thyra_PreconditionerFactoryHelpers.hpp"
00032 #include "Thyra_DefaultInverseLinearOp.hpp"
00033 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00034 #include "Thyra_EpetraThyraWrappers.hpp"
00035 #include "Thyra_EpetraLinearOp.hpp"
00036 #include "Thyra_get_Epetra_Operator.hpp"
00037 #include "Thyra_TestingTools.hpp"
00038 #include "EpetraExt_CrsMatrixIn.h"
00039 #include "Epetra_CrsMatrix.h"
00040 #include "Teuchos_GlobalMPISession.hpp"
00041 #include "Teuchos_VerboseObject.hpp"
00042 #include "Teuchos_XMLParameterListHelpers.hpp"
00043 #include "Teuchos_CommandLineProcessor.hpp"
00044 #include "Teuchos_StandardCatchMacros.hpp"
00045 #include "Teuchos_TimeMonitor.hpp"
00046 #ifdef HAVE_MPI
00047 #  include "Epetra_MpiComm.h"
00048 #else
00049 #  include "Epetra_SerialComm.h"
00050 #endif
00051 
00052 
00062 namespace {
00063 
00064 
00065 // Read a Epetra_CrsMatrix from a Matrix market file
00066 Teuchos::RCP<Epetra_CrsMatrix>
00067 readEpetraCrsMatrixFromMatrixMarket(
00068   const std::string fileName, const Epetra_Comm &comm
00069   )
00070 {
00071   Epetra_CrsMatrix *A = 0;
00072   TEUCHOS_TEST_FOR_EXCEPTION(
00073     0!=EpetraExt::MatrixMarketFileToCrsMatrix( fileName.c_str(), comm, A ),
00074     std::runtime_error,
00075     "Error, could not read matrix file \""<<fileName<<"\"!"
00076     );
00077   return Teuchos::rcp(A);
00078 }
00079 
00080 
00081 // Read an Epetra_CrsMatrix in as a wrapped Thyra::EpetraLinearOp object
00082 Teuchos::RCP<const Thyra::LinearOpBase<double> >
00083 readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
00084   const std::string fileName, const Epetra_Comm &comm,
00085   const std::string &label
00086   )
00087 {
00088   return Thyra::epetraLinearOp(
00089     readEpetraCrsMatrixFromMatrixMarket(fileName,comm),label
00090     );
00091 }
00092 
00093 
00094 } // namespace
00095 
00096 
00097 int main(int argc, char* argv[])
00098 {
00099 
00100   using Teuchos::describe;
00101   using Teuchos::rcp;
00102   using Teuchos::rcp_dynamic_cast;
00103   using Teuchos::rcp_const_cast;
00104   using Teuchos::RCP;
00105   using Teuchos::CommandLineProcessor;
00106   using Teuchos::ParameterList;
00107   using Teuchos::sublist;
00108   typedef ParameterList::PrintOptions PLPrintOptions;
00109   using Thyra::inverse;
00110   using Thyra::initializePreconditionedOp;
00111   using Thyra::initializeOp;
00112   using Thyra::unspecifiedPrec;
00113   using Thyra::solve;
00114   typedef RCP<const Thyra::LinearOpBase<double> > LinearOpPtr;
00115   typedef RCP<Thyra::VectorBase<double> > VectorPtr;
00116   
00117   bool success = true;
00118   bool verbose = true;
00119 
00120   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00121 
00122   Teuchos::RCP<Teuchos::FancyOStream>
00123     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00124 
00125   try {
00126 
00127     //
00128     // Read in options from the command line
00129     //
00130 
00131     const int numVerbLevels = 6;
00132     Teuchos::EVerbosityLevel
00133       verbLevelValues[] =
00134       {
00135         Teuchos::VERB_DEFAULT, Teuchos::VERB_NONE,
00136         Teuchos::VERB_LOW, Teuchos::VERB_MEDIUM,
00137         Teuchos::VERB_HIGH, Teuchos::VERB_EXTREME
00138       };
00139     const char*
00140       verbLevelNames[] =
00141       { "default", "none", "low", "medium", "high", "extreme" };
00142 
00143     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_MEDIUM;
00144     std::string paramsFile = "";
00145     std::string extraParamsFile = "";
00146     std::string baseDir = ".";
00147     bool useP1Prec = true;
00148     bool invertP1 = false;
00149     bool showParams = false;
00150     double solveTol = 1e-8;
00151 
00152     CommandLineProcessor clp(false); // Don't throw exceptions
00153 
00154     clp.setOption( "verb-level", &verbLevel,
00155       numVerbLevels, verbLevelValues, verbLevelNames,
00156       "Verbosity level used for all objects."
00157       );
00158     clp.setOption( "params-file", &paramsFile,
00159       "File containing parameters in XML format.",
00160       true // Required
00161       );
00162     clp.setOption( "extra-params-file", &extraParamsFile,
00163       "File containing extra parameters in XML format."
00164       );
00165     clp.setOption( "base-dir", &baseDir,
00166       "Base directory for all of the files."
00167       );
00168     clp.setOption( "use-P1-prec", "use-algebraic-prec", &useP1Prec,
00169       "Use the physics-based preconditioner for P2 based on P1 or just an algebraic preconditioner"
00170       );
00171     clp.setOption( "invert-P1", "prec-P1-only", &invertP1,
00172       "In the physics-based preconditioner, use P1's preconditioner or fully invert P1."
00173       );
00174     clp.setOption( "solve-tol", &solveTol,
00175       "Tolerance for the solution to determine success or failure!"
00176       );
00177     clp.setOption( "show-params", "no-show-params", &showParams,
00178       "Show the parameter list or not."
00179       );
00180     clp.setDocString(
00181       "This example program shows an attempt to create a physics-based\n"
00182       "preconditioner using the building blocks of Stratimikos and Thyra\n"
00183       "implicit linear operators.  The idea is to use the linear operator\n"
00184       "for first-order Lagrange finite elements as the preconditioner for the\n"
00185       "operator using second-order Lagrange finite elements.  The details of the\n"
00186       "PDE being represented, the mesh being used, or the boundary conditions are\n"
00187       "beyond the scope of this example.\n"
00188       "\n"
00189       "The matrices used in this example are:\n"
00190       "\n"
00191       "  P2: The discretized PDE matrix using second-order Lagrange finite elements.\n"
00192       "  P1: The discretized PDE matrix using first-order Lagrange finite elements.\n"
00193       "  M22: The mass matrix for the second-order Lagrange finite-element basis functions.\n"
00194       "  M11: The mass matrix for the first-order Lagrange finite-element basis functions.\n"
00195       "  M21: A rectangular matrix that uses mixed first- and second-order basis functions\n"
00196       "       to map from P2 space to P1 space.\n"
00197       "  M12: A rectangular matrix that uses mixed first- and second-order basis functions\n"
00198       "       to map from P1 space to P2 space.\n"
00199       "\n"
00200       "The above matrices are read from Matrix Market *.mtx files in the directory given by\n"
00201       "the --base-dir command-line option.\n"
00202       "\n"
00203       "The preconditioner operator created in this example program is:\n"
00204       "\n"
00205       "   precP2Op = (inv(M22) * M12) * prec(P1) * (inv(M11) * M21)\n"
00206       "\n"
00207       "where prec(P1) is either the algebraic preconditioner for P1 (--prec-P1-only)\n"
00208       "or is a full solve for P1 (--invert-P1).\n"
00209       "\n"
00210       "We use Stratimikos to specify the linear solvers and/or algebraic\n"
00211       "preconditioners and we use the Thyra implicit operators to build the\n"
00212       "implicitly multiplied linear operators associated with the preconditioner.\n"
00213       "\n"
00214       "Warning!  This physics-based preconditioner is singular and can not\n"
00215       "be used to solve the linear system given a random right-hand side.  However,\n"
00216       "singular or not, this example shows how one can use Thyra/Stratimikos to quickly\n"
00217       "try out these types of preconditioning ideas (even if they do not work).\n"
00218       );
00219 
00220     // Note: Use --help on the command line to see the above documentation
00221 
00222     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00223     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00224 
00225 
00226     //
00227     *out << "\nA) Reading in Epetra_CrsMatrix objects for P1, P2, M11, M12, M21, and M22 ...\n";
00228     //
00229     
00230 #ifdef HAVE_MPI
00231     Epetra_MpiComm comm(MPI_COMM_WORLD);
00232 #else
00233     Epetra_SerialComm comm;
00234 #endif
00235 
00236     LinearOpPtr P1=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
00237       baseDir+"/P1.mtx",comm,"P1");
00238     *out << "\nP1 = " << describe(*P1,verbLevel) << "\n";
00239     LinearOpPtr P2= readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
00240       baseDir+"/P2.mtx",comm,"P2");
00241     *out << "\nP2 = " << describe(*P2,verbLevel) << "\n";
00242     LinearOpPtr M11=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
00243       baseDir+"/M11.mtx",comm,"M11");
00244     *out << "\nM11 = " << describe(*M11,verbLevel) << "\n";
00245     LinearOpPtr M22=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
00246       baseDir+"/M22.mtx",comm,"M22");
00247     *out << "\nM22 = " << describe(*M22,verbLevel) << "\n";
00248     LinearOpPtr M12=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
00249       baseDir+"/M12.mtx",comm,"M12");
00250     *out << "\nM12 = " << describe(*M12,verbLevel) << "\n";
00251     LinearOpPtr M21=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
00252       baseDir+"/M21.mtx",comm,"M21");
00253     *out << "\nM21 = " << describe(*M21,verbLevel) << "\n";
00254 
00255     // ToDo: Replace the above functions with a general Thyra strategy object
00256     // to do the reading
00257     
00258     //
00259     *out << "\nB) Get the preconditioner and/or linear solver strategies to invert M11, M22, P1, and P2 ...\n";
00260     //
00261 
00262     //
00263     // Get separate parameter sublists for each square operator separately
00264     // that specify the type of linear solver and/or preconditioner to use.
00265     //
00266 
00267     RCP<ParameterList> paramList =
00268       Teuchos::getParametersFromXmlFile( baseDir+"/"+paramsFile );
00269     if (extraParamsFile.length()) {
00270       Teuchos::updateParametersFromXmlFile( baseDir+"/"+extraParamsFile, paramList.ptr() );
00271     }
00272     if (showParams) {
00273       *out << "\nRead in parameter list:\n\n";
00274       paramList->print(*out,PLPrintOptions().indent(2).showTypes(true));
00275     }
00276     
00277     Stratimikos::DefaultLinearSolverBuilder M11_linsolve_strategy_builder;
00278     M11_linsolve_strategy_builder.setParameterList(
00279       sublist(paramList,"M11 Solver",true) );
00280 
00281     Stratimikos::DefaultLinearSolverBuilder M22_linsolve_strategy_builder;
00282     M22_linsolve_strategy_builder.setParameterList(
00283       sublist(paramList,"M22 Solver",true) );
00284 
00285     Stratimikos::DefaultLinearSolverBuilder P1_linsolve_strategy_builder;
00286     P1_linsolve_strategy_builder.setParameterList(
00287       sublist(paramList,"P1 Solver",true) );
00288 
00289     Stratimikos::DefaultLinearSolverBuilder P2_linsolve_strategy_builder;
00290     P2_linsolve_strategy_builder.setParameterList(
00291       sublist(paramList,"P2 Solver",true) );
00292 
00293     // 
00294     // Create the linear solver and/or preconditioner strategies
00295     // (i.e. factories)
00296     //
00297 
00298     // For M11 and M22, we want full linear solver factories with embedded
00299     // algebraic preconditioner factories.
00300 
00301     RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > M11_linsolve_strategy
00302       = createLinearSolveStrategy(M11_linsolve_strategy_builder);
00303       
00304     RCP<Thyra::LinearOpWithSolveFactoryBase<double> > M22_linsolve_strategy
00305       = createLinearSolveStrategy(M22_linsolve_strategy_builder);
00306       
00307     // For P1, we only want its preconditioner factory.
00308 
00309     RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P1_linsolve_strategy;
00310     RCP<const Thyra::PreconditionerFactoryBase<double> > P1_prec_strategy;
00311     if(invertP1)
00312       P1_linsolve_strategy
00313         = createLinearSolveStrategy(P1_linsolve_strategy_builder);
00314     else
00315       P1_prec_strategy
00316         = createPreconditioningStrategy(P1_linsolve_strategy_builder);
00317 
00318     // For P2, we only want a linear solver factory.  We will supply the
00319     // preconditioner ourselves (that is the whole point of this example).
00320 
00321     RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P2_linsolve_strategy
00322       = createLinearSolveStrategy(P2_linsolve_strategy_builder);
00323 
00324     //
00325     *out << "\nC) Create the physics-based preconditioner! ...\n";
00326     //
00327 
00328     *out << "\nCreating inv(M11) ...\n";
00329     LinearOpPtr invM11 = inverse(*M11_linsolve_strategy, M11);
00330     *out << "\ninvM11 = " << describe(*invM11,verbLevel) << "\n";
00331 
00332     *out << "\nCreating inv(M22) ...\n";
00333     LinearOpPtr invM22 = inverse(*M22_linsolve_strategy, M22);
00334     *out << "\ninvM22 = " << describe(*invM22,verbLevel) << "\n";
00335 
00336     *out << "\nCreating prec(P1) ...\n";
00337     LinearOpPtr invP1;
00338     if(invertP1) {
00339       *out << "\nCreating prec(P1) as a full solver ...\n";
00340       invP1 = inverse(*P1_linsolve_strategy, P1);
00341     }
00342     else {
00343       *out << "\nCreating prec(P1) as just an algebraic preconditioner ...\n";
00344       RCP<Thyra::PreconditionerBase<double> >
00345         precP1 = prec(*P1_prec_strategy,P1);
00346       *out << "\nprecP1 = " << describe(*precP1,verbLevel) << "\n";
00347       invP1 = precP1->getUnspecifiedPrecOp();
00348     }
00349     rcp_const_cast<Thyra::LinearOpBase<double> >(
00350       invP1)->setObjectLabel("invP1"); // Cast to set label ...
00351     *out << "\ninvP1 = " << describe(*invP1,verbLevel) << "\n";
00352 
00353     LinearOpPtr P2ToP1 = multiply( invM11, M21 );
00354     *out << "\nP2ToP1 = " << describe(*P2ToP1,verbLevel) << "\n";
00355 
00356     LinearOpPtr P1ToP2 = multiply( invM22, M12 );
00357     *out << "\nP1ToP2 = " << describe(*P1ToP2,verbLevel) << "\n";
00358 
00359     LinearOpPtr precP2Op = multiply( P1ToP2, invP1, P2ToP1 );
00360     *out << "\nprecP2Op = " << describe(*precP2Op,verbLevel) << "\n";
00361       
00362     //
00363     *out << "\nD) Setup the solver for P2 ...\n";
00364     //
00365 
00366     RCP<Thyra::LinearOpWithSolveBase<double> >
00367       P2_lows = P2_linsolve_strategy->createOp();
00368     if(useP1Prec) {
00369       *out << "\nCreating the solver P2 using the specialized precP2Op\n";
00370       initializePreconditionedOp<double>( *P2_linsolve_strategy, P2,
00371         unspecifiedPrec(precP2Op), P2_lows.ptr());
00372     }
00373     else {
00374       *out << "\nCreating the solver P2 using algebraic preconditioner\n";
00375       initializeOp(*P2_linsolve_strategy, P2, P2_lows.ptr());
00376     }
00377     *out << "\nP2_lows = " << describe(*P2_lows, verbLevel) << "\n";
00378     
00379     //
00380     *out << "\nE) Solve P2 for a random RHS ...\n";
00381     //
00382 
00383     VectorPtr x = createMember(P2->domain());
00384     VectorPtr b = createMember(P2->range());
00385     Thyra::randomize(-1.0, +1.0, b.ptr());
00386     Thyra::assign(x.ptr(), 0.0); // Must give an initial guess!
00387 
00388     Thyra::SolveStatus<double>
00389       solveStatus = solve<double>( *P2_lows, Thyra::NOTRANS, *b, x.ptr() );
00390 
00391     *out << "\nSolve status:\n" << solveStatus;
00392 
00393     *out << "\nSolution ||x|| = " << Thyra::norm(*x) << "\n";
00394 
00395     if(showParams) {
00396       *out << "\nParameter list after use:\n\n";
00397       paramList->print(*out, PLPrintOptions().indent(2).showTypes(true));
00398     }
00399 
00400     //
00401     *out << "\nF) Checking the error in the solution of r=b-P2*x ...\n";
00402     //
00403     
00404     VectorPtr P2x = Thyra::createMember(b->space());
00405     Thyra::apply( *P2, Thyra::NOTRANS, *x, P2x.ptr() );
00406     VectorPtr r = Thyra::createMember(b->space());
00407     Thyra::V_VmV<double>(r.ptr(), *b, *P2x);
00408     
00409     double
00410       P2x_nrm = Thyra::norm(*P2x),
00411       r_nrm = Thyra::norm(*r),
00412       b_nrm = Thyra::norm(*b),
00413       r_nrm_over_b_nrm = r_nrm / b_nrm;
00414 
00415     bool result = r_nrm_over_b_nrm <= solveTol;
00416     if(!result) success = false;
00417     
00418     *out
00419       << "\n||P2*x|| = " << P2x_nrm << "\n";
00420     
00421     *out
00422       << "\n||P2*x-b||/||b|| = " << r_nrm << "/" << b_nrm
00423       << " = " << r_nrm_over_b_nrm << " <= " << solveTol
00424       << " : " << Thyra::passfail(result) << "\n";
00425     
00426   }
00427   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success)
00428     
00429   Teuchos::TimeMonitor::summarize(*out<<"\n");
00430   
00431   if (verbose) {
00432     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00433     else         *out << "\nOh no! At least one of the tests failed!\n";
00434   }
00435 
00436   return ( success ? 0 : 1 );
00437 
00438 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines