50 #ifndef BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 51 #define BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 71 #include "Teuchos_BLAS.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 );
221 initializeTFQMR(empty);
239 state.solnUpdate = solnUpdate_;
257 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms )
const;
279 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
280 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
294 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
295 const Teuchos::RCP<OutputManager<ScalarType> > om_;
296 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
306 std::vector<ScalarType> alpha_, rho_, rho_old_;
307 std::vector<MagnitudeType> tau_, cs_, theta_;
325 Teuchos::RCP<MV> U_, AU_;
326 Teuchos::RCP<MV> Rtilde_;
329 Teuchos::RCP<MV> solnUpdate_;
340 template <
class ScalarType,
class MV,
class OP>
344 Teuchos::ParameterList ¶ms
357 template <
class ScalarType,
class MV,
class OP>
358 Teuchos::RCP<const MV>
361 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
364 if ((
int) normvec->size() < numRHS_)
365 normvec->resize( numRHS_ );
368 for (
int i=0; i<numRHS_; i++) {
369 (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( iter_ + one )*tau_[i];
373 return Teuchos::null;
378 template <
class ScalarType,
class MV,
class OP>
382 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
383 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
386 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
388 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
389 "Belos::PseudoBlockTFQMRIter::initialize(): linear problem does not specify multivectors to clone from.");
405 alpha_.resize( numRHS_ );
406 rho_.resize( numRHS_ );
407 rho_old_.resize( numRHS_ );
408 tau_.resize( numRHS_ );
409 cs_.resize( numRHS_ );
410 theta_.resize( numRHS_ );
415 std::string errstr(
"Belos::PseudoBlockTFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
418 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
419 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
420 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
422 if (newstate.
R != Teuchos::null) {
425 std::invalid_argument, errstr );
428 if (newstate.
R != R_) {
444 lp_->apply( *U_, *V_ );
449 for (
int i=0; i<numRHS_; i++)
456 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
457 "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
467 template <
class ScalarType,
class MV,
class OP>
473 if (initialized_ ==
false) {
478 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
479 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
480 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
481 std::vector< ScalarType > eta(numRHS_,STzero), beta(numRHS_,STzero);
482 std::vector<int> index(1);
487 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
492 while (stest_->checkStatus(
this) !=
Passed) {
501 for (
int i=0; i<numRHS_; i++)
502 alpha_[i] = rho_old_[i]/alpha_[i];
509 for (
int i=0; i<numRHS_; ++i) {
529 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta[i], *D_i, *D_i );
551 lp_->apply( *U_, *AU_ );
560 for (
int i=0; i<numRHS_; ++i) {
561 theta_[i] /= tau_[i];
563 cs_[i] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
564 tau_[i] *= theta_[i]*cs_[i];
565 eta[i] = cs_[i]*cs_[i]*alpha_[i];
572 for (
int i=0; i<numRHS_; ++i) {
576 MVT::MvAddMv( STone, *update_i, eta[i], *D_i, *update_i );
587 for (
int i=0; i<numRHS_; ++i) {
588 beta[i] = rho_[i]/rho_old_[i];
589 rho_old_[i] = rho_[i];
609 lp_->apply( *U_, *AU_ );
612 for (
int i=0; i<numRHS_; ++i) {
628 #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.
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
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.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
Declaration of basic traits for the multivector type.
Teuchos::RCP< const MV > R
The current residual basis.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
int getNumIters() const
Get the current iteration count.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
Pure virtual base class which describes the basic interface to the linear solver iteration.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
PseudoBlockTFQMRIterInitFailure(const std::string &what_arg)
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
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()
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
PseudoBlockTFQMRIterInitFailure is thrown when the PseudoBlockTFQMRIter object is unable to generate ...
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ScalarTraits< ScalarType > SCT
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setBlockSize(int blockSize)
Set the blocksize.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Teuchos::RCP< const MV > D
Teuchos::RCP< const MV > U
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
PseudoBlockTFQMRIterateFailure(const std::string &what_arg)
virtual ~PseudoBlockTFQMRIter()
Belos::PseudoBlockTFQMRIter destructor.
Teuchos::RCP< const MV > W
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
Class which defines basic traits for the operator type.
void iterate()
This method performs block TFQMR iterations until the status test indicates the need to stop or an er...
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...
Teuchos::RCP< const MV > V
bool isInitialized()
States whether the solver has been initialized or not.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.