351 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
352 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
353 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
356 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Rtilde == Teuchos::null,std::invalid_argument,
357 "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
360 int numRHS = MVT::GetNumberVecs(*newstate.
Rtilde);
365 if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
368 if ( Teuchos::is_null(D_) )
369 D_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
370 MVT::MvInit( *D_, STzero );
373 if ( Teuchos::is_null(solnUpdate_) )
374 solnUpdate_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
375 MVT::MvInit( *solnUpdate_, STzero );
378 if (newstate.
Rtilde != Rtilde_)
379 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
380 W_ = MVT::CloneCopy( *Rtilde_ );
381 U_ = MVT::CloneCopy( *Rtilde_ );
382 V_ = MVT::Clone( *Rtilde_, numRHS_ );
386 lp_->apply( *U_, *V_ );
387 AU_ = MVT::CloneCopy( *V_ );
390 alpha_.resize( numRHS_, STone );
391 eta_.resize( numRHS_, STzero );
392 rho_.resize( numRHS_ );
393 rho_old_.resize( numRHS_ );
394 tau_.resize( numRHS_ );
395 theta_.resize( numRHS_, MTzero );
397 MVT::MvNorm( *Rtilde_, tau_ );
398 MVT::MvDot( *Rtilde_, *Rtilde_, rho_ );
403 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
404 W_ = MVT::CloneCopy( *newstate.
W );
405 U_ = MVT::CloneCopy( *newstate.
U );
406 AU_ = MVT::CloneCopy( *newstate.
AU );
407 D_ = MVT::CloneCopy( *newstate.
D );
408 V_ = MVT::CloneCopy( *newstate.
V );
412 solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
413 MVT::MvInit( *solnUpdate_, STzero );
416 alpha_ = newstate.
alpha;
420 theta_ = newstate.
theta;
436 if (initialized_ ==
false) {
441 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
442 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
443 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
444 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
445 std::vector< ScalarType > beta(numRHS_,STzero);
446 std::vector<int> index(1);
454 while (stest_->checkStatus(
this) !=
Passed) {
456 for (
int iIter=0; iIter<2; iIter++)
464 MVT::MvDot( *V_, *Rtilde_, alpha_ );
465 for (
int i=0; i<numRHS_; i++) {
466 rho_old_[i] = rho_[i];
467 alpha_[i] = rho_old_[i]/alpha_[i];
475 for (
int i=0; i<numRHS_; ++i) {
484 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
485 Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
486 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
493 Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
494 Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
495 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
506 Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
507 Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
508 MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
517 lp_->apply( *U_, *AU_ );
524 MVT::MvNorm( *W_, theta_ );
526 bool breakdown=
false;
527 for (
int i=0; i<numRHS_; ++i) {
528 theta_[i] /= tau_[i];
530 MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
531 tau_[i] *= theta_[i]*cs;
532 if ( tau_[i] == MTzero ) {
535 eta_[i] = cs*cs*alpha_[i];
542 for (
int i=0; i<numRHS_; ++i) {
544 Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
545 Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
546 MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
563 MVT::MvDot( *W_, *Rtilde_, rho_ );
565 for (
int i=0; i<numRHS_; ++i) {
566 beta[i] = rho_[i]/rho_old_[i];
575 Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
576 Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
577 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i );
580 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
581 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
582 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
586 lp_->apply( *U_, *AU_ );
589 for (
int i=0; i<numRHS_; ++i) {
591 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
592 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
593 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
PseudoBlockTFQMRIter(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)
Belos::PseudoBlockTFQMRIter constructor.