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>
96 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
100 Teuchos::RCP<const MV>
W;
101 Teuchos::RCP<const MV>
U;
102 Teuchos::RCP<const MV>
AU;
104 Teuchos::RCP<const MV>
D;
105 Teuchos::RCP<const MV>
V;
111 Rtilde(Teuchos::null),
D(Teuchos::null),
V(Teuchos::null)
147 template <
class ScalarType,
class MV,
class OP>
155 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
165 Teuchos::ParameterList ¶ms );
249 state.
alpha = alpha_;
253 state.
theta = theta_;
294 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
295 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
309 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
310 const Teuchos::RCP<OutputManager<ScalarType> > om_;
311 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
321 std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
322 std::vector<MagnitudeType> tau_, theta_;
339 Teuchos::RCP<MV> U_, AU_;
340 Teuchos::RCP<MV> Rtilde_;
343 Teuchos::RCP<MV> solnUpdate_;
354 template <
class ScalarType,
class MV,
class OP>
358 Teuchos::ParameterList ¶ms
371 template <
class ScalarType,
class MV,
class OP>
372 Teuchos::RCP<const MV>
375 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
378 if ((
int) normvec->size() < numRHS_)
379 normvec->resize( numRHS_ );
382 for (
int i=0; i<numRHS_; i++) {
383 (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
387 return Teuchos::null;
392 template <
class ScalarType,
class MV,
class OP>
396 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
397 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
398 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
401 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
Rtilde == Teuchos::null,std::invalid_argument,
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_) )
413 if ( Teuchos::is_null(D_) )
414 D_ = MVT::Clone( *newstate.
Rtilde, numRHS_ );
415 MVT::MvInit( *D_, STzero );
418 if ( Teuchos::is_null(solnUpdate_) )
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) {
486 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
487 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
488 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
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) {
528 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
529 Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
530 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
537 Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
538 Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
539 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
550 Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
551 Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
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];
573 MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
574 tau_[i] *= theta_[i]*cs;
575 eta_[i] = cs*cs*alpha_[i];
582 for (
int i=0; i<numRHS_; ++i) {
584 Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
585 Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
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];
607 Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
608 Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
609 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i );
612 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
613 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
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) {
623 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
624 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
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.
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
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.