50 #ifndef BELOS_TFQMR_ITER_HPP 51 #define BELOS_TFQMR_ITER_HPP 70 #include "Teuchos_SerialDenseMatrix.hpp" 71 #include "Teuchos_SerialDenseVector.hpp" 72 #include "Teuchos_ScalarTraits.hpp" 73 #include "Teuchos_ParameterList.hpp" 74 #include "Teuchos_TimeMonitor.hpp" 93 template <
class ScalarType,
class MV>
97 Teuchos::RCP<const MV>
R;
98 Teuchos::RCP<const MV>
W;
99 Teuchos::RCP<const MV>
U;
101 Teuchos::RCP<const MV>
D;
102 Teuchos::RCP<const MV>
V;
105 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
141 template <
class ScalarType,
class MV,
class OP>
149 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
159 Teuchos::ParameterList ¶ms );
228 state.solnUpdate = solnUpdate_;
268 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
269 "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
289 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
290 const Teuchos::RCP<OutputManager<ScalarType> > om_;
291 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
299 std::vector<ScalarType> alpha_, rho_, rho_old_;
300 std::vector<MagnitudeType> tau_, cs_, theta_;
313 bool stateStorageInitialized_;
323 Teuchos::RCP<MV> U_, AU_;
324 Teuchos::RCP<MV> Rtilde_;
327 Teuchos::RCP<MV> solnUpdate_;
337 template <
class ScalarType,
class MV,
class OP>
341 Teuchos::ParameterList &
353 stateStorageInitialized_(false),
360 template <
class ScalarType,
class MV,
class OP>
361 Teuchos::RCP<const MV>
364 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
366 (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
368 return Teuchos::null;
374 template <
class ScalarType,
class MV,
class OP>
377 if (!stateStorageInitialized_) {
380 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
381 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
382 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
383 stateStorageInitialized_ =
false;
390 if (R_ == Teuchos::null) {
392 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
393 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
394 "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
395 R_ = MVT::Clone( *tmp, 1 );
396 D_ = MVT::Clone( *tmp, 1 );
397 V_ = MVT::Clone( *tmp, 1 );
398 solnUpdate_ = MVT::Clone( *tmp, 1 );
402 stateStorageInitialized_ =
true;
409 template <
class ScalarType,
class MV,
class OP>
413 if (!stateStorageInitialized_)
416 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
417 "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
421 std::string errstr(
"Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
424 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
426 if (newstate.
R != Teuchos::null) {
428 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
429 std::invalid_argument, errstr );
430 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
431 std::invalid_argument, errstr );
434 if (newstate.
R != R_) {
436 MVT::Assign( *newstate.
R, *R_ );
442 W_ = MVT::CloneCopy( *R_ );
443 U_ = MVT::CloneCopy( *R_ );
444 Rtilde_ = MVT::CloneCopy( *R_ );
446 MVT::MvInit( *solnUpdate_ );
450 lp_->apply( *U_, *V_ );
451 AU_ = MVT::CloneCopy( *V_ );
456 MVT::MvNorm( *R_, tau_ );
457 MVT::MvDot( *R_, *Rtilde_, rho_old_ );
461 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
462 "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
472 template <
class ScalarType,
class MV,
class OP>
478 if (initialized_ ==
false) {
483 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
484 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
485 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
486 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
487 ScalarType eta = STzero, beta = STzero;
492 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
496 "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
502 while (stest_->checkStatus(
this) !=
Passed) {
504 for (
int iIter=0; iIter<2; iIter++)
512 MVT::MvDot( *V_, *Rtilde_, alpha_ );
513 alpha_[0] = rho_old_[0]/alpha_[0];
521 MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
528 MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
539 MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
542 lp_->apply( *U_, *AU_ );
549 MVT::MvNorm( *W_, theta_ );
550 theta_[0] /= tau_[0];
552 cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
553 tau_[0] *= theta_[0]*cs_[0];
554 eta = cs_[0]*cs_[0]*alpha_[0];
561 MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
566 if ( tau_[0] == MTzero ) {
576 MVT::MvDot( *W_, *Rtilde_, rho_ );
577 beta = rho_[0]/rho_old_[0];
578 rho_old_[0] = rho_[0];
585 MVT::MvAddMv( STone, *W_, beta, *U_, *U_ );
588 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
591 lp_->apply( *U_, *AU_ );
594 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
607 #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...