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_;
275 Teuchos::RCP<MV> AZ_;
285 Teuchos::RCP<MV> R2_;
291 Teuchos::RCP<MV> AP_;
297 template<
class ScalarType,
class MV,
class OP>
302 Teuchos::ParameterList ¶ms ):
306 convTest_(convTester),
308 stateStorageInitialized_(false),
311 foldConvergenceDetectionIntoAllreduce_ = params.get<
bool>(
"Fold Convergence Detection Into Allreduce",
false);
316 template <
class ScalarType,
class MV,
class OP>
319 if (!stateStorageInitialized_) {
322 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
323 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
324 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
325 stateStorageInitialized_ =
false;
332 if (R_ == Teuchos::null) {
334 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
335 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
336 "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
337 S_ = MVT::Clone( *tmp, 2 );
338 if (foldConvergenceDetectionIntoAllreduce_) {
339 T_ = MVT::Clone( *tmp, 2 );
341 std::vector<int> index(1,0);
342 Z_ = MVT::CloneViewNonConst( *T_, index );
344 R2_ = MVT::CloneViewNonConst( *T_, index );
347 Z_ = MVT::Clone( *tmp, 1 );
348 P_ = MVT::Clone( *tmp, 1 );
349 AP_ = MVT::Clone( *tmp, 1 );
352 std::vector<int> index(1,0);
353 R_ = MVT::CloneViewNonConst( *S_, index );
355 AZ_ = MVT::CloneViewNonConst( *S_, index );
359 stateStorageInitialized_ =
true;
367 template <
class ScalarType,
class MV,
class OP>
371 if (!stateStorageInitialized_)
374 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
375 "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
379 std::string errstr(
"Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
381 if (newstate.
R != Teuchos::null) {
383 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
384 std::invalid_argument, errstr );
385 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
386 std::invalid_argument, errstr );
389 if (newstate.
R != R_) {
391 MVT::Assign( *newstate.
R, *R_ );
397 if ( lp_->getLeftPrec() != Teuchos::null ) {
398 lp_->applyLeftPrec( *R_, *Z_ );
399 if ( lp_->getRightPrec() != Teuchos::null ) {
400 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
401 lp_->applyRightPrec( *Z_, *tmp );
402 MVT::Assign( *tmp, *Z_ );
405 else if ( lp_->getRightPrec() != Teuchos::null ) {
406 lp_->applyRightPrec( *R_, *Z_ );
409 MVT::Assign( *R_, *Z_ );
411 MVT::Assign( *Z_, *P_ );
414 lp_->applyOp( *Z_, *AZ_ );
417 MVT::Assign( *AZ_, *AP_ );
421 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
422 "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
432 template <
class ScalarType,
class MV,
class OP>
433 Teuchos::RCP<const MV>
436 (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHz_));
437 return Teuchos::null;
438 }
else if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
439 (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHr_));
440 return Teuchos::null;
448 template <
class ScalarType,
class MV,
class OP>
454 if (initialized_ ==
false) {
459 Teuchos::SerialDenseMatrix<int,ScalarType> sHz( 2, 1 );
460 Teuchos::SerialDenseMatrix<int,ScalarType> sHt( 2, 2 );
461 ScalarType rHz_old, alpha, beta, delta;
464 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
465 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
468 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
471 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
472 "Belos::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
474 if (foldConvergenceDetectionIntoAllreduce_) {
476 MVT::Assign( *R_, *R2_ );
477 MVT::MvTransMv( one, *S_, *T_, sHt );
483 MVT::MvTransMv( one, *S_, *Z_, sHz );
487 if ((Teuchos::ScalarTraits<ScalarType>::magnitude(delta) < Teuchos::ScalarTraits<ScalarType>::eps()) &&
488 (stest_->checkStatus(
this) ==
Passed))
490 alpha = rHz_ / delta;
494 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
499 if (foldConvergenceDetectionIntoAllreduce_) {
507 MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
511 MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
515 if ( lp_->getLeftPrec() != Teuchos::null ) {
516 lp_->applyLeftPrec( *R_, *Z_ );
517 if ( lp_->getRightPrec() != Teuchos::null ) {
518 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
519 lp_->applyRightPrec( *Z_, *tmp );
520 MVT::Assign( *tmp, *Z_ );
523 else if ( lp_->getRightPrec() != Teuchos::null ) {
524 lp_->applyRightPrec( *R_, *Z_ );
527 MVT::Assign( *R_, *Z_ );
531 lp_->applyOp( *Z_, *AZ_ );
534 MVT::Assign( *R_, *R2_ );
535 MVT::MvTransMv( one, *S_, *T_, sHt );
545 if (stest_->checkStatus(
this) ==
Passed) {
551 beta = rHz_ / rHz_old;
552 alpha = rHz_ / (delta - (beta*rHz_ / alpha));
556 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
560 MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
566 MVT::MvAddMv( one, *AZ_, beta, *AP_, *AP_ );
577 MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
581 MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
585 if (stest_->checkStatus(
this) ==
Passed) {
593 if ( lp_->getLeftPrec() != Teuchos::null ) {
594 lp_->applyLeftPrec( *R_, *Z_ );
595 if ( lp_->getRightPrec() != Teuchos::null ) {
596 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
597 lp_->applyRightPrec( *Z_, *tmp );
601 else if ( lp_->getRightPrec() != Teuchos::null ) {
602 lp_->applyRightPrec( *R_, *Z_ );
609 lp_->applyOp( *Z_, *AZ_ );
612 MVT::MvTransMv( one, *S_, *Z_, sHz );
619 beta = rHz_ / rHz_old;
620 alpha = rHz_ / (delta - (beta*rHz_ / alpha));
624 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
628 MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
634 MVT::MvAddMv( one, *AZ_, beta, *AP_, *AP_ );
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)