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.");
219 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
220 const Teuchos::RCP<OutputManager<ScalarType> > om_;
221 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
234 bool stateStorageInitialized_;
252 Teuchos::RCP<MV> AP_;
258 template<
class ScalarType,
class MV,
class OP>
262 Teuchos::ParameterList ¶ms ):
267 stateStorageInitialized_(false),
274 template <
class ScalarType,
class MV,
class OP>
277 if (!stateStorageInitialized_) {
280 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
281 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
282 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
283 stateStorageInitialized_ =
false;
290 if (R_ == Teuchos::null) {
292 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
293 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
294 "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
295 R_ = MVT::Clone( *tmp, 1 );
296 Z_ = MVT::Clone( *tmp, 1 );
297 P_ = MVT::Clone( *tmp, 1 );
298 AP_ = MVT::Clone( *tmp, 1 );
302 stateStorageInitialized_ =
true;
310 template <
class ScalarType,
class MV,
class OP>
314 if (!stateStorageInitialized_)
317 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
318 "Belos::CGIter::initialize(): Cannot initialize state storage!");
322 std::string errstr(
"Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
325 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
326 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
328 if (newstate.
R != Teuchos::null) {
330 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
331 std::invalid_argument, errstr );
332 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
333 std::invalid_argument, errstr );
336 if (newstate.
R != R_) {
338 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
344 if ( lp_->getLeftPrec() != Teuchos::null ) {
345 lp_->applyLeftPrec( *R_, *Z_ );
346 if ( lp_->getRightPrec() != Teuchos::null ) {
347 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
348 lp_->applyRightPrec( *Z_, *tmp );
352 else if ( lp_->getRightPrec() != Teuchos::null ) {
353 lp_->applyRightPrec( *R_, *Z_ );
358 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
362 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
363 "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
373 template <
class ScalarType,
class MV,
class OP>
379 if (initialized_ ==
false) {
384 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
385 Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
386 Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
389 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
390 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
393 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
396 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
CGIterateFailure,
397 "Belos::CGIter::iterate(): current linear system has more than one vector!" );
400 MVT::MvTransMv( one, *R_, *Z_, rHz );
405 while (stest_->checkStatus(
this) !=
Passed) {
411 lp_->applyOp( *P_, *AP_ );
414 MVT::MvTransMv( one, *P_, *AP_, pAp );
415 alpha(0,0) = rHz(0,0) / pAp(0,0);
418 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero,
CGIterateFailure,
419 "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
423 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
424 lp_->updateSolution();
428 rHz_old(0,0) = rHz(0,0);
432 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
437 if ( lp_->getLeftPrec() != Teuchos::null ) {
438 lp_->applyLeftPrec( *R_, *Z_ );
439 if ( lp_->getRightPrec() != Teuchos::null ) {
440 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
441 lp_->applyRightPrec( *Z_, *tmp );
445 else if ( lp_->getRightPrec() != Teuchos::null ) {
446 lp_->applyRightPrec( *R_, *Z_ );
452 MVT::MvTransMv( one, *R_, *Z_, rHz );
454 beta(0,0) = rHz(0,0) / rHz_old(0,0);
456 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.
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.
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.
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< 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
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.