42 #ifndef BELOS_BICGSTAB_ITER_HPP
43 #define BELOS_BICGSTAB_ITER_HPP
86 template <
class ScalarType,
class MV>
104 P(Teuchos::null),
V(Teuchos::null)
112 template<
class ScalarType,
class MV,
class OP>
207 state.
alpha = alpha_;
208 state.
omega = omega_;
252 "Belos::BiCGStabIter::setBlockSize(): Cannot use a block size that is not one.");
262 void axpy(
const ScalarType alpha,
const MV & A,
263 const std::vector<ScalarType> beta,
const MV& B, MV& mv,
bool minus=
false);
307 std::vector<ScalarType> rho_old_, alpha_, omega_;
312 template<
class ScalarType,
class MV,
class OP>
330 template <
class ScalarType,
class MV,
class OP>
337 "Belos::BiCGStabIter::initialize(): Cannot initialize state storage!");
343 int numRHS = MVT::GetNumberVecs(*tmp);
349 R_ = MVT::Clone( *tmp, numRHS_ );
350 Rhat_ = MVT::Clone( *tmp, numRHS_ );
351 P_ = MVT::Clone( *tmp, numRHS_ );
352 V_ = MVT::Clone( *tmp, numRHS_ );
354 rho_old_.resize(numRHS_);
355 alpha_.resize(numRHS_);
356 omega_.resize(numRHS_);
364 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
367 const ScalarType one = SCT::one();
372 std::invalid_argument, errstr );
374 std::invalid_argument, errstr );
377 if (newstate.
R != R_) {
379 MVT::Assign(*newstate.
R, *R_);
383 lp_->computeCurrResVec(R_.get());
389 MVT::Assign(*newstate.
Rhat, *Rhat_);
393 MVT::Assign(*R_, *Rhat_);
399 MVT::Assign(*newstate.
V, *V_);
409 MVT::Assign(*newstate.
P, *P_);
417 if (newstate.
rho_old.size () ==
static_cast<size_t> (numRHS_)) {
423 rho_old_.assign(numRHS_,one);
427 if (newstate.
alpha.size() ==
static_cast<size_t> (numRHS_)) {
429 alpha_ = newstate.
alpha;
433 alpha_.assign(numRHS_,one);
437 if (newstate.
omega.size() ==
static_cast<size_t> (numRHS_)) {
439 omega_ = newstate.
omega;
443 omega_.assign(numRHS_,one);
450 "Belos::BiCGStabIter::initialize(): BiCGStabStateIterState does not have initial residual.");
460 template <
class ScalarType,
class MV,
class OP>
468 if (initialized_ ==
false) {
474 std::vector<ScalarType> rho_new( numRHS_ ), beta( numRHS_ );
475 std::vector<ScalarType> rhatV( numRHS_ ), tT( numRHS_ ), tS( numRHS_ );
478 const ScalarType one = SCT::one();
481 RCP<MV> leftPrecVec, leftPrecVec2;
484 S = MVT::Clone( *R_, numRHS_ );
485 T = MVT::Clone( *R_, numRHS_ );
486 if (lp_->isLeftPrec() || lp_->isRightPrec()) {
487 Y = MVT::Clone( *R_, numRHS_ );
488 Z = MVT::Clone( *R_, numRHS_ );
501 while (stest_->checkStatus(
this) !=
Passed && !breakdown_) {
507 MVT::MvDot(*R_,*Rhat_,rho_new);
511 for(i=0; i<numRHS_; i++) {
514 if (SCT::magnitude(rho_new[i]) < MT::sfmin())
517 beta[i] = (rho_new[i] / rho_old_[i]) * (alpha_[i] / omega_[i]);
523 axpy(one, *P_, omega_, *V_, *P_,
true);
524 axpy(one, *R_, beta, *P_, *P_);
528 if(lp_->isLeftPrec()) {
529 if(lp_->isRightPrec()) {
530 if(leftPrecVec == Teuchos::null) {
531 leftPrecVec = MVT::Clone( *R_, numRHS_ );
533 lp_->applyLeftPrec(*P_,*leftPrecVec);
534 lp_->applyRightPrec(*leftPrecVec,*Y);
537 lp_->applyLeftPrec(*P_,*Y);
540 else if(lp_->isRightPrec()) {
541 lp_->applyRightPrec(*P_,*Y);
545 lp_->applyOp(*Y,*V_);
548 MVT::MvDot(*V_,*Rhat_,rhatV);
549 for(i=0; i<numRHS_; i++) {
550 if (SCT::magnitude(rhatV[i]) < MT::sfmin())
556 alpha_[i] = rho_new[i] / rhatV[i];
560 axpy(one, *R_, alpha_, *V_, *S,
true);
563 if(lp_->isLeftPrec()) {
564 if(lp_->isRightPrec()) {
565 if(leftPrecVec == Teuchos::null) {
566 leftPrecVec = MVT::Clone( *R_, numRHS_ );
568 lp_->applyLeftPrec(*S,*leftPrecVec);
569 lp_->applyRightPrec(*leftPrecVec,*Z);
572 lp_->applyLeftPrec(*S,*Z);
575 else if(lp_->isRightPrec()) {
576 lp_->applyRightPrec(*S,*Z);
583 if(lp_->isLeftPrec()) {
584 if(leftPrecVec == Teuchos::null) {
585 leftPrecVec = MVT::Clone( *R_, numRHS_ );
587 if(leftPrecVec2 == Teuchos::null) {
588 leftPrecVec2 = MVT::Clone( *R_, numRHS_ );
590 lp_->applyLeftPrec(*T,*leftPrecVec2);
591 MVT::MvDot(*leftPrecVec2,*leftPrecVec2,tT);
592 MVT::MvDot(*leftPrecVec,*leftPrecVec2,tS);
595 MVT::MvDot(*T,*T,tT);
596 MVT::MvDot(*S,*T,tS);
598 for(i=0; i<numRHS_; i++) {
599 if (SCT::magnitude(tT[i]) < MT::sfmin())
601 omega_[i] = SCT::zero();
605 omega_[i] = tS[i] / tT[i];
609 axpy(one, *X, alpha_, *Y, *X);
610 axpy(one, *X, omega_, *Z, *X);
613 axpy(one, *S, omega_, *T, *R_,
true);
623 template <
class ScalarType,
class MV,
class OP>
625 const std::vector<ScalarType> beta,
const MV& B, MV& mv,
bool minus)
629 std::vector<int> index(1);
631 for(
int i=0; i<numRHS_; i++) {
633 A1 = MVT::CloneView(A,index);
634 B1 = MVT::CloneView(B,index);
635 mv1 = MVT::CloneViewNonConst(mv,index);
637 MVT::MvAddMv(alpha,*A1,-beta[i],*B1,*mv1);
640 MVT::MvAddMv(alpha,*A1,beta[i],*B1,*mv1);
Collection of types and exceptions used within the Belos solvers.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the pseudo-block BiCGStab iteration, where the basic BiCGStab algorithm is perf...
Class which manages the output and verbosity of the Belos solvers.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > R
The current residual.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > Rhat
The initial residual.
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)
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Structure to contain pointers to BiCGStabIteration state variables.
Declaration of basic traits for the multivector type.
void initializeBiCGStab(BiCGStabIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
void setBlockSize(int blockSize)
Set the blocksize.
void resetNumIters(int iter=0)
Reset the iteration count.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
std::vector< ScalarType > alpha
Traits class which defines basic operations on multivectors.
std::vector< ScalarType > omega
BiCGStabIter(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)
BiCGStabIter constructor with linear problem, solver utilities, and parameter list of solver options...
A linear system to solve, and its associated information.
Teuchos::RCP< const MV > V
A * M * the first decent direction vector.
Class which describes the linear problem to be solved by the iterative solver.
virtual ~BiCGStabIter()
Destructor.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< const MV > P
The first decent direction vector.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void iterate()
This method performs BiCGStab iterations on each linear system until the status test indicates the ne...
std::vector< ScalarType > rho_old
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
SCT::magnitudeType MagnitudeType
MultiVecTraits< ScalarType, MV > MVT
bool breakdownDetected()
Has breakdown been detected in any linear system.
Belos header file which uses auto-configuration information to include necessary C++ headers...
BiCGStabIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...