10 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP
11 #define BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP
22 #include "BelosPseudoBlockCGIter.hpp"
41 template<
class Storage,
class MV,
class OP>
43 virtual public CGIteration<Sacado::MP::Vector<Storage>, MV, OP> {
51 typedef MultiVecTraits<ScalarType,MV>
MVT;
52 typedef OperatorTraits<ScalarType,MV,OP>
OPT;
67 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
131 auto state =
Teuchos::rcp(
new PseudoBlockCGIterationState<ScalarType,MV>());
140 auto s = Teuchos::rcp_dynamic_cast<PseudoBlockCGIterationState<ScalarType,MV> >(state,
true);
174 const LinearProblem<ScalarType,MV,OP>&
getProblem()
const {
return *lp_; }
182 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
199 if (static_cast<size_type> (iter_) >= diag_.size ()) {
202 return diag_ (0, iter_);
214 if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
217 return offdiag_ (0, iter_);
277 template<
class StorageType,
class MV,
class OP>
280 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
281 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
289 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness",
true) ),
290 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
291 breakDownTol_(params.get(
"Ensemble Breakdown Tolerance", 0.0))
298 template <
class StorageType,
class MV,
class OP>
299 void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
306 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
312 int numRHS = MVT::GetNumberVecs(*tmp);
316 if (!Teuchos::rcp_dynamic_cast<PseudoBlockCGIterationState<ScalarType,MV> >(newstate,
true)->matches(tmp, numRHS_))
317 newstate->initialize(tmp, numRHS_);
321 if(numEntriesForCondEst_ > 0) {
322 diag_.resize(numEntriesForCondEst_);
323 offdiag_.resize(numEntriesForCondEst_-1);
326 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
331 std::invalid_argument, errstr );
333 std::invalid_argument, errstr );
338 MVT::Assign( *R_0, *R_ );
345 lp_->applyLeftPrec( *R_, *Z_ );
348 lp_->applyRightPrec( *Z_, *tmp1 );
353 lp_->applyRightPrec( *R_, *Z_ );
356 MVT::Assign( *R_, *Z_ );
358 MVT::Assign( *Z_, *P_ );
368 template <
class StorageType,
class MV,
class OP>
369 void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
375 if (initialized_ ==
false) {
381 std::vector<int> index(1);
382 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
390 ScalarType pAp_old=one, beta_old=one ,rHz_old2=one;
396 MVT::MvDot( *R_, *Z_, rHz );
398 if ( assertPositiveDefiniteness_ )
399 for (i=0; i<numRHS_; ++i)
402 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
407 while (stest_->checkStatus(
this) != Passed) {
413 lp_->applyOp( *P_, *AP_ );
416 MVT::MvDot( *P_, *AP_, pAp );
418 for (i=0; i<numRHS_; ++i) {
419 if ( assertPositiveDefiniteness_ )
423 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
425 const int ensemble_size = pAp[i].size();
426 for (
int k=0; k<ensemble_size; ++k) {
427 if (SVT::magnitude(pAp[i].coeff(k)) <= breakDownTol_)
428 alpha(i,i).fastAccessCoeff(k) = 0.0;
430 alpha(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / pAp[i].coeff(k);
437 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
438 lp_->updateSolution();
442 for (i=0; i<numRHS_; ++i) {
448 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
454 lp_->applyLeftPrec( *R_, *Z_ );
457 lp_->applyRightPrec( *Z_, *tmp );
462 lp_->applyRightPrec( *R_, *Z_ );
468 MVT::MvDot( *R_, *Z_, rHz );
469 if ( assertPositiveDefiniteness_ )
470 for (i=0; i<numRHS_; ++i)
473 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
476 for (i=0; i<numRHS_; ++i) {
477 const int ensemble_size = rHz[i].size();
478 for (
int k=0; k<ensemble_size; ++k) {
479 if (SVT::magnitude(rHz[i].coeff(k)) <= breakDownTol_)
480 beta(i,i).fastAccessCoeff(k) = 0.0;
482 beta(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / rHz_old[i].coeff(k);
487 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
499 rHz_old2 = rHz_old[0];
500 beta_old = beta(0,0);
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
bool assertPositiveDefiniteness_
Teuchos::ScalarTraits< typename Storage::value_type > SVT
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(T a)
void setState(Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > state)
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
OperatorTraits< ScalarType, MV, OP > OPT
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
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
virtual ~PseudoBlockCGIter()
Destructor.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
Teuchos::ScalarTraits< ScalarType > SCT
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setBlockSize(int blockSize)
Set the blocksize.
int getNumIters() const
Get the current iteration count.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< CGIterationStateBase< ScalarType, MV > > getState() const
Get the current state of the linear solver.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
MultiVecTraits< ScalarType, MV > MVT
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int numEntriesForCondEst_
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
SCT::magnitudeType MagnitudeType
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Sacado::MP::Vector< Storage > ScalarType
bool isInitialized()
States whether the solver has been initialized or not.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
SVT::magnitudeType breakDownTol_
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::ArrayRCP< MagnitudeType > offdiag_
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_