42 #ifndef BELOS_CG_ITER_HPP
43 #define BELOS_CG_ITER_HPP
78 template<
class ScalarType,
class MV,
class OP>
182 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
184 return Teuchos::null;
208 "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
217 if (numEntriesForCondEst_) doCondEst_=val;
226 if (static_cast<size_type> (iter_) >= diag_.
size ()) {
229 return diag_ (0, iter_);
241 if (static_cast<size_type> (iter_) >= offdiag_.
size ()) {
244 return offdiag_ (0, iter_);
277 bool stateStorageInitialized_;
283 bool foldConvergenceDetectionIntoAllreduce_;
289 bool assertPositiveDefiniteness_;
293 int numEntriesForCondEst_;
319 template<
class ScalarType,
class MV,
class OP>
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_) {
349 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
350 stateStorageInitialized_ =
false;
357 if (R_ == Teuchos::null) {
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_)
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) {
405 std::invalid_argument, errstr );
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 ) {
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_ );
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);
467 ScalarType pAp_old = one, beta_old = one, rHz_old2 = one;
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 ) {
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_ );
556 rHz_old2 = rHz_old[0];
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::ScalarTraits< ScalarType > SCT
Belos'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.
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)
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.
An implementation of StatusTestResNorm using a family of residual norms.
MultiVecTraits< ScalarType, MV > MVT
A pure virtual class for defining the status tests for the Belos iterative solvers.
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...
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.
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
int getNumIters() const
Get the current iteration count.
Traits class which defines basic operations on multivectors.
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero...
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
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
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.
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.