10 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_HPP
11 #define BELOS_PSEUDO_BLOCK_CG_ITER_HPP
55 template <
class ScalarType,
class MV>
69 this->
R = MVT::Clone( *tmp, _numVectors );
70 this->
Z = MVT::Clone( *tmp, _numVectors );
71 this->
P = MVT::Clone( *tmp, _numVectors );
72 this->
AP = MVT::Clone(*tmp, _numVectors );
82 template<
class ScalarType,
class MV,
class OP>
220 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
315 template<
class ScalarType,
class MV,
class OP>
326 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness", true) ),
327 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
335 template <
class ScalarType,
class MV,
class OP>
342 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
348 int numRHS = MVT::GetNumberVecs(*tmp);
354 newstate->initialize(tmp, numRHS_);
358 if(numEntriesForCondEst_ > 0) {
359 diag_.resize(numEntriesForCondEst_);
360 offdiag_.resize(numEntriesForCondEst_-1);
363 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
368 std::invalid_argument, errstr );
370 std::invalid_argument, errstr );
375 MVT::Assign( *R_0, *R_ );
382 lp_->applyLeftPrec( *R_, *Z_ );
385 lp_->applyRightPrec( *Z_, *tmp1 );
390 lp_->applyRightPrec( *R_, *Z_ );
393 MVT::Assign( *R_, *Z_ );
395 MVT::Assign( *Z_, *P_ );
405 template <
class ScalarType,
class MV,
class OP>
417 std::vector<int> index(1);
418 std::vector<ScalarType> rHz( numRHS_ );
419 std::vector<ScalarType> rHz_old( numRHS_ );
420 std::vector<ScalarType> pAp( numRHS_ );
421 std::vector<ScalarType> beta( numRHS_ );
432 MVT::MvDot( *R_, *Z_, rHz );
434 if ( assertPositiveDefiniteness_ )
435 for (i=0; i<numRHS_; ++i)
438 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
443 while (stest_->checkStatus(
this) !=
Passed) {
449 lp_->applyOp( *P_, *AP_ );
452 MVT::MvDot( *P_, *AP_, pAp );
454 for (i=0; i<numRHS_; ++i) {
455 if ( assertPositiveDefiniteness_ )
459 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
461 alpha(i,i) = rHz[i] / pAp[i];
467 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
468 lp_->updateSolution();
472 for (i=0; i<numRHS_; ++i) {
478 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
484 lp_->applyLeftPrec( *R_, *Z_ );
487 lp_->applyRightPrec( *Z_, *tmp );
492 lp_->applyRightPrec( *R_, *Z_ );
498 MVT::MvDot( *R_, *Z_, rHz );
499 if ( assertPositiveDefiniteness_ )
500 for (i=0; i<numRHS_; ++i)
503 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
506 for (i=0; i<numRHS_; ++i) {
507 beta[i] = rHz[i] / rHz_old[i];
511 MVT::MvAddMv( one, *Z_i, beta[i], *P_i, *P_i );
515 if (doCondEst_ && (iter_ - 1) < diag_.size()) {
523 rHz_old2_ = rHz_old[0];
void initialize(Teuchos::RCP< const MV > tmp, int _numVectors)
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...
PseudoBlockCGIterationState()=default
Class which manages the output and verbosity of the Belos solvers.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
virtual ~PseudoBlockCGIter()=default
Destructor.
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_
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::ArrayRCP< MagnitudeType > offdiag_
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setState(Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > state)
static magnitudeType real(T a)
Declaration of basic traits for the multivector type.
typename SCT::magnitudeType MagnitudeType
void setBlockSize(int blockSize)
Set the blocksize.
bool matches(Teuchos::RCP< const MV > tmp, int _numVectors=1) const
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< MV > P
The current decent direction vector.
Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > getState() const
Get the current state of the linear solver.
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...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Structure to contain pointers to PseudoBlockCGIteration state variables.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
virtual void initialize(Teuchos::RCP< const MV > tmp, int _numVectors)
bool assertPositiveDefiniteness_
A linear system to solve, and its associated information.
void initializeCG(Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > newstate, Teuchos::RCP< MV > R_0)
Initialize the solver to an iterate, providing a complete state.
int getNumIters() const
Get the current iteration count.
virtual ~PseudoBlockCGIterationState()=default
Class which describes the linear problem to be solved by the iterative solver.
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< 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::RCP< MV > AP
The matrix A applied to current decent direction vector.
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.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< MV > R
The current residual.
Belos header file which uses auto-configuration information to include necessary C++ headers...
const Teuchos::RCP< OutputManager< ScalarType > > om_
virtual bool matches(Teuchos::RCP< const MV > tmp, int _numVectors=1) const
PseudoBlockCGIterationState(Teuchos::RCP< const MV > tmp)
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ArrayRCP< MagnitudeType > diag_