Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosLSQRStatusTest.hpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ************************************************************************
4//
5// Belos: Block Linear Solvers Package
6// Copyright 2004 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42*/
43
44#ifndef BELOS_LSQR_STATUS_TEST_HPP
45#define BELOS_LSQR_STATUS_TEST_HPP
46
52#include "BelosStatusTest.hpp"
53#include "BelosLSQRIter.hpp"
54
61namespace Belos {
62
63
64template <class ScalarType, class MV, class OP>
65class LSQRStatusTest: public Belos::StatusTest<ScalarType,MV,OP> {
66
67public:
68
69 // Convenience typedefs
70 typedef Teuchos::ScalarTraits<ScalarType> SCT;
71 typedef typename SCT::magnitudeType MagnitudeType;
73
75
76
78
83 LSQRStatusTest( MagnitudeType condMax = 0.0,
84 int term_iter_max = 1,
85 MagnitudeType rel_rhs_err = 0.0,
86 MagnitudeType rel_mat_err = 0.0 );
87
89 virtual ~LSQRStatusTest();
91
93
94
96
99
102
104
106
107
109 void reset();
110
113 condMax_ = condMax;
114 rcondMin_ = (condMax > 0) ? (Teuchos::ScalarTraits< MagnitudeType >::one() / condMax) : Teuchos::ScalarTraits< MagnitudeType >::eps();
115 return(0);}
116
117 int setTermIterMax(int term_iter_max) {
118 term_iter_max_ = term_iter_max;
119 if (term_iter_max_ < 1)
120 term_iter_max_ = 1;
121 return(0);}
122
123 int setRelRhsErr(MagnitudeType rel_rhs_err) {
124 rel_rhs_err_ = rel_rhs_err;
125 return(0);}
126
127 int setRelMatErr(MagnitudeType rel_mat_err) {
128 rel_mat_err_ = rel_mat_err;
129 return(0);}
130
132
134
135
138
140 int getTermIterMax() const {return(term_iter_max_);}
141
144
147
150
153
155 int getTermIter() const { return term_iter_; }
156
159
163
164
166
167
169 void print(std::ostream& os, int indent = 0) const;
170
172 void printStatus(std::ostream& os, Belos::StatusType type) const;
173
175
178
183
186
188 std::string description() const
189 {
190 std::ostringstream oss;
191 oss << "LSQRStatusTest<>: [ limit of condition number = " << condMax_ << " ]";
192 return oss.str();
193 }
195
196private:
197
199
200
203
206
209
212
215
218
219 // term_iter_ records the number of consecutive "successful" iterations.
220 // convergence requires that term_iter_max consecutive iterates satisfy the other convergence tests
222
223 // condition number of the operator
225
226 // Frobenius norm of the operator
228
229 // residual norm for the linear system
231
232 // least squares residual, operator^Transpose * residual
234
236
237};
238
239template <class ScalarType, class MV, class OP>
241LSQRStatusTest (MagnitudeType condMax /* = 0 */,
242 int term_iter_max /* = 1 */,
243 MagnitudeType rel_rhs_err /* = 0 */,
244 MagnitudeType rel_mat_err /* = 0 */)
245 : condMax_(condMax),
246 term_iter_max_ (term_iter_max),
247 rel_rhs_err_ (rel_rhs_err),
248 rel_mat_err_ (rel_mat_err),
249 rcondMin_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
250 status_ (Belos::Undefined),
251 term_iter_ (0),
252 matCondNum_ ( Teuchos::ScalarTraits<MagnitudeType>::one() ),
253 matNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
254 resNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
255 matResNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() )
256{}
257
258template <class ScalarType, class MV, class OP>
261
262template <class ScalarType, class MV, class OP>
267
268template <class ScalarType, class MV, class OP>
270{
271 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
272 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
273 if (condMax_ > MTzero )
274 {
275 rcondMin_ = MTone / condMax_;
276 }
277 else
278 {
279 rcondMin_ = Teuchos::ScalarTraits< MagnitudeType >::eps();
280 }
281
282 bool termIterFlag = false;
283 LSQRIter<ScalarType,MV,OP>* solver = dynamic_cast< LSQRIter<ScalarType,MV,OP>* > (iSolver);
284 TEUCHOS_ASSERT(solver != NULL);
286 //
287 // LSQR solves a least squares problem. A converged preconditioned residual norm
288 // suffices for convergence, but is not necessary. LSQR sometimes returns a larger
289 // relative residual norm than what would have been returned by a linear solver.
290 // This section evaluates three stopping criteria. In the Solver Manager, this test
291 // is combined with a generic number of iteration test.
292 // If the linear system includes a preconditioner, then the least squares problem
293 // is solved for the preconditioned linear system. Preconditioning changes the least
294 // squares problem (in the sense of changing the norms), and the solution depends
295 // on the preconditioner in this sense.
296 // In the context of Linear Least Squares problems, preconditioning refers
297 // to the regularization matrix. Here the regularization matrix is always a scalar
298 // multiple of the identity (standard form least squres).
299 // The "loss of accuracy" concept is not yet implemented here, becuase it is unclear
300 // what this means for linear least squares. LSQR solves an inconsistent system
301 // in a least-squares sense. "Loss of accuracy" would correspond to
302 // the difference between the preconditioned residual and the unpreconditioned residual.
303 //
304
305 std::cout << " X " << state.sol_norm
306 << " b-AX " << state.resid_norm
307 << " Atr " << state.mat_resid_norm
308 << " A " << state.frob_mat_norm
309 << " cond " << state.mat_cond_num
310 << " relResNorm " << state.resid_norm/state.bnorm
311 << " LS " << state.mat_resid_norm /( state.resid_norm * state.frob_mat_norm )
312 << std::endl;
313
314 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
315 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
316 ScalarType stop_crit_1 = zero; // b = 0, done
317 if( state.bnorm > zero )
318 {
319 stop_crit_1 = state.resid_norm / state.bnorm;
320 }
321 ScalarType stop_crit_2 = zero;
322 if( state.frob_mat_norm > zero && state.resid_norm > zero )
323 {
324 stop_crit_2 = (state.resid_norm > zero) ? state.mat_resid_norm / (state.frob_mat_norm * state.resid_norm) : zero;
325 }
326 else
327 {
328 if( state.resid_norm == zero )
329 {
330 stop_crit_2 = zero;
331 }
332 else
333 {
334 stop_crit_2 = one; // Initial mat_norm always vanishes
335 }
336 }
337 ScalarType stop_crit_3 = one / state.mat_cond_num;
338 ScalarType resid_tol = rel_rhs_err_ + rel_mat_err_ * state.frob_mat_norm * state.sol_norm / state.bnorm;
339 ScalarType resid_tol_mach = Teuchos::ScalarTraits< MagnitudeType >::eps() + Teuchos::ScalarTraits< MagnitudeType >::eps() * state.frob_mat_norm * state.sol_norm / state.bnorm;
340
341 // The expected use case for our users is that the linear system will almost
342 // always be compatible, but occasionally may not be. However, some users
343 // may use LSQR for more general cases. This is why we include the full
344 // suite of tests, for both compatible and incompatible systems.
345 //
346 // Users will have to be educated that sometimes they will get an answer X
347 // that does _not_ satisfy the linear system AX=B, but _does_ satisfy the
348 // corresponding least-squares problem. Perhaps the solution manager should
349 // provide them with a way to find out.
350
351 // stop_crit_1 is for compatible linear systems.
352 // stop_crit_2 is for incompatible linear systems.
353 // stop_crit_3 is for either compatible or incompatible linear systems.
354
355 // Have we met any of the stopping criteria?
356 if (stop_crit_1 <= resid_tol || stop_crit_2 <= rel_mat_err_ || stop_crit_3 <= rcondMin_ || stop_crit_1 <= resid_tol_mach || stop_crit_2 <= Teuchos::ScalarTraits< MagnitudeType >::eps() || stop_crit_3 <= Teuchos::ScalarTraits< MagnitudeType >::eps()) {
357 termIterFlag = true;
358
359 if (stop_crit_1 <= resid_tol )
360 std::cout << "Conv: stop_crit_1 " << stop_crit_1 << " resid_tol " << resid_tol << std::endl;
361
362 if (stop_crit_1 <= resid_tol_mach )
363 std::cout << "Conv: stop_crit_1 " << stop_crit_1 << " resid_tol_mach " << resid_tol_mach << std::endl;
364
365 if (stop_crit_2 <= rel_mat_err_ )
366 std::cout << "Conv: stop_crit_2 " << stop_crit_2 << " rel_mat_err " << rel_mat_err_ << std::endl;
367
368 if (stop_crit_2 <= Teuchos::ScalarTraits< MagnitudeType >::eps() )
369 std::cout << "Conv: stop_crit_2 " << stop_crit_2 << " eps " << Teuchos::ScalarTraits< MagnitudeType >::eps() << std::endl;
370
371 if (stop_crit_3 <= rcondMin_ )
372 std::cout << "Conv: stop_crit_3 " << stop_crit_3 << " rcondMin_ " << rcondMin_ << std::endl;
373
374 if (stop_crit_3 <= Teuchos::ScalarTraits< MagnitudeType >::eps() )
375 std::cout << "Conv: stop_crit_3 " << stop_crit_3 << " eps " << Teuchos::ScalarTraits< MagnitudeType >::eps() << std::endl;
376 }
377
378 // update number of consecutive successful iterations
379 if (!termIterFlag) {
380 term_iter_ = 0;
381 } else {
382 term_iter_++;
383 }
384 status_ = (term_iter_ < term_iter_max_) ? Belos::Failed : Belos::Passed;
385
386 matCondNum_ = state.mat_cond_num; // information that defined convergence
387 matNorm_ = state.frob_mat_norm; // in accessible variables
388 resNorm_ = state.resid_norm;
389 matResNorm_ = state.mat_resid_norm;
390
391 return status_;
392}
393
394template <class ScalarType, class MV, class OP>
395void LSQRStatusTest<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
396{
397 for (int j = 0; j < indent; j++)
398 os << ' ';
399 printStatus(os, status_);
400 os << "limit of condition number = " << condMax_ << std::endl;
401 os << "limit of condition number = " << condMax_ << std::endl;
402}
403
404template <class ScalarType, class MV, class OP>
406{
407 os << std::left << std::setw(13) << std::setfill('.');
408 switch (type) {
409 case Belos::Passed:
410 os << "Passed";
411 break;
412 case Belos::Failed:
413 os << "Failed";
414 break;
415 case Belos::Undefined:
416 default:
417 os << "Undefined";
418 break;
419 }
420 os << std::left << std::setfill(' ');
421 return;
422}
423
424} // end Belos namespace
425
426
427#endif /* BELOS_LSQR_STATUS_TEST_HPP */
Belos concrete class that iterates LSQR.
Pure virtual base class for defining the status testing capabilities of Belos.
Implementation of the LSQR iteration.
LSQRIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
int setRelMatErr(MagnitudeType rel_mat_err)
int setRelRhsErr(MagnitudeType rel_rhs_err)
void print(std::ostream &os, int indent=0) const
Output formatted description of stopping test to output stream.
Belos::MultiVecTraits< ScalarType, MV > MVT
Teuchos::ScalarTraits< ScalarType > SCT
int getTermIter() const
!Returns the current number of successful iterations from the most recent StatusTest call.
Belos::StatusType status_
Status.
MagnitudeType getMatNorm() const
Returns the value of the observed (Frobenius) norm of A.
LSQRStatusTest(MagnitudeType condMax=0.0, int term_iter_max=1, MagnitudeType rel_rhs_err=0.0, MagnitudeType rel_mat_err=0.0)
Constructor.
int setCondLim(MagnitudeType condMax)
Set the tolerances.
MagnitudeType condMax_
Upper limit of condition number of Abar.
void printStatus(std::ostream &os, Belos::StatusType type) const
Print message for each status specific to this stopping test.
Belos::StatusType firstCallCheckStatusSetup(Belos::Iteration< ScalarType, MV, OP > *iSolver)
Called in checkStatus exactly once, on the first call to checkStatus.
int getTermIterMax() const
Returns the number of successful convergent iterations required set in the constructor.
void reset()
Resets the status test to the initial internal state.
MagnitudeType rcondMin_
One of the tolerances defining convergence, the reciprocal of condMax_ or, if that is zero,...
MagnitudeType getMatErr() const
Returns the value of the estimate of the relative error in the data defining A set in the constructor...
SCT::magnitudeType MagnitudeType
MagnitudeType rel_mat_err_
Error in data defining A.
MagnitudeType getResidNorm() const
Returns the value of the observed norm of the residual r = b-Ax.
MagnitudeType getCondMaxLim() const
Returns the value of the upper limit of the condition number of Abar set in the constructor.
Belos::StatusType getStatus() const
Return the result of the most recent CheckStatus call.
MagnitudeType getRelRhsErr() const
Returns the value of the estimate of the relative error in the data defining b set in the constructor...
MagnitudeType getLSResidNorm() const
Returns the value of the observed norm of the Least Squares residual A^T r.
int setTermIterMax(int term_iter_max)
Belos::StatusType checkStatus(Belos::Iteration< ScalarType, MV, OP > *iSolver)
Check convergence status of the iterative solver: Unconverged, Converged, Failed.
std::string description() const
Method to return description of the maximum iteration status test
virtual ~LSQRStatusTest()
Destructor.
MagnitudeType getMatCondNum() const
Returns the value of the observed condition number of Abar.
MagnitudeType rel_rhs_err_
Error in data defining b.
int term_iter_max_
How many iterations in a row a passing test for convergence is required.
Traits class which defines basic operations on multivectors.
A pure virtual class for defining the status tests for the Belos iterative solvers.
StatusType
Whether the StatusTest wants iteration to stop.
Structure to contain pointers to LSQRIteration state variables, ...
ScalarType sol_norm
An estimate of the norm of the solution.
ScalarType frob_mat_norm
An approximation to the Frobenius norm of A.
ScalarType mat_resid_norm
An estimate of the norm of A^T*resid.
ScalarType bnorm
The norm of the RHS vector b.
ScalarType resid_norm
The current residual norm.
ScalarType mat_cond_num
An approximation to the condition number of A.