Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
48 #include "Thyra_EpetraLinearOp.hpp"
50 #include "Thyra_TestingTools.hpp"
51 #include "Thyra_LinearOpTester.hpp"
52 #include "EpetraExt_CrsMatrixIn.h"
53 #include "Epetra_CrsMatrix.h"
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
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  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 
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;
149  verbLevelValues[] =
150  {
154  };
155  const char*
156  verbLevelNames[] =
157  { "default", "none", "low", "medium", "high", "extreme" };
158 
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 
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[])
#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)
Ptr< T > ptr() const
ParameterList::PrintOptions PLPrintOptions
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...