50 #ifndef BELOS_TFQMR_ITER_HPP 51 #define BELOS_TFQMR_ITER_HPP 71 #include "Teuchos_BLAS.hpp" 72 #include "Teuchos_SerialDenseMatrix.hpp" 73 #include "Teuchos_SerialDenseVector.hpp" 74 #include "Teuchos_ScalarTraits.hpp" 75 #include "Teuchos_ParameterList.hpp" 76 #include "Teuchos_TimeMonitor.hpp" 95 template <
class ScalarType,
class MV>
99 Teuchos::RCP<const MV>
R;
100 Teuchos::RCP<const MV>
W;
101 Teuchos::RCP<const MV>
U;
103 Teuchos::RCP<const MV>
D;
104 Teuchos::RCP<const MV>
V;
107 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
143 template <
class ScalarType,
class MV,
class OP>
151 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
161 Teuchos::ParameterList ¶ms );
241 state.solnUpdate = solnUpdate_;
281 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
282 "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
302 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
303 const Teuchos::RCP<OutputManager<ScalarType> > om_;
304 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
312 std::vector<ScalarType> alpha_, rho_, rho_old_;
313 std::vector<MagnitudeType> tau_, cs_, theta_;
326 bool stateStorageInitialized_;
336 Teuchos::RCP<MV> U_, AU_;
337 Teuchos::RCP<MV> Rtilde_;
340 Teuchos::RCP<MV> solnUpdate_;
350 template <
class ScalarType,
class MV,
class OP>
354 Teuchos::ParameterList ¶ms
366 stateStorageInitialized_(false),
373 template <
class ScalarType,
class MV,
class OP>
374 Teuchos::RCP<const MV>
377 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
379 (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
381 return Teuchos::null;
387 template <
class ScalarType,
class MV,
class OP>
390 if (!stateStorageInitialized_) {
393 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
394 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
395 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
396 stateStorageInitialized_ =
false;
403 if (R_ == Teuchos::null) {
405 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
406 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
407 "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
408 R_ = MVT::Clone( *tmp, 1 );
409 D_ = MVT::Clone( *tmp, 1 );
410 V_ = MVT::Clone( *tmp, 1 );
411 solnUpdate_ = MVT::Clone( *tmp, 1 );
415 stateStorageInitialized_ =
true;
422 template <
class ScalarType,
class MV,
class OP>
426 if (!stateStorageInitialized_)
429 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
430 "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
434 std::string errstr(
"Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
437 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
438 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
439 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
441 if (newstate.
R != Teuchos::null) {
443 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
R) != MVT::GetGlobalLength(*R_),
444 std::invalid_argument, errstr );
445 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
R) != 1,
446 std::invalid_argument, errstr );
449 if (newstate.
R != R_) {
451 MVT::MvAddMv( one, *newstate.
R, STzero, *newstate.
R, *R_ );
457 W_ = MVT::CloneCopy( *R_ );
458 U_ = MVT::CloneCopy( *R_ );
459 Rtilde_ = MVT::CloneCopy( *R_ );
461 MVT::MvInit( *solnUpdate_ );
465 lp_->apply( *U_, *V_ );
466 AU_ = MVT::CloneCopy( *V_ );
471 MVT::MvNorm( *R_, tau_ );
472 MVT::MvDot( *R_, *Rtilde_, rho_old_ );
476 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
477 "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
487 template <
class ScalarType,
class MV,
class OP>
493 if (initialized_ ==
false) {
498 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
499 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
500 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
501 ScalarType eta = STzero, beta = STzero;
506 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
510 "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
516 while (stest_->checkStatus(
this) !=
Passed) {
518 for (
int iIter=0; iIter<2; iIter++)
526 MVT::MvDot( *V_, *Rtilde_, alpha_ );
527 alpha_[0] = rho_old_[0]/alpha_[0];
535 MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
542 MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
553 MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
556 lp_->apply( *U_, *AU_ );
563 MVT::MvNorm( *W_, theta_ );
564 theta_[0] /= tau_[0];
566 cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
567 tau_[0] *= theta_[0]*cs_[0];
568 eta = cs_[0]*cs_[0]*alpha_[0];
576 MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
584 MVT::MvDot( *W_, *Rtilde_, rho_ );
585 beta = rho_[0]/rho_old_[0];
586 rho_old_[0] = rho_[0];
593 MVT::MvAddMv( STone, *W_, beta, *U_, *U_ );
596 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
599 lp_->apply( *U_, *AU_ );
602 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );
615 #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 block TFQMR iterations until the status test indicates the need to stop or an er...
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...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...