50 #ifndef BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 51 #define BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP 93 template <
class ScalarType,
class MV>
147 template <
class ScalarType,
class MV,
class OP>
249 state.
alpha = alpha_;
253 state.
theta = theta_;
295 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
321 std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
322 std::vector<MagnitudeType> tau_, theta_;
354 template <
class ScalarType,
class MV,
class OP>
371 template <
class ScalarType,
class MV,
class OP>
378 if ((
int) normvec->size() < numRHS_)
379 normvec->resize( numRHS_ );
382 for (
int i=0; i<numRHS_; i++) {
387 return Teuchos::null;
392 template <
class ScalarType,
class MV,
class OP>
402 "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
405 int numRHS = MVT::GetNumberVecs(*newstate.
Rtilde);
410 if (
Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
414 D_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
415 MVT::MvInit( *D_, STzero );
419 solnUpdate_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
420 MVT::MvInit( *solnUpdate_, STzero );
423 if (newstate.
Rtilde != Rtilde_)
424 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
425 W_ = MVT::CloneCopy( *Rtilde_ );
426 U_ = MVT::CloneCopy( *Rtilde_ );
427 V_ = MVT::Clone( *Rtilde_, numRHS_ );
431 lp_->apply( *U_, *V_ );
432 AU_ = MVT::CloneCopy( *V_ );
435 alpha_.resize( numRHS_, STone );
436 eta_.resize( numRHS_, STzero );
437 rho_.resize( numRHS_ );
438 rho_old_.resize( numRHS_ );
439 tau_.resize( numRHS_ );
440 theta_.resize( numRHS_, MTzero );
442 MVT::MvNorm( *Rtilde_, tau_ );
443 MVT::MvDot( *Rtilde_, *Rtilde_, rho_ );
448 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
449 W_ = MVT::CloneCopy( *newstate.
W );
450 U_ = MVT::CloneCopy( *newstate.
U );
451 AU_ = MVT::CloneCopy( *newstate.
AU );
452 D_ = MVT::CloneCopy( *newstate.
D );
453 V_ = MVT::CloneCopy( *newstate.
V );
457 solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
458 MVT::MvInit( *solnUpdate_, STzero );
461 alpha_ = newstate.
alpha;
465 theta_ = newstate.
theta;
475 template <
class ScalarType,
class MV,
class OP>
481 if (initialized_ ==
false) {
489 std::vector< ScalarType > beta(numRHS_,STzero);
490 std::vector<int> index(1);
498 while (stest_->checkStatus(
this) !=
Passed) {
500 for (
int iIter=0; iIter<2; iIter++)
508 MVT::MvDot( *V_, *Rtilde_, alpha_ );
509 for (
int i=0; i<numRHS_; i++) {
510 rho_old_[i] = rho_[i];
511 alpha_[i] = rho_old_[i]/alpha_[i];
519 for (
int i=0; i<numRHS_; ++i) {
530 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
539 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
552 MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
561 lp_->apply( *U_, *AU_ );
568 MVT::MvNorm( *W_, theta_ );
570 for (
int i=0; i<numRHS_; ++i) {
571 theta_[i] /= tau_[i];
574 tau_[i] *= theta_[i]*cs;
575 eta_[i] = cs*cs*alpha_[i];
582 for (
int i=0; i<numRHS_; ++i) {
586 MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
595 MVT::MvDot( *W_, *Rtilde_, rho_ );
597 for (
int i=0; i<numRHS_; ++i) {
598 beta[i] = rho_[i]/rho_old_[i];
609 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i );
614 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
618 lp_->apply( *U_, *AU_ );
621 for (
int i=0; i<numRHS_; ++i) {
625 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
638 #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.
int getNumIters() const
Get the current iteration count.
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.
bool is_null(const std::shared_ptr< T > &p)
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< const MV > AU
Pure virtual base class which describes the basic interface to the linear solver iteration.
std::vector< ScalarType > alpha
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
PseudoBlockTFQMRIterInitFailure(const std::string &what_arg)
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()
PseudoBlockTFQMRIterInitFailure is thrown when the PseudoBlockTFQMRIter object is unable to generate ...
Teuchos::ScalarTraits< ScalarType > SCT
A linear system to solve, and its associated information.
std::vector< MagnitudeType > theta
Class which describes the linear problem to be solved by the iterative solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::RCP< const MV > D
std::vector< ScalarType > eta
Teuchos::RCP< const MV > U
std::vector< MagnitudeType > tau
PseudoBlockTFQMRIterateFailure(const std::string &what_arg)
virtual ~PseudoBlockTFQMRIter()
Belos::PseudoBlockTFQMRIter destructor.
Teuchos::RCP< const MV > W
The current residual basis.
std::vector< ScalarType > rho
Class which defines basic traits for the operator type.
Teuchos::ScalarTraits< ScalarType > SCT
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...
SCT::magnitudeType MagnitudeType
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > V
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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.