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)
126 template <
class ScalarType,
class MV,
class OP>
134 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
144 Teuchos::ParameterList ¶ms );
213 state.solnUpdate = solnUpdate_;
253 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
254 "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
274 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
275 const Teuchos::RCP<OutputManager<ScalarType> > om_;
276 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
284 std::vector<ScalarType> alpha_, rho_, rho_old_;
285 std::vector<MagnitudeType> tau_, cs_, theta_;
298 bool stateStorageInitialized_;
308 Teuchos::RCP<MV> U_, AU_;
309 Teuchos::RCP<MV> Rtilde_;
312 Teuchos::RCP<MV> solnUpdate_;
322 template <
class ScalarType,
class MV,
class OP>
326 Teuchos::ParameterList &
338 stateStorageInitialized_(false),
345 template <
class ScalarType,
class MV,
class OP>
346 Teuchos::RCP<const MV>
349 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
351 (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
353 return Teuchos::null;
359 template <
class ScalarType,
class MV,
class OP>
362 if (!stateStorageInitialized_) {
365 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
366 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
367 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
368 stateStorageInitialized_ =
false;
375 if (R_ == Teuchos::null) {
377 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
378 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
379 "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
380 R_ = MVT::Clone( *tmp, 1 );
381 D_ = MVT::Clone( *tmp, 1 );
382 V_ = MVT::Clone( *tmp, 1 );
383 solnUpdate_ = MVT::Clone( *tmp, 1 );
387 stateStorageInitialized_ =
true;
394 template <
class ScalarType,
class MV,
class OP>
398 if (!stateStorageInitialized_)
401 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
402 "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
406 std::string errstr(
"Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
409 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
411 if (newstate.
R != Teuchos::null) {
413 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
414 std::invalid_argument, errstr );
415 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
416 std::invalid_argument, errstr );
419 if (newstate.
R != R_) {
421 MVT::Assign( *newstate.
R, *R_ );
427 W_ = MVT::CloneCopy( *R_ );
428 U_ = MVT::CloneCopy( *R_ );
429 Rtilde_ = MVT::CloneCopy( *R_ );
431 MVT::MvInit( *solnUpdate_ );
435 lp_->apply( *U_, *V_ );
436 AU_ = MVT::CloneCopy( *V_ );
441 MVT::MvNorm( *R_, tau_ );
442 MVT::MvDot( *R_, *Rtilde_, rho_old_ );
446 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
447 "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
457 template <
class ScalarType,
class MV,
class OP>
463 if (initialized_ ==
false) {
468 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
469 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
470 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
471 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
472 ScalarType eta = STzero, beta = STzero;
477 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
481 "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
487 while (stest_->checkStatus(
this) !=
Passed) {
489 for (
int iIter=0; iIter<2; iIter++)
497 MVT::MvDot( *V_, *Rtilde_, alpha_ );
498 alpha_[0] = rho_old_[0]/alpha_[0];
506 MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
513 MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
524 MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
527 lp_->apply( *U_, *AU_ );
534 MVT::MvNorm( *W_, theta_ );
535 theta_[0] /= tau_[0];
537 cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
538 tau_[0] *= theta_[0]*cs_[0];
539 eta = cs_[0]*cs_[0]*alpha_[0];
546 MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
551 if ( tau_[0] == MTzero ) {
561 MVT::MvDot( *W_, *Rtilde_, rho_ );
562 beta = rho_[0]/rho_old_[0];
563 rho_old_[0] = rho_[0];
570 MVT::MvAddMv( STone, *W_, beta, *U_, *U_ );
573 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
576 lp_->apply( *U_, *AU_ );
579 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which describes the basic interface to the linear solver iteration.
Class which describes the linear problem to be solved by the iterative solver.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A pure virtual class for defining the status tests for the Belos iterative solvers.
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
Teuchos::ScalarTraits< ScalarType > SCT
void resetNumIters(int iter=0)
Reset the iteration count.
virtual ~TFQMRIter()
Belos::TFQMRIter destructor.
void setBlockSize(int blockSize)
Set the blocksize.
bool isInitialized()
States whether the solver has been initialized or not.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
void iterate()
This method performs TFQMR iterations until the status test indicates the need to stop or an error oc...
SCT::magnitudeType MagnitudeType
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
TFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
MultiVecTraits< ScalarType, MV > MVT
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.
int getNumIters() const
Get the current iteration count.
void initializeTFQMR(const TFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
TFQMRIterateFailure is thrown when the TFQMRIter object is unable to compute the next iterate in the ...
TFQMRIterateFailure(const std::string &what_arg)
Structure to contain pointers to TFQMRIter state variables.
Teuchos::RCP< const MV > W
Teuchos::RCP< const MV > V
Teuchos::RCP< const MV > Rtilde
Teuchos::RCP< const MV > R
The current residual basis.
Teuchos::RCP< const MV > D
Teuchos::RCP< const MV > U