42 #ifndef BELOS_BICGSTAB_ITER_HPP 43 #define BELOS_BICGSTAB_ITER_HPP 60 #include "Teuchos_SerialDenseMatrix.hpp" 61 #include "Teuchos_SerialDenseVector.hpp" 62 #include "Teuchos_ScalarTraits.hpp" 63 #include "Teuchos_ParameterList.hpp" 64 #include "Teuchos_TimeMonitor.hpp" 86 template <
class ScalarType,
class MV>
90 Teuchos::RCP<const MV>
R;
93 Teuchos::RCP<const MV>
Rhat;
96 Teuchos::RCP<const MV>
P;
99 Teuchos::RCP<const MV>
V;
104 P(Teuchos::null),
V(Teuchos::null)
112 template<
class ScalarType,
class MV,
class OP>
122 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
124 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
137 Teuchos::ParameterList ¶ms );
207 state.
alpha = alpha_;
208 state.
omega = omega_;
251 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
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);
268 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
269 const Teuchos::RCP<OutputManager<ScalarType> > om_;
270 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
296 Teuchos::RCP<MV> Rhat_;
307 std::vector<ScalarType> rho_old_, alpha_, omega_;
312 template<
class ScalarType,
class MV,
class OP>
316 Teuchos::ParameterList & ):
330 template <
class ScalarType,
class MV,
class OP>
334 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
335 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
336 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
337 "Belos::BiCGStabIter::initialize(): Cannot initialize state storage!");
340 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
343 int numRHS = MVT::GetNumberVecs(*tmp);
348 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
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();
369 if (!Teuchos::is_null(newstate.
R)) {
371 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
372 std::invalid_argument, errstr );
373 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != numRHS_,
374 std::invalid_argument, errstr );
377 if (newstate.
R != R_) {
379 MVT::Assign(*newstate.
R, *R_);
383 lp_->computeCurrResVec(R_.get());
387 if (!Teuchos::is_null(newstate.
Rhat) && newstate.
Rhat != Rhat_) {
389 MVT::Assign(*newstate.
Rhat, *Rhat_);
393 MVT::Assign(*R_, *Rhat_);
397 if (!Teuchos::is_null(newstate.
V) && newstate.
V != V_) {
399 MVT::Assign(*newstate.
V, *V_);
407 if (!Teuchos::is_null(newstate.
P) && newstate.
P != P_) {
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);
449 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.
R),std::invalid_argument,
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_ );
496 Teuchos::RCP<MV> X = lp_->getCurrLHSVec();
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(*T,*S,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)
627 Teuchos::RCP<const MV> A1, B1;
628 Teuchos::RCP<MV> mv1;
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.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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...
Pure virtual base class for defining the status testing capabilities of Belos.
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.
int getNumIters() const
Get the current iteration count.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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 > P
The first decent direction vector.
BiCGStabIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
SCT::magnitudeType MagnitudeType
MultiVecTraits< ScalarType, MV > MVT
bool breakdownDetected()
Has breakdown been detected in any linear system.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Belos header file which uses auto-configuration information to include necessary C++ headers...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...