42 #ifndef BELOS_CG_SINGLE_RED_ITER_HPP
43 #define BELOS_CG_SINGLE_RED_ITER_HPP
78 template<
class ScalarType,
class MV,
class OP>
202 "Belos::CGSingleRedIter::setBlockSize(): Cannot use a block size that is not one.");
252 bool stateStorageInitialized_;
258 bool foldConvergenceDetectionIntoAllreduce_;
294 template<
class ScalarType,
class MV,
class OP>
303 convTest_(convTester),
305 stateStorageInitialized_(false),
308 foldConvergenceDetectionIntoAllreduce_ = params.
get<
bool>(
"Fold Convergence Detection Into Allreduce",
false);
313 template <
class ScalarType,
class MV,
class OP>
316 if (!stateStorageInitialized_) {
321 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
322 stateStorageInitialized_ =
false;
329 if (R_ == Teuchos::null) {
333 "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
336 W_ = MVT::Clone( *tmp, 3 );
337 std::vector<int> index2(2,0);
338 std::vector<int> index(1,0);
343 S_ = MVT::CloneViewNonConst( *W_, index2 );
348 U_ = MVT::CloneViewNonConst( *W_, index2 );
351 R_ = MVT::CloneViewNonConst( *W_, index );
353 AZ_ = MVT::CloneViewNonConst( *W_, index );
355 Z_ = MVT::CloneViewNonConst( *W_, index );
360 T_ = MVT::CloneViewNonConst( *W_, index2 );
363 V_ = MVT::Clone( *tmp, 2 );
365 AP_ = MVT::CloneViewNonConst( *V_, index );
367 P_ = MVT::CloneViewNonConst( *V_, index );
372 stateStorageInitialized_ =
true;
380 template <
class ScalarType,
class MV,
class OP>
384 if (!stateStorageInitialized_)
388 "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
392 std::string errstr(
"Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
394 if (newstate.
R != Teuchos::null) {
397 std::invalid_argument, errstr );
399 std::invalid_argument, errstr );
402 if (newstate.
R != R_) {
404 MVT::Assign( *newstate.
R, *R_ );
410 if ( lp_->getLeftPrec() != Teuchos::null ) {
411 lp_->applyLeftPrec( *R_, *Z_ );
412 if ( lp_->getRightPrec() != Teuchos::null ) {
414 lp_->applyRightPrec( *Z_, *tmp );
415 MVT::Assign( *tmp, *Z_ );
418 else if ( lp_->getRightPrec() != Teuchos::null ) {
419 lp_->applyRightPrec( *R_, *Z_ );
422 MVT::Assign( *R_, *Z_ );
426 lp_->applyOp( *Z_, *AZ_ );
430 MVT::Assign( *U_, *V_);
435 "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
445 template <
class ScalarType,
class MV,
class OP>
450 return Teuchos::null;
451 }
else if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
453 return Teuchos::null;
461 template <
class ScalarType,
class MV,
class OP>
467 if (initialized_ ==
false) {
474 ScalarType rHz_old, alpha, beta, delta;
485 "Belos::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
487 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
489 MVT::MvTransMv( one, *S_, *T_, sHt );
495 MVT::MvTransMv( one, *S_, *Z_, sHz );
500 (stest_->checkStatus(
this) ==
Passed))
502 alpha = rHz_ / delta;
506 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
511 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
519 MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
523 MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
527 if ( lp_->getLeftPrec() != Teuchos::null ) {
528 lp_->applyLeftPrec( *R_, *Z_ );
529 if ( lp_->getRightPrec() != Teuchos::null ) {
531 lp_->applyRightPrec( *tmp, *Z_ );
534 else if ( lp_->getRightPrec() != Teuchos::null ) {
535 lp_->applyRightPrec( *R_, *Z_ );
538 MVT::Assign( *R_, *Z_ );
542 lp_->applyOp( *Z_, *AZ_ );
545 MVT::MvTransMv( one, *S_, *T_, sHt );
558 if (stest_->checkStatus(
this) ==
Passed) {
562 beta = rHz_ / rHz_old;
563 alpha = rHz_ / (delta - (beta*rHz_ / alpha));
567 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
574 MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
585 MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
589 MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
596 if (stest_->checkStatus(
this) ==
Passed) {
602 if ( lp_->getLeftPrec() != Teuchos::null ) {
603 lp_->applyLeftPrec( *R_, *Z_ );
604 if ( lp_->getRightPrec() != Teuchos::null ) {
606 lp_->applyRightPrec( *tmp, *Z_ );
609 else if ( lp_->getRightPrec() != Teuchos::null ) {
610 lp_->applyRightPrec( *R_, *Z_ );
613 MVT::Assign( *R_, *Z_ );
617 lp_->applyOp( *Z_, *AZ_ );
620 MVT::MvTransMv( one, *S_, *Z_, sHz );
627 beta = rHz_ / rHz_old;
628 alpha = rHz_ / (delta - (beta*rHz_ / alpha));
632 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
639 MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
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'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 'alpha = p^H*A*P' 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 ¶ms)
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)