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 );
190 initializeBiCGStab(empty);
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 ¶ms ):
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 );
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)) {
362 std::invalid_argument, errstr );
364 std::invalid_argument, errstr );
367 if (newstate.
R != R_) {
373 lp_->computeCurrResVec(R_.get());
377 if (!Teuchos::is_null(newstate.
Rhat) && newstate.
Rhat != Rhat_) {
387 if (!Teuchos::is_null(newstate.
V) && newstate.
V != V_) {
397 if (!Teuchos::is_null(newstate.
P) && 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;
476 if (lp_->isLeftPrec() || lp_->isRightPrec()) {
486 Teuchos::RCP<MV> X = lp_->getCurrLHSVec();
491 while (stest_->checkStatus(
this) !=
Passed) {
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) {
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_);
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) {
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) {
584 if(leftPrecVec2 == Teuchos::null) {
587 lp_->applyLeftPrec(*T,*leftPrecVec2);
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++) {
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
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
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...
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.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
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.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
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
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void Assign(const MV &A, MV &mv)
mv := A
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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...
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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
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 ...