Belos  Version of the Day
BelosGmresPolyOp.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_GMRESPOLYOP_HPP
43 #define BELOS_GMRESPOLYOP_HPP
44 
49 #include "BelosOperatorTraits.hpp"
50 #include "BelosMultiVecTraits.hpp"
51 #include "BelosLinearProblem.hpp"
52 #include "BelosConfigDefs.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_SerialDenseMatrix.hpp"
55 #include "Teuchos_SerialDenseVector.hpp"
56 
68 namespace Belos {
69 
70  template <class ScalarType, class MV, class OP>
71  class GmresPolyOp {
72  public:
73 
75 
76 
79 
81  GmresPolyOp( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem_in,
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
85  )
86  : problem_(problem_in), LP_(problem_in->getLeftPrec()), RP_(problem_in->getRightPrec()), H_(hess), y_(comb), r0_(scal)
87  {
88  dim_ = y_->numRows();
89  }
90 
92  virtual ~GmresPolyOp() {};
94 
96 
97 
103  void Apply ( const MV& x, MV& y, ETrans trans=NOTRANS );
104 
105  private:
106 
107  typedef MultiVecTraits<ScalarType,MV> MVT;
108  typedef Teuchos::ScalarTraits<ScalarType> SCT ;
109 
110  int dim_;
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_;
116  };
117 
118  template <class ScalarType, class MV, class OP>
119  void GmresPolyOp<ScalarType, MV, OP>::Apply( const MV& x, MV& y, ETrans trans )
120  {
121  // Initialize vector storage.
122  if (V_ == Teuchos::null) {
123  V_ = MVT::Clone( x, dim_ );
124  if (LP_ != Teuchos::null) {
125  wL_ = MVT::Clone( y, 1 );
126  }
127  if (RP_ != Teuchos::null) {
128  wR_ = MVT::Clone( y, 1 );
129  }
130  }
131  //
132  // Apply polynomial to x.
133  //
134  int n = MVT::GetNumberVecs( x );
135  std::vector<int> idxi(1), idxi2, idxj(1);
136 
137  // Select vector x[j].
138  for (int j=0; j<n; ++j) {
139 
140  idxi[0] = 0;
141  idxj[0] = 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 ); // Left precondition x into the first vector of V
147  } else {
148  MVT::SetBlock( *x_view, idxi, *V_ ); // Set x as the first vector of V
149  }
150 
151  for (int i=0; i<dim_-1; ++i) {
152 
153  // Get views into the current and next vectors
154  idxi2.resize(i+1);
155  for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
156  Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
157  // the tricks below with wR_ and wL_ (potentially set to v_curr and v_next) unfortunately imply that
158  // v_curr and v_next must be non-const views.
159  Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
160  idxi[0] = i+1;
161  Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
162 
163  //---------------------------------------------
164  // Apply operator to next vector
165  //---------------------------------------------
166  // 1) Apply right preconditioner, if we have one.
167  if (RP_ != Teuchos::null) {
168  problem_->applyRightPrec( *v_curr, *wR_ );
169  } else {
170  wR_ = v_curr;
171  }
172  // 2) Check for right preconditioner, if none exists, point at the next vector.
173  if (LP_ == Teuchos::null) {
174  wL_ = v_next;
175  }
176  // 3) Apply operator A.
177  problem_->applyOp( *wR_, *wL_ );
178  // 4) Apply left preconditioner, if we have one.
179  if (LP_ != Teuchos::null) {
180  problem_->applyLeftPrec( *wL_, *v_next );
181  }
182 
183  // Compute A*v_curr - v_prev*H(1:i,i)
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 );
186 
187  // Scale by H(i+1,i)
188  MVT::MvScale( *v_next, SCT::one()/(*H_)(i+1,i) );
189  }
190 
191  // Compute output y = V*y_./r0_
192  if (RP_ != Teuchos::null) {
193  MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *wR_ );
194  problem_->applyRightPrec( *wR_, *y_view );
195  }
196  else {
197  MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *y_view );
198  }
199  } // (int j=0; j<n; ++j)
200  } // end Apply()
201 
203  //
204  // Implementation of the Belos::OperatorTraits for Belos::GmresPolyOp
205  //
207 
211  template <class ScalarType, class MV, class OP>
212  class OperatorTraits < ScalarType, MV, GmresPolyOp<ScalarType,MV,OP> >
213  {
214  public:
215 
217  static void Apply ( const GmresPolyOp<ScalarType,MV,OP>& Op,
218  const MV& x, MV& y,
219  ETrans trans=NOTRANS )
220  { Op.Apply( x, y, trans ); }
221 
222  };
223 
224 } // end Belos namespace
225 
226 #endif
227 
228 // end of file BelosGmresPolyOp.hpp
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.
Definition: BelosTypes.hpp:81
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&#39;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...

Generated for Belos by doxygen 1.8.14