50 #ifndef BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
51 #define BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
91 template <
class ScalarType,
class MV>
109 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
113 template <
class ScalarType,
class MV,
class OP>
204 state.
alpha = alpha_;
208 state.
theta = theta_;
250 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
276 std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
277 std::vector<MagnitudeType> tau_, theta_;
309 template <
class ScalarType,
class MV,
class OP>
326 template <
class ScalarType,
class MV,
class OP>
333 if ((
int) normvec->size() < numRHS_)
334 normvec->resize( numRHS_ );
337 for (
int i=0; i<numRHS_; i++) {
342 return Teuchos::null;
347 template <
class ScalarType,
class MV,
class OP>
357 "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
360 int numRHS = MVT::GetNumberVecs(*newstate.
Rtilde);
365 if (
Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
369 D_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
370 MVT::MvInit( *D_, STzero );
374 solnUpdate_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
375 MVT::MvInit( *solnUpdate_, STzero );
378 if (newstate.
Rtilde != Rtilde_)
379 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
380 W_ = MVT::CloneCopy( *Rtilde_ );
381 U_ = MVT::CloneCopy( *Rtilde_ );
382 V_ = MVT::Clone( *Rtilde_, numRHS_ );
386 lp_->apply( *U_, *V_ );
387 AU_ = MVT::CloneCopy( *V_ );
390 alpha_.resize( numRHS_, STone );
391 eta_.resize( numRHS_, STzero );
392 rho_.resize( numRHS_ );
393 rho_old_.resize( numRHS_ );
394 tau_.resize( numRHS_ );
395 theta_.resize( numRHS_, MTzero );
397 MVT::MvNorm( *Rtilde_, tau_ );
398 MVT::MvDot( *Rtilde_, *Rtilde_, rho_ );
403 Rtilde_ = MVT::CloneCopy( *newstate.
Rtilde );
404 W_ = MVT::CloneCopy( *newstate.
W );
405 U_ = MVT::CloneCopy( *newstate.
U );
406 AU_ = MVT::CloneCopy( *newstate.
AU );
407 D_ = MVT::CloneCopy( *newstate.
D );
408 V_ = MVT::CloneCopy( *newstate.
V );
412 solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
413 MVT::MvInit( *solnUpdate_, STzero );
416 alpha_ = newstate.
alpha;
420 theta_ = newstate.
theta;
430 template <
class ScalarType,
class MV,
class OP>
436 if (initialized_ ==
false) {
445 std::vector< ScalarType > beta(numRHS_,STzero);
446 std::vector<int> index(1);
454 while (stest_->checkStatus(
this) !=
Passed) {
456 for (
int iIter=0; iIter<2; iIter++)
464 MVT::MvDot( *V_, *Rtilde_, alpha_ );
465 for (
int i=0; i<numRHS_; i++) {
466 rho_old_[i] = rho_[i];
467 alpha_[i] = rho_old_[i]/alpha_[i];
475 for (
int i=0; i<numRHS_; ++i) {
486 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
495 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
508 MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
517 lp_->apply( *U_, *AU_ );
524 MVT::MvNorm( *W_, theta_ );
526 bool breakdown=
false;
527 for (
int i=0; i<numRHS_; ++i) {
528 theta_[i] /= tau_[i];
531 tau_[i] *= theta_[i]*cs;
532 if ( tau_[i] == MTzero ) {
535 eta_[i] = cs*cs*alpha_[i];
542 for (
int i=0; i<numRHS_; ++i) {
546 MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
563 MVT::MvDot( *W_, *Rtilde_, rho_ );
565 for (
int i=0; i<numRHS_; ++i) {
566 beta[i] = rho_[i]/rho_old_[i];
577 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i );
582 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
586 lp_->apply( *U_, *AU_ );
589 for (
int i=0; i<numRHS_; ++i) {
593 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
606 #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.
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.
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.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > AU
int getNumIters() const
Get the current iteration count.
Pure virtual base class which describes the basic interface to the linear solver iteration.
std::vector< ScalarType > alpha
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.
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()
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.
std::vector< MagnitudeType > theta
Class which describes the linear problem to be solved by the iterative solver.
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
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 pseudo-block TFQMR iterations until the status test indicates the need to stop o...
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
Teuchos::RCP< const MV > V
bool isInitialized()
States whether the solver has been initialized or not.
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.