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.");
215 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
216 if (static_cast<size_type> (
iter_) >=
diag_.size ()) {
230 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
251 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
252 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
253 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
300 template<
class ScalarType,
class MV,
class OP>
304 Teuchos::ParameterList ¶ms ):
309 stateStorageInitialized_(false),
311 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
312 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) )
318 template <
class ScalarType,
class MV,
class OP>
321 if (!stateStorageInitialized_) {
324 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
325 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
326 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
327 stateStorageInitialized_ =
false;
334 if (R_ == Teuchos::null) {
336 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
337 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
338 "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
339 R_ = MVT::Clone( *tmp, 1 );
340 Z_ = MVT::Clone( *tmp, 1 );
341 P_ = MVT::Clone( *tmp, 1 );
342 AP_ = MVT::Clone( *tmp, 1 );
346 if(numEntriesForCondEst_ > 0) {
347 diag_.resize(numEntriesForCondEst_);
348 offdiag_.resize(numEntriesForCondEst_-1);
352 stateStorageInitialized_ =
true;
360 template <
class ScalarType,
class MV,
class OP>
364 if (!stateStorageInitialized_)
367 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
368 "Belos::CGIter::initialize(): Cannot initialize state storage!");
372 std::string errstr(
"Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
375 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
376 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
378 if (newstate.
R != Teuchos::null) {
380 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
381 std::invalid_argument, errstr );
382 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
383 std::invalid_argument, errstr );
386 if (newstate.
R != R_) {
388 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
394 if ( lp_->getLeftPrec() != Teuchos::null ) {
395 lp_->applyLeftPrec( *R_, *Z_ );
396 if ( lp_->getRightPrec() != Teuchos::null ) {
397 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
398 lp_->applyRightPrec( *Z_, *tmp );
402 else if ( lp_->getRightPrec() != Teuchos::null ) {
403 lp_->applyRightPrec( *R_, *Z_ );
408 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
412 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
413 "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
423 template <
class ScalarType,
class MV,
class OP>
429 if (initialized_ ==
false) {
434 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
435 Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
436 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
439 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
440 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
443 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
446 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
447 "Belos::CGIter::iterate(): current linear system has more than one vector!" );
450 MVT::MvTransMv( one, *R_, *Z_, rHz );
455 while (stest_->checkStatus(
this) !=
Passed) {
461 lp_->applyOp( *P_, *AP_ );
464 MVT::MvTransMv( one, *P_, *AP_, pAp );
465 alpha(0,0) = rHz(0,0) / pAp(0,0);
468 if(assertPositiveDefiniteness_)
469 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero,
CGIterateFailure,
470 "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
474 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
475 lp_->updateSolution();
479 rHz_old(0,0) = rHz(0,0);
483 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
488 if ( lp_->getLeftPrec() != Teuchos::null ) {
489 lp_->applyLeftPrec( *R_, *Z_ );
490 if ( lp_->getRightPrec() != Teuchos::null ) {
491 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
492 lp_->applyRightPrec( *Z_, *tmp );
496 else if ( lp_->getRightPrec() != Teuchos::null ) {
497 lp_->applyRightPrec( *R_, *Z_ );
503 MVT::MvTransMv( one, *R_, *Z_, rHz );
505 beta(0,0) = rHz(0,0) / rHz_old(0,0);
507 MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
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.
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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.