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 );
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_;
276 Teuchos::ArrayRCP<MagnitudeType>
diag_, offdiag_;
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>
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>
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);
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::ScalarTraits< ScalarType > SCT
void resetNumIters(int iter=0)
Reset the iteration count.
int getNumIters() const
Get the current iteration count.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Sacado::MP::Vector< Storage > ScalarType
Teuchos::ScalarTraits< typename Storage::value_type > SVT
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...
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.
MultiVecTraits< ScalarType, MV > MVT
const Teuchos::RCP< OutputManager< ScalarType > > om_
bool assertPositiveDefiniteness_
SCT::magnitudeType MagnitudeType
virtual ~PseudoBlockCGIter()
Destructor.
Teuchos::ArrayRCP< MagnitudeType > diag_
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
int numEntriesForCondEst_
bool isInitialized()
States whether the solver has been initialized or not.
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::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
SVT::magnitudeType breakDownTol_
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)