Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockCGIter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_HPP
11 #define BELOS_PSEUDO_BLOCK_CG_ITER_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 #include "BelosCGIteration.hpp"
20 
21 #include "BelosLinearProblem.hpp"
22 #include "BelosMatOrthoManager.hpp"
23 #include "BelosOutputManager.hpp"
24 #include "BelosStatusTest.hpp"
25 #include "BelosOperatorTraits.hpp"
26 #include "BelosMultiVecTraits.hpp"
27 
28 #include "Teuchos_Assert.hpp"
31 #include "Teuchos_ScalarTraits.hpp"
33 #include "Teuchos_TimeMonitor.hpp"
34 
46 namespace Belos {
47 
49 
50 
55  template <class ScalarType, class MV>
56  class PseudoBlockCGIterationState : public CGIterationStateBase<ScalarType, MV> {
57 
58  public:
59  PseudoBlockCGIterationState() = default;
60 
62  initialize(tmp);
63  }
64 
65  virtual ~PseudoBlockCGIterationState() = default;
66 
67  void initialize(Teuchos::RCP<const MV> tmp, int _numVectors) {
69  this->R = MVT::Clone( *tmp, _numVectors );
70  this->Z = MVT::Clone( *tmp, _numVectors );
71  this->P = MVT::Clone( *tmp, _numVectors );
72  this->AP = MVT::Clone(*tmp, _numVectors );
73 
75  }
76 
77  bool matches(Teuchos::RCP<const MV> tmp, int _numVectors=1) const {
78  return CGIterationStateBase<ScalarType, MV>::matches(tmp, _numVectors);
79  }
80 };
81 
82  template<class ScalarType, class MV, class OP>
83  class PseudoBlockCGIter : virtual public CGIteration<ScalarType,MV,OP> {
84 
85  public:
86 
87  //
88  // Convenience typedefs
89  //
94 
96 
97 
104  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
106  Teuchos::ParameterList &params );
107 
109  virtual ~PseudoBlockCGIter() = default;
111 
112 
114 
115 
129  void iterate();
130 
152 
156  void initialize()
157  {
159  }
160 
170  state->R = R_;
171  state->P = P_;
172  state->AP = AP_;
173  state->Z = Z_;
174  return state;
175  }
176 
178  auto s = Teuchos::rcp_dynamic_cast<PseudoBlockCGIterationState<ScalarType,MV> >(state, true);
179  R_ = s->R;
180  Z_ = s->Z;
181  P_ = s->P;
182  AP_ = s->AP;
183  }
184 
186 
187 
189 
190 
192  int getNumIters() const { return iter_; }
193 
195  void resetNumIters( int iter = 0 ) { iter_ = iter; }
196 
199  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
200 
202 
205 
207 
209 
210 
212  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
213 
215  int getBlockSize() const { return 1; }
216 
218  void setBlockSize(int blockSize) {
219  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
220  "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
221  }
222 
224  bool isInitialized() { return initialized_; }
225 
227 
229  void setDoCondEst(bool val) {
230  if (numEntriesForCondEst_ != 0) doCondEst_=val;
231  }
232 
235  // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
236  // getDiag() didn't actually throw for me in that case, but why
237  // not be cautious?
238  using size_type = typename Teuchos::ArrayView<MagnitudeType>::size_type;
239  if (static_cast<size_type> (iter_) >= diag_.size ()) {
240  return diag_ ();
241  } else {
242  return diag_ (0, iter_);
243  }
244  }
245 
248  // NOTE (mfh 30 Jul 2015) The implementation as I found it
249  // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
250  // debug mode) when the maximum number of iterations has been
251  // reached, because iter_ == offdiag_.size() in that case. The
252  // new logic fixes this.
253  using size_type = typename Teuchos::ArrayView<MagnitudeType>::size_type;
254  if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
255  return offdiag_ ();
256  } else {
257  return offdiag_ (0, iter_);
258  }
259  }
260 
261  private:
262 
263  //
264  // Classes inputed through constructor that define the linear problem to be solved.
265  //
269 
270  //
271  // Algorithmic parameters
272  //
273  // numRHS_ is the current number of linear systems being solved.
274  int numRHS_;
275 
276  //
277  // Current solver state
278  //
279  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
280  // is capable of running; _initialize is controlled by the initialize() member method
281  // For the implications of the state of initialized_, please see documentation for initialize()
283 
284  // Current number of iterations performed.
285  int iter_;
286 
287  // Assert that the matrix is positive definite
289 
290  // Tridiagonal system for condition estimation (if needed)
292  ScalarType pAp_old_, beta_old_, rHz_old2_; // Put scalars here so that estimate is correct for multiple RHS, when deflation occurs.
295 
296  //
297  // State Storage
298  //
299  // Residual
301  //
302  // Preconditioned residual
304  //
305  // Direction vector
307  //
308  // Operator applied to direction vector
310 
311  };
312 
314  // Constructor.
315  template<class ScalarType, class MV, class OP>
317  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
319  Teuchos::ParameterList &params ):
320  lp_(problem),
321  om_(printer),
322  stest_(tester),
323  numRHS_(0),
324  initialized_(false),
325  iter_(0),
326  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
327  numEntriesForCondEst_(params.get("Max Size For Condest",0) ),
328  doCondEst_(false)
329  {
330  }
331 
332 
334  // Initialize this iteration object
335  template <class ScalarType, class MV, class OP>
337 
338  // Check if there is any mltivector to clone from.
339  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
340  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
341  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
342  "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
343 
344  // Get the multivector that is not null.
345  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
346 
347  // Get the number of right-hand sides we're solving for now.
348  int numRHS = MVT::GetNumberVecs(*tmp);
349  numRHS_ = numRHS;
350 
351  // Initialize the state storage if it isn't already.
352  TEUCHOS_ASSERT(!newstate.is_null());
353  if (!Teuchos::rcp_dynamic_cast<PseudoBlockCGIterationState<ScalarType,MV> >(newstate, true)->matches(tmp, numRHS_))
354  newstate->initialize(tmp, numRHS_);
355  setState(newstate);
356 
357  // Tracking information for condition number estimation
358  if(numEntriesForCondEst_ > 0) {
359  diag_.resize(numEntriesForCondEst_);
360  offdiag_.resize(numEntriesForCondEst_-1);
361  }
362 
363  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
364 
365  {
366 
367  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*R_0) != MVT::GetGlobalLength(*R_),
368  std::invalid_argument, errstr );
369  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*R_0) != numRHS_,
370  std::invalid_argument, errstr );
371 
372  // Copy basis vectors from newstate into V
373  if (R_0 != R_) {
374  // copy over the initial residual (unpreconditioned).
375  MVT::Assign( *R_0, *R_ );
376  }
377 
378  // Compute initial direction vectors
379  // Initially, they are set to the preconditioned residuals
380  //
381  if ( lp_->getLeftPrec() != Teuchos::null ) {
382  lp_->applyLeftPrec( *R_, *Z_ );
383  if ( lp_->getRightPrec() != Teuchos::null ) {
384  Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
385  lp_->applyRightPrec( *Z_, *tmp1 );
386  Z_ = tmp1;
387  }
388  }
389  else if ( lp_->getRightPrec() != Teuchos::null ) {
390  lp_->applyRightPrec( *R_, *Z_ );
391  }
392  else {
393  MVT::Assign( *R_, *Z_ );
394  }
395  MVT::Assign( *Z_, *P_ );
396  }
397 
398  // The solver is initialized
399  initialized_ = true;
400  }
401 
402 
404  // Iterate until the status test informs us we should stop.
405  template <class ScalarType, class MV, class OP>
407  {
408  //
409  // Allocate/initialize data structures
410  //
411  if (!initialized_) {
412  initialize();
413  }
414 
415  // Allocate memory for scalars.
416  int i=0;
417  std::vector<int> index(1);
418  std::vector<ScalarType> rHz( numRHS_ );
419  std::vector<ScalarType> rHz_old( numRHS_ );
420  std::vector<ScalarType> pAp( numRHS_ );
421  std::vector<ScalarType> beta( numRHS_ );
422  Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ );
423 
424  // Create convenience variables for zero and one.
425  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
427 
428  // Get the current solution std::vector.
429  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
430 
431  // Compute first <r,z> a.k.a. rHz
432  MVT::MvDot( *R_, *Z_, rHz );
433 
434  if ( assertPositiveDefiniteness_ )
435  for (i=0; i<numRHS_; ++i)
436  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
438  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
439 
441  // Iterate until the status test tells us to stop.
442  //
443  while (stest_->checkStatus(this) != Passed) {
444 
445  // Increment the iteration
446  iter_++;
447 
448  // Multiply the current direction std::vector by A and store in AP_
449  lp_->applyOp( *P_, *AP_ );
450 
451  // Compute alpha := <R_,Z_> / <P_,AP_>
452  MVT::MvDot( *P_, *AP_, pAp );
453 
454  for (i=0; i<numRHS_; ++i) {
455  if ( assertPositiveDefiniteness_ )
456  // Check that pAp[i] is a positive number!
457  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
459  "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
460 
461  alpha(i,i) = rHz[i] / pAp[i];
462  }
463 
464  //
465  // Update the solution std::vector x := x + alpha * P_
466  //
467  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
468  lp_->updateSolution();// what does this do?
469  //
470  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
471  //
472  for (i=0; i<numRHS_; ++i) {
473  rHz_old[i] = rHz[i];
474  }
475  //
476  // Compute the new residual R_ := R_ - alpha * AP_
477  //
478  MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
479  //
480  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
481  // and the new direction std::vector p.
482  //
483  if ( lp_->getLeftPrec() != Teuchos::null ) {
484  lp_->applyLeftPrec( *R_, *Z_ );
485  if ( lp_->getRightPrec() != Teuchos::null ) {
486  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
487  lp_->applyRightPrec( *Z_, *tmp );
488  Z_ = tmp;
489  }
490  }
491  else if ( lp_->getRightPrec() != Teuchos::null ) {
492  lp_->applyRightPrec( *R_, *Z_ );
493  }
494  else {
495  Z_ = R_;
496  }
497  //
498  MVT::MvDot( *R_, *Z_, rHz );
499  if ( assertPositiveDefiniteness_ )
500  for (i=0; i<numRHS_; ++i)
501  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
503  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
504  //
505  // Update the search directions.
506  for (i=0; i<numRHS_; ++i) {
507  beta[i] = rHz[i] / rHz_old[i];
508  index[0] = i;
509  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
510  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
511  MVT::MvAddMv( one, *Z_i, beta[i], *P_i, *P_i );
512  }
513 
514  // Condition estimate (if needed)
515  if (doCondEst_ && (iter_ - 1) < diag_.size()) {
516  if (iter_ > 1) {
517  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old_ * beta_old_ * pAp_old_ + pAp[0]) / rHz_old[0]);
518  offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old_ * pAp_old_ / (sqrt( rHz_old[0] * rHz_old2_)));
519  }
520  else {
521  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
522  }
523  rHz_old2_ = rHz_old[0];
524  beta_old_ = beta[0];
525  pAp_old_ = pAp[0];
526  }
527 
528 
529  //
530  } // end while (sTest_->checkStatus(this) != Passed)
531  }
532 
533 } // end Belos namespace
534 
535 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
void initialize(Teuchos::RCP< const MV > tmp, int _numVectors)
Collection of types and exceptions used within the Belos solvers.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
virtual ~PseudoBlockCGIter()=default
Destructor.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
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.
Teuchos::ArrayRCP< MagnitudeType > offdiag_
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setState(Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > state)
static magnitudeType real(T a)
Declaration of basic traits for the multivector type.
typename SCT::magnitudeType MagnitudeType
void setBlockSize(int blockSize)
Set the blocksize.
bool matches(Teuchos::RCP< const MV > tmp, int _numVectors=1) const
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.
Teuchos::RCP< MV > P
The current decent direction vector.
Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > getState() const
Get the current state of the linear solver.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
CGPositiveDefiniteFailure is thrown when the the CG &#39;alpha = p^H*A*P&#39; value is less than zero...
void resetNumIters(int iter=0)
Reset the iteration count.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Structure to contain pointers to PseudoBlockCGIteration state variables.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
virtual void initialize(Teuchos::RCP< const MV > tmp, int _numVectors)
A linear system to solve, and its associated information.
void initializeCG(Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > newstate, Teuchos::RCP< MV > R_0)
Initialize the solver to an iterate, providing a complete state.
int getNumIters() const
Get the current iteration count.
virtual ~PseudoBlockCGIterationState()=default
Class which describes the linear problem to be solved by the iterative solver.
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< MV > Z
The current preconditioned residual.
PseudoBlockCGIter(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)
PseudoBlockCGIter constructor with linear problem, solver utilities, and parameter list of solver opt...
TransListIter iter
Teuchos::RCP< MV > AP
The matrix A applied to current decent direction vector.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< MV > R
The current residual.
Belos header file which uses auto-configuration information to include necessary C++ headers...
const Teuchos::RCP< OutputManager< ScalarType > > om_
virtual bool matches(Teuchos::RCP< const MV > tmp, int _numVectors=1) const
PseudoBlockCGIterationState(Teuchos::RCP< const MV > tmp)
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ArrayRCP< MagnitudeType > diag_