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 // Belos: Block Linear Solvers Package
5 //
6 // Copyright 2004-2016 NTESS and the Belos contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 */
11 
12 #ifndef BELOS_LSQR_STATUS_TEST_HPP
13 #define BELOS_LSQR_STATUS_TEST_HPP
14 
20 #include "BelosStatusTest.hpp"
21 #include "BelosLSQRIter.hpp"
22 
29 namespace Belos {
30 
31 
32 template <class ScalarType, class MV, class OP>
33 class LSQRStatusTest: public Belos::StatusTest<ScalarType,MV,OP> {
34 
35 public:
36 
37  // Convenience typedefs
41 
43 
44 
46 
51  LSQRStatusTest( MagnitudeType condMax = 0.0,
52  int term_iter_max = 1,
53  MagnitudeType rel_rhs_err = 0.0,
54  MagnitudeType rel_mat_err = 0.0 );
55 
57  virtual ~LSQRStatusTest();
59 
61 
62 
64 
67 
69  Belos::StatusType getStatus() const {return(status_);}
70 
72 
74 
75 
77  void reset();
78 
80  int setCondLim(MagnitudeType condMax) {
81  condMax_ = condMax;
83  return(0);}
84 
85  int setTermIterMax(int term_iter_max) {
86  term_iter_max_ = term_iter_max;
87  if (term_iter_max_ < 1)
88  term_iter_max_ = 1;
89  return(0);}
90 
91  int setRelRhsErr(MagnitudeType rel_rhs_err) {
92  rel_rhs_err_ = rel_rhs_err;
93  return(0);}
94 
95  int setRelMatErr(MagnitudeType rel_mat_err) {
96  rel_mat_err_ = rel_mat_err;
97  return(0);}
98 
100 
102 
103 
105  MagnitudeType getCondMaxLim() const {return(condMax_);}
106 
108  int getTermIterMax() const {return(term_iter_max_);}
109 
111  MagnitudeType getRelRhsErr() const {return(rel_rhs_err_);}
112 
114  MagnitudeType getMatErr() const {return(rel_mat_err_);}
115 
117  MagnitudeType getMatCondNum() const {return(matCondNum_);}
118 
120  MagnitudeType getMatNorm() const {return(matNorm_);}
121 
123  int getTermIter() const { return term_iter_; }
124 
126  MagnitudeType getResidNorm() const {return(resNorm_);}
127 
129  MagnitudeType getLSResidNorm() const {return(matResNorm_);}
131 
132 
134 
135 
137  void print(std::ostream& os, int indent = 0) const;
138 
140  void printStatus(std::ostream& os, Belos::StatusType type) const;
141 
143 
146 
151 
154 
156  std::string description() const
157  {
158  std::ostringstream oss;
159  oss << "LSQRStatusTest<>: [ limit of condition number = " << condMax_ << " ]";
160  return oss.str();
161  }
163 
164 private:
165 
167 
168 
170  MagnitudeType condMax_;
171 
173  int term_iter_max_;
174 
176  MagnitudeType rel_rhs_err_;
177 
179  MagnitudeType rel_mat_err_;
180 
182  MagnitudeType rcondMin_;
183 
185  Belos::StatusType status_;
186 
187  // term_iter_ records the number of consecutive "successful" iterations.
188  // convergence requires that term_iter_max consecutive iterates satisfy the other convergence tests
189  int term_iter_;
190 
191  // condition number of the operator
192  MagnitudeType matCondNum_;
193 
194  // Frobenius norm of the operator
195  MagnitudeType matNorm_;
196 
197  // residual norm for the linear system
198  MagnitudeType resNorm_;
199 
200  // least squares residual, operator^Transpose * residual
201  MagnitudeType matResNorm_;
202 
204 
205 };
206 
207 template <class ScalarType, class MV, class OP>
209 LSQRStatusTest (MagnitudeType condMax /* = 0 */,
210  int term_iter_max /* = 1 */,
211  MagnitudeType rel_rhs_err /* = 0 */,
212  MagnitudeType rel_mat_err /* = 0 */)
213  : condMax_(condMax),
214  term_iter_max_ (term_iter_max),
215  rel_rhs_err_ (rel_rhs_err),
216  rel_mat_err_ (rel_mat_err),
217  rcondMin_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
218  status_ (Belos::Undefined),
219  term_iter_ (0),
220  matCondNum_ ( Teuchos::ScalarTraits<MagnitudeType>::one() ),
221  matNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
222  resNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
223  matResNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() )
224 {}
225 
226 template <class ScalarType, class MV, class OP>
228 {}
229 
230 template <class ScalarType, class MV, class OP>
232 {
233  status_ = Belos::Undefined;
234 }
235 
236 template <class ScalarType, class MV, class OP>
238 {
241  if (condMax_ > MTzero )
242  {
243  rcondMin_ = MTone / condMax_;
244  }
245  else
246  {
248  }
249 
250  bool termIterFlag = false;
251  LSQRIter<ScalarType,MV,OP>* solver = dynamic_cast< LSQRIter<ScalarType,MV,OP>* > (iSolver);
252  TEUCHOS_ASSERT(solver != NULL);
254  //
255  // LSQR solves a least squares problem. A converged preconditioned residual norm
256  // suffices for convergence, but is not necessary. LSQR sometimes returns a larger
257  // relative residual norm than what would have been returned by a linear solver.
258  // This section evaluates three stopping criteria. In the Solver Manager, this test
259  // is combined with a generic number of iteration test.
260  // If the linear system includes a preconditioner, then the least squares problem
261  // is solved for the preconditioned linear system. Preconditioning changes the least
262  // squares problem (in the sense of changing the norms), and the solution depends
263  // on the preconditioner in this sense.
264  // In the context of Linear Least Squares problems, preconditioning refers
265  // to the regularization matrix. Here the regularization matrix is always a scalar
266  // multiple of the identity (standard form least squres).
267  // The "loss of accuracy" concept is not yet implemented here, becuase it is unclear
268  // what this means for linear least squares. LSQR solves an inconsistent system
269  // in a least-squares sense. "Loss of accuracy" would correspond to
270  // the difference between the preconditioned residual and the unpreconditioned residual.
271  //
272 
273  std::cout << " X " << state.sol_norm
274  << " b-AX " << state.resid_norm
275  << " Atr " << state.mat_resid_norm
276  << " A " << state.frob_mat_norm
277  << " cond " << state.mat_cond_num
278  << " relResNorm " << state.resid_norm/state.bnorm
279  << " LS " << state.mat_resid_norm /( state.resid_norm * state.frob_mat_norm )
280  << std::endl;
281 
283  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
284  ScalarType stop_crit_1 = zero; // b = 0, done
285  if( state.bnorm > zero )
286  {
287  stop_crit_1 = state.resid_norm / state.bnorm;
288  }
289  ScalarType stop_crit_2 = zero;
290  if( state.frob_mat_norm > zero && state.resid_norm > zero )
291  {
292  stop_crit_2 = (state.resid_norm > zero) ? state.mat_resid_norm / (state.frob_mat_norm * state.resid_norm) : zero;
293  }
294  else
295  {
296  if( state.resid_norm == zero )
297  {
298  stop_crit_2 = zero;
299  }
300  else
301  {
302  stop_crit_2 = one; // Initial mat_norm always vanishes
303  }
304  }
305  ScalarType stop_crit_3 = one / state.mat_cond_num;
306  ScalarType resid_tol = rel_rhs_err_ + rel_mat_err_ * state.frob_mat_norm * state.sol_norm / state.bnorm;
308 
309  // The expected use case for our users is that the linear system will almost
310  // always be compatible, but occasionally may not be. However, some users
311  // may use LSQR for more general cases. This is why we include the full
312  // suite of tests, for both compatible and incompatible systems.
313  //
314  // Users will have to be educated that sometimes they will get an answer X
315  // that does _not_ satisfy the linear system AX=B, but _does_ satisfy the
316  // corresponding least-squares problem. Perhaps the solution manager should
317  // provide them with a way to find out.
318 
319  // stop_crit_1 is for compatible linear systems.
320  // stop_crit_2 is for incompatible linear systems.
321  // stop_crit_3 is for either compatible or incompatible linear systems.
322 
323  // Have we met any of the stopping criteria?
324  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()) {
325  termIterFlag = true;
326 
327  if (stop_crit_1 <= resid_tol )
328  std::cout << "Conv: stop_crit_1 " << stop_crit_1 << " resid_tol " << resid_tol << std::endl;
329 
330  if (stop_crit_1 <= resid_tol_mach )
331  std::cout << "Conv: stop_crit_1 " << stop_crit_1 << " resid_tol_mach " << resid_tol_mach << std::endl;
332 
333  if (stop_crit_2 <= rel_mat_err_ )
334  std::cout << "Conv: stop_crit_2 " << stop_crit_2 << " rel_mat_err " << rel_mat_err_ << std::endl;
335 
336  if (stop_crit_2 <= Teuchos::ScalarTraits< MagnitudeType >::eps() )
337  std::cout << "Conv: stop_crit_2 " << stop_crit_2 << " eps " << Teuchos::ScalarTraits< MagnitudeType >::eps() << std::endl;
338 
339  if (stop_crit_3 <= rcondMin_ )
340  std::cout << "Conv: stop_crit_3 " << stop_crit_3 << " rcondMin_ " << rcondMin_ << std::endl;
341 
342  if (stop_crit_3 <= Teuchos::ScalarTraits< MagnitudeType >::eps() )
343  std::cout << "Conv: stop_crit_3 " << stop_crit_3 << " eps " << Teuchos::ScalarTraits< MagnitudeType >::eps() << std::endl;
344  }
345 
346  // update number of consecutive successful iterations
347  if (!termIterFlag) {
348  term_iter_ = 0;
349  } else {
350  term_iter_++;
351  }
352  status_ = (term_iter_ < term_iter_max_) ? Belos::Failed : Belos::Passed;
353 
354  matCondNum_ = state.mat_cond_num; // information that defined convergence
355  matNorm_ = state.frob_mat_norm; // in accessible variables
356  resNorm_ = state.resid_norm;
357  matResNorm_ = state.mat_resid_norm;
358 
359  return status_;
360 }
361 
362 template <class ScalarType, class MV, class OP>
363 void LSQRStatusTest<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
364 {
365  for (int j = 0; j < indent; j++)
366  os << ' ';
367  printStatus(os, status_);
368  os << "limit of condition number = " << condMax_ << std::endl;
369  os << "limit of condition number = " << condMax_ << std::endl;
370 }
371 
372 template <class ScalarType, class MV, class OP>
374 {
375  os << std::left << std::setw(13) << std::setfill('.');
376  switch (type) {
377  case Belos::Passed:
378  os << "Passed";
379  break;
380  case Belos::Failed:
381  os << "Failed";
382  break;
383  case Belos::Undefined:
384  default:
385  os << "Undefined";
386  break;
387  }
388  os << std::left << std::setfill(' ');
389  return;
390 }
391 
392 } // end Belos namespace
393 
394 
395 #endif /* BELOS_LSQR_STATUS_TEST_HPP */
Belos::StatusType checkStatus(Belos::Iteration< ScalarType, MV, OP > *iSolver)
Check convergence status of the iterative solver: Unconverged, Converged, Failed. ...
Teuchos::ScalarTraits< ScalarType > SCT
void print(std::ostream &os, int indent=0) const
Output formatted description of stopping test to output stream.
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.
static magnitudeType eps()
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.
void printStatus(std::ostream &os, Belos::StatusType type) const
Print message for each status specific to this stopping test.
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:157
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.
A Belos::StatusTest class for specifying convergence of LSQR. The outer status tests passes if an inn...
SCT::magnitudeType MagnitudeType
ScalarType mat_resid_norm
An estimate of the norm of A^T*resid.
LSQRStatusTest(MagnitudeType condMax=0.0, int term_iter_max=1, MagnitudeType rel_rhs_err=0.0, MagnitudeType rel_mat_err=0.0)
Constructor.
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.
void reset()
Resets the status test to the initial internal state.
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...
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.
ScalarType bnorm
The norm of the RHS vector b.
#define TEUCHOS_ASSERT(assertion_test)
int setRelRhsErr(MagnitudeType rel_rhs_err)
virtual ~LSQRStatusTest()
Destructor.
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 Mon Jul 15 2024 09:24:24 for Belos by doxygen 1.8.5