42 #ifndef BELOS_CG_ITER_HPP 43 #define BELOS_CG_ITER_HPP 59 #include "Teuchos_SerialDenseMatrix.hpp" 60 #include "Teuchos_SerialDenseVector.hpp" 61 #include "Teuchos_ScalarTraits.hpp" 62 #include "Teuchos_ParameterList.hpp" 63 #include "Teuchos_TimeMonitor.hpp" 77 template<
class ScalarType,
class MV,
class OP>
87 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
101 Teuchos::ParameterList ¶ms );
199 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
200 "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
217 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
218 if (static_cast<size_type> (
iter_) >=
diag_.size ()) {
232 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
253 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
254 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
255 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
302 template<
class ScalarType,
class MV,
class OP>
306 Teuchos::ParameterList ¶ms ):
311 stateStorageInitialized_(false),
313 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
314 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
321 template <
class ScalarType,
class MV,
class OP>
324 if (!stateStorageInitialized_) {
327 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
328 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
329 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
330 stateStorageInitialized_ =
false;
337 if (R_ == Teuchos::null) {
339 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
340 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
341 "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
342 R_ = MVT::Clone( *tmp, 1 );
343 Z_ = MVT::Clone( *tmp, 1 );
344 P_ = MVT::Clone( *tmp, 1 );
345 AP_ = MVT::Clone( *tmp, 1 );
349 if(numEntriesForCondEst_ > 0) {
350 diag_.resize(numEntriesForCondEst_);
351 offdiag_.resize(numEntriesForCondEst_-1);
355 stateStorageInitialized_ =
true;
363 template <
class ScalarType,
class MV,
class OP>
367 if (!stateStorageInitialized_)
370 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
371 "Belos::CGIter::initialize(): Cannot initialize state storage!");
375 std::string errstr(
"Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
378 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
379 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
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::MvAddMv( one, *newstate.
R, zero, *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 );
405 else if ( lp_->getRightPrec() != Teuchos::null ) {
406 lp_->applyRightPrec( *R_, *Z_ );
411 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
415 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
416 "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
426 template <
class ScalarType,
class MV,
class OP>
432 if (initialized_ ==
false) {
437 std::vector<ScalarType> alpha(1);
438 std::vector<ScalarType> beta(1);
439 std::vector<ScalarType> rHz(1), rHz_old(1), pAp(1);
442 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
443 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
446 ScalarType pAp_old = one, beta_old = one, rHz_old2 = one;
449 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
452 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
453 "Belos::CGIter::iterate(): current linear system has more than one vector!" );
456 MVT::MvDot( *R_, *Z_, rHz );
461 while (stest_->checkStatus(
this) !=
Passed) {
467 lp_->applyOp( *P_, *AP_ );
470 MVT::MvDot( *P_, *AP_, pAp );
471 alpha[0] = rHz[0] / pAp[0];
474 if(assertPositiveDefiniteness_)
476 "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
480 MVT::MvAddMv( one, *cur_soln_vec, alpha[0], *P_, *cur_soln_vec );
481 lp_->updateSolution();
489 MVT::MvAddMv( one, *R_, -alpha[0], *AP_, *R_ );
494 if ( lp_->getLeftPrec() != Teuchos::null ) {
495 lp_->applyLeftPrec( *R_, *Z_ );
496 if ( lp_->getRightPrec() != Teuchos::null ) {
497 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
498 lp_->applyRightPrec( *Z_, *tmp );
502 else if ( lp_->getRightPrec() != Teuchos::null ) {
503 lp_->applyRightPrec( *R_, *Z_ );
509 MVT::MvDot( *R_, *Z_, rHz );
511 beta[0] = rHz[0] / rHz_old[0];
513 MVT::MvAddMv( one, *Z_, beta[0], *P_, *P_ );
518 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
519 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
522 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
524 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.
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::ArrayRCP< MagnitudeType > diag_
Teuchos::ScalarTraits< ScalarType > SCT
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Class which manages the output and verbosity of the Belos solvers.
bool assertPositiveDefiniteness_
void setStateSize()
Method for initalizing the state storage needed by CG.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Teuchos::ArrayRCP< MagnitudeType > offdiag_
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.
MultiVecTraits< ScalarType, MV > MVT
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
CGIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, Teuchos::ParameterList ¶ms)
CGIter constructor with linear problem, solver utilities, and parameter list of solver options...
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::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
bool stateStorageInitialized_
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.
SCT::magnitudeType MagnitudeType
int numEntriesForCondEst_
void resetNumIters(int iter=0)
Reset the iteration count.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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.