42 #ifndef BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP 43 #define BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP 54 #include "BelosPseudoBlockCGIter.hpp" 73 template<
class Storage,
class MV,
class OP>
75 virtual public CGIteration<Sacado::MP::Vector<Storage>, MV, OP> {
83 typedef MultiVecTraits<ScalarType,MV>
MVT;
84 typedef OperatorTraits<ScalarType,MV,OP>
OPT;
85 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
87 typedef Teuchos::ScalarTraits<typename Storage::value_type>
SVT;
98 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
99 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
100 Teuchos::ParameterList ¶ms );
145 void initializeCG(CGIterationState<ScalarType,MV>& newstate);
152 CGIterationState<ScalarType,MV> empty;
164 CGIterationState<ScalarType,MV> state;
199 const LinearProblem<ScalarType,MV,OP>&
getProblem()
const {
return *lp_; }
206 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
207 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
223 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
224 if (static_cast<size_type> (iter_) >= diag_.size ()) {
227 return diag_ (0, iter_);
238 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
239 if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
242 return offdiag_ (0, iter_);
251 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
252 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
253 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
302 template<
class StorageType,
class MV,
class OP>
304 PseudoBlockCGIter(
const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
305 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
306 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
307 Teuchos::ParameterList ¶ms ):
314 assertPositiveDefiniteness_( params.get(
"Assert Positive Definiteness",
true) ),
315 numEntriesForCondEst_(params.get(
"Max Size For Condest",0) ),
316 breakDownTol_(params.get(
"Ensemble Breakdown Tolerance", 0.0))
323 template <
class StorageType,
class MV,
class OP>
324 void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
325 initializeCG(CGIterationState<ScalarType,MV>& newstate)
328 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
329 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
330 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
331 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
334 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
337 int numRHS = MVT::GetNumberVecs(*tmp);
342 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
343 R_ = MVT::Clone( *tmp, numRHS_ );
344 Z_ = MVT::Clone( *tmp, numRHS_ );
345 P_ = MVT::Clone( *tmp, numRHS_ );
346 AP_ = MVT::Clone( *tmp, numRHS_ );
350 if(numEntriesForCondEst_ > 0) {
351 diag_.resize(numEntriesForCondEst_);
352 offdiag_.resize(numEntriesForCondEst_-1);
357 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
360 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
361 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
363 if (!Teuchos::is_null(newstate.R)) {
365 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
366 std::invalid_argument, errstr );
367 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
368 std::invalid_argument, errstr );
371 if (newstate.R != R_) {
373 MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
379 if ( lp_->getLeftPrec() != Teuchos::null ) {
380 lp_->applyLeftPrec( *R_, *Z_ );
381 if ( lp_->getRightPrec() != Teuchos::null ) {
382 Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
383 lp_->applyRightPrec( *Z_, *tmp1 );
387 else if ( lp_->getRightPrec() != Teuchos::null ) {
388 lp_->applyRightPrec( *R_, *Z_ );
393 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
397 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
398 "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
408 template <
class StorageType,
class MV,
class OP>
409 void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
415 if (initialized_ ==
false) {
421 std::vector<int> index(1);
422 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
423 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
426 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
427 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
430 ScalarType pAp_old=one, beta_old=one ,rHz_old2=one;
433 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
436 MVT::MvDot( *R_, *Z_, rHz );
438 if ( assertPositiveDefiniteness_ )
439 for (i=0; i<numRHS_; ++i)
440 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
442 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
447 while (stest_->checkStatus(
this) != Passed) {
453 lp_->applyOp( *P_, *AP_ );
456 MVT::MvDot( *P_, *AP_, pAp );
458 for (i=0; i<numRHS_; ++i) {
459 if ( assertPositiveDefiniteness_ )
461 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
463 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
465 const int ensemble_size = pAp[i].size();
466 for (
int k=0; k<ensemble_size; ++k) {
467 if (SVT::magnitude(pAp[i].coeff(k)) <= breakDownTol_)
468 alpha(i,i).fastAccessCoeff(k) = 0.0;
470 alpha(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / pAp[i].coeff(k);
477 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
478 lp_->updateSolution();
482 for (i=0; i<numRHS_; ++i) {
488 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
493 if ( lp_->getLeftPrec() != Teuchos::null ) {
494 lp_->applyLeftPrec( *R_, *Z_ );
495 if ( lp_->getRightPrec() != Teuchos::null ) {
496 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
497 lp_->applyRightPrec( *Z_, *tmp );
501 else if ( lp_->getRightPrec() != Teuchos::null ) {
502 lp_->applyRightPrec( *R_, *Z_ );
508 MVT::MvDot( *R_, *Z_, rHz );
509 if ( assertPositiveDefiniteness_ )
510 for (i=0; i<numRHS_; ++i)
511 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
513 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
516 for (i=0; i<numRHS_; ++i) {
517 const int ensemble_size = rHz[i].size();
518 for (
int k=0; k<ensemble_size; ++k) {
519 if (SVT::magnitude(rHz[i].coeff(k)) <= breakDownTol_)
520 beta(i,i).fastAccessCoeff(k) = 0.0;
522 beta(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / rHz_old[i].coeff(k);
525 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
526 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
527 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
533 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
534 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (
sqrt( rHz_old[0] * rHz_old2)));
537 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
539 rHz_old2 = rHz_old[0];
540 beta_old = beta(0,0);
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
bool assertPositiveDefiniteness_
Teuchos::ScalarTraits< typename Storage::value_type > SVT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
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.
virtual ~PseudoBlockCGIter()
Destructor.
int getNumIters() const
Get the current iteration count.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void setBlockSize(int blockSize)
Set the blocksize.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
MultiVecTraits< ScalarType, MV > MVT
int numEntriesForCondEst_
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
SCT::magnitudeType MagnitudeType
CGIterationState< 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.
Sacado::MP::Vector< Storage > ScalarType
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
bool isInitialized()
States whether the solver has been initialized or not.
SVT::magnitudeType breakDownTol_
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::ArrayRCP< MagnitudeType > offdiag_
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_