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" 78 template<
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!" );
506 MVT::MvAddMv( one, *cur_soln_vec, alpha[0], *P_, *cur_soln_vec );
507 lp_->updateSolution();
515 MVT::MvAddMv( one, *R_, -alpha[0], *AP_, *R_ );
520 if ( lp_->getLeftPrec() != Teuchos::null ) {
521 lp_->applyLeftPrec( *R_, *Z_ );
522 if ( lp_->getRightPrec() != Teuchos::null ) {
523 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_);
524 lp_->applyRightPrec( *tmp, *Z_ );
527 else if ( lp_->getRightPrec() != Teuchos::null ) {
528 lp_->applyRightPrec( *R_, *Z_ );
531 MVT::Assign( *R_, *Z_ );
534 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() ==
Belos::TwoNorm) {
535 MVT::MvTransMv( one, *R_, *S_, rHs );
539 MVT::MvDot( *R_, *Z_, rHz );
541 beta[0] = rHz[0] / rHz_old[0];
543 MVT::MvAddMv( one, *Z_, beta[0], *P_, *P_ );
548 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
549 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
552 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
554 rHz_old2 = rHz_old[0];
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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...
Pure virtual base class for defining the status testing capabilities of Belos.
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.
bool isInitialized()
States whether the solver has been initialized or not.
OperatorTraits< ScalarType, MV, OP > OPT
Traits class which defines basic operations on multivectors.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
int getNumIters() const
Get the current iteration count.
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.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
SCT::magnitudeType MagnitudeType
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
void resetNumIters(int iter=0)
Reset the iteration count.
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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.