Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosCGIter.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_CG_ITER_HPP
43 #define BELOS_CG_ITER_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 #include "BelosCGIteration.hpp"
52 
53 #include "BelosLinearProblem.hpp"
54 #include "BelosOutputManager.hpp"
55 #include "BelosStatusTest.hpp"
56 #include "BelosOperatorTraits.hpp"
57 #include "BelosMultiVecTraits.hpp"
58 
61 #include "Teuchos_ScalarTraits.hpp"
63 #include "Teuchos_TimeMonitor.hpp"
64 
75 namespace Belos {
76 
77 template<class ScalarType, class MV, class OP>
78 class CGIter : virtual public CGIteration<ScalarType,MV,OP> {
79 
80  public:
81 
82  //
83  // Convenience typedefs
84  //
89 
91 
92 
99  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
101  Teuchos::ParameterList &params );
102 
104  virtual ~CGIter() {};
106 
107 
109 
110 
123  void iterate();
124 
140 
144  void initialize()
145  {
147  initializeCG(empty);
148  }
149 
158  state.R = R_;
159  state.P = P_;
160  state.AP = AP_;
161  state.Z = Z_;
162  return state;
163  }
164 
166 
167 
169 
170 
172  int getNumIters() const { return iter_; }
173 
175  void resetNumIters( int iter = 0 ) { iter_ = iter; }
176 
179  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
180 
182 
184  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
185 
187 
189 
190 
192  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
193 
195  int getBlockSize() const { return 1; }
196 
198  void setBlockSize(int blockSize) {
199  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
200  "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
201  }
202 
204  bool isInitialized() { return initialized_; }
205 
206 
208  void setDoCondEst(bool val) {
209  if (numEntriesForCondEst_) doCondEst_=val;
210  }
211 
214  // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
215  // getDiag() didn't actually throw for me in that case, but why
216  // not be cautious?
217  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
218  if (static_cast<size_type> (iter_) >= diag_.size ()) {
219  return diag_ ();
220  } else {
221  return diag_ (0, iter_);
222  }
223  }
224 
227  // NOTE (mfh 30 Jul 2015) The implementation as I found it
228  // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
229  // debug mode) when the maximum number of iterations has been
230  // reached, because iter_ == offdiag_.size() in that case. The
231  // new logic fixes this.
232  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
233  if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
234  return offdiag_ ();
235  } else {
236  return offdiag_ (0, iter_);
237  }
238  }
239 
241 
242  private:
243 
244  //
245  // Internal methods
246  //
248  void setStateSize();
249 
250  //
251  // Classes inputed through constructor that define the linear problem to be solved.
252  //
256 
257  //
258  // Current solver state
259  //
260  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
261  // is capable of running; _initialize is controlled by the initialize() member method
262  // For the implications of the state of initialized_, please see documentation for initialize()
263  bool initialized_;
264 
265  // stateStorageInitialized_ specifies that the state storage has been initialized.
266  // This initialization may be postponed if the linear problem was generated without
267  // the right-hand side or solution vectors.
268  bool stateStorageInitialized_;
269 
270  // Current number of iterations performed.
271  int iter_;
272 
273  // Assert that the matrix is positive definite
274  bool assertPositiveDefiniteness_;
275 
276  // Tridiagonal system for condition estimation (if needed)
277  Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
278  int numEntriesForCondEst_;
279  bool doCondEst_;
280 
281 
282 
283  //
284  // State Storage
285  //
286  // Residual
287  Teuchos::RCP<MV> R_;
288  //
289  // Preconditioned residual
290  Teuchos::RCP<MV> Z_;
291  //
292  // Direction vector
293  Teuchos::RCP<MV> P_;
294  //
295  // Operator applied to direction vector
296  Teuchos::RCP<MV> AP_;
297 
298 };
299 
301  // Constructor.
302  template<class ScalarType, class MV, class OP>
304  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
306  Teuchos::ParameterList &params ):
307  lp_(problem),
308  om_(printer),
309  stest_(tester),
310  initialized_(false),
311  stateStorageInitialized_(false),
312  iter_(0),
313  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
314  numEntriesForCondEst_(params.get("Max Size For Condest",0) ),
315  doCondEst_(false)
316  {
317  }
318 
320  // Setup the state storage.
321  template <class ScalarType, class MV, class OP>
323  {
324  if (!stateStorageInitialized_) {
325 
326  // Check if there is any multivector to clone from.
327  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
328  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
329  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
330  stateStorageInitialized_ = false;
331  return;
332  }
333  else {
334 
335  // Initialize the state storage
336  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
337  if (R_ == Teuchos::null) {
338  // Get the multivector that is not null.
339  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
340  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
341  "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
342  R_ = MVT::Clone( *tmp, 1 );
343  Z_ = MVT::Clone( *tmp, 1 );
344  P_ = MVT::Clone( *tmp, 1 );
345  AP_ = MVT::Clone( *tmp, 1 );
346  }
347 
348  // Tracking information for condition number estimation
349  if(numEntriesForCondEst_ > 0) {
350  diag_.resize(numEntriesForCondEst_);
351  offdiag_.resize(numEntriesForCondEst_-1);
352  }
353 
354  // State storage has now been initialized.
355  stateStorageInitialized_ = true;
356  }
357  }
358  }
359 
360 
362  // Initialize this iteration object
363  template <class ScalarType, class MV, class OP>
365  {
366  // Initialize the state storage if it isn't already.
367  if (!stateStorageInitialized_)
368  setStateSize();
369 
370  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
371  "Belos::CGIter::initialize(): Cannot initialize state storage!");
372 
373  // NOTE: In CGIter R_, the initial residual, is required!!!
374  //
375  std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
376 
377  // Create convenience variables for zero and one.
378  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
379  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
380 
381  if (newstate.R != Teuchos::null) {
382 
383  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
384  std::invalid_argument, errstr );
385  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
386  std::invalid_argument, errstr );
387 
388  // Copy basis vectors from newstate into V
389  if (newstate.R != R_) {
390  // copy over the initial residual (unpreconditioned).
391  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
392  }
393 
394  // Compute initial direction vectors
395  // Initially, they are set to the preconditioned residuals
396  //
397  if ( lp_->getLeftPrec() != Teuchos::null ) {
398  lp_->applyLeftPrec( *R_, *Z_ );
399  if ( lp_->getRightPrec() != Teuchos::null ) {
400  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
401  lp_->applyRightPrec( *Z_, *tmp );
402  Z_ = tmp;
403  }
404  }
405  else if ( lp_->getRightPrec() != Teuchos::null ) {
406  lp_->applyRightPrec( *R_, *Z_ );
407  }
408  else {
409  Z_ = R_;
410  }
411  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
412  }
413  else {
414 
415  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
416  "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
417  }
418 
419  // The solver is initialized
420  initialized_ = true;
421  }
422 
423 
425  // Iterate until the status test informs us we should stop.
426  template <class ScalarType, class MV, class OP>
428  {
429  //
430  // Allocate/initialize data structures
431  //
432  if (initialized_ == false) {
433  initialize();
434  }
435 
436  // Allocate memory for scalars.
439  Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
440 
441  // Create convenience variables for zero and one.
442  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
444 
445  // Scalars for condition estimation (if needed) - These will always use entry zero, for convenience
446  ScalarType pAp_old = one, beta_old = one, rHz_old2 = one;
447 
448  // Get the current solution vector.
449  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
450 
451  // Check that the current solution vector only has one column.
452  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
453  "Belos::CGIter::iterate(): current linear system has more than one vector!" );
454 
455  // Compute first <r,z> a.k.a. rHz
456  MVT::MvTransMv( one, *R_, *Z_, rHz );
457 
459  // Iterate until the status test tells us to stop.
460  //
461  while (stest_->checkStatus(this) != Passed) {
462 
463  // Increment the iteration
464  iter_++;
465 
466  // Multiply the current direction vector by A and store in AP_
467  lp_->applyOp( *P_, *AP_ );
468 
469  // Compute alpha := <R_,Z_> / <P_,AP_>
470  MVT::MvTransMv( one, *P_, *AP_, pAp );
471  alpha(0,0) = rHz(0,0) / pAp(0,0);
472 
473  // Check that alpha is a positive number!
474  if(assertPositiveDefiniteness_)
475  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero, CGIterateFailure,
476  "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
477  //
478  // Update the solution vector x := x + alpha * P_
479  //
480  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
481  lp_->updateSolution();
482  //
483  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
484  //
485  rHz_old(0,0) = rHz(0,0);
486  //
487  // Compute the new residual R_ := R_ - alpha * AP_
488  //
489  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
490  //
491  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
492  // and the new direction vector p.
493  //
494  if ( lp_->getLeftPrec() != Teuchos::null ) {
495  lp_->applyLeftPrec( *R_, *Z_ );
496  if ( lp_->getRightPrec() != Teuchos::null ) {
497  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
498  lp_->applyRightPrec( *Z_, *tmp );
499  Z_ = tmp;
500  }
501  }
502  else if ( lp_->getRightPrec() != Teuchos::null ) {
503  lp_->applyRightPrec( *R_, *Z_ );
504  }
505  else {
506  Z_ = R_;
507  }
508  //
509  MVT::MvTransMv( one, *R_, *Z_, rHz );
510  //
511  beta(0,0) = rHz(0,0) / rHz_old(0,0);
512  //
513  MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
514 
515  // Condition estimate (if needed)
516  if (doCondEst_) {
517  if (iter_ > 1) {
518  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp(0,0)) / rHz_old(0,0));
519  offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old(0,0) * rHz_old2)));
520  }
521  else {
522  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp(0,0) / rHz_old(0,0));
523  }
524  rHz_old2 = rHz_old(0,0);
525  beta_old = beta(0,0);
526  pAp_old = pAp(0,0);
527  }
528 
529  } // end while (sTest_->checkStatus(this) != Passed)
530  }
531 
532 } // end Belos namespace
533 
534 #endif /* BELOS_CG_ITER_HPP */
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::ScalarTraits< ScalarType > SCT
Definition: BelosCGIter.hpp:87
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:78
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(T a)
Declaration of basic traits for the multivector type.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
MultiVecTraits< ScalarType, MV > MVT
Definition: BelosCGIter.hpp:85
size_type size() const
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
CGIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, Teuchos::ParameterList &params)
CGIter constructor with linear problem, solver utilities, and parameter list of solver options...
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
bool isInitialized()
States whether the solver has been initialized or not.
OperatorTraits< ScalarType, MV, OP > OPT
Definition: BelosCGIter.hpp:86
int getNumIters() const
Get the current iteration count.
Traits class which defines basic operations on multivectors.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
virtual ~CGIter()
Destructor.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
SCT::magnitudeType MagnitudeType
Definition: BelosCGIter.hpp:88
void resetNumIters(int iter=0)
Reset the iteration count.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Class which defines basic traits for the operator type.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const MV > P
The current decent direction vector.
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.

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