Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
75namespace {
76
77
78// Read a Epetra_CrsMatrix from a Matrix market file
79Teuchos::RCP<Epetra_CrsMatrix>
80readEpetraCrsMatrixFromMatrixMarket(
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
95Teuchos::RCP<const Thyra::LinearOpBase<double> >
96readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
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
110int 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");
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[])
LinearSolverBuilder< double > DefaultLinearSolverBuilder