42#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
44#include "Teuchos_GlobalMPISession.hpp"
45#include "Teuchos_VerboseObject.hpp"
46#include "Teuchos_XMLParameterListHelpers.hpp"
47#include "Teuchos_CommandLineProcessor.hpp"
48#include "Teuchos_StandardCatchMacros.hpp"
50# include "Epetra_MpiComm.h"
52# include "Epetra_SerialComm.h"
62 double norm[1] = { -1.0 };
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;
105 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
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";
169 RCP<Epetra_CrsMatrix> epetra_A;
170 RCP<Epetra_Vector> 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";
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 );
int Norm2(double *Result) const
void readEpetraLinearSystem(const std::string &fileName, const Epetra_Comm &comm, Teuchos::RCP< Epetra_CrsMatrix > *A=NULL, Teuchos::RCP< Epetra_Map > *map=NULL, Teuchos::RCP< Epetra_Vector > *x=NULL, Teuchos::RCP< Epetra_Vector > *b=NULL, Teuchos::RCP< Epetra_Vector > *xExact=NULL)
#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
#include "Teuchos_GlobalMPISession.hpp"
#include "Teuchos_VerboseObject.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_StandardCatchMacros.hpp"
#ifdef HAVE_MPI
# include "Epetra_MpiComm.h"
#else
# include "Epetra_SerialComm.h"
#endif
namespace {
{
double norm[1] = { -1.0 };
v.Norm2(&norm[0]);
return norm[0];
}
}
int main(int argc, char* argv[])
{
using Teuchos::rcp;
using Teuchos::RCP;
using Teuchos::CommandLineProcessor;
typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
bool success = true;
bool verbose = true;
Teuchos::GlobalMPISession mpiSession(&argc, &argv);
Teuchos::RCP<Teuchos::FancyOStream>
out = Teuchos::VerboseObjectBase::getDefaultOStream();
try {
std::string matrixFile = "";
std::string extraParamsFile = "";
double tol = 1e-5;
bool onlyPrintOptions = false;
bool printXmlFormat = false;
bool showDoc = true;
Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
CommandLineProcessor clp(false);
linearSolverBuilder.setupCLP(&clp);
clp.setOption( "matrix-file", &matrixFile
,"Defines the matrix and perhaps the RHS and LHS for a linear system to be solved." );
clp.setOption( "tol", &tol, "Tolerance to check against the scaled residual norm." );
clp.setOption( "extra-params-file", &extraParamsFile, "File containing extra parameters in XML format.");
clp.setOption( "only-print-options", "continue-after-printing-options", &onlyPrintOptions
,"Only print options and stop or continue on" );
clp.setOption( "print-xml-format", "print-readable-format", &printXmlFormat
,"Print the valid options in XML format or in readable format." );
clp.setOption( "show-doc", "hide-doc", &showDoc
,"Print the valid options with or without documentation." );
clp.setDocString(
"Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.\n"
"\n"
"To print out just the valid options use --matrix-file=\"\" --only-print-options with --print-xml-format"
" or --print-readable-format.\n"
"\n"
"To solve a linear system read from a file use --matrix-file=\"SomeFile.mtx\""
" with options read from an XML file using --linear-solver-params-file=\"SomeFile.xml\"\n"
);
CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
if(onlyPrintOptions) {
if(printXmlFormat)
Teuchos::writeParameterListToXmlOStream(
*linearSolverBuilder.getValidParameters()
,*out
);
else
linearSolverBuilder.getValidParameters()->print(
*out,PLPrintOptions().indent(2).showTypes(true).showDoc(showDoc)
);
return 0;
}
*out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
#ifdef HAVE_MPI
#else
#endif
RCP<Epetra_CrsMatrix> epetra_A;
RCP<Epetra_Vector> epetra_x, epetra_b;
if(!epetra_b.get()) {
*out << "\nThe RHS b was not read in so generate a new random vector ...\n";
epetra_b->Random();
}
if(!epetra_x.get()) {
*out << "\nThe LHS x was not read in so generate a new zero vector ...\n";
epetra_x = rcp(
new Epetra_Vector(epetra_A->OperatorDomainMap()));
epetra_x->PutScalar(0.0);
}
*out << "\nPrinting statistics of the Epetra linear system ...\n";
*out
<< "\n Epetra_CrsMatrix epetra_A of dimension "
<< epetra_A->OperatorRangeMap().NumGlobalElements()
<< " x " << epetra_A->OperatorDomainMap().NumGlobalElements()
<< "\n ||epetraA||inf = " << epetra_A->NormInf()
<< "\n ||epetra_b||2 = " << epetraNorm2(*epetra_b)
<< "\n ||epetra_x||2 = " << epetraNorm2(*epetra_x)
<< "\n";
RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp( epetra_A );
RCP<Thyra::VectorBase<double> > x = Thyra::create_Vector( epetra_x, A->domain() );
RCP<const Thyra::VectorBase<double> > b = Thyra::create_Vector( epetra_b, A->range() );
linearSolverBuilder.readParameters(out.get());
if(extraParamsFile.length()) {
Teuchos::updateParametersFromXmlFile( "./"+extraParamsFile,
linearSolverBuilder.getNonconstParameterList().ptr() );
}
RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
linearSolverBuilder.createLinearSolveStrategy("");
lowsFactory->setOStream(out);
lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
RCP<Thyra::LinearOpWithSolveBase<double> > lows =
Thyra::linearOpWithSolve(*lowsFactory, A);
Thyra::SolveStatus<double> status =
Thyra::solve<double>(*lows, Thyra::NOTRANS, *b, x.ptr());
*out << "\nSolve status:\n" << status;
linearSolverBuilder.writeParamsFile(*lowsFactory);
x = Teuchos::null;
*out
<< "\nSolution ||epetra_x||2 = " << epetraNorm2(*epetra_x) << "\n";
*out << "\nTesting the solution error ||b-A*x||/||b|| computed through the Epetra objects ...\n";
{
epetra_A->Apply(*epetra_x,epetra_A_x);
epetra_r.Update(-1.0,epetra_A_x,1.0);
}
const double
nrm_r = epetraNorm2(epetra_r),
nrm_b = epetraNorm2(*epetra_b),
rel_err = ( nrm_r / nrm_b );
const bool
passed = (rel_err <= tol);
*out
<< "||b-A*x||/||b|| = " << nrm_r << "/" << nrm_b << " = " << rel_err
<< " < tol = " << tol << " ? " << ( passed ? "passed" : "failed" ) << "\n";
if(!passed) success = false;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success)
if (verbose) {
if(success) *out << "\nCongratulations! All of the tests checked out!\n";
else *out << "\nOh no! At least one of the tests failed!\n";
}
return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
}