Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MixedOrderPhysicsBasedPreconditioner.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stratimikos: Thyra-based strategies for linear solvers
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
43 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
44 #include "Thyra_PreconditionerFactoryHelpers.hpp"
45 #include "Thyra_DefaultInverseLinearOp.hpp"
46 #include "Thyra_DefaultMultipliedLinearOp.hpp"
48 #include "Thyra_EpetraLinearOp.hpp"
50 #include "Thyra_TestingTools.hpp"
51 #include "EpetraExt_CrsMatrixIn.h"
52 #include "Epetra_CrsMatrix.h"
58 #include "Teuchos_TimeMonitor.hpp"
59 #ifdef HAVE_MPI
60 # include "Epetra_MpiComm.h"
61 #else
62 # include "Epetra_SerialComm.h"
63 #endif
64 
65 
75 namespace {
76 
77 
78 // Read a Epetra_CrsMatrix from a Matrix market file
80 readEpetraCrsMatrixFromMatrixMarket(
81  const std::string fileName, const Epetra_Comm &comm
82  )
83 {
84  Epetra_CrsMatrix *A = 0;
86  0!=EpetraExt::MatrixMarketFileToCrsMatrix( fileName.c_str(), comm, A ),
87  std::runtime_error,
88  "Error, could not read matrix file \""<<fileName<<"\"!"
89  );
90  return Teuchos::rcp(A);
91 }
92 
93 
94 // Read an Epetra_CrsMatrix in as a wrapped Thyra::EpetraLinearOp object
96 readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
97  const std::string fileName, const Epetra_Comm &comm,
98  const std::string &label
99  )
100 {
101  return Thyra::epetraLinearOp(
102  readEpetraCrsMatrixFromMatrixMarket(fileName,comm),label
103  );
104 }
105 
106 
107 } // namespace
108 
109 
110 int main(int argc, char* argv[])
111 {
112 
113  using Teuchos::describe;
114  using Teuchos::rcp;
115  using Teuchos::rcp_dynamic_cast;
116  using Teuchos::rcp_const_cast;
117  using Teuchos::RCP;
120  using Teuchos::sublist;
121  typedef ParameterList::PrintOptions PLPrintOptions;
122  using Thyra::inverse;
123  using Thyra::initializePreconditionedOp;
124  using Thyra::initializeOp;
125  using Thyra::unspecifiedPrec;
126  using Thyra::solve;
127  typedef RCP<const Thyra::LinearOpBase<double> > LinearOpPtr;
128  typedef RCP<Thyra::VectorBase<double> > VectorPtr;
129 
130  bool success = true;
131  bool verbose = true;
132 
133  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
134 
136  out = Teuchos::VerboseObjectBase::getDefaultOStream();
137 
138  try {
139 
140  //
141  // Read in options from the command line
142  //
143 
144  const int numVerbLevels = 6;
146  verbLevelValues[] =
147  {
151  };
152  const char*
153  verbLevelNames[] =
154  { "default", "none", "low", "medium", "high", "extreme" };
155 
157  std::string paramsFile = "";
158  std::string extraParamsFile = "";
159  std::string baseDir = ".";
160  bool useP1Prec = true;
161  bool invertP1 = false;
162  bool showParams = false;
163  double solveTol = 1e-8;
164 
165  CommandLineProcessor clp(false); // Don't throw exceptions
166 
167  clp.setOption( "verb-level", &verbLevel,
168  numVerbLevels, verbLevelValues, verbLevelNames,
169  "Verbosity level used for all objects."
170  );
171  clp.setOption( "params-file", &paramsFile,
172  "File containing parameters in XML format.",
173  true // Required
174  );
175  clp.setOption( "extra-params-file", &extraParamsFile,
176  "File containing extra parameters in XML format."
177  );
178  clp.setOption( "base-dir", &baseDir,
179  "Base directory for all of the files."
180  );
181  clp.setOption( "use-P1-prec", "use-algebraic-prec", &useP1Prec,
182  "Use the physics-based preconditioner for P2 based on P1 or just an algebraic preconditioner"
183  );
184  clp.setOption( "invert-P1", "prec-P1-only", &invertP1,
185  "In the physics-based preconditioner, use P1's preconditioner or fully invert P1."
186  );
187  clp.setOption( "solve-tol", &solveTol,
188  "Tolerance for the solution to determine success or failure!"
189  );
190  clp.setOption( "show-params", "no-show-params", &showParams,
191  "Show the parameter list or not."
192  );
193  clp.setDocString(
194  "This example program shows an attempt to create a physics-based\n"
195  "preconditioner using the building blocks of Stratimikos and Thyra\n"
196  "implicit linear operators. The idea is to use the linear operator\n"
197  "for first-order Lagrange finite elements as the preconditioner for the\n"
198  "operator using second-order Lagrange finite elements. The details of the\n"
199  "PDE being represented, the mesh being used, or the boundary conditions are\n"
200  "beyond the scope of this example.\n"
201  "\n"
202  "The matrices used in this example are:\n"
203  "\n"
204  " P2: The discretized PDE matrix using second-order Lagrange finite elements.\n"
205  " P1: The discretized PDE matrix using first-order Lagrange finite elements.\n"
206  " M22: The mass matrix for the second-order Lagrange finite-element basis functions.\n"
207  " M11: The mass matrix for the first-order Lagrange finite-element basis functions.\n"
208  " M21: A rectangular matrix that uses mixed first- and second-order basis functions\n"
209  " to map from P2 space to P1 space.\n"
210  " M12: A rectangular matrix that uses mixed first- and second-order basis functions\n"
211  " to map from P1 space to P2 space.\n"
212  "\n"
213  "The above matrices are read from Matrix Market *.mtx files in the directory given by\n"
214  "the --base-dir command-line option.\n"
215  "\n"
216  "The preconditioner operator created in this example program is:\n"
217  "\n"
218  " precP2Op = (inv(M22) * M12) * prec(P1) * (inv(M11) * M21)\n"
219  "\n"
220  "where prec(P1) is either the algebraic preconditioner for P1 (--prec-P1-only)\n"
221  "or is a full solve for P1 (--invert-P1).\n"
222  "\n"
223  "We use Stratimikos to specify the linear solvers and/or algebraic\n"
224  "preconditioners and we use the Thyra implicit operators to build the\n"
225  "implicitly multiplied linear operators associated with the preconditioner.\n"
226  "\n"
227  "Warning! This physics-based preconditioner is singular and can not\n"
228  "be used to solve the linear system given a random right-hand side. However,\n"
229  "singular or not, this example shows how one can use Thyra/Stratimikos to quickly\n"
230  "try out these types of preconditioning ideas (even if they do not work).\n"
231  );
232 
233  // Note: Use --help on the command line to see the above documentation
234 
235  CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
236  if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
237 
238 
239  //
240  *out << "\nA) Reading in Epetra_CrsMatrix objects for P1, P2, M11, M12, M21, and M22 ...\n";
241  //
242 
243 #ifdef HAVE_MPI
244  Epetra_MpiComm comm(MPI_COMM_WORLD);
245 #else
246  Epetra_SerialComm comm;
247 #endif
248 
249  LinearOpPtr P1=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
250  baseDir+"/P1.mtx",comm,"P1");
251  *out << "\nP1 = " << describe(*P1,verbLevel) << "\n";
252  LinearOpPtr P2= readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
253  baseDir+"/P2.mtx",comm,"P2");
254  *out << "\nP2 = " << describe(*P2,verbLevel) << "\n";
255  LinearOpPtr M11=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
256  baseDir+"/M11.mtx",comm,"M11");
257  *out << "\nM11 = " << describe(*M11,verbLevel) << "\n";
258  LinearOpPtr M22=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
259  baseDir+"/M22.mtx",comm,"M22");
260  *out << "\nM22 = " << describe(*M22,verbLevel) << "\n";
261  LinearOpPtr M12=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
262  baseDir+"/M12.mtx",comm,"M12");
263  *out << "\nM12 = " << describe(*M12,verbLevel) << "\n";
264  LinearOpPtr M21=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
265  baseDir+"/M21.mtx",comm,"M21");
266  *out << "\nM21 = " << describe(*M21,verbLevel) << "\n";
267 
268  // ToDo: Replace the above functions with a general Thyra strategy object
269  // to do the reading
270 
271  //
272  *out << "\nB) Get the preconditioner and/or linear solver strategies to invert M11, M22, P1, and P2 ...\n";
273  //
274 
275  //
276  // Get separate parameter sublists for each square operator separately
277  // that specify the type of linear solver and/or preconditioner to use.
278  //
279 
280  RCP<ParameterList> paramList =
281  Teuchos::getParametersFromXmlFile( baseDir+"/"+paramsFile );
282  if (extraParamsFile.length()) {
283  Teuchos::updateParametersFromXmlFile( baseDir+"/"+extraParamsFile, paramList.ptr() );
284  }
285  if (showParams) {
286  *out << "\nRead in parameter list:\n\n";
287  paramList->print(*out,PLPrintOptions().indent(2).showTypes(true));
288  }
289 
290  Stratimikos::DefaultLinearSolverBuilder M11_linsolve_strategy_builder;
291  M11_linsolve_strategy_builder.setParameterList(
292  sublist(paramList,"M11 Solver",true) );
293 
294  Stratimikos::DefaultLinearSolverBuilder M22_linsolve_strategy_builder;
295  M22_linsolve_strategy_builder.setParameterList(
296  sublist(paramList,"M22 Solver",true) );
297 
298  Stratimikos::DefaultLinearSolverBuilder P1_linsolve_strategy_builder;
299  P1_linsolve_strategy_builder.setParameterList(
300  sublist(paramList,"P1 Solver",true) );
301 
302  Stratimikos::DefaultLinearSolverBuilder P2_linsolve_strategy_builder;
303  P2_linsolve_strategy_builder.setParameterList(
304  sublist(paramList,"P2 Solver",true) );
305 
306  //
307  // Create the linear solver and/or preconditioner strategies
308  // (i.e. factories)
309  //
310 
311  // For M11 and M22, we want full linear solver factories with embedded
312  // algebraic preconditioner factories.
313 
314  RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > M11_linsolve_strategy
315  = createLinearSolveStrategy(M11_linsolve_strategy_builder);
316 
317  RCP<Thyra::LinearOpWithSolveFactoryBase<double> > M22_linsolve_strategy
318  = createLinearSolveStrategy(M22_linsolve_strategy_builder);
319 
320  // For P1, we only want its preconditioner factory.
321 
322  RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P1_linsolve_strategy;
323  RCP<const Thyra::PreconditionerFactoryBase<double> > P1_prec_strategy;
324  if(invertP1)
325  P1_linsolve_strategy
326  = createLinearSolveStrategy(P1_linsolve_strategy_builder);
327  else
328  P1_prec_strategy
329  = createPreconditioningStrategy(P1_linsolve_strategy_builder);
330 
331  // For P2, we only want a linear solver factory. We will supply the
332  // preconditioner ourselves (that is the whole point of this example).
333 
334  RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P2_linsolve_strategy
335  = createLinearSolveStrategy(P2_linsolve_strategy_builder);
336 
337  //
338  *out << "\nC) Create the physics-based preconditioner! ...\n";
339  //
340 
341  *out << "\nCreating inv(M11) ...\n";
342  LinearOpPtr invM11 = inverse(*M11_linsolve_strategy, M11);
343  *out << "\ninvM11 = " << describe(*invM11,verbLevel) << "\n";
344 
345  *out << "\nCreating inv(M22) ...\n";
346  LinearOpPtr invM22 = inverse(*M22_linsolve_strategy, M22);
347  *out << "\ninvM22 = " << describe(*invM22,verbLevel) << "\n";
348 
349  *out << "\nCreating prec(P1) ...\n";
350  LinearOpPtr invP1;
351  if(invertP1) {
352  *out << "\nCreating prec(P1) as a full solver ...\n";
353  invP1 = inverse(*P1_linsolve_strategy, P1);
354  }
355  else {
356  *out << "\nCreating prec(P1) as just an algebraic preconditioner ...\n";
357  RCP<Thyra::PreconditionerBase<double> >
358  precP1 = prec(*P1_prec_strategy,P1);
359  *out << "\nprecP1 = " << describe(*precP1,verbLevel) << "\n";
360  invP1 = precP1->getUnspecifiedPrecOp();
361  }
362  rcp_const_cast<Thyra::LinearOpBase<double> >(
363  invP1)->setObjectLabel("invP1"); // Cast to set label ...
364  *out << "\ninvP1 = " << describe(*invP1,verbLevel) << "\n";
365 
366  LinearOpPtr P2ToP1 = multiply( invM11, M21 );
367  *out << "\nP2ToP1 = " << describe(*P2ToP1,verbLevel) << "\n";
368 
369  LinearOpPtr P1ToP2 = multiply( invM22, M12 );
370  *out << "\nP1ToP2 = " << describe(*P1ToP2,verbLevel) << "\n";
371 
372  LinearOpPtr precP2Op = multiply( P1ToP2, invP1, P2ToP1 );
373  *out << "\nprecP2Op = " << describe(*precP2Op,verbLevel) << "\n";
374 
375  //
376  *out << "\nD) Setup the solver for P2 ...\n";
377  //
378 
379  RCP<Thyra::LinearOpWithSolveBase<double> >
380  P2_lows = P2_linsolve_strategy->createOp();
381  if(useP1Prec) {
382  *out << "\nCreating the solver P2 using the specialized precP2Op\n";
383  initializePreconditionedOp<double>( *P2_linsolve_strategy, P2,
384  unspecifiedPrec(precP2Op), P2_lows.ptr());
385  }
386  else {
387  *out << "\nCreating the solver P2 using algebraic preconditioner\n";
388  initializeOp(*P2_linsolve_strategy, P2, P2_lows.ptr());
389  }
390  *out << "\nP2_lows = " << describe(*P2_lows, verbLevel) << "\n";
391 
392  //
393  *out << "\nE) Solve P2 for a random RHS ...\n";
394  //
395 
396  VectorPtr x = createMember(P2->domain());
397  VectorPtr b = createMember(P2->range());
398  Thyra::randomize(-1.0, +1.0, b.ptr());
399  Thyra::assign(x.ptr(), 0.0); // Must give an initial guess!
400 
401  Thyra::SolveStatus<double>
402  solveStatus = solve<double>( *P2_lows, Thyra::NOTRANS, *b, x.ptr() );
403 
404  *out << "\nSolve status:\n" << solveStatus;
405 
406  *out << "\nSolution ||x|| = " << Thyra::norm(*x) << "\n";
407 
408  if(showParams) {
409  *out << "\nParameter list after use:\n\n";
410  paramList->print(*out, PLPrintOptions().indent(2).showTypes(true));
411  }
412 
413  //
414  *out << "\nF) Checking the error in the solution of r=b-P2*x ...\n";
415  //
416 
417  VectorPtr P2x = Thyra::createMember(b->space());
418  Thyra::apply( *P2, Thyra::NOTRANS, *x, P2x.ptr() );
419  VectorPtr r = Thyra::createMember(b->space());
420  Thyra::V_VmV<double>(r.ptr(), *b, *P2x);
421 
422  double
423  P2x_nrm = Thyra::norm(*P2x),
424  r_nrm = Thyra::norm(*r),
425  b_nrm = Thyra::norm(*b),
426  r_nrm_over_b_nrm = r_nrm / b_nrm;
427 
428  bool result = r_nrm_over_b_nrm <= solveTol;
429  if(!result) success = false;
430 
431  *out
432  << "\n||P2*x|| = " << P2x_nrm << "\n";
433 
434  *out
435  << "\n||P2*x-b||/||b|| = " << r_nrm << "/" << b_nrm
436  << " = " << r_nrm_over_b_nrm << " <= " << solveTol
437  << " : " << Thyra::passfail(result) << "\n";
438 
440  }
441  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success)
442 
443  if (verbose) {
444  if(success) *out << "\nCongratulations! All of the tests checked out!\n";
445  else *out << "\nOh no! At least one of the tests failed!\n";
446  }
447 
448  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
449 
450 }
int main(int argc, char *argv[])
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setParameterList(RCP< ParameterList > const &paramList)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
ParameterList::PrintOptions PLPrintOptions
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...