Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
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[])
LinearSolverBuilder< double > DefaultLinearSolverBuilder

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 "EpetraExt_readEpetraLinearSystem.h"
#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;
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 {
//
// 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
Epetra_SerialComm comm;
#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::SolveStatus<double> status =
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 );
}