42 #ifndef BELOS_CG_SINGLE_RED_ITER_HPP 43 #define BELOS_CG_SINGLE_RED_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::CGSingleRedIter::setBlockSize(): Cannot use a block size that is not one.");
211 Teuchos::ArrayView<MagnitudeType> temp;
217 Teuchos::ArrayView<MagnitudeType> temp;
234 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
235 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
236 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
280 template<
class ScalarType,
class MV,
class OP>
284 Teuchos::ParameterList ¶ms ):
289 stateStorageInitialized_(false),
296 template <
class ScalarType,
class MV,
class OP>
299 if (!stateStorageInitialized_) {
302 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
303 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
304 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
305 stateStorageInitialized_ =
false;
312 if (R_ == Teuchos::null) {
314 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
315 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
316 "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
317 S_ = MVT::Clone( *tmp, 2 );
318 Z_ = MVT::Clone( *tmp, 1 );
319 P_ = MVT::Clone( *tmp, 1 );
320 AP_ = MVT::Clone( *tmp, 1 );
323 std::vector<int> index(1,0);
324 R_ = MVT::CloneViewNonConst( *S_, index );
326 AZ_ = MVT::CloneViewNonConst( *S_, index );
330 stateStorageInitialized_ =
true;
338 template <
class ScalarType,
class MV,
class OP>
342 if (!stateStorageInitialized_)
345 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
346 "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
350 std::string errstr(
"Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
352 if (newstate.
R != Teuchos::null) {
354 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
355 std::invalid_argument, errstr );
356 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
357 std::invalid_argument, errstr );
360 if (newstate.
R != R_) {
362 MVT::Assign( *newstate.
R, *R_ );
368 if ( lp_->getLeftPrec() != Teuchos::null ) {
369 lp_->applyLeftPrec( *R_, *Z_ );
370 if ( lp_->getRightPrec() != Teuchos::null ) {
371 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
372 lp_->applyRightPrec( *Z_, *tmp );
376 else if ( lp_->getRightPrec() != Teuchos::null ) {
377 lp_->applyRightPrec( *R_, *Z_ );
382 MVT::Assign( *Z_, *P_ );
385 lp_->applyOp( *Z_, *AZ_ );
388 MVT::Assign( *AZ_, *AP_ );
392 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
393 "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
403 template <
class ScalarType,
class MV,
class OP>
409 if (initialized_ ==
false) {
414 Teuchos::SerialDenseMatrix<int,ScalarType> sHz( 2, 1 );
415 ScalarType rHz, rHz_old, alpha, beta, delta;
418 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
419 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
422 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
425 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
426 "Belos::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
429 MVT::MvTransMv( one, *S_, *Z_, sHz );
436 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
445 MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
446 lp_->updateSolution();
450 MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
454 if (stest_->checkStatus(
this) ==
Passed) {
462 if ( lp_->getLeftPrec() != Teuchos::null ) {
463 lp_->applyLeftPrec( *R_, *Z_ );
464 if ( lp_->getRightPrec() != Teuchos::null ) {
465 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
466 lp_->applyRightPrec( *Z_, *tmp );
470 else if ( lp_->getRightPrec() != Teuchos::null ) {
471 lp_->applyRightPrec( *R_, *Z_ );
478 lp_->applyOp( *Z_, *AZ_ );
481 MVT::MvTransMv( one, *S_, *Z_, sHz );
488 beta = rHz / rHz_old;
489 alpha = rHz / (delta - (beta*rHz / alpha));
493 "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
497 MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
503 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.
void setStateSize()
Method for initalizing the state storage needed by CG.
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Declaration of basic traits for the multivector type.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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.
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, Teuchos::ParameterList ¶ms)
CGSingleRedIter constructor with linear problem, solver utilities, and parameter list of solver optio...
const Teuchos::RCP< OutputManager< ScalarType > > om_
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
bool stateStorageInitialized_
Class which defines basic traits for the operator type.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
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.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation (NOT_IMPLEMENTED)