50 #ifndef BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 51 #define BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 70 #include "Teuchos_ScalarTraits.hpp" 71 #include "Teuchos_ParameterList.hpp" 72 #include "Teuchos_TimeMonitor.hpp" 91 template <
class ScalarType,
class MV>
94 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
98 Teuchos::RCP<const MV>
W;
99 Teuchos::RCP<const MV>
U;
100 Teuchos::RCP<const MV>
AU;
102 Teuchos::RCP<const MV>
D;
103 Teuchos::RCP<const MV>
V;
109 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
145 template <
class ScalarType,
class MV,
class OP>
153 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
163 Teuchos::ParameterList ¶ms );
236 state.
alpha = alpha_;
240 state.
theta = theta_;
281 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
282 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
296 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
297 const Teuchos::RCP<OutputManager<ScalarType> > om_;
298 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
308 std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
309 std::vector<MagnitudeType> tau_, theta_;
326 Teuchos::RCP<MV> U_, AU_;
327 Teuchos::RCP<MV> Rtilde_;
330 Teuchos::RCP<MV> solnUpdate_;
341 template <
class ScalarType,
class MV,
class OP>
345 Teuchos::ParameterList &
358 template <
class ScalarType,
class MV,
class OP>
359 Teuchos::RCP<const MV>
362 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
365 if ((
int) normvec->size() < numRHS_)
366 normvec->resize( numRHS_ );
369 for (
int i=0; i<numRHS_; i++) {
370 (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
374 return Teuchos::null;
379 template <
class ScalarType,
class MV,
class OP>
383 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
384 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
385 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
388 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Rtilde == Teuchos::null,std::invalid_argument,
389 "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
392 int numRHS = MVT::GetNumberVecs(*newstate.
Rtilde);
397 if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
400 if ( Teuchos::is_null(D_) )
401 D_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
402 MVT::MvInit( *D_, STzero );
405 if ( Teuchos::is_null(solnUpdate_) )
406 solnUpdate_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
407 MVT::MvInit( *solnUpdate_, STzero );
410 if (newstate.
Rtilde != Rtilde_)
411 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
412 W_ = MVT::CloneCopy( *Rtilde_ );
413 U_ = MVT::CloneCopy( *Rtilde_ );
414 V_ = MVT::Clone( *Rtilde_, numRHS_ );
418 lp_->apply( *U_, *V_ );
419 AU_ = MVT::CloneCopy( *V_ );
422 alpha_.resize( numRHS_, STone );
423 eta_.resize( numRHS_, STzero );
424 rho_.resize( numRHS_ );
425 rho_old_.resize( numRHS_ );
426 tau_.resize( numRHS_ );
427 theta_.resize( numRHS_, MTzero );
429 MVT::MvNorm( *Rtilde_, tau_ );
430 MVT::MvDot( *Rtilde_, *Rtilde_, rho_ );
435 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
436 W_ = MVT::CloneCopy( *newstate.
W );
437 U_ = MVT::CloneCopy( *newstate.
U );
438 AU_ = MVT::CloneCopy( *newstate.
AU );
439 D_ = MVT::CloneCopy( *newstate.
D );
440 V_ = MVT::CloneCopy( *newstate.
V );
444 solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
445 MVT::MvInit( *solnUpdate_, STzero );
448 alpha_ = newstate.
alpha;
452 theta_ = newstate.
theta;
462 template <
class ScalarType,
class MV,
class OP>
468 if (initialized_ ==
false) {
473 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
474 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
475 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
476 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
477 std::vector< ScalarType > beta(numRHS_,STzero);
478 std::vector<int> index(1);
486 while (stest_->checkStatus(
this) !=
Passed) {
488 for (
int iIter=0; iIter<2; iIter++)
496 MVT::MvDot( *V_, *Rtilde_, alpha_ );
497 for (
int i=0; i<numRHS_; i++) {
498 rho_old_[i] = rho_[i];
499 alpha_[i] = rho_old_[i]/alpha_[i];
507 for (
int i=0; i<numRHS_; ++i) {
516 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
517 Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
518 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
525 Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
526 Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
527 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
538 Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
539 Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
540 MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
549 lp_->apply( *U_, *AU_ );
556 MVT::MvNorm( *W_, theta_ );
558 bool breakdown=
false;
559 for (
int i=0; i<numRHS_; ++i) {
560 theta_[i] /= tau_[i];
562 MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
563 tau_[i] *= theta_[i]*cs;
564 if ( tau_[i] == MTzero ) {
567 eta_[i] = cs*cs*alpha_[i];
574 for (
int i=0; i<numRHS_; ++i) {
576 Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
577 Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
578 MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
595 MVT::MvDot( *W_, *Rtilde_, rho_ );
597 for (
int i=0; i<numRHS_; ++i) {
598 beta[i] = rho_[i]/rho_old_[i];
607 Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
608 Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
609 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i );
612 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
613 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
614 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
618 lp_->apply( *U_, *AU_ );
621 for (
int i=0; i<numRHS_; ++i) {
623 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
624 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
625 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
638 #endif // BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 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...
OperatorTraits< ScalarType, MV, OP > OPT
Class which manages the output and verbosity of the Belos solvers.
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< const MV > Rtilde
SCT::magnitudeType MagnitudeType
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
PseudoBlockTFQMRIterateFailure is thrown when the PseudoBlockTFQMRIter object is unable to compute th...
Structure to contain pointers to PseudoBlockTFQMRIter state variables.
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
Declaration of basic traits for the multivector type.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< const MV > AU
Pure virtual base class which describes the basic interface to the linear solver iteration.
std::vector< ScalarType > alpha
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
PseudoBlockTFQMRIterInitFailure(const std::string &what_arg)
Traits class which defines basic operations on multivectors.
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.
PseudoBlockTFQMRIterState()
PseudoBlockTFQMRIterInitFailure is thrown when the PseudoBlockTFQMRIter object is unable to generate ...
Teuchos::ScalarTraits< ScalarType > SCT
A linear system to solve, and its associated information.
std::vector< MagnitudeType > theta
Class which describes the linear problem to be solved by the iterative solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::RCP< const MV > D
std::vector< ScalarType > eta
Teuchos::RCP< const MV > U
std::vector< MagnitudeType > tau
PseudoBlockTFQMRIterateFailure(const std::string &what_arg)
virtual ~PseudoBlockTFQMRIter()
Belos::PseudoBlockTFQMRIter destructor.
Teuchos::RCP< const MV > W
The current residual basis.
std::vector< ScalarType > rho
Class which defines basic traits for the operator type.
Teuchos::ScalarTraits< ScalarType > SCT
void iterate()
This method performs pseudo-block TFQMR iterations until the status test indicates the need to stop o...
Parent class to all Belos exceptions.
void resetNumIters(int iter=0)
Reset the iteration count.
Belos header file which uses auto-configuration information to include necessary C++ headers...
SCT::magnitudeType MagnitudeType
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > V
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.