42 #ifndef BELOS_GMRESPOLYOP_HPP 43 #define BELOS_GMRESPOLYOP_HPP 53 #include "Teuchos_RCP.hpp" 54 #include "Teuchos_SerialDenseMatrix.hpp" 55 #include "Teuchos_SerialDenseVector.hpp" 70 template <
class ScalarType,
class MV,
class OP>
82 const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> >& hess,
83 const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> >& comb,
84 const Teuchos::RCP<Teuchos::SerialDenseVector<int, ScalarType> >& scal
86 : problem_(problem_in), LP_(problem_in->getLeftPrec()), RP_(problem_in->getRightPrec()), H_(hess), y_(comb), r0_(scal)
108 typedef Teuchos::ScalarTraits<ScalarType> SCT ;
111 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
112 Teuchos::RCP<const OP> LP_, RP_;
113 Teuchos::RCP<MV> V_, wL_, wR_;
114 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_, y_;
115 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > r0_;
118 template <
class ScalarType,
class MV,
class OP>
122 if (V_ == Teuchos::null) {
123 V_ = MVT::Clone( x, dim_ );
124 if (LP_ != Teuchos::null) {
125 wL_ = MVT::Clone( y, 1 );
127 if (RP_ != Teuchos::null) {
128 wR_ = MVT::Clone( y, 1 );
134 int n = MVT::GetNumberVecs( x );
135 std::vector<int> idxi(1), idxi2, idxj(1);
138 for (
int j=0; j<n; ++j) {
142 Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
143 Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
144 if (LP_ != Teuchos::null) {
145 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
146 problem_->applyLeftPrec( *x_view, *v_curr );
148 MVT::SetBlock( *x_view, idxi, *V_ );
151 for (
int i=0; i<dim_-1; ++i) {
155 for (
int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
156 Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
159 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
161 Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
167 if (RP_ != Teuchos::null) {
168 problem_->applyRightPrec( *v_curr, *wR_ );
173 if (LP_ == Teuchos::null) {
177 problem_->applyOp( *wR_, *wL_ );
179 if (LP_ != Teuchos::null) {
180 problem_->applyLeftPrec( *wL_, *v_next );
184 Teuchos::SerialDenseMatrix<int,ScalarType> h(Teuchos::View,*H_,i+1,1,0,i);
185 MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
188 MVT::MvScale( *v_next, SCT::one()/(*H_)(i+1,i) );
192 if (RP_ != Teuchos::null) {
193 MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *wR_ );
194 problem_->applyRightPrec( *wR_, *y_view );
197 MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *y_view );
211 template <
class ScalarType,
class MV,
class OP>
220 { Op.
Apply( x, y, trans ); }
static void Apply(const GmresPolyOp< ScalarType, MV, OP > &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
ETrans
Whether to apply the (conjugate) transpose of an operator.
Traits class which defines basic operations on multivectors.
GmresPolyOp()
Default constructor.
void Apply(const MV &x, MV &y, ETrans trans=NOTRANS)
This routine takes the MV x and applies the polynomial operator phi(OP) to it resulting in the MV y...
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in, const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > &hess, const Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > &comb, const Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > &scal)
Basic contstructor.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Belos's class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
virtual ~GmresPolyOp()
Destructor.
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...