Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
61 namespace Belos {
62 
63 
64 template <class ScalarType, class MV, class OP>
65 class LSQRStatusTest: public Belos::StatusTest<ScalarType,MV,OP> {
66 
67 public:
68 
69  // Convenience typedefs
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 
101  Belos::StatusType getStatus() const {return(status_);}
102 
104 
106 
107 
109  void reset();
110 
112  int setCondLim(MagnitudeType condMax) {
113  condMax_ = condMax;
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 
137  MagnitudeType getCondMaxLim() const {return(condMax_);}
138 
140  int getTermIterMax() const {return(term_iter_max_);}
141 
143  MagnitudeType getRelRhsErr() const {return(rel_rhs_err_);}
144 
146  MagnitudeType getMatErr() const {return(rel_mat_err_);}
147 
149  MagnitudeType getMatCondNum() const {return(matCondNum_);}
150 
152  MagnitudeType getMatNorm() const {return(matNorm_);}
153 
155  int getTermIter() const { return term_iter_; }
156 
158  MagnitudeType getResidNorm() const {return(resNorm_);}
159 
161  MagnitudeType getLSResidNorm() const {return(matResNorm_);}
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 
196 private:
197 
199 
200 
202  MagnitudeType condMax_;
203 
205  int term_iter_max_;
206 
208  MagnitudeType rel_rhs_err_;
209 
211  MagnitudeType rel_mat_err_;
212 
214  MagnitudeType rcondMin_;
215 
217  Belos::StatusType status_;
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
221  int term_iter_;
222 
223  // condition number of the operator
224  MagnitudeType matCondNum_;
225 
226  // Frobenius norm of the operator
227  MagnitudeType matNorm_;
228 
229  // residual norm for the linear system
230  MagnitudeType resNorm_;
231 
232  // least squares residual, operator^Transpose * residual
233  MagnitudeType matResNorm_;
234 
236 
237 };
238 
239 template <class ScalarType, class MV, class OP>
241 LSQRStatusTest (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 
258 template <class ScalarType, class MV, class OP>
260 {}
261 
262 template <class ScalarType, class MV, class OP>
264 {
265  status_ = Belos::Undefined;
266 }
267 
268 template <class ScalarType, class MV, class OP>
270 {
273  if (condMax_ > MTzero )
274  {
275  rcondMin_ = MTone / condMax_;
276  }
277  else
278  {
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 
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;
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 
394 template <class ScalarType, class MV, class OP>
395 void 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 
404 template <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 */
Teuchos::ScalarTraits< ScalarType > SCT
Belos concrete class that iterates LSQR.
MagnitudeType getMatCondNum() const
Returns the value of the observed condition number of Abar.
MagnitudeType getMatErr() const
Returns the value of the estimate of the relative error in the data defining A set in the constructor...
int setTermIterMax(int term_iter_max)
ScalarType resid_norm
The current residual norm.
void printStatus(std::ostream &os, Belos::StatusType type) const
Print message for each status specific to this stopping test.
static magnitudeType eps()
void print(std::ostream &os, int indent=0) const
Output formatted description of stopping test to output stream.
MagnitudeType getResidNorm() const
Returns the value of the observed norm of the residual r = b-Ax.
Pure virtual base class for defining the status testing capabilities of Belos.
MagnitudeType getMatNorm() const
Returns the value of the observed (Frobenius) norm of A.
ScalarType sol_norm
An estimate of the norm of the solution.
ScalarType mat_cond_num
An approximation to the condition number of A.
A pure virtual class for defining the status tests for the Belos iterative solvers.
StatusType
Whether the StatusTest wants iteration to stop.
Definition: BelosTypes.hpp:189
std::string description() const
Method to return description of the maximum iteration status test.
Traits class which defines basic operations on multivectors.
Belos::StatusType getStatus() const
Return the result of the most recent CheckStatus call.
ScalarType frob_mat_norm
An approximation to the Frobenius norm of A.
SCT::magnitudeType MagnitudeType
ScalarType mat_resid_norm
An estimate of the norm of A^T*resid.
Belos::StatusType checkStatus(Belos::Iteration< ScalarType, MV, OP > *iSolver)
Check convergence status of the iterative solver: Unconverged, Converged, Failed. ...
Belos::MultiVecTraits< ScalarType, MV > MVT
Belos::StatusType firstCallCheckStatusSetup(Belos::Iteration< ScalarType, MV, OP > *iSolver)
Called in checkStatus exactly once, on the first call to checkStatus.
int setCondLim(MagnitudeType condMax)
Set the tolerances.
MagnitudeType getRelRhsErr() const
Returns the value of the estimate of the relative error in the data defining b set in the constructor...
int getTermIter() const
!Returns the current number of successful iterations from the most recent StatusTest call...
LSQRStatusTest(MagnitudeType condMax=0.0, int term_iter_max=1, MagnitudeType rel_rhs_err=0.0, MagnitudeType rel_mat_err=0.0)
Constructor.
int getTermIterMax() const
Returns the number of successful convergent iterations required set in the constructor.
int setRelMatErr(MagnitudeType rel_mat_err)
MagnitudeType getCondMaxLim() const
Returns the value of the upper limit of the condition number of Abar set in the constructor.
LSQRIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Implementation of the LSQR iteration.
void reset()
Resets the status test to the initial internal state.
ScalarType bnorm
The norm of the RHS vector b.
virtual ~LSQRStatusTest()
Destructor.
#define TEUCHOS_ASSERT(assertion_test)
int setRelRhsErr(MagnitudeType rel_rhs_err)
Structure to contain pointers to LSQRIteration state variables, ...
MagnitudeType getLSResidNorm() const
Returns the value of the observed norm of the Least Squares residual A^T r.

Generated on Fri Aug 14 2020 10:48:33 for Belos by doxygen 1.8.5