50#ifndef BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
51#define BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
70#include "Teuchos_ScalarTraits.hpp"
71#include "Teuchos_ParameterList.hpp"
72#include "Teuchos_TimeMonitor.hpp"
91 template <
class ScalarType,
class MV>
94 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
98 Teuchos::RCP<const MV>
W;
99 Teuchos::RCP<const MV>
U;
100 Teuchos::RCP<const MV>
AU;
102 Teuchos::RCP<const MV>
D;
103 Teuchos::RCP<const MV>
V;
109 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
113 template <
class ScalarType,
class MV,
class OP>
121 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
131 Teuchos::ParameterList ¶ms );
249 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
250 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
264 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
265 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
266 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
309 template <
class ScalarType,
class MV,
class OP>
313 Teuchos::ParameterList &
326 template <
class ScalarType,
class MV,
class OP>
327 Teuchos::RCP<const MV>
330 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
333 if ((
int) normvec->size() < numRHS_)
334 normvec->resize( numRHS_ );
337 for (
int i=0; i<numRHS_; i++) {
338 (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
342 return Teuchos::null;
347 template <
class ScalarType,
class MV,
class OP>
351 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
352 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
353 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
356 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Rtilde == Teuchos::null,std::invalid_argument,
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_) )
368 if ( Teuchos::is_null(D_) )
369 D_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
370 MVT::MvInit( *D_, STzero );
373 if ( Teuchos::is_null(solnUpdate_) )
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) {
441 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
442 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
443 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
444 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
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) {
484 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
485 Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
486 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
493 Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
494 Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
495 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
506 Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
507 Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
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];
530 MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[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) {
544 Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
545 Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
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];
575 Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
576 Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
577 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i );
580 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
581 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
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) {
591 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
592 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
593 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i );
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.
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...
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
std::vector< ScalarType > rho_
int getNumIters() const
Get the current iteration count.
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.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
std::vector< ScalarType > rho_old_
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
bool isInitialized()
States whether the solver has been initialized or not.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
virtual ~PseudoBlockTFQMRIter()
Belos::PseudoBlockTFQMRIter destructor.
std::vector< MagnitudeType > theta_
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
std::vector< MagnitudeType > tau_
const Teuchos::RCP< OutputManager< ScalarType > > om_
MultiVecTraits< ScalarType, MV > MVT
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
std::vector< ScalarType > eta_
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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.
Teuchos::RCP< MV > solnUpdate_
Teuchos::ScalarTraits< ScalarType > SCT
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< MV > Rtilde_
SCT::magnitudeType MagnitudeType
void setBlockSize(int blockSize)
Set the blocksize.
std::vector< ScalarType > alpha_
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
OperatorTraits< ScalarType, MV, OP > OPT
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockTFQMRIter state variables.
std::vector< MagnitudeType > theta
std::vector< MagnitudeType > tau
SCT::magnitudeType MagnitudeType
Teuchos::RCP< const MV > AU
std::vector< ScalarType > rho
Teuchos::RCP< const MV > D
std::vector< ScalarType > alpha
Teuchos::RCP< const MV > Rtilde
Teuchos::RCP< const MV > W
The current residual basis.
Teuchos::RCP< const MV > U
Teuchos::RCP< const MV > V
PseudoBlockTFQMRIterState()
Teuchos::ScalarTraits< ScalarType > SCT
std::vector< ScalarType > eta