Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosCGSingleRedIter.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_CG_SINGLE_RED_ITER_HPP
11 #define BELOS_CG_SINGLE_RED_ITER_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 #include "BelosCGIteration.hpp"
20 
21 #include "BelosLinearProblem.hpp"
22 #include "BelosOutputManager.hpp"
23 #include "BelosStatusTest.hpp"
25 #include "BelosOperatorTraits.hpp"
26 #include "BelosMultiVecTraits.hpp"
27 
30 #include "Teuchos_ScalarTraits.hpp"
32 #include "Teuchos_TimeMonitor.hpp"
33 
44 namespace Belos {
45 
46 template<class ScalarType, class MV, class OP>
47 class CGSingleRedIter : virtual public CGIteration<ScalarType,MV,OP> {
48 
49  public:
50 
51  //
52  // Convenience typedefs
53  //
58 
60 
61 
68  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
71  Teuchos::ParameterList &params );
72 
74  virtual ~CGSingleRedIter() {};
76 
77 
79 
80 
93  void iterate();
94 
110 
114  void initialize()
115  {
117  initializeCG(empty);
118  }
119 
128  state.R = R_;
129  state.P = P_;
130  state.AP = AP_;
131  state.Z = Z_;
132  return state;
133  }
134 
136 
137 
139 
140 
142  int getNumIters() const { return iter_; }
143 
145  void resetNumIters( int iter = 0 ) { iter_ = iter; }
146 
149  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
150 
152 
154  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
155 
157 
159 
160 
162  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
163 
165  int getBlockSize() const { return 1; }
166 
168  void setBlockSize(int blockSize) {
169  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
170  "Belos::CGSingleRedIter::setBlockSize(): Cannot use a block size that is not one.");
171  }
172 
174  bool isInitialized() { return initialized_; }
175 
177  void setDoCondEst(bool /* val */){/*ignored*/}
178 
182  return temp;
183  }
184 
188  return temp;
189  }
190 
192 
193  private:
194 
195  //
196  // Internal methods
197  //
199  void setStateSize();
200 
201  //
202  // Classes inputed through constructor that define the linear problem to be solved.
203  //
208 
209  //
210  // Current solver state
211  //
212  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
213  // is capable of running; _initialize is controlled by the initialize() member method
214  // For the implications of the state of initialized_, please see documentation for initialize()
215  bool initialized_;
216 
217  // stateStorageInitialized_ specifies that the state storage has been initialized.
218  // This initialization may be postponed if the linear problem was generated without
219  // the right-hand side or solution vectors.
220  bool stateStorageInitialized_;
221 
222  // Current number of iterations performed.
223  int iter_;
224 
225  // Should the allreduce for convergence detection be merged with the inner products?
226  bool foldConvergenceDetectionIntoAllreduce_;
227 
228  // <r,z>
229  ScalarType rHz_;
230  // <r,r>
231  ScalarType rHr_;
232 
233  //
234  // State Storage
235  //
236  // Residual
237  Teuchos::RCP<MV> R_;
238  // Preconditioned residual
239  Teuchos::RCP<MV> Z_;
240  // Operator applied to preconditioned residual
241  Teuchos::RCP<MV> AZ_;
242  // Direction vector
243  Teuchos::RCP<MV> P_;
244  // Operator applied to direction vector
245  Teuchos::RCP<MV> AP_;
246  //
247  // W_ = (R_, AZ_, Z_)
248  Teuchos::RCP<MV> W_;
249  // S_ = (R_, AZ_)
250  Teuchos::RCP<MV> S_;
251  // T_ = (Z_, R_)
252  Teuchos::RCP<MV> T_;
253  // U_ = (AZ_, Z_)
254  Teuchos::RCP<MV> U_;
255  // V_ = (AP_, P_)
256  Teuchos::RCP<MV> V_;
257 
258 };
259 
261  // Constructor.
262  template<class ScalarType, class MV, class OP>
264  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
267  Teuchos::ParameterList &params ):
268  lp_(problem),
269  om_(printer),
270  stest_(tester),
271  convTest_(convTester),
272  initialized_(false),
273  stateStorageInitialized_(false),
274  iter_(0)
275  {
276  foldConvergenceDetectionIntoAllreduce_ = params.get<bool>("Fold Convergence Detection Into Allreduce",false);
277  }
278 
280  // Setup the state storage.
281  template <class ScalarType, class MV, class OP>
283  {
284  if (!stateStorageInitialized_) {
285 
286  // Check if there is any multivector to clone from.
287  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
288  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
289  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
290  stateStorageInitialized_ = false;
291  return;
292  }
293  else {
294 
295  // Initialize the state storage
296  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
297  if (R_ == Teuchos::null) {
298  // Get the multivector that is not null.
299  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
300  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
301  "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
302 
303  // W_ = (R_, AZ_, Z_)
304  W_ = MVT::Clone( *tmp, 3 );
305  std::vector<int> index2(2,0);
306  std::vector<int> index(1,0);
307 
308  // S_ = (R_, AZ_)
309  index2[0] = 0;
310  index2[1] = 1;
311  S_ = MVT::CloneViewNonConst( *W_, index2 );
312 
313  // U_ = (AZ_, Z_)
314  index2[0] = 1;
315  index2[1] = 2;
316  U_ = MVT::CloneViewNonConst( *W_, index2 );
317 
318  index[0] = 0;
319  R_ = MVT::CloneViewNonConst( *W_, index );
320  index[0] = 1;
321  AZ_ = MVT::CloneViewNonConst( *W_, index );
322  index[0] = 2;
323  Z_ = MVT::CloneViewNonConst( *W_, index );
324 
325  // T_ = (R_, Z_)
326  index2[0] = 0;
327  index2[1] = 2;
328  T_ = MVT::CloneViewNonConst( *W_, index2 );
329 
330  // V_ = (AP_, P_)
331  V_ = MVT::Clone( *tmp, 2 );
332  index[0] = 0;
333  AP_ = MVT::CloneViewNonConst( *V_, index );
334  index[0] = 1;
335  P_ = MVT::CloneViewNonConst( *V_, index );
336 
337  }
338 
339  // State storage has now been initialized.
340  stateStorageInitialized_ = true;
341  }
342  }
343  }
344 
345 
347  // Initialize this iteration object
348  template <class ScalarType, class MV, class OP>
350  {
351  // Initialize the state storage if it isn't already.
352  if (!stateStorageInitialized_)
353  setStateSize();
354 
355  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
356  "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
357 
358  // NOTE: In CGSingleRedIter R_, the initial residual, is required!!!
359  //
360  std::string errstr("Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
361 
362  if (newstate.R != Teuchos::null) {
363 
364  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
365  std::invalid_argument, errstr );
366  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
367  std::invalid_argument, errstr );
368 
369  // Copy basis vectors from newstate into V
370  if (newstate.R != R_) {
371  // copy over the initial residual (unpreconditioned).
372  MVT::Assign( *newstate.R, *R_ );
373  }
374 
375  // Compute initial direction vectors
376  // Initially, they are set to the preconditioned residuals
377  //
378  if ( lp_->getLeftPrec() != Teuchos::null ) {
379  lp_->applyLeftPrec( *R_, *Z_ );
380  if ( lp_->getRightPrec() != Teuchos::null ) {
381  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
382  lp_->applyRightPrec( *Z_, *tmp );
383  MVT::Assign( *tmp, *Z_ );
384  }
385  }
386  else if ( lp_->getRightPrec() != Teuchos::null ) {
387  lp_->applyRightPrec( *R_, *Z_ );
388  }
389  else {
390  MVT::Assign( *R_, *Z_ );
391  }
392 
393  // Multiply the current preconditioned residual vector by A and store in AZ_
394  lp_->applyOp( *Z_, *AZ_ );
395 
396  // P_ := Z_
397  // Logically, AP_ := AZ_
398  MVT::Assign( *U_, *V_);
399  }
400  else {
401 
402  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
403  "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
404  }
405 
406  // The solver is initialized
407  initialized_ = true;
408  }
409 
410 
412  // Get the norms of the residuals native to the solver.
413  template <class ScalarType, class MV, class OP>
415  CGSingleRedIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const {
416  if (convTest_->getResNormType() == Belos::PreconditionerNorm) {
417  (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHz_));
418  return Teuchos::null;
419  } else if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
420  (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHr_));
421  return Teuchos::null;
422  } else
423  return R_;
424  }
425 
426 
428  // Iterate until the status test informs us we should stop.
429  template <class ScalarType, class MV, class OP>
431  {
432  //
433  // Allocate/initialize data structures
434  //
435  if (initialized_ == false) {
436  initialize();
437  }
438 
439  // Allocate memory for scalars.
442  ScalarType rHz_old, alpha, beta, delta;
443 
444  // Create convenience variables for zero and one.
445  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::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::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
454 
455  if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
456  // Compute first <S_,T_> a.k.a. <R_,Z_>, <AZ_,Z_> and <R_,R_> combined (also computes unneeded <AZ_,R_>)
457  MVT::MvTransMv( one, *S_, *T_, sHt );
458  rHz_ = sHt(0,1);
459  delta = sHt(1,1);
460  rHr_ = sHt(0,0);
461  } else {
462  // Compute first <s,z> a.k.a. <r,z> and <Az,z> combined
463  MVT::MvTransMv( one, *S_, *Z_, sHz );
464  rHz_ = sHz(0,0);
465  delta = sHz(1,0);
466  }
468  (stest_->checkStatus(this) == Passed))
469  return;
470  alpha = rHz_ / delta;
471 
472  // Check that alpha is a positive number!
473  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGPositiveDefiniteFailure,
474  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
475 
477  // Iterate until the status test tells us to stop.
478  //
479  if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
481  // Iterate until the status test tells us to stop.
482  //
483  while (1) {
484 
485  // Update the solution vector x := x + alpha * P_
486  //
487  MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
488  //
489  // Compute the new residual R_ := R_ - alpha * AP_
490  //
491  MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
492  //
493  // Apply preconditioner to new residual to update Z_
494  //
495  if ( lp_->getLeftPrec() != Teuchos::null ) {
496  lp_->applyLeftPrec( *R_, *Z_ );
497  if ( lp_->getRightPrec() != Teuchos::null ) {
498  Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
499  lp_->applyRightPrec( *tmp, *Z_ );
500  }
501  }
502  else if ( lp_->getRightPrec() != Teuchos::null ) {
503  lp_->applyRightPrec( *R_, *Z_ );
504  }
505  else {
506  MVT::Assign( *R_, *Z_ );
507  }
508  //
509  // Multiply the current preconditioned residual vector by A and store in AZ_
510  lp_->applyOp( *Z_, *AZ_ );
511  //
512  // Compute <S_,T_> a.k.a. <R_,Z_>, <AZ_,Z_> and <R_,R_> combined (also computes unneeded <AZ_,R_>)
513  MVT::MvTransMv( one, *S_, *T_, sHt );
514  //
515  // Update scalars.
516  rHz_old = rHz_;
517  rHz_ = sHt(0,1);
518  delta = sHt(1,1);
519  rHr_ = sHt(0,0);
520 
521  // Increment the iteration
522  iter_++;
523  //
524  // Check the status test, now that the solution and residual have been updated
525  //
526  if (stest_->checkStatus(this) == Passed) {
527  break;
528  }
529  //
530  beta = rHz_ / rHz_old;
531  alpha = rHz_ / (delta - (beta*rHz_ / alpha));
532  //
533  // Check that alpha is a positive number!
534  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGPositiveDefiniteFailure,
535  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
536 
537  //
538  // Update the direction vector P_ := Z_ + beta * P_
539  // Update AP_ through recurrence relation AP_ := AZ_ + beta * AP_
540  // Hence: V_ = (AP_, P_) := (AZ_, Z_) + beta (AP_, P_) = U_ + beta * V_
541  //
542  MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
543 
544  } // end while (1)
545  } else {
547  // Iterate until the status test tells us to stop.
548  //
549  while (1) {
550 
551  // Update the solution vector x := x + alpha * P_
552  //
553  MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
554  //
555  // Compute the new residual R_ := R_ - alpha * AP_
556  //
557  MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
558 
559  // Increment the iteration
560  iter_++;
561  //
562  // Check the status test, now that the solution and residual have been updated
563  //
564  if (stest_->checkStatus(this) == Passed) {
565  break;
566  }
567  //
568  // Apply preconditioner to new residual to update Z_
569  //
570  if ( lp_->getLeftPrec() != Teuchos::null ) {
571  lp_->applyLeftPrec( *R_, *Z_ );
572  if ( lp_->getRightPrec() != Teuchos::null ) {
573  Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
574  lp_->applyRightPrec( *tmp, *Z_ );
575  }
576  }
577  else if ( lp_->getRightPrec() != Teuchos::null ) {
578  lp_->applyRightPrec( *R_, *Z_ );
579  }
580  else {
581  MVT::Assign( *R_, *Z_ );
582  }
583  //
584  // Multiply the current preconditioned residual vector by A and store in AZ_
585  lp_->applyOp( *Z_, *AZ_ );
586  //
587  // Compute <S_,Z_> a.k.a. <R_,Z_> and <AZ_,Z_> combined
588  MVT::MvTransMv( one, *S_, *Z_, sHz );
589  //
590  // Update scalars.
591  rHz_old = rHz_;
592  rHz_ = sHz(0,0);
593  delta = sHz(1,0);
594  //
595  beta = rHz_ / rHz_old;
596  alpha = rHz_ / (delta - (beta*rHz_ / alpha));
597  //
598  // Check that alpha is a positive number!
599  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGPositiveDefiniteFailure,
600  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
601 
602  //
603  // Update the direction vector P_ := Z_ + beta * P_
604  // Update AP_ through recurrence relation AP_ := AZ_ + beta * AP_
605  // Hence: V_ = (AP_, P_) := (AZ_, Z_) + beta (AP_, P_) = U_ + beta * V_
606  //
607  MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
608 
609  } // end while (1)
610  }
611  }
612 
613 } // end Belos namespace
614 
615 #endif /* BELOS_CG_SINGLE_RED_ITER_HPP */
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
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.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
int getNumIters() const
Get the current iteration count.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
T & get(const std::string &name, T def_value)
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
An implementation of StatusTestResNorm using a family of residual norms.
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
MultiVecTraits< ScalarType, MV > MVT
Traits class which defines basic operations on multivectors.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
CGPositiveDefiniteFailure is thrown when the the CG &#39;alpha = p^H*A*P&#39; value is less than zero...
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation (NOT_IMPLEMENTED)
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
bool isInitialized()
States whether the solver has been initialized or not.
virtual ~CGSingleRedIter()
Destructor.
SCT::magnitudeType MagnitudeType
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
CGSingleRedIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > &convTester, Teuchos::ParameterList &params)
CGSingleRedIter constructor with linear problem, solver utilities, and parameter list of solver optio...
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Class which defines basic traits for the operator type.
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
Teuchos::ScalarTraits< ScalarType > SCT
Belos header file which uses auto-configuration information to include necessary C++ headers...
void resetNumIters(int iter=0)
Reset the iteration count.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > P
The current decent direction vector.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void setDoCondEst(bool)
Sets whether or not to store the diagonal for condition estimation.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation (NOT_IMPLEMENTED)

Generated on Mon Jul 15 2024 09:24:23 for Belos by doxygen 1.8.5