71int main(
int argc,
char* argv[])
76 using Teuchos::CommandLineProcessor;
77 typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
82 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
84 Teuchos::RCP<Teuchos::FancyOStream>
85 out = Teuchos::VerboseObjectBase::getDefaultOStream();
98 std::string matrixFile =
"";
99 std::string extraParamsFile =
"";
101 bool onlyPrintOptions =
false;
102 bool printXmlFormat =
false;
107 CommandLineProcessor clp(
false);
110 linearSolverBuilder.setupCLP(&clp);
112 clp.setOption(
"matrix-file", &matrixFile
113 ,
"Defines the matrix and perhaps the RHS and LHS for a linear system to be solved." );
114 clp.setOption(
"tol", &tol,
"Tolerance to check against the scaled residual norm." );
115 clp.setOption(
"extra-params-file", &extraParamsFile,
"File containing extra parameters in XML format.");
116 clp.setOption(
"only-print-options",
"continue-after-printing-options", &onlyPrintOptions
117 ,
"Only print options and stop or continue on" );
118 clp.setOption(
"print-xml-format",
"print-readable-format", &printXmlFormat
119 ,
"Print the valid options in XML format or in readable format." );
120 clp.setOption(
"show-doc",
"hide-doc", &showDoc
121 ,
"Print the valid options with or without documentation." );
124 "Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.\n"
126 "To print out just the valid options use --matrix-file=\"\" --only-print-options with --print-xml-format"
127 " or --print-readable-format.\n"
129 "To solve a linear system read from a file use --matrix-file=\"SomeFile.mtx\""
130 " with options read from an XML file using --linear-solver-params-file=\"SomeFile.xml\"\n"
133 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
134 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
return parse_return;
140 if(onlyPrintOptions) {
142 Teuchos::writeParameterListToXmlOStream(
143 *linearSolverBuilder.getValidParameters()
147 linearSolverBuilder.getValidParameters()->print(
148 *out,PLPrintOptions().indent(2).showTypes(
true).showDoc(showDoc)
162 *out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile<<
"\' ...\n";
165 Epetra_MpiComm comm(MPI_COMM_WORLD);
167 Epetra_SerialComm comm;
169 RCP<Epetra_CrsMatrix> epetra_A;
170 RCP<Epetra_Vector> epetra_x, epetra_b;
171 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, &epetra_x, &epetra_b );
173 if(!epetra_b.get()) {
174 *out <<
"\nThe RHS b was not read in so generate a new random vector ...\n";
175 epetra_b = rcp(
new Epetra_Vector(epetra_A->OperatorRangeMap()));
179 if(!epetra_x.get()) {
180 *out <<
"\nThe LHS x was not read in so generate a new zero vector ...\n";
181 epetra_x = rcp(
new Epetra_Vector(epetra_A->OperatorDomainMap()));
182 epetra_x->PutScalar(0.0);
185 *out <<
"\nPrinting statistics of the Epetra linear system ...\n";
188 <<
"\n Epetra_CrsMatrix epetra_A of dimension "
189 << epetra_A->OperatorRangeMap().NumGlobalElements()
190 <<
" x " << epetra_A->OperatorDomainMap().NumGlobalElements()
191 <<
"\n ||epetraA||inf = " << epetra_A->NormInf()
192 <<
"\n ||epetra_b||2 = " << epetraNorm2(*epetra_b)
193 <<
"\n ||epetra_x||2 = " << epetraNorm2(*epetra_x)
206 RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp( epetra_A );
207 RCP<Thyra::VectorBase<double> > x = Thyra::create_Vector( epetra_x, A->domain() );
208 RCP<const Thyra::VectorBase<double> > b = Thyra::create_Vector( epetra_b, A->range() );
226 linearSolverBuilder.readParameters(out.get());
229 if(extraParamsFile.length()) {
230 Teuchos::updateParametersFromXmlFile(
"./"+extraParamsFile,
231 linearSolverBuilder.getNonconstParameterList().ptr() );
236 RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
237 linearSolverBuilder.createLinearSolveStrategy(
"");
240 lowsFactory->setOStream(out);
241 lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
244 RCP<Thyra::LinearOpWithSolveBase<double> > lows =
245 Thyra::linearOpWithSolve(*lowsFactory, A);
248 Thyra::SolveStatus<double> status =
249 Thyra::solve<double>(*lows, Thyra::NOTRANS, *b, x.ptr());
250 *out <<
"\nSolve status:\n" << status;
253 linearSolverBuilder.writeParamsFile(*lowsFactory);
272 <<
"\nSolution ||epetra_x||2 = " << epetraNorm2(*epetra_x) <<
"\n";
274 *out <<
"\nTesting the solution error ||b-A*x||/||b|| computed through the Epetra objects ...\n";
277 Epetra_Vector epetra_r(*epetra_b);
279 Epetra_Vector epetra_A_x(epetra_A->OperatorRangeMap());
280 epetra_A->Apply(*epetra_x,epetra_A_x);
281 epetra_r.Update(-1.0,epetra_A_x,1.0);
285 nrm_r = epetraNorm2(epetra_r),
286 nrm_b = epetraNorm2(*epetra_b),
287 rel_err = ( nrm_r / nrm_b );
289 passed = (rel_err <= tol);
292 <<
"||b-A*x||/||b|| = " << nrm_r <<
"/" << nrm_b <<
" = " << rel_err
293 <<
" < tol = " << tol <<
" ? " << ( passed ?
"passed" :
"failed" ) <<
"\n";
295 if(!passed) success =
false;
298 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success)
301 if(success) *out <<
"\nCongratulations! All of the tests checked out!\n";
302 else *out <<
"\nOh no! At least one of the tests failed!\n";
305 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );