Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Belos_PseudoBlockCGIter_MP_Vector.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP
11 #define BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP
12 
22 #include "BelosPseudoBlockCGIter.hpp"
23 
39 namespace Belos {
40 
41  template<class Storage, class MV, class OP>
42  class PseudoBlockCGIter<Sacado::MP::Vector<Storage>, MV, OP> :
43  virtual public CGIteration<Sacado::MP::Vector<Storage>, MV, OP> {
44 
45  public:
46 
47  //
48  // Convenience typedefs
49  //
51  typedef MultiVecTraits<ScalarType,MV> MVT;
52  typedef OperatorTraits<ScalarType,MV,OP> OPT;
56 
58 
59 
65  PseudoBlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
66  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
67  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
68  Teuchos::ParameterList &params );
69 
71  virtual ~PseudoBlockCGIter() {};
73 
74 
76 
77 
91  void iterate();
92 
113  void initializeCG(Teuchos::RCP<CGIterationStateBase<ScalarType,MV> > newstate, Teuchos::RCP<MV> R_0);
114 
118  void initialize()
119  {
120  initializeCG(Teuchos::null, Teuchos::null);
121  }
122 
131  auto state = Teuchos::rcp(new PseudoBlockCGIterationState<ScalarType,MV>());
132  state->R = R_;
133  state->P = P_;
134  state->AP = AP_;
135  state->Z = Z_;
136  return state;
137  }
138 
139  void setState(Teuchos::RCP<CGIterationStateBase<ScalarType,MV> > state) {
140  auto s = Teuchos::rcp_dynamic_cast<PseudoBlockCGIterationState<ScalarType,MV> >(state, true);
141  R_ = s->R;
142  Z_ = s->Z;
143  P_ = s->P;
144  AP_ = s->AP;
145  }
146 
148 
149 
151 
152 
154  int getNumIters() const { return iter_; }
155 
157  void resetNumIters( int iter = 0 ) { iter_ = iter; }
158 
161  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
162 
164 
167 
169 
171 
172 
174  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
175 
177  int getBlockSize() const { return 1; }
178 
180  void setBlockSize(int blockSize) {
181  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
182  "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
183  }
184 
186  bool isInitialized() { return initialized_; }
187 
189 
191  void setDoCondEst(bool val){doCondEst_=val;}
192 
195  // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
196  // getDiag() didn't actually throw for me in that case, but why
197  // not be cautious?
198  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
199  if (static_cast<size_type> (iter_) >= diag_.size ()) {
200  return diag_ ();
201  } else {
202  return diag_ (0, iter_);
203  }
204  }
205 
208  // NOTE (mfh 30 Jul 2015) The implementation as I found it
209  // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
210  // debug mode) when the maximum number of iterations has been
211  // reached, because iter_ == offdiag_.size() in that case. The
212  // new logic fixes this.
213  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
214  if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
215  return offdiag_ ();
216  } else {
217  return offdiag_ (0, iter_);
218  }
219  }
220 
221  private:
222 
223  //
224  // Classes inputed through constructor that define the linear problem to be solved.
225  //
229 
230  //
231  // Algorithmic parameters
232  //
233  // numRHS_ is the current number of linear systems being solved.
234  int numRHS_;
235 
236  //
237  // Current solver state
238  //
239  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
240  // is capable of running; _initialize is controlled by the initialize() member method
241  // For the implications of the state of initialized_, please see documentation for initialize()
243 
244  // Current number of iterations performed.
245  int iter_;
246 
247  // Assert that the matrix is positive definite
249 
250  // Tridiagonal system for condition estimation (if needed)
254 
255  //
256  // State Storage
257  //
258  // Residual
260  //
261  // Preconditioned residual
263  //
264  // Direction vector
266  //
267  // Operator applied to direction vector
269 
270  // Tolerance for ensemble breakdown
272 
273  };
274 
276  // Constructor.
277  template<class StorageType, class MV, class OP>
279  PseudoBlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
280  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
281  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
282  Teuchos::ParameterList &params ):
283  lp_(problem),
284  om_(printer),
285  stest_(tester),
286  numRHS_(0),
287  initialized_(false),
288  iter_(0),
289  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
290  numEntriesForCondEst_(params.get("Max Size For Condest",0) ),
291  breakDownTol_(params.get("Ensemble Breakdown Tolerance", 0.0))
292  {
293  }
294 
295 
297  // Initialize this iteration object
298  template <class StorageType, class MV, class OP>
299  void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
300  initializeCG(Teuchos::RCP<CGIterationStateBase<ScalarType,MV> > newstate, Teuchos::RCP<MV> R_0) {
301 
302  // Check if there is any mltivector to clone from.
303  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
304  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
305  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
306  "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
307 
308  // Get the multivector that is not null.
309  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
310 
311  // Get the number of right-hand sides we're solving for now.
312  int numRHS = MVT::GetNumberVecs(*tmp);
313  numRHS_ = numRHS;
314 
315  // Initialize the state storage if it isn't already.
316  if (!Teuchos::rcp_dynamic_cast<PseudoBlockCGIterationState<ScalarType,MV> >(newstate, true)->matches(tmp, numRHS_))
317  newstate->initialize(tmp, numRHS_);
318  setState(newstate);
319 
320  // Tracking information for condition number estimation
321  if(numEntriesForCondEst_ > 0) {
322  diag_.resize(numEntriesForCondEst_);
323  offdiag_.resize(numEntriesForCondEst_-1);
324  }
325 
326  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
327 
328  {
329 
330  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*R_0) != MVT::GetGlobalLength(*R_),
331  std::invalid_argument, errstr );
332  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*R_0) != numRHS_,
333  std::invalid_argument, errstr );
334 
335  // Copy basis vectors from newstate into V
336  if (R_0 != R_) {
337  // copy over the initial residual (unpreconditioned).
338  MVT::Assign( *R_0, *R_ );
339  }
340 
341  // Compute initial direction vectors
342  // Initially, they are set to the preconditioned residuals
343  //
344  if ( lp_->getLeftPrec() != Teuchos::null ) {
345  lp_->applyLeftPrec( *R_, *Z_ );
346  if ( lp_->getRightPrec() != Teuchos::null ) {
347  Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
348  lp_->applyRightPrec( *Z_, *tmp1 );
349  Z_ = tmp1;
350  }
351  }
352  else if ( lp_->getRightPrec() != Teuchos::null ) {
353  lp_->applyRightPrec( *R_, *Z_ );
354  }
355  else {
356  MVT::Assign( *R_, *Z_ );
357  }
358  MVT::Assign( *Z_, *P_ );
359  }
360 
361  // The solver is initialized
362  initialized_ = true;
363  }
364 
365 
367  // Iterate until the status test informs us we should stop.
368  template <class StorageType, class MV, class OP>
369  void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
370  iterate()
371  {
372  //
373  // Allocate/initialize data structures
374  //
375  if (initialized_ == false) {
376  initialize();
377  }
378 
379  // Allocate memory for scalars.
380  int i=0;
381  std::vector<int> index(1);
382  std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
383  Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
384 
385  // Create convenience variables for zero and one.
386  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
387  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
388 
389  // Scalars for condition estimation (if needed) - These will always use entry zero, for convenience
390  ScalarType pAp_old=one, beta_old=one ,rHz_old2=one;
391 
392  // Get the current solution std::vector.
393  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
394 
395  // Compute first <r,z> a.k.a. rHz
396  MVT::MvDot( *R_, *Z_, rHz );
397 
398  if ( assertPositiveDefiniteness_ )
399  for (i=0; i<numRHS_; ++i)
400  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
401  CGIterateFailure,
402  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
403 
405  // Iterate until the status test tells us to stop.
406  //
407  while (stest_->checkStatus(this) != Passed) {
408 
409  // Increment the iteration
410  iter_++;
411 
412  // Multiply the current direction std::vector by A and store in AP_
413  lp_->applyOp( *P_, *AP_ );
414 
415  // Compute alpha := <R_,Z_> / <P_,AP_>
416  MVT::MvDot( *P_, *AP_, pAp );
417 
418  for (i=0; i<numRHS_; ++i) {
419  if ( assertPositiveDefiniteness_ )
420  // Check that pAp[i] is a positive number!
421  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
422  CGIterateFailure,
423  "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
424 
425  const int ensemble_size = pAp[i].size();
426  for (int k=0; k<ensemble_size; ++k) {
427  if (SVT::magnitude(pAp[i].coeff(k)) <= breakDownTol_)
428  alpha(i,i).fastAccessCoeff(k) = 0.0;
429  else
430  alpha(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / pAp[i].coeff(k);
431  }
432  }
433 
434  //
435  // Update the solution std::vector x := x + alpha * P_
436  //
437  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
438  lp_->updateSolution();
439  //
440  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
441  //
442  for (i=0; i<numRHS_; ++i) {
443  rHz_old[i] = rHz[i];
444  }
445  //
446  // Compute the new residual R_ := R_ - alpha * AP_
447  //
448  MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
449  //
450  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
451  // and the new direction std::vector p.
452  //
453  if ( lp_->getLeftPrec() != Teuchos::null ) {
454  lp_->applyLeftPrec( *R_, *Z_ );
455  if ( lp_->getRightPrec() != Teuchos::null ) {
456  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
457  lp_->applyRightPrec( *Z_, *tmp );
458  Z_ = tmp;
459  }
460  }
461  else if ( lp_->getRightPrec() != Teuchos::null ) {
462  lp_->applyRightPrec( *R_, *Z_ );
463  }
464  else {
465  Z_ = R_;
466  }
467  //
468  MVT::MvDot( *R_, *Z_, rHz );
469  if ( assertPositiveDefiniteness_ )
470  for (i=0; i<numRHS_; ++i)
471  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
472  CGIterateFailure,
473  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
474  //
475  // Update the search directions.
476  for (i=0; i<numRHS_; ++i) {
477  const int ensemble_size = rHz[i].size();
478  for (int k=0; k<ensemble_size; ++k) {
479  if (SVT::magnitude(rHz[i].coeff(k)) <= breakDownTol_)
480  beta(i,i).fastAccessCoeff(k) = 0.0;
481  else
482  beta(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / rHz_old[i].coeff(k);
483  }
484  index[0] = i;
485  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
486  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
487  MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
488  }
489 
490  // Condition estimate (if needed)
491  if(doCondEst_ > 0) {
492  if(iter_ > 1 ) {
493  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
494  offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
495  }
496  else {
497  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
498  }
499  rHz_old2 = rHz_old[0];
500  beta_old = beta(0,0);
501  pAp_old = pAp[0];
502  }
503 
504 
505  //
506  } // end while (sTest_->checkStatus(this) != Passed)
507  }
508 
509 } // end Belos namespace
510 
511 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(T a)
void setState(Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > state)
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > getState() const
Get the current state of the linear solver.
expr val()
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.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
bool isInitialized()
States whether the solver has been initialized or not.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.