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 
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)
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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.
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Belos header file which uses auto-configuration information to include necessary C++ headers...

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11