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.");
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);
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);
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.