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"
61# include "Epetra_MpiComm.h"
63# include "Epetra_SerialComm.h"
79Teuchos::RCP<Epetra_CrsMatrix>
80readEpetraCrsMatrixFromMatrixMarket(
81 const std::string fileName,
const Epetra_Comm &comm
84 Epetra_CrsMatrix *A = 0;
85 TEUCHOS_TEST_FOR_EXCEPTION(
86 0!=EpetraExt::MatrixMarketFileToCrsMatrix( fileName.c_str(), comm, A ),
88 "Error, could not read matrix file \""<<fileName<<
"\"!"
90 return Teuchos::rcp(A);
95Teuchos::RCP<const Thyra::LinearOpBase<double> >
96readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
97 const std::string fileName,
const Epetra_Comm &comm,
98 const std::string &label
101 return Thyra::epetraLinearOp(
102 readEpetraCrsMatrixFromMatrixMarket(fileName,comm),label
110int main(
int argc,
char* argv[])
113 using Teuchos::describe;
115 using Teuchos::rcp_dynamic_cast;
116 using Teuchos::rcp_const_cast;
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;
128 typedef RCP<const Thyra::LinearOpBase<double> > LinearOpPtr;
129 typedef RCP<Thyra::VectorBase<double> > VectorPtr;
134 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
136 Teuchos::RCP<Teuchos::FancyOStream>
137 out = Teuchos::VerboseObjectBase::getDefaultOStream();
145 CommandLineProcessor clp(
false);
147 const int numVerbLevels = 6;
148 Teuchos::EVerbosityLevel
151 Teuchos::VERB_DEFAULT, Teuchos::VERB_NONE,
152 Teuchos::VERB_LOW, Teuchos::VERB_MEDIUM,
153 Teuchos::VERB_HIGH, Teuchos::VERB_EXTREME
157 {
"default",
"none",
"low",
"medium",
"high",
"extreme" };
159 Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_MEDIUM;
160 clp.setOption(
"verb-level", &verbLevel,
161 numVerbLevels, verbLevelValues, verbLevelNames,
162 "Verbosity level used for all objects."
165 std::string matrixFile =
".";
166 clp.setOption(
"matrix-file", &matrixFile,
170 std::string paramListFile =
"";
171 clp.setOption(
"param-list-file", ¶mListFile,
172 "Parameter list for preconditioner and solver blocks."
175 bool showParams =
false;
176 clp.setOption(
"show-params",
"no-show-params", &showParams,
177 "Show the parameter list or not."
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."
185 double solveTol = 1e-8;
186 clp.setOption(
"solve-tol", &solveTol,
187 "Tolerance for the solution to determine success or failure!"
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"
197 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
198 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
return parse_return;
202 *out <<
"\nA) Reading in the matrix ...\n";
206 Epetra_MpiComm comm(MPI_COMM_WORLD);
208 Epetra_SerialComm comm;
211 const LinearOpPtr A = readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
212 matrixFile, comm,
"A");
213 *out <<
"\nA = " << describe(*A,verbLevel) <<
"\n";
215 const RCP<ParameterList> paramList = getParametersFromXmlFile(paramListFile);
217 *out <<
"\nRead in parameter list:\n\n";
218 paramList->print(*out, PLPrintOptions().indent(2).showTypes(
true));
222 *out <<
"\nB) Get the preconditioner as a forward solver\n";
225 const RCP<ParameterList> precParamList = sublist(paramList,
"Preconditioner Solver");
228 const RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > precSolverStrategy
229 = createLinearSolveStrategy(precSolverBuilder);
232 const LinearOpPtr A_inv_prec = inverse<double>(*precSolverStrategy, A,
233 Thyra::SUPPORT_SOLVE_FORWARD_ONLY,
235 Thyra::IGNORE_SOLVE_FAILURE
237 *out <<
"\nA_inv_prec = " << describe(*A_inv_prec, verbLevel) <<
"\n";
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; }
248 *out <<
"\nC) Create the forward solver using the created preconditioner ...\n";
251 const RCP<ParameterList> fwdSolverParamList = sublist(paramList,
"Forward Solver");
254 const RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > fwdSolverSolverStrategy
255 = createLinearSolveStrategy(fwdSolverSolverBuilder);
257 const RCP<Thyra::LinearOpWithSolveBase<double> >
258 A_lows = fwdSolverSolverStrategy->createOp();
260 initializePreconditionedOp<double>( *fwdSolverSolverStrategy, A,
261 unspecifiedPrec(A_inv_prec), A_lows.ptr());
263 *out <<
"\nA_lows = " << describe(*A_lows, verbLevel) <<
"\n";
266 *out <<
"\nD) Solve the linear system for a random RHS ...\n";
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);
274 Thyra::SolveStatus<double>
275 solveStatus = solve<double>( *A_lows, Thyra::NOTRANS, *b, x.ptr() );
277 *out <<
"\nSolve status:\n" << solveStatus;
279 *out <<
"\nSolution ||x|| = " << Thyra::norm(*x) <<
"\n";
282 *out <<
"\nParameter list after use:\n\n";
283 paramList->print(*out, PLPrintOptions().indent(2).showTypes(
true));
287 *out <<
"\nF) Checking the error in the solution of r=b-A*x ...\n";
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);
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;
301 bool resid_tol_check = ( r_nrm_over_b_nrm <= solveTol );
302 if(!resid_tol_check) success =
false;
305 <<
"\n||A*x|| = " << Ax_nrm <<
"\n";
308 <<
"\n||A*x-b||/||b|| = " << r_nrm <<
"/" << b_nrm
309 <<
" = " << r_nrm_over_b_nrm <<
" <= " << solveTol
310 <<
" : " << Thyra::passfail(resid_tol_check) <<
"\n";
312 Teuchos::TimeMonitor::summarize(*out<<
"\n");
314 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success)
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";
321 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
int main(int argc, char *argv[])
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...
void setParameterList(RCP< ParameterList > const ¶mList)