Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
49 #ifdef HAVE_MPI
50 # include "Epetra_MpiComm.h"
51 #else
52 # include "Epetra_SerialComm.h"
53 #endif
54 
55 
56 namespace {
57 
58 
59 // Helper function to compute a single norm for a vector
60 double 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 
71 int main(int argc, char* argv[])
72 {
73 
74  using Teuchos::rcp;
75  using Teuchos::RCP;
78 
79  bool success = true;
80  bool verbose = true;
81 
82  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
83 
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[])
void setupCLP(Teuchos::CommandLineProcessor *clp)
Setup the command-line processor to read in the needed data to extra the parameters from...
T * get() const
void readParameters(std::ostream *out)
Force the parameters to be read from a file and/or an extra XML string.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
RCP< Thyra::LinearOpWithSolveFactoryBase< double > > createLinearSolveStrategy(const std::string &linearSolveStrategyName) const
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space=Teuchos::null)
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
void writeParamsFile(const Thyra::LinearOpWithSolveFactoryBase< double > &lowsFactory, const std::string &outputXmlFileName="") const
Write the parameters list for a LinearOpWithSolveFactoryBase object to a file after the parameters ar...
ParameterList::PrintOptions PLPrintOptions
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...