Stratimikos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
simple_stratimikos_example.cpp

Simple example program for how the Stratimikos::DefaultLinearSolverBuilder class is used to take a parameter list (e.g. read from a file) and use it so solve a linear system read in as Epetra objects from a file.

Outline

Example program source (with line numbers)

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 
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"
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;
77  typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
78 
79  bool success = true;
80  bool verbose = true;
81 
82  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
83 
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)
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 }
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< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space=Teuchos::null)
static RCP< FancyOStream > getDefaultOStream()
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)
RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > createLinearSolveStrategy(const std::string &linearSolveStrategyName) const
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
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...
void setupCLP(Teuchos::CommandLineProcessor *clp)
Setup the command-line processor to read in the needed data to extra the parameters from...
RCP< const ParameterList > getValidParameters() const

Command-line options for the driver program

Teuchos::GlobalMPISession::GlobalMPISession(): started serial run

Echoing the command-line:

/home/rabartl/PROJECTS/Trilinos.base/BUILDS/CMAKE/SERIAL_DEBUG/packages/stratimikos/example/Stratimikos_simple_stratimikos_example.exe --echo-command-line --help 

Usage: /home/rabartl/PROJECTS/Trilinos.base/BUILDS/CMAKE/SERIAL_DEBUG/packages/stratimikos/example/Stratimikos_simple_stratimikos_example.exe [options]
  options:
  --help                                     Prints this help message
  --pause-for-debugging                      Pauses for user input to allow attaching a debugger
  --echo-command-line                        Echo the command-line but continue as normal
  --linear-solver-params-file        string  Name of an XML file containing parameters for linear solver options to be appended first.
                                             (default: --linear-solver-params-file="")
  --extra-linear-solver-params       string  An XML string containing linear solver parameters to be appended second.
                                             (default: --extra-linear-solver-params="")
  --linear-solver-params-used-file   string  Name of an XML file that can be written with the parameter list after it has been used on completion of this program.
                                             (default: --linear-solver-params-used-file="")
  --matrix-file                      string  Defines the matrix and perhaps the RHS and LHS for a linear system to be solved.
                                             (default: --matrix-file="")
  --tol                              double  Tolerance to check against the scaled residual norm.
                                             (default: --tol=1e-05)
  --only-print-options               bool    Only print options and stop or continue on
  --continue-after-printing-options          (default: --continue-after-printing-options)
  --print-xml-format                 bool    Print the valid options in XML format or in readable format.
  --print-readable-format                    (default: --print-readable-format)
  --show-doc                         bool    Print the valid options with or without documentation.
  --hide-doc                                 (default: --show-doc)

DETAILED DOCUMENTATION:

Simple example for the use of the Stratimikos facade Stratimikos::DefaultLinearSolverBuilder.

To print out just the valid options use --matrix-file="" --only-print-options with --print-xml-format or --print-readable-format.

To solve a linear system from a system read from a file use --matrix-file="SomeFile.mtx" with options read from an XML file using --linear-solver-params-file="SomeFile.xml"


Example input XML file

Here is an example input XML file, called amesos.klu.xml, that specifies that the KLU solver from Amesos be used to solve the linear system:

<ParameterList>
  <Parameter isUsed="true" name="Linear Solver Type" type="string" value="Amesos"/>
  <ParameterList name="Linear Solver Types">
    <ParameterList name="Amesos">
      <ParameterList name="Amesos Settings"/>
      <Parameter name="Solver Type" type="string" value="Klu"/>
    </ParameterList>
  </ParameterList>
  <Parameter name="Preconditioner Type" type="string" value="None"/>
</ParameterList>

Example program output

Sample output looks like when using the above amesos.klu.xml file as input with a test system:

Teuchos::GlobalMPISession::GlobalMPISession(): started serial run

Echoing the command-line:

/home/rabartl/PROJECTS/Trilinos.base/BUILDS/CMAKE/SERIAL_DEBUG/packages/stratimikos/example/Stratimikos_simple_stratimikos_example.exe --echo-command-line --matrix-file=FourByFour.mtx --linear-solver-params-file=amesos.klu.xml 


Reading linear system in Epetra format from the file 'FourByFour.mtx' ...

Printing statistics of the Epetra linear system ...

  Epetra_CrsMatrix epetra_A of dimension 4 x 4
  ||epetraA||inf = 1.643
  ||epetra_b||2 = 1.52593
  ||epetra_x||2 = 0

Reading parameters from XML file "amesos.klu.xml" ...
 
 Solving block system using Amesos solver Amesos_Klu ...
 
 
 Total solve time = 0 sec

Solve status:
 solveStatus = SOLVE_STATUS_CONVERGED
 achievedTol = unknownTolerance()
 message:
  Solver Amesos_Klu converged!
 extraParameters: NONE

Solution ||epetra_x||2 = 1.24984

Testing the solution error ||b-A*x||/||b|| computed through the Epetra objects ...
||b-A*x||/||b|| = 1.11022e-16/1.52593 = 7.27571e-17 < tol = 1e-05 ? passed

Congratulations! All of the tests checked out!

Example program source (without line numbers)

Here is the main driver program itself without line numbers:

// @HEADER
// ***********************************************************************
//
// Stratimikos: Thyra-based strategies for linear solvers
// Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
//
// ***********************************************************************
// @HEADER
#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 {
// Helper function to compute a single norm for a vector
double epetraNorm2( const Epetra_Vector &v )
{
double norm[1] = { -1.0 };
v.Norm2(&norm[0]);
return norm[0];
}
} // namespace
int main(int argc, char* argv[])
{
using Teuchos::rcp;
using Teuchos::RCP;
typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
bool success = true;
bool verbose = true;
Teuchos::GlobalMPISession mpiSession(&argc, &argv);
try {
//
// A) Program setup code
//
//
// Read options from command-line
//
std::string matrixFile = "";
std::string extraParamsFile = "";
double tol = 1e-5;
bool onlyPrintOptions = false;
bool printXmlFormat = false;
bool showDoc = true;
CommandLineProcessor clp(false); // Don't throw exceptions
// Set up command-line options for the linear solver that will be used!
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;
//
// Print out the valid options and stop if asked
//
if(onlyPrintOptions) {
if(printXmlFormat)
Teuchos::writeParameterListToXmlOStream(
*linearSolverBuilder.getValidParameters()
,*out
);
else
linearSolverBuilder.getValidParameters()->print(
*out,PLPrintOptions().indent(2).showTypes(true).showDoc(showDoc)
);
return 0;
}
//
// B) Epetra-specific code that sets up the linear system to be solved
//
// While the below code reads in the Epetra objects from a file, you can
// setup the Epetra objects any way you would like. Note that this next
// set of code as nothing to do with Thyra at all, and it should not.
//
*out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
#ifdef HAVE_MPI
Epetra_MpiComm comm(MPI_COMM_WORLD);
#else
#endif
RCP<Epetra_CrsMatrix> epetra_A;
RCP<Epetra_Vector> epetra_x, epetra_b;
EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, &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 = rcp(new Epetra_Vector(epetra_A->OperatorRangeMap()));
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); // Initial guess is critical!
}
*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";
//
// C) The "Glue" code that takes Epetra objects and wraps them as Thyra
// objects
//
// This next set of code wraps the Epetra objects that define the linear
// system to be solved as Thyra objects so that they can be passed to the
// linear solver.
//
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() );
// Note that above Thyra is only interacted with in the most trival of
// ways. For most users, Thyra will only be seen as a thin wrapper that
// they need know little about in order to wrap their objects in order to
// pass them to Thyra-enabled solvers.
//
// D) Thyra-specific code for solving the linear system
//
// Note that this code has no mention of any concrete implementation and
// therefore can be used in any use case.
//
// Reading in the solver parameters from the parameters file and/or from
// the command line. This was setup by the command-line options
// set by the setupCLP(...) function above.
linearSolverBuilder.readParameters(out.get());
// Augment parameters if appropriate
if(extraParamsFile.length()) {
Teuchos::updateParametersFromXmlFile( "./"+extraParamsFile,
linearSolverBuilder.getNonconstParameterList().ptr() );
}
// Create a linear solver factory given information read from the
// parameter list.
RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
linearSolverBuilder.createLinearSolveStrategy("");
// Setup output stream and the verbosity level
lowsFactory->setOStream(out);
lowsFactory->setVerbLevel(Teuchos::VERB_LOW);
// Create a linear solver based on the forward operator A
RCP<Thyra::LinearOpWithSolveBase<double> > lows =
Thyra::linearOpWithSolve(*lowsFactory, A);
// Solve the linear system (note: the initial guess in 'x' is critical)
Thyra::solve<double>(*lows, Thyra::NOTRANS, *b, x.ptr());
*out << "\nSolve status:\n" << status;
// Write the linear solver parameters after they were read
linearSolverBuilder.writeParamsFile(*lowsFactory);
//
// E) Post process the solution and check the error
//
// Note that the below code is based only on the Epetra objects themselves
// and does not in any way depend or interact with any Thyra-based
// objects. The point is that most users of Thyra can largely gloss over
// the fact that Thyra is really being used for anything.
//
// Wipe out the Thyra wrapper for x to guarantee that the solution will be
// written back to epetra_x! At the time of this writting this is not
// really needed but the behavior may change at some point so this is a
// good idea.
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";
// r = b - A*x
Epetra_Vector epetra_r(*epetra_b);
{
Epetra_Vector epetra_A_x(epetra_A->OperatorRangeMap());
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 );
}

Generated on Mon Jul 4 2022 11:26:35 for Stratimikos by doxygen 1.8.5