Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
simple_stratimikos_example.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 "EpetraExt_readEpetraLinearSystem.h"
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"
49#ifdef HAVE_MPI
50# include "Epetra_MpiComm.h"
51#else
52# include "Epetra_SerialComm.h"
53#endif
54
55
56namespace {
57
58
59// Helper function to compute a single norm for a vector
60double epetraNorm2( const Epetra_Vector &v )
61{
62 double norm[1] = { -1.0 };
63 v.Norm2(&norm[0]);
64 return norm[0];
65}
66
67
68} // namespace
69
70
71int main(int argc, char* argv[])
72{
73
74 using Teuchos::rcp;
75 using Teuchos::RCP;
76 using Teuchos::CommandLineProcessor;
77 typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
78
79 bool success = true;
80 bool verbose = true;
81
82 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
83
84 Teuchos::RCP<Teuchos::FancyOStream>
85 out = Teuchos::VerboseObjectBase::getDefaultOStream();
86
87 try {
88
89
90 //
91 // A) Program setup code
92 //
93
94 //
95 // Read options from command-line
96 //
97
98 std::string matrixFile = "";
99 std::string extraParamsFile = "";
100 double tol = 1e-5;
101 bool onlyPrintOptions = false;
102 bool printXmlFormat = false;
103 bool showDoc = true;
104
105 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
106
107 CommandLineProcessor clp(false); // Don't throw exceptions
108
109 // Set up command-line options for the linear solver that will be used!
110 linearSolverBuilder.setupCLP(&clp);
111
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." );
122
123 clp.setDocString(
124 "Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.\n"
125 "\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"
128 "\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"
131 );
132
133 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
134 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
135
136 //
137 // Print out the valid options and stop if asked
138 //
139
140 if(onlyPrintOptions) {
141 if(printXmlFormat)
142 Teuchos::writeParameterListToXmlOStream(
143 *linearSolverBuilder.getValidParameters()
144 ,*out
145 );
146 else
147 linearSolverBuilder.getValidParameters()->print(
148 *out,PLPrintOptions().indent(2).showTypes(true).showDoc(showDoc)
149 );
150 return 0;
151 }
152
153
154 //
155 // B) Epetra-specific code that sets up the linear system to be solved
156 //
157 // While the below code reads in the Epetra objects from a file, you can
158 // setup the Epetra objects any way you would like. Note that this next
159 // set of code as nothing to do with Thyra at all, and it should not.
160 //
161
162 *out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
163
164#ifdef HAVE_MPI
165 Epetra_MpiComm comm(MPI_COMM_WORLD);
166#else
167 Epetra_SerialComm comm;
168#endif
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 );
172
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()));
176 epetra_b->Random();
177 }
178
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); // Initial guess is critical!
183 }
184
185 *out << "\nPrinting statistics of the Epetra linear system ...\n";
186
187 *out
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)
194 << "\n";
195
196
197 //
198 // C) The "Glue" code that takes Epetra objects and wraps them as Thyra
199 // objects
200 //
201 // This next set of code wraps the Epetra objects that define the linear
202 // system to be solved as Thyra objects so that they can be passed to the
203 // linear solver.
204 //
205
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() );
209
210 // Note that above Thyra is only interacted with in the most trival of
211 // ways. For most users, Thyra will only be seen as a thin wrapper that
212 // they need know little about in order to wrap their objects in order to
213 // pass them to Thyra-enabled solvers.
214
215
216 //
217 // D) Thyra-specific code for solving the linear system
218 //
219 // Note that this code has no mention of any concrete implementation and
220 // therefore can be used in any use case.
221 //
222
223 // Reading in the solver parameters from the parameters file and/or from
224 // the command line. This was setup by the command-line options
225 // set by the setupCLP(...) function above.
226 linearSolverBuilder.readParameters(out.get());
227
228 // Augment parameters if appropriate
229 if(extraParamsFile.length()) {
230 Teuchos::updateParametersFromXmlFile( "./"+extraParamsFile,
231 linearSolverBuilder.getNonconstParameterList().ptr() );
232 }
233
234 // Create a linear solver factory given information read from the
235 // parameter list.
236 RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
237 linearSolverBuilder.createLinearSolveStrategy("");
238
239 // Setup output stream and the verbosity level
240 lowsFactory->setOStream(out);
241 lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
242
243 // Create a linear solver based on the forward operator A
244 RCP<Thyra::LinearOpWithSolveBase<double> > lows =
245 Thyra::linearOpWithSolve(*lowsFactory, A);
246
247 // Solve the linear system (note: the initial guess in 'x' is critical)
248 Thyra::SolveStatus<double> status =
249 Thyra::solve<double>(*lows, Thyra::NOTRANS, *b, x.ptr());
250 *out << "\nSolve status:\n" << status;
251
252 // Write the linear solver parameters after they were read
253 linearSolverBuilder.writeParamsFile(*lowsFactory);
254
255
256 //
257 // E) Post process the solution and check the error
258 //
259 // Note that the below code is based only on the Epetra objects themselves
260 // and does not in any way depend or interact with any Thyra-based
261 // objects. The point is that most users of Thyra can largely gloss over
262 // the fact that Thyra is really being used for anything.
263 //
264
265 // Wipe out the Thyra wrapper for x to guarantee that the solution will be
266 // written back to epetra_x! At the time of this writting this is not
267 // really needed but the behavior may change at some point so this is a
268 // good idea.
269 x = Teuchos::null;
270
271 *out
272 << "\nSolution ||epetra_x||2 = " << epetraNorm2(*epetra_x) << "\n";
273
274 *out << "\nTesting the solution error ||b-A*x||/||b|| computed through the Epetra objects ...\n";
275
276 // r = b - A*x
277 Epetra_Vector epetra_r(*epetra_b);
278 {
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);
282 }
283
284 const double
285 nrm_r = epetraNorm2(epetra_r),
286 nrm_b = epetraNorm2(*epetra_b),
287 rel_err = ( nrm_r / nrm_b );
288 const bool
289 passed = (rel_err <= tol);
290
291 *out
292 << "||b-A*x||/||b|| = " << nrm_r << "/" << nrm_b << " = " << rel_err
293 << " < tol = " << tol << " ? " << ( passed ? "passed" : "failed" ) << "\n";
294
295 if(!passed) success = false;
296
297 }
298 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success)
299
300 if (verbose) {
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";
303 }
304
305 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
306
307}
int main(int argc, char *argv[])
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...
RCP< const ParameterList > getValidParameters() const
RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > createLinearSolveStrategy(const std::string &linearSolveStrategyName) const
void setupCLP(Teuchos::CommandLineProcessor *clp)
Setup the command-line processor to read in the needed data to extra the parameters from.
void readParameters(std::ostream *out)
Force the parameters to be read from a file and/or an extra XML string.
void writeParamsFile(const Thyra::LinearOpWithSolveFactoryBase< Scalar > &lowsFactory, const std::string &outputXmlFileName="") const
Write the parameters list for a LinearOpWithSolveFactoryBase object to a file after the parameters ar...