Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
ForwardSolverAsPreconditioner.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"
47 #include "Thyra_EpetraThyraWrappers.hpp"
48 #include "Thyra_EpetraLinearOp.hpp"
49 #include "Thyra_get_Epetra_Operator.hpp"
50 #include "Thyra_TestingTools.hpp"
51 #include "Thyra_LinearOpTester.hpp"
52 #include "EpetraExt_CrsMatrixIn.h"
53 #include "Epetra_CrsMatrix.h"
54 #include "Teuchos_GlobalMPISession.hpp"
55 #include "Teuchos_VerboseObject.hpp"
56 #include "Teuchos_XMLParameterListHelpers.hpp"
57 #include "Teuchos_CommandLineProcessor.hpp"
58 #include "Teuchos_StandardCatchMacros.hpp"
59 #include "Teuchos_TimeMonitor.hpp"
60 #ifdef HAVE_MPI
61 # include "Epetra_MpiComm.h"
62 #else
63 # include "Epetra_SerialComm.h"
64 #endif
65 
66 
75 namespace {
76 
77 
78 // Read a Epetra_CrsMatrix from a Matrix market file
79 Teuchos::RCP<Epetra_CrsMatrix>
80 readEpetraCrsMatrixFromMatrixMarket(
81  const std::string fileName, const Epetra_Comm &comm
82  )
83 {
84  Epetra_CrsMatrix *A = 0;
85  TEUCHOS_TEST_FOR_EXCEPTION(
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
95 Teuchos::RCP<const Thyra::LinearOpBase<double> >
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;
118  using Teuchos::CommandLineProcessor;
119  using Teuchos::ParameterList;
120  using Teuchos::sublist;
121  using Teuchos::getParametersFromXmlFile;
122  typedef ParameterList::PrintOptions PLPrintOptions;
123  using Thyra::inverse;
124  using Thyra::initializePreconditionedOp;
125  using Thyra::initializeOp;
126  using Thyra::unspecifiedPrec;
127  using Thyra::solve;
128  typedef RCP<const Thyra::LinearOpBase<double> > LinearOpPtr;
129  typedef RCP<Thyra::VectorBase<double> > VectorPtr;
130 
131  bool success = true;
132  bool verbose = true;
133 
134  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
135 
136  Teuchos::RCP<Teuchos::FancyOStream>
137  out = Teuchos::VerboseObjectBase::getDefaultOStream();
138 
139  try {
140 
141  //
142  // Read in options from the command line
143  //
144 
145  CommandLineProcessor clp(false); // Don't throw exceptions
146 
147  const int numVerbLevels = 6;
148  Teuchos::EVerbosityLevel
149  verbLevelValues[] =
150  {
151  Teuchos::VERB_DEFAULT, Teuchos::VERB_NONE,
152  Teuchos::VERB_LOW, Teuchos::VERB_MEDIUM,
153  Teuchos::VERB_HIGH, Teuchos::VERB_EXTREME
154  };
155  const char*
156  verbLevelNames[] =
157  { "default", "none", "low", "medium", "high", "extreme" };
158 
159  Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_MEDIUM;
160  clp.setOption( "verb-level", &verbLevel,
161  numVerbLevels, verbLevelValues, verbLevelNames,
162  "Verbosity level used for all objects."
163  );
164 
165  std::string matrixFile = ".";
166  clp.setOption( "matrix-file", &matrixFile,
167  "Matrix file."
168  );
169 
170  std::string paramListFile = "";
171  clp.setOption( "param-list-file", &paramListFile,
172  "Parameter list for preconditioner and solver blocks."
173  );
174 
175  bool showParams = false;
176  clp.setOption( "show-params", "no-show-params", &showParams,
177  "Show the parameter list or not."
178  );
179 
180  bool testPrecIsLinearOp = true;
181  clp.setOption( "test-prec-is-linear-op", "test-prec-is-linear-op", &testPrecIsLinearOp,
182  "Test if the preconditioner is a linear operator or not."
183  );
184 
185  double solveTol = 1e-8;
186  clp.setOption( "solve-tol", &solveTol,
187  "Tolerance for the solution to determine success or failure!"
188  );
189 
190  clp.setDocString(
191  "This example program shows how to use one linear solver (e.g. AztecOO)\n"
192  "as a preconditioner for another iterative solver (e.g. Belos).\n"
193  );
194 
195  // Note: Use --help on the command line to see the above documentation
196 
197  CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
198  if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
199 
200 
201  //
202  *out << "\nA) Reading in the matrix ...\n";
203  //
204 
205 #ifdef HAVE_MPI
206  Epetra_MpiComm comm(MPI_COMM_WORLD);
207 #else
208  Epetra_SerialComm comm;
209 #endif
210 
211  const LinearOpPtr A = readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
212  matrixFile, comm, "A");
213  *out << "\nA = " << describe(*A,verbLevel) << "\n";
214 
215  const RCP<ParameterList> paramList = getParametersFromXmlFile(paramListFile);
216  if (showParams) {
217  *out << "\nRead in parameter list:\n\n";
218  paramList->print(*out, PLPrintOptions().indent(2).showTypes(true));
219  }
220 
221  //
222  *out << "\nB) Get the preconditioner as a forward solver\n";
223  //
224 
225  const RCP<ParameterList> precParamList = sublist(paramList, "Preconditioner Solver");
226  Stratimikos::DefaultLinearSolverBuilder precSolverBuilder;
227  precSolverBuilder.setParameterList(precParamList);
228  const RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > precSolverStrategy
229  = createLinearSolveStrategy(precSolverBuilder);
230  //precSolverStrategy->setVerbLevel(verbLevel);
231 
232  const LinearOpPtr A_inv_prec = inverse<double>(*precSolverStrategy, A,
233  Thyra::SUPPORT_SOLVE_FORWARD_ONLY,
234  Teuchos::null, // Use internal solve criteria
235  Thyra::IGNORE_SOLVE_FAILURE // Ignore solve failures since this is just a prec
236  );
237  *out << "\nA_inv_prec = " << describe(*A_inv_prec, verbLevel) << "\n";
238 
239  if (testPrecIsLinearOp) {
240  *out << "\nTest that the preconditioner A_inv_prec is indeed a linear operator.\n";
241  Thyra::LinearOpTester<double> linearOpTester;
242  linearOpTester.check_adjoint(false);
243  const bool linearOpCheck = linearOpTester.check(*A_inv_prec, out.ptr());
244  if (!linearOpCheck) { success = false; }
245  }
246 
247  //
248  *out << "\nC) Create the forward solver using the created preconditioner ...\n";
249  //
250 
251  const RCP<ParameterList> fwdSolverParamList = sublist(paramList, "Forward Solver");
252  Stratimikos::DefaultLinearSolverBuilder fwdSolverSolverBuilder;
253  fwdSolverSolverBuilder.setParameterList(fwdSolverParamList);
254  const RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > fwdSolverSolverStrategy
255  = createLinearSolveStrategy(fwdSolverSolverBuilder);
256 
257  const RCP<Thyra::LinearOpWithSolveBase<double> >
258  A_lows = fwdSolverSolverStrategy->createOp();
259 
260  initializePreconditionedOp<double>( *fwdSolverSolverStrategy, A,
261  unspecifiedPrec(A_inv_prec), A_lows.ptr());
262  //A_lows->setVerbLevel(verbLevel);
263  *out << "\nA_lows = " << describe(*A_lows, verbLevel) << "\n";
264 
265  //
266  *out << "\nD) Solve the linear system for a random RHS ...\n";
267  //
268 
269  VectorPtr x = createMember(A->domain());
270  VectorPtr b = createMember(A->range());
271  Thyra::randomize(-1.0, +1.0, b.ptr());
272  Thyra::assign(x.ptr(), 0.0); // Must give an initial guess!
273 
274  Thyra::SolveStatus<double>
275  solveStatus = solve<double>( *A_lows, Thyra::NOTRANS, *b, x.ptr() );
276 
277  *out << "\nSolve status:\n" << solveStatus;
278 
279  *out << "\nSolution ||x|| = " << Thyra::norm(*x) << "\n";
280 
281  if(showParams) {
282  *out << "\nParameter list after use:\n\n";
283  paramList->print(*out, PLPrintOptions().indent(2).showTypes(true));
284  }
285 
286  //
287  *out << "\nF) Checking the error in the solution of r=b-A*x ...\n";
288  //
289 
290  VectorPtr Ax = Thyra::createMember(b->space());
291  Thyra::apply( *A, Thyra::NOTRANS, *x, Ax.ptr() );
292  VectorPtr r = Thyra::createMember(b->space());
293  Thyra::V_VmV<double>(r.ptr(), *b, *Ax);
294 
295  double
296  Ax_nrm = Thyra::norm(*Ax),
297  r_nrm = Thyra::norm(*r),
298  b_nrm = Thyra::norm(*b),
299  r_nrm_over_b_nrm = r_nrm / b_nrm;
300 
301  bool resid_tol_check = ( r_nrm_over_b_nrm <= solveTol );
302  if(!resid_tol_check) success = false;
303 
304  *out
305  << "\n||A*x|| = " << Ax_nrm << "\n";
306 
307  *out
308  << "\n||A*x-b||/||b|| = " << r_nrm << "/" << b_nrm
309  << " = " << r_nrm_over_b_nrm << " <= " << solveTol
310  << " : " << Thyra::passfail(resid_tol_check) << "\n";
311 
312  Teuchos::TimeMonitor::summarize(*out<<"\n");
313  }
314  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success)
315 
316  if (verbose) {
317  if(success) *out << "\nCongratulations! All of the tests checked out!\n";
318  else *out << "\nOh no! At least one of the tests failed!\n";
319  }
320 
321  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
322 }
int main(int argc, char *argv[])
void setParameterList(RCP< ParameterList > const &paramList)
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...