42 #ifndef BELOS_BICGSTAB_ITER_HPP 43 #define BELOS_BICGSTAB_ITER_HPP 60 #include "Teuchos_BLAS.hpp" 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Teuchos_SerialDenseVector.hpp" 63 #include "Teuchos_ScalarTraits.hpp" 64 #include "Teuchos_ParameterList.hpp" 65 #include "Teuchos_TimeMonitor.hpp" 87 template <
class ScalarType,
class MV>
91 Teuchos::RCP<const MV>
R;
94 Teuchos::RCP<const MV>
Rhat;
97 Teuchos::RCP<const MV>
P;
100 Teuchos::RCP<const MV>
V;
105 P(Teuchos::null),
V(Teuchos::null)
113 template<
class ScalarType,
class MV,
class OP>
123 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
137 Teuchos::ParameterList ¶ms );
207 state.
alpha = alpha_;
208 state.
omega = omega_;
248 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
249 "Belos::BiCGStabIter::setBlockSize(): Cannot use a block size that is not one.");
259 void axpy(
const ScalarType alpha,
const MV & A,
260 const std::vector<ScalarType> beta,
const MV& B, MV& mv,
bool minus=
false);
265 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
266 const Teuchos::RCP<OutputManager<ScalarType> > om_;
267 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
290 Teuchos::RCP<MV> Rhat_;
301 std::vector<ScalarType> rho_old_, alpha_, omega_;
306 template<
class ScalarType,
class MV,
class OP>
310 Teuchos::ParameterList & ):
323 template <
class ScalarType,
class MV,
class OP>
327 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
328 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
329 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
330 "Belos::BiCGStabIter::initialize(): Cannot initialize state storage!");
333 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
336 int numRHS = MVT::GetNumberVecs(*tmp);
341 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
342 R_ = MVT::Clone( *tmp, numRHS_ );
343 Rhat_ = MVT::Clone( *tmp, numRHS_ );
344 P_ = MVT::Clone( *tmp, numRHS_ );
345 V_ = MVT::Clone( *tmp, numRHS_ );
347 rho_old_.resize(numRHS_);
348 alpha_.resize(numRHS_);
349 omega_.resize(numRHS_);
354 std::string errstr(
"Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
357 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
359 if (!Teuchos::is_null(newstate.
R)) {
361 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
362 std::invalid_argument, errstr );
363 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != numRHS_,
364 std::invalid_argument, errstr );
367 if (newstate.
R != R_) {
369 MVT::Assign(*newstate.
R, *R_);
373 lp_->computeCurrResVec(R_.get());
377 if (!Teuchos::is_null(newstate.
Rhat) && newstate.
Rhat != Rhat_) {
379 MVT::Assign(*newstate.
Rhat, *Rhat_);
383 MVT::Assign(*R_, *Rhat_);
387 if (!Teuchos::is_null(newstate.
V) && newstate.
V != V_) {
389 MVT::Assign(*newstate.
V, *V_);
397 if (!Teuchos::is_null(newstate.
P) && newstate.
P != P_) {
399 MVT::Assign(*newstate.
P, *P_);
407 if (newstate.
rho_old.size () ==
static_cast<size_t> (numRHS_)) {
413 rho_old_.assign(numRHS_,one);
417 if (newstate.
alpha.size() ==
static_cast<size_t> (numRHS_)) {
419 alpha_ = newstate.
alpha;
423 alpha_.assign(numRHS_,one);
427 if (newstate.
omega.size() ==
static_cast<size_t> (numRHS_)) {
429 omega_ = newstate.
omega;
433 omega_.assign(numRHS_,one);
439 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.
R),std::invalid_argument,
440 "Belos::BiCGStabIter::initialize(): BiCGStabStateIterState does not have initial residual.");
450 template <
class ScalarType,
class MV,
class OP>
458 if (initialized_ ==
false) {
464 std::vector<ScalarType> rho_new( numRHS_ ), beta( numRHS_ );
465 std::vector<ScalarType> rhatV( numRHS_ ), tT( numRHS_ ), tS( numRHS_ );
468 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
471 RCP<MV> leftPrecVec, leftPrecVec2;
474 S = MVT::Clone( *R_, numRHS_ );
475 T = MVT::Clone( *R_, numRHS_ );
476 if (lp_->isLeftPrec() || lp_->isRightPrec()) {
477 Y = MVT::Clone( *R_, numRHS_ );
478 Z = MVT::Clone( *R_, numRHS_ );
486 Teuchos::RCP<MV> X = lp_->getCurrLHSVec();
491 while (stest_->checkStatus(
this) !=
Passed) {
503 MVT::MvDot(*R_,*Rhat_,rho_new);
511 for(i=0; i<numRHS_; i++) {
512 beta[i] = (rho_new[i] / rho_old_[i]) * (alpha_[i] / omega_[i]);
519 axpy(one, *P_, omega_, *V_, *P_,
true);
520 axpy(one, *R_, beta, *P_, *P_);
524 if(lp_->isLeftPrec()) {
525 if(lp_->isRightPrec()) {
526 if(leftPrecVec == Teuchos::null) {
527 leftPrecVec = MVT::Clone( *R_, numRHS_ );
529 lp_->applyLeftPrec(*P_,*leftPrecVec);
530 lp_->applyRightPrec(*leftPrecVec,*Y);
533 lp_->applyLeftPrec(*P_,*Y);
536 else if(lp_->isRightPrec()) {
537 lp_->applyRightPrec(*P_,*Y);
541 lp_->applyOp(*Y,*V_);
544 MVT::MvDot(*V_,*Rhat_,rhatV);
545 for(i=0; i<numRHS_; i++) {
546 alpha_[i] = rho_new[i] / rhatV[i];
554 axpy(one, *R_, alpha_, *V_, *S,
true);
557 if(lp_->isLeftPrec()) {
558 if(lp_->isRightPrec()) {
559 if(leftPrecVec == Teuchos::null) {
560 leftPrecVec = MVT::Clone( *R_, numRHS_ );
562 lp_->applyLeftPrec(*S,*leftPrecVec);
563 lp_->applyRightPrec(*leftPrecVec,*Z);
566 lp_->applyLeftPrec(*S,*Z);
569 else if(lp_->isRightPrec()) {
570 lp_->applyRightPrec(*S,*Z);
580 if(lp_->isLeftPrec()) {
581 if(leftPrecVec == Teuchos::null) {
582 leftPrecVec = MVT::Clone( *R_, numRHS_ );
584 if(leftPrecVec2 == Teuchos::null) {
585 leftPrecVec2 = MVT::Clone( *R_, numRHS_ );
587 lp_->applyLeftPrec(*T,*leftPrecVec2);
588 MVT::MvDot(*leftPrecVec2,*leftPrecVec2,tT);
589 MVT::MvDot(*leftPrecVec,*leftPrecVec2,tS);
592 MVT::MvDot(*T,*T,tT);
593 MVT::MvDot(*T,*S,tS);
595 for(i=0; i<numRHS_; i++) {
596 omega_[i] = tS[i] / tT[i];
604 axpy(one, *X, alpha_, *Y, *X);
605 axpy(one, *X, omega_, *Z, *X);
608 axpy(one, *S, omega_, *T, *R_,
true);
618 template <
class ScalarType,
class MV,
class OP>
620 const std::vector<ScalarType> beta,
const MV& B, MV& mv,
bool minus)
622 Teuchos::RCP<const MV> A1, B1;
623 Teuchos::RCP<MV> mv1;
624 std::vector<int> index(1);
626 for(
int i=0; i<numRHS_; i++) {
628 A1 = MVT::CloneView(A,index);
629 B1 = MVT::CloneView(B,index);
630 mv1 = MVT::CloneViewNonConst(mv,index);
632 MVT::MvAddMv(alpha,*A1,-beta[i],*B1,*mv1);
635 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::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
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 ...