42 #ifndef BELOS_CG_SINGLE_RED_ITER_HPP 43 #define BELOS_CG_SINGLE_RED_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" 78 template<
class ScalarType,
class MV,
class OP>
88 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
103 Teuchos::ParameterList ¶ms );
201 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
202 "Belos::CGSingleRedIter::setBlockSize(): Cannot use a block size that is not one.");
213 Teuchos::ArrayView<MagnitudeType> temp;
219 Teuchos::ArrayView<MagnitudeType> temp;
236 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
237 const Teuchos::RCP<OutputManager<ScalarType> > om_;
238 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
239 const Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
252 bool stateStorageInitialized_;
258 bool foldConvergenceDetectionIntoAllreduce_;
273 Teuchos::RCP<MV> AZ_;
277 Teuchos::RCP<MV> AP_;
294 template<
class ScalarType,
class MV,
class OP>
299 Teuchos::ParameterList ¶ms ):
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_) {
319 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
320 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
321 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
322 stateStorageInitialized_ =
false;
329 if (R_ == Teuchos::null) {
331 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
332 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
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_)
387 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
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) {
396 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
397 std::invalid_argument, errstr );
398 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
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 ) {
413 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
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_);
434 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
435 "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
445 template <
class ScalarType,
class MV,
class OP>
446 Teuchos::RCP<const MV>
449 (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHz_));
450 return Teuchos::null;
451 }
else if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
452 (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHr_));
453 return Teuchos::null;
461 template <
class ScalarType,
class MV,
class OP>
467 if (initialized_ ==
false) {
472 Teuchos::SerialDenseMatrix<int,ScalarType> sHz( 2, 1 );
473 Teuchos::SerialDenseMatrix<int,ScalarType> sHt( 2, 2 );
474 ScalarType rHz_old, alpha, beta, delta;
477 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
478 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
481 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
484 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
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 );
499 if ((Teuchos::ScalarTraits<ScalarType>::magnitude(delta) < Teuchos::ScalarTraits<ScalarType>::eps()) &&
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 ) {
530 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
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 ) {
605 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
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.
int getNumIters() const
Get the current iteration count.
Collection of types and exceptions used within the Belos solvers.
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.
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Pure virtual base class for defining the status testing capabilities of Belos.
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.
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.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation (NOT_IMPLEMENTED)
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
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.
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)