42 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_HPP
43 #define BELOS_PSEUDO_BLOCK_CG_ITER_HPP
79 template<
class ScalarType,
class MV,
class OP>
210 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
305 template<
class ScalarType,
class MV,
class OP>
316 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
317 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
325 template <
class ScalarType,
class MV,
class OP>
332 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
338 int numRHS = MVT::GetNumberVecs(*tmp);
344 R_ = MVT::Clone( *tmp, numRHS_ );
345 Z_ = MVT::Clone( *tmp, numRHS_ );
346 P_ = MVT::Clone( *tmp, numRHS_ );
347 AP_ = MVT::Clone( *tmp, numRHS_ );
351 if(numEntriesForCondEst_ > 0) {
352 diag_.resize(numEntriesForCondEst_);
353 offdiag_.resize(numEntriesForCondEst_-1);
358 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
367 std::invalid_argument, errstr );
369 std::invalid_argument, errstr );
372 if (newstate.
R != R_) {
374 MVT::MvAddMv( one, *newstate.
R, zero, *newstate.
R, *R_ );
381 lp_->applyLeftPrec( *R_, *Z_ );
384 lp_->applyRightPrec( *Z_, *tmp1 );
389 lp_->applyRightPrec( *R_, *Z_ );
394 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
399 "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
409 template <
class ScalarType,
class MV,
class OP>
415 if (initialized_ ==
false) {
421 std::vector<int> index(1);
422 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ ), beta( numRHS_ );
433 MVT::MvDot( *R_, *Z_, rHz );
435 if ( assertPositiveDefiniteness_ )
436 for (i=0; i<numRHS_; ++i)
439 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
444 while (stest_->checkStatus(
this) !=
Passed) {
450 lp_->applyOp( *P_, *AP_ );
453 MVT::MvDot( *P_, *AP_, pAp );
455 for (i=0; i<numRHS_; ++i) {
456 if ( assertPositiveDefiniteness_ )
460 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
462 alpha(i,i) = rHz[i] / pAp[i];
468 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
469 lp_->updateSolution();
473 for (i=0; i<numRHS_; ++i) {
479 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
485 lp_->applyLeftPrec( *R_, *Z_ );
488 lp_->applyRightPrec( *Z_, *tmp );
493 lp_->applyRightPrec( *R_, *Z_ );
499 MVT::MvDot( *R_, *Z_, rHz );
500 if ( assertPositiveDefiniteness_ )
501 for (i=0; i<numRHS_; ++i)
504 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
507 for (i=0; i<numRHS_; ++i) {
508 beta[i] = rHz[i] / rHz_old[i];
512 MVT::MvAddMv( one, *Z_i, beta[i], *P_i, *P_i );
516 if (doCondEst_ && (iter_ - 1) < diag_.size()) {
524 rHz_old2_ = rHz_old[0];
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
int numEntriesForCondEst_
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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.
virtual ~PseudoBlockCGIter()
Destructor.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
bool is_null(const std::shared_ptr< T > &p)
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::ArrayRCP< MagnitudeType > offdiag_
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(T a)
Declaration of basic traits for the multivector type.
void setBlockSize(int blockSize)
Set the blocksize.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
OperatorTraits< ScalarType, MV, OP > OPT
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero...
void resetNumIters(int iter=0)
Reset the iteration count.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
SCT::magnitudeType MagnitudeType
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
MultiVecTraits< ScalarType, MV > MVT
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
bool assertPositiveDefiniteness_
A linear system to solve, and its associated information.
int getNumIters() const
Get the current iteration count.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::ScalarTraits< ScalarType > SCT
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
PseudoBlockCGIter(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)
PseudoBlockCGIter constructor with linear problem, solver utilities, and parameter list of solver opt...
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...
const Teuchos::RCP< OutputManager< ScalarType > > om_
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ArrayRCP< MagnitudeType > diag_