Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
stratimikos_example.cpp
Go to the documentation of this file.
1//
2//
3// @HEADER
4// ***********************************************************************
5//
6// Stratimikos: Thyra-based strategies for linear solvers
7// Copyright (2006) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// This library is free software; you can redistribute it and/or modify
13// it under the terms of the GNU Lesser General Public License as
14// published by the Free Software Foundation; either version 2.1 of the
15// License, or (at your option) any later version.
16//
17// This library is distributed in the hope that it will be useful, but
18// WITHOUT ANY WARRANTY; without even the implied warranty of
19// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20// Lesser General Public License for more details.
21//
22// You should have received a copy of the GNU Lesser General Public
23// License along with this library; if not, write to the Free Software
24// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
25// USA
26// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
27//
28// ***********************************************************************
29// @HEADER
30
31#include "Stratimikos_Config.h"
32#include "Teuchos_ConfigDefs.hpp"
33#include "Teuchos_FancyOStream.hpp"
35#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
36#include "Thyra_EpetraLinearOp.hpp"
37// #include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
38#include "EpetraExt_readEpetraLinearSystem.h"
39#include "Teuchos_ParameterList.hpp"
40
41namespace Teuchos { class ParameterList; }
42
43#ifdef HAVE_MPI
44# include "Epetra_MpiComm.h"
45#else
46# include "Epetra_SerialComm.h"
47#endif
48
50 Teuchos::ParameterList *paramList_inout
51 ,const bool dumpAll
52 ,Teuchos::FancyOStream *out
53 )
54{
55
56 using Teuchos::rcp;
57 using Teuchos::RCP;
58 using Teuchos::OSTab;
59 using Teuchos::ParameterList;
60 using Teuchos::getParameter;
61 bool result, success = true;
62
63 try {
64
65 TEUCHOS_TEST_FOR_EXCEPT(!paramList_inout);
66
67 RCP<ParameterList>
68 paramList = rcp(paramList_inout,false);
69
70#if 0
71 if(out) {
72 *out << "\nEchoing input parameters ...\n";
73 paramList->print(*out,1,true,false);
74 }
75#endif
76
77 // Create list of valid parameter sublists
78 Teuchos::ParameterList validParamList("TestSingleStratimikosSolver");
79 validParamList.set("Matrix File","fileName");
80 validParamList.sublist("Linear Solver Builder");
81
82 Teuchos::ParameterList &StratimikosList = paramList_inout->sublist("Linear Solver Builder");
83
84#if 0
85
86 this fails because validParamList has only been set one level deep - see 5-10 lines above
87
88 if(out) *out << "\nValidating top-level input parameters ...\n";
89 StratimikosList.validateParameters(validParamList.sublist("Linear Solver Builder"),0);
90#endif
91
92#if 0
93 paramList->validateParameters(validParamList,0);
94
95 if(out) *out << "\nValidating top-level input parameters ...\n";
96 paramList->validateParameters(StratimikosList,0);
97#endif
98
99 const std::string
100 &matrixFile = getParameter<std::string>(*paramList,"Matrix File");
101
102 if(out) *out << "\nReading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
103
104#ifdef HAVE_MPI
105 Epetra_MpiComm comm(MPI_COMM_WORLD);
106#else
107 Epetra_SerialComm comm;
108#endif
109 Teuchos::RCP<Epetra_CrsMatrix> epetra_A;
110 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
111
112 const int num_vecs = 1 ;
113 const Epetra_Map& A_domainmap = epetra_A->DomainMap() ;
114 Epetra_MultiVector X( A_domainmap, num_vecs ) ;
115 const Epetra_Map& A_rangemap = epetra_A->RangeMap() ;
116 Epetra_MultiVector B( A_rangemap, num_vecs ) ;
117 B.Random();
118
119 int status = simpleStratimikosSolve( *epetra_A, B, &X, &StratimikosList );
120 assert( status==0 ) ;
121
122
123 int nrows =epetra_A->NumGlobalRows();
124 Epetra_MultiVector Residual( B ) ;
125 epetra_A->Apply( X, Residual );
126 Residual.Update( 1.0, B, -1.0 );
127 double ResidualNorm;
128 double BNorm;
129 Residual.Norm2( &ResidualNorm );
130 B.Norm2( &BNorm );
131
132 const Teuchos::ParameterList& LST_PL = StratimikosList.sublist("Linear Solver Types");
133 const Teuchos::ParameterList& Amesos_PL = LST_PL.sublist("Amesos");
134 std::string SolverType = Amesos_PL.get<std::string>("Solver Type");
135 assert( SolverType == "Lapack" ) ; // for now - remove later
136 const Teuchos::ParameterList& AmesosSettings_PL = Amesos_PL.sublist("Amesos Settings");
137 const Teuchos::ParameterList& SolverType_PL = AmesosSettings_PL.sublist(SolverType);
138 double AddToDiag = -13.0;
139 if ( SolverType == "Lapack" ) {
140 AddToDiag = SolverType_PL.get<double>("AddToDiag");
141 }
142 assert( AddToDiag >= 0.0 ) ;
143 double MaxError = 1e-15*BNorm + 10.0 * AddToDiag ;
144 double MinError = 0.02 * ( 1e-15*BNorm + AddToDiag );
145 success = ResidualNorm < nrows * MaxError &&
146 ResidualNorm > MinError ;
147
148#if 0
149 std::cout << __FILE__ << "::" << __LINE__
150 << " ResidualNorm = " << ResidualNorm
151 << " MinError = " << MinError
152 << " MaxError = " << MaxError << std::endl ;
153 std::cout << " B = " ;
154 B.Print( std::cout ) ;
155 std::cout << " X = " ;
156 X.Print( std::cout ) ;
157 std::cout << " epetra_A = " ;
158 epetra_A->Print( std::cout ) ;
159 std::cout << " success = " << success << " ResidualNorm = " << ResidualNorm << std::endl ;
160#endif
161
162 if(false && out) {
163 *out << "\nPrinting the parameter list (showing what was used) ...\n";
164 paramList->print(*out,1,true,true);
165 }
166
167 }
168 catch( const std::exception &excpt ) {
169 std::cerr << "*** Caught standard exception : " << excpt.what() << std::endl;
170 success = false;
171 }
172
173 return success;
174
175}
176
177#include "Teuchos_GlobalMPISession.hpp"
178#include "Teuchos_VerboseObject.hpp"
179#include "Teuchos_CommandLineProcessor.hpp"
180#include "Teuchos_XMLParameterListHelpers.hpp"
181#include "Teuchos_StandardCatchMacros.hpp"
182
183int main(int argc, char* argv[])
184{
185
186 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
187
188 using Teuchos::CommandLineProcessor;
189
190 bool success = true;
191 bool verbose = true;
192
193 Teuchos::RCP<Teuchos::FancyOStream>
194 out = Teuchos::VerboseObjectBase::getDefaultOStream();
195
196 try {
197
198 //
199 // Read options from command-line
200 //
201
202 std::string inputFile = "stratimikos_amesos_lapack.xml";
203 std::string extraParams = "";
204 bool dumpAll = false;
205
206 CommandLineProcessor clp(false); // Don't throw exceptions
207 std::string default_inputFile ;
208 sprintf( &default_inputFile[0], "Input file [%s]", &inputFile[0] );
209 clp.setOption( "input-file", &inputFile, &default_inputFile[0], false );
210 clp.setOption( "extra-params", &extraParams, "Extra parameters overriding the parameters read in from --input-file");
211 clp.setDocString(
212 "Testing program for Trilinos (and non-Trilinos) linear solvers access through Thyra."
213 );
214
215 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
216 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
217
218 Teuchos::ParameterList paramList;
219 if(verbose) *out << "\nReading parameters from XML file \""<<inputFile<<"\" ...\n";
220 Teuchos::updateParametersFromXmlFile(inputFile,&paramList);
221 if(extraParams.length()) {
222 if(verbose) *out << "\nAppending extra parameters from the XML string \""<<extraParams<<"\" ...\n";
223 Teuchos::updateParametersFromXmlString(extraParams,&paramList);
224 }
225
226 success
228 &paramList,dumpAll,verbose?&*out:0
229 );
230
231 }
232 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success)
233
234 if (verbose) {
235 if(success) *out << "\nCongratulations! All of the tests checked out!\n";
236 else *out << "\nOh no! At least one of the tests failed!\n";
237 }
238
239 return ( success ? 0 : 1 );
240}
static bool verbose
Definition Amesos.cpp:67
int simpleStratimikosSolve(Epetra_CrsMatrix const &epetra_A, Epetra_MultiVector const &epetra_B, Epetra_MultiVector *epetra_X, Teuchos::ParameterList *paramList)
int main(int argc, char *argv[])
bool TestSingleStratimikosSolver(Teuchos::ParameterList *paramList_inout, const bool dumpAll, Teuchos::FancyOStream *out)