50 #ifndef BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 51 #define BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 70 #include "Teuchos_BLAS.hpp" 71 #include "Teuchos_ScalarTraits.hpp" 72 #include "Teuchos_ParameterList.hpp" 73 #include "Teuchos_TimeMonitor.hpp" 92 template <
class ScalarType,
class MV>
95 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
99 Teuchos::RCP<const MV>
W;
100 Teuchos::RCP<const MV>
U;
101 Teuchos::RCP<const MV>
AU;
103 Teuchos::RCP<const MV>
D;
104 Teuchos::RCP<const MV>
V;
110 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
146 template <
class ScalarType,
class MV,
class OP>
154 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
164 Teuchos::ParameterList ¶ms );
237 state.
alpha = alpha_;
241 state.
theta = theta_;
282 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
283 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
297 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
298 const Teuchos::RCP<OutputManager<ScalarType> > om_;
299 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
309 std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
310 std::vector<MagnitudeType> tau_, theta_;
327 Teuchos::RCP<MV> U_, AU_;
328 Teuchos::RCP<MV> Rtilde_;
331 Teuchos::RCP<MV> solnUpdate_;
342 template <
class ScalarType,
class MV,
class OP>
346 Teuchos::ParameterList ¶ms
359 template <
class ScalarType,
class MV,
class OP>
360 Teuchos::RCP<const MV>
363 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
366 if ((
int) normvec->size() < numRHS_)
367 normvec->resize( numRHS_ );
370 for (
int i=0; i<numRHS_; i++) {
371 (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
375 return Teuchos::null;
380 template <
class ScalarType,
class MV,
class OP>
384 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
385 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
386 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
389 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Rtilde == Teuchos::null,std::invalid_argument,
390 "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
393 int numRHS = MVT::GetNumberVecs(*newstate.
Rtilde);
398 if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
401 if ( Teuchos::is_null(D_) )
402 D_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
403 MVT::MvInit( *D_, STzero );
406 if ( Teuchos::is_null(solnUpdate_) )
407 solnUpdate_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
408 MVT::MvInit( *solnUpdate_, STzero );
411 if (newstate.
Rtilde != Rtilde_)
412 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
413 W_ = MVT::CloneCopy( *Rtilde_ );
414 U_ = MVT::CloneCopy( *Rtilde_ );
415 V_ = MVT::Clone( *Rtilde_, numRHS_ );
419 lp_->apply( *U_, *V_ );
420 AU_ = MVT::CloneCopy( *V_ );
423 alpha_.resize( numRHS_, STone );
424 eta_.resize( numRHS_, STzero );
425 rho_.resize( numRHS_ );
426 rho_old_.resize( numRHS_ );
427 tau_.resize( numRHS_ );
428 theta_.resize( numRHS_, MTzero );
430 MVT::MvNorm( *Rtilde_, tau_ );
431 MVT::MvDot( *Rtilde_, *Rtilde_, rho_ );
436 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
437 W_ = MVT::CloneCopy( *newstate.
W );
438 U_ = MVT::CloneCopy( *newstate.
U );
439 AU_ = MVT::CloneCopy( *newstate.
AU );
440 D_ = MVT::CloneCopy( *newstate.
D );
441 V_ = MVT::CloneCopy( *newstate.
V );
445 solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
446 MVT::MvInit( *solnUpdate_, STzero );
449 alpha_ = newstate.
alpha;
453 theta_ = newstate.
theta;
463 template <
class ScalarType,
class MV,
class OP>
469 if (initialized_ ==
false) {
474 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
475 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
476 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
477 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
478 std::vector< ScalarType > beta(numRHS_,STzero);
479 std::vector<int> index(1);
487 while (stest_->checkStatus(
this) !=
Passed) {
489 for (
int iIter=0; iIter<2; iIter++)
497 MVT::MvDot( *V_, *Rtilde_, alpha_ );
498 for (
int i=0; i<numRHS_; i++) {
499 rho_old_[i] = rho_[i];
500 alpha_[i] = rho_old_[i]/alpha_[i];
508 for (
int i=0; i<numRHS_; ++i) {
517 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
518 Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
519 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
526 Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
527 Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
528 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
539 Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
540 Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
541 MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
550 lp_->apply( *U_, *AU_ );
557 MVT::MvNorm( *W_, theta_ );
559 bool breakdown=
false;
560 for (
int i=0; i<numRHS_; ++i) {
561 theta_[i] /= tau_[i];
563 MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
564 tau_[i] *= theta_[i]*cs;
565 if ( tau_[i] == MTzero ) {
568 eta_[i] = cs*cs*alpha_[i];
575 for (
int i=0; i<numRHS_; ++i) {
577 Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
578 Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
579 MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
596 MVT::MvDot( *W_, *Rtilde_, rho_ );
598 for (
int i=0; i<numRHS_; ++i) {
599 beta[i] = rho_[i]/rho_old_[i];
608 Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
609 Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
610 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i );
613 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
614 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
615 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
619 lp_->apply( *U_, *AU_ );
622 for (
int i=0; i<numRHS_; ++i) {
624 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
625 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
626 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
639 #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.