MixedOrderPhysicsBasedPreconditioner.cpp

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   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 );
00271     if(showParams) {
00272       *out << "\nRead in parameter list:\n\n";
00273       paramList->print(*out,PLPrintOptions().indent(2).showTypes(true));
00274     }
00275     
00276     Stratimikos::DefaultLinearSolverBuilder M11_linsolve_strategy_builder;
00277     M11_linsolve_strategy_builder.setParameterList(
00278       sublist(paramList,"M11 Solver",true) );
00279 
00280     Stratimikos::DefaultLinearSolverBuilder M22_linsolve_strategy_builder;
00281     M22_linsolve_strategy_builder.setParameterList(
00282       sublist(paramList,"M22 Solver",true) );
00283 
00284     Stratimikos::DefaultLinearSolverBuilder P1_linsolve_strategy_builder;
00285     P1_linsolve_strategy_builder.setParameterList(
00286       sublist(paramList,"P1 Solver",true) );
00287 
00288     Stratimikos::DefaultLinearSolverBuilder P2_linsolve_strategy_builder;
00289     P2_linsolve_strategy_builder.setParameterList(
00290       sublist(paramList,"P2 Solver",true) );
00291 
00292     // 
00293     // Create the linear solver and/or preconditioner strategies
00294     // (i.e. factories)
00295     //
00296 
00297     // For M11 and M22, we want full linear solver factories with embedded
00298     // algebraic preconditioner factories.
00299 
00300     RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > M11_linsolve_strategy
00301       = createLinearSolveStrategy(M11_linsolve_strategy_builder);
00302       
00303     RCP<Thyra::LinearOpWithSolveFactoryBase<double> > M22_linsolve_strategy
00304       = createLinearSolveStrategy(M22_linsolve_strategy_builder);
00305       
00306     // For P1, we only want its preconditioner factory.
00307 
00308     RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P1_linsolve_strategy;
00309     RCP<const Thyra::PreconditionerFactoryBase<double> > P1_prec_strategy;
00310     if(invertP1)
00311       P1_linsolve_strategy
00312         = createLinearSolveStrategy(P1_linsolve_strategy_builder);
00313     else
00314       P1_prec_strategy
00315         = createPreconditioningStrategy(P1_linsolve_strategy_builder);
00316 
00317     // For P2, we only want a linear solver factory.  We will supply the
00318     // preconditioner ourselves (that is the whole point of this example).
00319 
00320     RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P2_linsolve_strategy
00321       = createLinearSolveStrategy(P2_linsolve_strategy_builder);
00322 
00323     //
00324     *out << "\nC) Create the physics-based preconditioner! ...\n";
00325     //
00326 
00327     *out << "\nCreating inv(M11) ...\n";
00328     LinearOpPtr invM11 = inverse(*M11_linsolve_strategy,M11);
00329     *out << "\ninvM11 = " << describe(*invM11,verbLevel) << "\n";
00330 
00331     *out << "\nCreating inv(M22) ...\n";
00332     LinearOpPtr invM22 = inverse(*M22_linsolve_strategy,M22);
00333     *out << "\ninvM22 = " << describe(*invM22,verbLevel) << "\n";
00334 
00335     *out << "\nCreating prec(P1) ...\n";
00336     LinearOpPtr invP1;
00337     if(invertP1) {
00338       *out << "\nCreating prec(P1) as a full solver ...\n";
00339       invP1 = inverse(*P1_linsolve_strategy,P1);
00340     }
00341     else {
00342       *out << "\nCreating prec(P1) as just an algebraic preconditioner ...\n";
00343       RCP<Thyra::PreconditionerBase<double> >
00344         precP1 = prec(*P1_prec_strategy,P1);
00345       *out << "\nprecP1 = " << describe(*precP1,verbLevel) << "\n";
00346       invP1 = precP1->getUnspecifiedPrecOp();
00347     }
00348     rcp_const_cast<Thyra::LinearOpBase<double> >(
00349       invP1)->setObjectLabel("invP1"); // Cast to set label ...
00350     *out << "\ninvP1 = " << describe(*invP1,verbLevel) << "\n";
00351 
00352     LinearOpPtr P2ToP1 = multiply( invM11, M21 );
00353     *out << "\nP2ToP1 = " << describe(*P2ToP1,verbLevel) << "\n";
00354 
00355     LinearOpPtr P1ToP2 = multiply( invM22, M12 );
00356     *out << "\nP1ToP2 = " << describe(*P1ToP2,verbLevel) << "\n";
00357 
00358     LinearOpPtr precP2Op = multiply( P1ToP2, invP1, P2ToP1 );
00359     *out << "\nprecP2Op = " << describe(*precP2Op,verbLevel) << "\n";
00360       
00361     //
00362     *out << "\nD) Setup the solver for P2 ...\n";
00363     //
00364 
00365     RCP<Thyra::LinearOpWithSolveBase<double> >
00366       P2_lows = P2_linsolve_strategy->createOp();
00367     if(useP1Prec) {
00368       *out << "\nCreating the solver P2 using the specialized precP2Op\n";
00369       initializePreconditionedOp<double>( *P2_linsolve_strategy, P2,
00370         unspecifiedPrec(precP2Op), &*P2_lows );
00371     }
00372     else {
00373       *out << "\nCreating the solver P2 using algebraic preconditioner\n";
00374       initializeOp( *P2_linsolve_strategy, P2, &*P2_lows );
00375     }
00376     *out << "\nP2_lows = " << describe(*P2_lows,verbLevel) << "\n";
00377     
00378     //
00379     *out << "\nE) Solve P2 for a random RHS ...\n";
00380     //
00381 
00382     VectorPtr x = createMember(P2->domain());
00383     VectorPtr b = createMember(P2->range());
00384     Thyra::randomize(-1.0,+1.0,&*b);
00385     Thyra::assign(&*x,0.0); // Must give an initial guess!
00386 
00387     Thyra::SolveStatus<double>
00388       solveStatus = solve( *P2_lows, Thyra::NOTRANS, *b, &*x );
00389 
00390     *out << "\nSolve status:\n" << solveStatus;
00391 
00392     *out << "\nSolution ||x|| = " << Thyra::norm(*x) << "\n";
00393 
00394     if(showParams) {
00395       *out << "\nParameter list after use:\n\n";
00396       paramList->print(*out,PLPrintOptions().indent(2).showTypes(true));
00397     }
00398 
00399     //
00400     *out << "\nF) Checking the error in the solution of r=b-P2*x ...\n";
00401     //
00402     
00403     VectorPtr P2x = Thyra::createMember(b->space());
00404     Thyra::apply( *P2, Thyra::NOTRANS, *x, &*P2x );
00405     VectorPtr r = Thyra::createMember(b->space());
00406     Thyra::V_VmV(&*r,*b,*P2x);
00407     
00408     double
00409       P2x_nrm = Thyra::norm(*P2x),
00410       r_nrm = Thyra::norm(*r),
00411       b_nrm = Thyra::norm(*b),
00412       r_nrm_over_b_nrm = r_nrm / b_nrm;
00413 
00414     bool result = r_nrm_over_b_nrm <= solveTol;
00415     if(!result) success = false;
00416     
00417     *out
00418       << "\n||P2*x|| = " << P2x_nrm << "\n";
00419     
00420     *out
00421       << "\n||P2*x-b||/||b|| = " << r_nrm << "/" << b_nrm
00422       << " = " << r_nrm_over_b_nrm << " <= " << solveTol
00423       << " : " << Thyra::passfail(result) << "\n";
00424     
00425   }
00426   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success)
00427     
00428   Teuchos::TimeMonitor::summarize(*out<<"\n");
00429   
00430   if (verbose) {
00431     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00432     else         *out << "\nOh no! At least one of the tests failed!\n";
00433   }
00434 
00435   return ( success ? 0 : 1 );
00436 
00437 }

Generated on Wed May 12 21:59:20 2010 for Stratimikos by  doxygen 1.4.7