Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Thyra_Belos_StatusTest_UnitTests.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Stratimikos: Thyra-based strategies for linear solvers
6// Copyright (2006) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44#include "Thyra_BelosLinearOpWithSolveFactory.hpp"
45#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
46#include "Thyra_MultiVectorStdOps.hpp"
47#include "Thyra_VectorBase.hpp"
48#include "Thyra_VectorStdOps.hpp"
49#include "Thyra_EpetraLinearOp.hpp"
50#include "EpetraExt_readEpetraLinearSystem.h"
51#include "Epetra_SerialComm.h"
52#include "Teuchos_XMLParameterListHelpers.hpp"
53#include "Teuchos_toString.hpp"
54
55#include "Teuchos_UnitTestHarness.hpp"
56
57
58namespace Thyra {
59
60
61//
62// Helper code
63//
64
65
66const std::string matrixFileName = "nos1.mtx";
67
68
69RCP<const LinearOpBase<double> > getFwdLinearOp()
70{
71 static RCP<const LinearOpBase<double> > fwdLinearOp;
72 if (is_null(fwdLinearOp)) {
73 Teuchos::RCP<Epetra_CrsMatrix> epetraCrsMatrix;
74 EpetraExt::readEpetraLinearSystem( matrixFileName, Epetra_SerialComm(), &epetraCrsMatrix );
75 fwdLinearOp = epetraLinearOp(epetraCrsMatrix);
76 }
77 return fwdLinearOp;
78}
79
80
82template<class Scalar>
83class MockNormInfReductionFunctional : public ReductionFunctional<Scalar> {
84protected:
85
88
90 virtual typename ScalarTraits<Scalar>::magnitudeType
91 reduceImpl(const VectorBase<Scalar> &v) const
92 { return norm_inf(v); }
93
95 virtual bool isCompatibleImpl( const VectorBase<Scalar> &v ) const
96 { return true; }
97
99
100};
101
102
107template<class Scalar>
108RCP<MockNormInfReductionFunctional<Scalar> >
110{
111 return Teuchos::rcp(new MockNormInfReductionFunctional<Scalar>());
112}
113
114
117template<class Scalar>
118class MockMaxNormInfEpsReductionFunctional : public ReductionFunctional<Scalar> {
119protected:
120
123
125 virtual typename ScalarTraits<Scalar>::magnitudeType
126 reduceImpl(const VectorBase<Scalar> &v) const
127 {
128 typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
129 return std::max(norm_inf(v), ScalarTraits<ScalarMag>::eps());
130 }
131
133 virtual bool isCompatibleImpl( const VectorBase<Scalar> &v ) const
134 { return true; }
135
137
138};
139
140
145template<class Scalar>
146RCP<MockMaxNormInfEpsReductionFunctional<Scalar> >
148{
149 return Teuchos::rcp(new MockMaxNormInfEpsReductionFunctional<Scalar>());
150}
151
152
153template<class Scalar>
155 const SolveCriteria<Scalar> &solveCriteria,
156 const Ptr<RCP<const VectorBase<Scalar> > > &x_out,
157 const Ptr<RCP<const VectorBase<Scalar> > > &r_out,
158 bool &success,
159 FancyOStream &out
160 )
161{
162
163 using Teuchos::describe; using Teuchos::optInArg; using Teuchos::rcpFromRef;
164 using Teuchos::toString;
165
166 typedef ScalarTraits<Scalar> ST;
167
168 // A) Set up the linear system
169
170 const RCP<const LinearOpBase<Scalar> > fwdOp = getFwdLinearOp();
171 out << "\nfwdOp = " << describe(*fwdOp, Teuchos::VERB_MEDIUM) << "\n";
172 const RCP<VectorBase<Scalar> > b = createMember(fwdOp->range());
173 V_S(b.ptr(), ST::one());
174 const RCP<VectorBase<Scalar> > x = createMember(fwdOp->domain());
175
176 // B) Print out the specialized SolveCriteria object
177
178 out << "\nsolveCriteria:\n" << solveCriteria;
179
180 // ToDo: Fill in the rest of the fields!
181
182 // C) Solve the system with the given SolveCriteria object
183
184 const int convergenceTestFrequency = 10;
185
186 const RCP<ParameterList> pl = Teuchos::getParametersFromXmlString(
187 "<ParameterList name=\"Belos\">"
188 " <Parameter name=\"Solver Type\" type=\"string\" value=\"Pseudo Block GMRES\"/>"
189 " <Parameter name=\"Convergence Test Frequency\" type=\"int\" value=\""+toString(convergenceTestFrequency)+"\"/>"
190 " <ParameterList name=\"Solver Types\">"
191 " <ParameterList name=\"Pseudo Block GMRES\">"
192 " <Parameter name=\"Block Size\" type=\"int\" value=\"1\"/>"
193 " <Parameter name=\"Convergence Tolerance\" type=\"double\" value=\"1e-13\"/>"
194 " <Parameter name=\"Output Frequency\" type=\"int\" value=\""+toString(convergenceTestFrequency)+"\"/>"
195 " <Parameter name=\"Show Maximum Residual Norm Only\" type=\"bool\" value=\"1\"/>"
196 " <Parameter name=\"Maximum Iterations\" type=\"int\" value=\"400\"/>"
197 " <Parameter name=\"Verbosity\" type=\"int\" value=\"1\"/>"
198 " </ParameterList>"
199 " </ParameterList>"
200 "</ParameterList>"
201 );
202 out << "\n\npl:\n" << *pl;
203
205 lowsFactory.setParameterList(pl);
206 lowsFactory.setOStream(rcpFromRef(out));
207 //lowsFactory.setVerbLevel(Teuchos::VERB_MEDIUM);
208 lowsFactory.setVerbLevel(Teuchos::VERB_HIGH);
209
210 // NOTE: To get Belos ouptut to be quite, turn down the Belos "Verbosity"
211 // option above. To get just the StatusTest VerboseObject output, turn up
212 // lowsFactory output level.
213
214
215 const RCP<LinearOpWithSolveBase<Scalar> > lows = linearOpWithSolve<Scalar>(
216 lowsFactory, fwdOp);
217
218 V_S(x.ptr(), ST::zero());
219 SolveStatus<Scalar> solveStatus = solve<Scalar>(*lows, NOTRANS, *b, x.ptr(),
220 optInArg(solveCriteria));
221 out << "\nsolveStatus:\n" << solveStatus;
222
223 TEST_COMPARE( solveStatus.achievedTol, <=, solveCriteria.requestedTol );
224
225 // D) Compute the actual residual and return x and r
226
227 const RCP<VectorBase<Scalar> > r = b->clone_v();
228 fwdOp->apply(NOTRANS, *x, r.ptr(), ST::one(), -ST::one());
229
230 *x_out = x;
231 *r_out = r;
232
233}
234
235
236
237//
238// GeneralSolveCriteriaBelosStatusTest Unit Tests
239//
240
241
243{
244
245 using Teuchos::outArg;
246
247 typedef double Scalar;
248 typedef ScalarTraits<Scalar> ST;
249 typedef ST::magnitudeType ScalarMag;
250
251 SolveCriteria<Scalar> solveCriteria;
252 solveCriteria.solveMeasureType.numerator = SOLVE_MEASURE_NORM_RESIDUAL;
253 solveCriteria.numeratorReductionFunc = createMockNormReductionFunctional<Scalar>();
254 solveCriteria.solveMeasureType.denominator = SOLVE_MEASURE_NORM_SOLUTION;
255 solveCriteria.denominatorReductionFunc = createMockMaxNormInfEpsReductionFunctional<Scalar>();
256 solveCriteria.requestedTol = 0.9;
257
258 RCP<const VectorBase<Scalar> > x, r;
259 runGeneralSolveCriteriaBelosStatusTestCase(solveCriteria, outArg(x), outArg(r),
260 success, out);
261
262 out << "\nChecking convergence ...\n\n";
263
264 const ScalarMag r_nrm_inf = norm_inf(*r);
265 const ScalarMag x_nrm_inf = norm_inf(*x);
266
267 out << "||r||inf = " << r_nrm_inf << "\n";
268 out << "||x||inf = " << x_nrm_inf << "\n";
269
270 TEST_COMPARE( r_nrm_inf / x_nrm_inf, <=, solveCriteria.requestedTol );
271
272}
273
274
276{
277
278 using Teuchos::outArg;
279
280 typedef double Scalar;
281 typedef ScalarTraits<Scalar> ST;
282 typedef ST::magnitudeType ScalarMag;
283
284 SolveCriteria<Scalar> solveCriteria;
285 solveCriteria.solveMeasureType.numerator = SOLVE_MEASURE_NORM_RESIDUAL;
286 solveCriteria.numeratorReductionFunc = createMockNormReductionFunctional<Scalar>();
287 solveCriteria.solveMeasureType.denominator = SOLVE_MEASURE_ONE;
288 solveCriteria.requestedTol = 0.9;
289
290 RCP<const VectorBase<Scalar> > x, r;
291 runGeneralSolveCriteriaBelosStatusTestCase(solveCriteria, outArg(x), outArg(r),
292 success, out);
293
294 out << "\nChecking convergence ...\n\n";
295
296 const ScalarMag r_nrm_inf = norm_inf(*r);
297 const ScalarMag x_nrm_inf = norm_inf(*x);
298
299 out << "||r||inf = " << r_nrm_inf << "\n";
300 out << "||x||inf = " << x_nrm_inf << "\n";
301
302 TEST_COMPARE( r_nrm_inf, <=, solveCriteria.requestedTol );
303
304}
305
306
307
308} // namespace Thyra
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Subclass of Belos::StatusTest that implements every possible form of SolveCriteria that exists by for...
Mock max(NormInf, eps) ReductionFunctional subclass used for unit testing.
virtual bool isCompatibleImpl(const VectorBase< Scalar > &v) const
virtual ScalarTraits< Scalar >::magnitudeType reduceImpl(const VectorBase< Scalar > &v) const
Mock NormInf ReductionFunctional subclass used for unit testing.
RCP< MockNormInfReductionFunctional< Scalar > > createMockNormReductionFunctional()
Non-member constructor.
RCP< MockMaxNormInfEpsReductionFunctional< Scalar > > createMockMaxNormInfEpsReductionFunctional()
Non-member constructor.
virtual ScalarTraits< Scalar >::magnitudeType reduceImpl(const VectorBase< Scalar > &v) const
virtual bool isCompatibleImpl(const VectorBase< Scalar > &v) const
NOTRANS
TEUCHOS_UNIT_TEST(GeneralSolveCriteriaBelosStatusTest, norm_inf_r_over_norm_inf_r0)
RCP< const LinearOpBase< double > > getFwdLinearOp()
void runGeneralSolveCriteriaBelosStatusTestCase(const SolveCriteria< Scalar > &solveCriteria, const Ptr< RCP< const VectorBase< Scalar > > > &x_out, const Ptr< RCP< const VectorBase< Scalar > > > &r_out, bool &success, FancyOStream &out)
const std::string matrixFileName