50 #ifndef BELOS_TFQMR_ITER_HPP 51 #define BELOS_TFQMR_ITER_HPP 70 #include "Teuchos_BLAS.hpp" 71 #include "Teuchos_SerialDenseMatrix.hpp" 72 #include "Teuchos_SerialDenseVector.hpp" 73 #include "Teuchos_ScalarTraits.hpp" 74 #include "Teuchos_ParameterList.hpp" 75 #include "Teuchos_TimeMonitor.hpp" 94 template <
class ScalarType,
class MV>
98 Teuchos::RCP<const MV>
R;
99 Teuchos::RCP<const MV>
W;
100 Teuchos::RCP<const MV>
U;
102 Teuchos::RCP<const MV>
D;
103 Teuchos::RCP<const MV>
V;
106 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
142 template <
class ScalarType,
class MV,
class OP>
150 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
160 Teuchos::ParameterList ¶ms );
229 state.solnUpdate = solnUpdate_;
269 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
270 "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
290 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
291 const Teuchos::RCP<OutputManager<ScalarType> > om_;
292 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
300 std::vector<ScalarType> alpha_, rho_, rho_old_;
301 std::vector<MagnitudeType> tau_, cs_, theta_;
314 bool stateStorageInitialized_;
324 Teuchos::RCP<MV> U_, AU_;
325 Teuchos::RCP<MV> Rtilde_;
328 Teuchos::RCP<MV> solnUpdate_;
338 template <
class ScalarType,
class MV,
class OP>
342 Teuchos::ParameterList ¶ms
354 stateStorageInitialized_(false),
361 template <
class ScalarType,
class MV,
class OP>
362 Teuchos::RCP<const MV>
365 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
367 (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
369 return Teuchos::null;
375 template <
class ScalarType,
class MV,
class OP>
378 if (!stateStorageInitialized_) {
381 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
382 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
383 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
384 stateStorageInitialized_ =
false;
391 if (R_ == Teuchos::null) {
393 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
394 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
395 "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
396 R_ = MVT::Clone( *tmp, 1 );
397 D_ = MVT::Clone( *tmp, 1 );
398 V_ = MVT::Clone( *tmp, 1 );
399 solnUpdate_ = MVT::Clone( *tmp, 1 );
403 stateStorageInitialized_ =
true;
410 template <
class ScalarType,
class MV,
class OP>
414 if (!stateStorageInitialized_)
417 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
418 "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
422 std::string errstr(
"Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
425 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
427 if (newstate.
R != Teuchos::null) {
429 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
430 std::invalid_argument, errstr );
431 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
432 std::invalid_argument, errstr );
435 if (newstate.
R != R_) {
437 MVT::Assign( *newstate.
R, *R_ );
443 W_ = MVT::CloneCopy( *R_ );
444 U_ = MVT::CloneCopy( *R_ );
445 Rtilde_ = MVT::CloneCopy( *R_ );
447 MVT::MvInit( *solnUpdate_ );
451 lp_->apply( *U_, *V_ );
452 AU_ = MVT::CloneCopy( *V_ );
457 MVT::MvNorm( *R_, tau_ );
458 MVT::MvDot( *R_, *Rtilde_, rho_old_ );
462 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
463 "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
473 template <
class ScalarType,
class MV,
class OP>
479 if (initialized_ ==
false) {
484 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
485 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
486 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
487 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
488 ScalarType eta = STzero, beta = STzero;
493 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
497 "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
503 while (stest_->checkStatus(
this) !=
Passed) {
505 for (
int iIter=0; iIter<2; iIter++)
513 MVT::MvDot( *V_, *Rtilde_, alpha_ );
514 alpha_[0] = rho_old_[0]/alpha_[0];
522 MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
529 MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
540 MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
543 lp_->apply( *U_, *AU_ );
550 MVT::MvNorm( *W_, theta_ );
551 theta_[0] /= tau_[0];
553 cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
554 tau_[0] *= theta_[0]*cs_[0];
555 eta = cs_[0]*cs_[0]*alpha_[0];
562 MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
567 if ( tau_[0] == MTzero ) {
577 MVT::MvDot( *W_, *Rtilde_, rho_ );
578 beta = rho_[0]/rho_old_[0];
579 rho_old_[0] = rho_[0];
586 MVT::MvAddMv( STone, *W_, beta, *U_, *U_ );
589 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
592 lp_->apply( *U_, *AU_ );
595 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
608 #endif // BELOS_TFQMR_ITER_HPP Teuchos::RCP< const MV > V
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Collection of types and exceptions used within the Belos solvers.
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
TFQMRIterateFailure(const std::string &what_arg)
void setBlockSize(int blockSize)
Set the blocksize.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Declaration of basic traits for the multivector type.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
TFQMRIter(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::TFQMRIter constructor.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
int getNumIters() const
Get the current iteration count.
TFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
virtual ~TFQMRIter()
Belos::TFQMRIter destructor.
Teuchos::RCP< const MV > R
The current residual basis.
Teuchos::RCP< const MV > U
void resetNumIters(int iter=0)
Reset the iteration count.
Pure virtual base class which describes the basic interface to the linear solver iteration.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
void initializeTFQMR(const TFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > Rtilde
void iterate()
This method performs TFQMR iterations until the status test indicates the need to stop or an error oc...
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< const MV > W
A linear system to solve, and its associated information.
Structure to contain pointers to TFQMRIter state variables.
Class which describes the linear problem to be solved by the iterative solver.
TFQMRIterInitFailure(const std::string &what_arg)
TFQMRIterInitFailure is thrown when the TFQMRIter object is unable to generate an initial iterate in ...
MultiVecTraits< ScalarType, MV > MVT
TFQMRIterateFailure is thrown when the TFQMRIter object is unable to compute the next iterate in the ...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Teuchos::RCP< const MV > D
Belos header file which uses auto-configuration information to include necessary C++ headers...