42#ifndef BELOS_CG_ITER_HPP
43#define BELOS_CG_ITER_HPP
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_ScalarTraits.hpp"
63#include "Teuchos_ParameterList.hpp"
64#include "Teuchos_TimeMonitor.hpp"
78template<
class ScalarType,
class MV,
class OP>
88 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
103 Teuchos::ParameterList ¶ms );
182 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
183 (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHr_));
184 return Teuchos::null;
207 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
208 "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
217 if (numEntriesForCondEst_) doCondEst_=val;
225 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
226 if (
static_cast<size_type
> (iter_) >= diag_.size ()) {
229 return diag_ (0, iter_);
240 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
241 if (
static_cast<size_type
> (iter_) >= offdiag_.size ()) {
244 return offdiag_ (0, iter_);
261 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
262 const Teuchos::RCP<OutputManager<ScalarType> > om_;
263 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
264 const Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
277 bool stateStorageInitialized_;
283 bool foldConvergenceDetectionIntoAllreduce_;
289 bool assertPositiveDefiniteness_;
292 Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
293 int numEntriesForCondEst_;
311 Teuchos::RCP<MV> AP_;
319 template<
class ScalarType,
class MV,
class OP>
324 Teuchos::ParameterList ¶ms ):
328 convTest_(convTester),
330 stateStorageInitialized_(false),
332 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
333 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
336 foldConvergenceDetectionIntoAllreduce_ = params.get<
bool>(
"Fold Convergence Detection Into Allreduce",
false);
341 template <
class ScalarType,
class MV,
class OP>
344 if (!stateStorageInitialized_) {
347 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
348 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
349 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
350 stateStorageInitialized_ =
false;
357 if (R_ == Teuchos::null) {
359 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
360 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
361 "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
362 S_ = MVT::Clone( *tmp, 2 );
363 std::vector<int> index(1,0);
365 R_ = MVT::CloneViewNonConst( *S_, index );
367 Z_ = MVT::CloneViewNonConst( *S_, index );
368 P_ = MVT::Clone( *tmp, 1 );
369 AP_ = MVT::Clone( *tmp, 1 );
374 if(numEntriesForCondEst_ > 0) {
375 diag_.resize(numEntriesForCondEst_);
376 offdiag_.resize(numEntriesForCondEst_-1);
380 stateStorageInitialized_ =
true;
388 template <
class ScalarType,
class MV,
class OP>
392 if (!stateStorageInitialized_)
395 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
396 "Belos::CGIter::initialize(): Cannot initialize state storage!");
400 std::string errstr(
"Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
402 if (newstate.
R != Teuchos::null) {
404 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
405 std::invalid_argument, errstr );
406 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
407 std::invalid_argument, errstr );
410 if (newstate.
R != R_) {
412 MVT::Assign( *newstate.
R, *R_ );
418 if ( lp_->getLeftPrec() != Teuchos::null ) {
419 lp_->applyLeftPrec( *R_, *Z_ );
420 if ( lp_->getRightPrec() != Teuchos::null ) {
421 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
422 lp_->applyRightPrec( *tmp, *Z_ );
425 else if ( lp_->getRightPrec() != Teuchos::null ) {
426 lp_->applyRightPrec( *R_, *Z_ );
429 MVT::Assign( *R_, *Z_ );
431 MVT::Assign( *Z_, *P_ );
435 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
436 "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
446 template <
class ScalarType,
class MV,
class OP>
452 if (initialized_ ==
false) {
457 std::vector<ScalarType> alpha(1);
458 std::vector<ScalarType> beta(1);
459 std::vector<ScalarType> rHz(1), rHz_old(1), pAp(1);
460 Teuchos::SerialDenseMatrix<int,ScalarType> rHs( 1, 2 );
463 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
464 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
467 ScalarType pAp_old = one, beta_old = one, rHz_old2 = one;
470 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
473 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
474 "Belos::CGIter::iterate(): current linear system has more than one vector!" );
477 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
478 MVT::MvTransMv( one, *R_, *S_, rHs );
482 MVT::MvDot( *R_, *Z_, rHz );
487 while (stest_->checkStatus(
this) !=
Passed) {
493 lp_->applyOp( *P_, *AP_ );
496 MVT::MvDot( *P_, *AP_, pAp );
497 alpha[0] = rHz[0] / pAp[0];
500 if(assertPositiveDefiniteness_) {
502 "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
508 MVT::MvAddMv( one, *cur_soln_vec, alpha[0], *P_, *cur_soln_vec );
509 lp_->updateSolution();
517 MVT::MvAddMv( one, *R_, -alpha[0], *AP_, *R_ );
522 if ( lp_->getLeftPrec() != Teuchos::null ) {
523 lp_->applyLeftPrec( *R_, *Z_ );
524 if ( lp_->getRightPrec() != Teuchos::null ) {
525 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_);
526 lp_->applyRightPrec( *tmp, *Z_ );
529 else if ( lp_->getRightPrec() != Teuchos::null ) {
530 lp_->applyRightPrec( *R_, *Z_ );
533 MVT::Assign( *R_, *Z_ );
536 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
537 MVT::MvTransMv( one, *R_, *S_, rHs );
541 MVT::MvDot( *R_, *Z_, rHz );
543 beta[0] = rHz[0] / rHz_old[0];
545 MVT::MvAddMv( one, *Z_, beta[0], *P_, *P_ );
550 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
551 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
554 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
556 rHz_old2 = rHz_old[0];
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
virtual ~CGIter()
Destructor.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Teuchos::ScalarTraits< ScalarType > SCT
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getNumIters() const
Get the current iteration count.
bool isInitialized()
States whether the solver has been initialized or not.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
void resetNumIters(int iter=0)
Reset the iteration count.
MultiVecTraits< ScalarType, MV > MVT
CGIter(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 ¶ms)
CGIter constructor with linear problem, solver utilities, and parameter list of solver options.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero,...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
An implementation of StatusTestResNorm using a family of residual norms.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< const MV > P
The current decent direction vector.
Teuchos::RCP< const MV > R
The current residual.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.