Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Thyra_GeneralSolveCriteriaBelosStatusTest_def.hpp
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
45#ifndef THYRA_GENERAL_SOLVE_CRITERIA_BELOS_STATUS_TEST_DEF_HPP
46#define THYRA_GENERAL_SOLVE_CRITERIA_BELOS_STATUS_TEST_DEF_HPP
47
48#include "Thyra_GeneralSolveCriteriaBelosStatusTest.hpp"
49#include "Thyra_VectorSpaceBase.hpp"
50#include "Thyra_MultiVectorBase.hpp"
51#include "Thyra_VectorBase.hpp"
52#include "Thyra_MultiVectorStdOps.hpp"
53#include "Thyra_VectorStdOps.hpp"
54#include "BelosThyraAdapter.hpp"
55#include "Teuchos_DebugDefaultAsserts.hpp"
56#include "Teuchos_VerbosityLevel.hpp"
57#include "Teuchos_as.hpp"
58
59
60namespace Thyra {
61
62
63// Constructors/initializers/accessors
64
65
66template<class Scalar>
68 :convergenceTestFrequency_(-1),
69 compute_x_(false),
70 compute_r_(false),
71 lastCurrIter_(-1),
72 lastRtnStatus_(Belos::Undefined)
73{
75}
76
77
78template<class Scalar>
80 const SolveCriteria<Scalar> &solveCriteria,
81 const int convergenceTestFrequency
82 )
83{
84
85 using Teuchos::as;
86 typedef ScalarTraits<ScalarMag> SMT;
87
88 // Make sure the solve criteria is valid
89
90 TEUCHOS_ASSERT_INEQUALITY(solveCriteria.requestedTol, >=, SMT::zero());
91 TEUCHOS_ASSERT_INEQUALITY(solveCriteria.solveMeasureType.numerator, !=, SOLVE_MEASURE_ONE);
92 TEUCHOS_ASSERT(nonnull(solveCriteria.numeratorReductionFunc) ||
93 nonnull(solveCriteria.denominatorReductionFunc) );
94
95 // Remember the solve criteria sett
96
97 solveCriteria_ = solveCriteria;
98 convergenceTestFrequency_ = convergenceTestFrequency;
99
100 // Determine what needs to be computed on each check
101
102 compute_r_ = solveCriteria.solveMeasureType.contains(SOLVE_MEASURE_NORM_RESIDUAL);
103
104 compute_x_ = (compute_r_ ||
105 solveCriteria.solveMeasureType.contains(SOLVE_MEASURE_NORM_SOLUTION));
106
107}
108
109
110template<class Scalar>
111ArrayView<const typename ScalarTraits<Scalar>::magnitudeType>
113{
114 return lastAchievedTol_;
115}
116
117
118// Overridden public functions from Belos::StatusTest
119
120
121template <class Scalar>
122Belos::StatusType
124 Belos::Iteration<Scalar,MV,OP> *iSolver
125 )
126{
127
128 using Teuchos::null;
129
130#ifdef THYRA_DEBUG
131 TEUCHOS_ASSERT(iSolver);
132#endif
133
134 const int currIter = iSolver->getNumIters();
135
136 if (currIter == 0 || currIter % convergenceTestFrequency_ != 0) {
137 // We will skip this check!
138 return Belos::Undefined;
139 }
140
141 const RCP<FancyOStream> out = this->getOStream();
142 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
143
144 const Belos::LinearProblem<Scalar,MV,OP>& lp = iSolver->getProblem();
145 const int numRhs = lp.getRHS()->domain()->dim();
146
147 // Compute the rhs norm if requested and not already computed
148 if (solveCriteria_.solveMeasureType.contains(SOLVE_MEASURE_NORM_RHS)) {
149 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||b||");
150 }
151
152 // Compute the initial residual norm if requested and not already computed
153 if (solveCriteria_.solveMeasureType.contains(SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
154 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||r0||");
155 }
156
157 // Compute X if requested
158 RCP<MultiVectorBase<Scalar> > X;
159 if (compute_x_) {
160 RCP<MV> X_update = iSolver->getCurrentUpdate();
161 X = lp.updateSolution(X_update);
162 }
163
164 // Compute R if requested
165 RCP<MultiVectorBase<Scalar> > R;
166 if (compute_r_) {
167 R = createMembers(lp.getOperator()->range(), X->domain());
168 lp.computeCurrResVec(&*R, &*X);
169 }
170
171 // Form numerators and denominaotors gN(vN)/gD(vD) for each system
172
173 lastNumerator_.resize(numRhs);
174 lastDenominator_.resize(numRhs);
175
176 for (int j = 0; j < numRhs; ++j) {
177 const RCP<const VectorBase<Scalar> > x_j = (nonnull(X) ? X->col(j) : null);
178 const RCP<const VectorBase<Scalar> > r_j = (nonnull(R) ? R->col(j) : null);
179 lastNumerator_[j] = computeReductionFunctional(
180 solveCriteria_.solveMeasureType.numerator,
181 solveCriteria_.numeratorReductionFunc.ptr(),
182 x_j.ptr(), r_j.ptr() );
183 lastDenominator_[j] = computeReductionFunctional(
184 solveCriteria_.solveMeasureType.denominator,
185 solveCriteria_.denominatorReductionFunc.ptr(),
186 x_j.ptr(), r_j.ptr() );
187 }
188
189 // Determine if convRatio <= requestedTol
190
191 bool systemsAreConverged = true;
192 lastAchievedTol_.resize(numRhs);
193
194 for (int j = 0; j < numRhs; ++j) {
195 const ScalarMag convRatio = lastNumerator_[j] / lastDenominator_[j];
196 lastAchievedTol_[j] = convRatio;
197 const bool sys_converged_j = (convRatio <= solveCriteria_.requestedTol);
198 if (includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
199 printRhsStatus(currIter, j, *out);
200 }
201 if (!sys_converged_j) {
202 systemsAreConverged = false;
203 }
204 }
205
206 lastRtnStatus_ = (systemsAreConverged ? Belos::Passed : Belos::Failed);
207 lastCurrIter_ = currIter;
208
209 return lastRtnStatus_;
210
211}
212
213template <class Scalar>
214Belos::StatusType
216{
217 return lastRtnStatus_;
218}
219
220
221template <class Scalar>
223{
224 r0_nrm_.clear();
225 b_nrm_.clear();
226 lastNumerator_.clear();
227 lastDenominator_.clear();
228 lastAchievedTol_.clear();
229 lastCurrIter_ = -1;
230 lastRtnStatus_ = Belos::Undefined;
231}
232
233
234template <class Scalar>
236 std::ostream& os, int indent
237 ) const
238{
239 const int numRhs = lastNumerator_.size();
240 for (int j = 0; j < numRhs; ++j) {
241 printRhsStatus(lastCurrIter_, j, os, indent);
242 }
243}
244
245
246// private
247
248
249template <class Scalar>
252 ESolveMeasureNormType measureType,
253 const Ptr<const ReductionFunctional<Scalar> > &reductFunc,
254 const Ptr<const VectorBase<Scalar> > &x,
255 const Ptr<const VectorBase<Scalar> > &r
256 ) const
257{
258 typedef ScalarTraits<ScalarMag> SMT;
259 ScalarMag rtn = -SMT::one();
260 Ptr<const VectorBase<Scalar> > v;
261 switch(measureType) {
262 case SOLVE_MEASURE_ONE:
263 rtn = SMT::one();
264 break;
265 case SOLVE_MEASURE_NORM_RESIDUAL:
266 v = r;
267 break;
268 case SOLVE_MEASURE_NORM_SOLUTION:
269 v = x;
270 break;
271 case SOLVE_MEASURE_NORM_INIT_RESIDUAL:
272 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||r0||!)")
273 case SOLVE_MEASURE_NORM_RHS:
274 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Handle ||b||!)");
275 TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT();
276 }
277 if (rtn >= SMT::zero()) {
278 // We already have what we need!
279 }
280 else if (nonnull(v) && rtn < SMT::zero()) {
281 if (nonnull(reductFunc)) {
282 rtn = reductFunc->reduce(*v);
283 }
284 else {
285 rtn = norm(*v);
286 }
287 }
288 TEUCHOS_IF_ELSE_DEBUG_ASSERT();
289 return rtn;
290}
291
292
293template <class Scalar>
294void
296 const int currIter, const int j, std::ostream &out,
297 int indent
298 ) const
299{
300 const ScalarMag convRatio = lastNumerator_[j] / lastDenominator_[j];
301 const bool sys_converged_j = (convRatio <= solveCriteria_.requestedTol);
302 for (int i = 0; i < indent; ++i) { out << " "; }
303 out
304 << "["<<currIter<<"] "
305 << "gN(vN("<<j<<"))/gD(vD("<<j<<")) = "
306 << lastNumerator_[j] << "/" << lastDenominator_[j] << " = "
307 << convRatio << " <= " << solveCriteria_.requestedTol << " : "
308 << (sys_converged_j ? " true" : "false")
309 << "\n";
310}
311
312
313} // namespace Thyra
314
315
316#endif // THYRA_GENERAL_SOLVE_CRITERIA_BELOS_STATUS_TEST_DEF_HPP
Thyra specializations of MultiVecTraits and OperatorTraits.
Subclass of Belos::StatusTest that implements every possible form of SolveCriteria that exists by for...
void setSolveCriteria(const SolveCriteria< Scalar > &solveCriteria, const int convergenceTestFrequency)