BelosGmresPolyOp.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_GMRESPOLYOP_HPP
00030 #define BELOS_GMRESPOLYOP_HPP
00031 
00036 #include "BelosOperatorTraits.hpp"
00037 #include "BelosMultiVecTraits.hpp"
00038 #include "BelosLinearProblem.hpp"
00039 #include "BelosConfigDefs.hpp"
00040 #include "Teuchos_RCP.hpp"
00041 #include "Teuchos_SerialDenseMatrix.hpp"
00042 #include "Teuchos_SerialDenseVector.hpp"
00043 
00055 namespace Belos {
00056   
00057   template <class ScalarType, class MV, class OP>
00058   class GmresPolyOp {
00059   public:
00060     
00062 
00063     
00065     GmresPolyOp() {}
00066     
00068     GmresPolyOp( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem_in, 
00069                  const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> >& hess,
00070                  const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> >& comb,
00071                  const Teuchos::RCP<Teuchos::SerialDenseVector<int, ScalarType> >& scal
00072                ) 
00073       : problem_(problem_in), LP_(problem_in->getLeftPrec()), RP_(problem_in->getRightPrec()), H_(hess), y_(comb), r0_(scal)
00074     {
00075       dim_ = y_->numRows();
00076     }
00077 
00079     virtual ~GmresPolyOp() {};
00081     
00083 
00084     
00090     void Apply ( const MV& x, MV& y, ETrans trans=NOTRANS );
00091 
00092     private:
00093 
00094     typedef MultiVecTraits<ScalarType,MV> MVT;
00095     typedef Teuchos::ScalarTraits<ScalarType> SCT ;
00096 
00097     int dim_;
00098     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00099     Teuchos::RCP<const OP> LP_, RP_;
00100     Teuchos::RCP<MV> V_, wL_, wR_;
00101     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_, y_;
00102     Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > r0_;
00103   };
00104   
00105   template <class ScalarType, class MV, class OP>
00106   void GmresPolyOp<ScalarType, MV, OP>::Apply( const MV& x, MV& y, ETrans trans ) 
00107   {
00108     // Initialize vector storage.
00109     if (V_ == Teuchos::null) {
00110       V_ = MVT::Clone( x, dim_ );
00111       if (LP_ != Teuchos::null) {
00112         wL_ = MVT::Clone( y, 1 );
00113       }
00114       if (RP_ != Teuchos::null) {
00115         wR_ = MVT::Clone( y, 1 );
00116       }
00117     }
00118     //
00119     // Apply polynomial to x.
00120     // 
00121     int n = MVT::GetNumberVecs( x );
00122     std::vector<int> idxi(1), idxi2, idxj(1);
00123 
00124     // Select vector x[j].
00125     for (int j=0; j<n; ++j) {
00126 
00127       idxi[0] = 0;
00128       idxj[0] = j;
00129       Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
00130       Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
00131       if (LP_ != Teuchos::null) {
00132         Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
00133         problem_->applyLeftPrec( *x_view, *v_curr ); // Left precondition x into the first vector of V
00134       } else {
00135         MVT::SetBlock( *x_view, idxi, *V_ );  // Set x as the first vector of V
00136       }
00137 
00138       for (int i=0; i<dim_-1; ++i) {
00139 
00140         // Get views into the current and next vectors
00141         idxi2.resize(i+1);
00142         for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
00143         Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
00144         // the tricks below with wR_ and wL_ (potentially set to v_curr and v_next) unfortunately imply that 
00145         // v_curr and v_next must be non-const views.
00146         Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
00147         idxi[0] = i+1;
00148         Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
00149 
00150         //---------------------------------------------
00151         // Apply operator to next vector
00152         //---------------------------------------------
00153         // 1) Apply right preconditioner, if we have one.
00154         if (RP_ != Teuchos::null) {
00155           problem_->applyRightPrec( *v_curr, *wR_ );
00156         } else {
00157           wR_ = v_curr;
00158         }
00159         // 2) Check for right preconditioner, if none exists, point at the next vector.
00160         if (LP_ == Teuchos::null) {
00161           wL_ = v_next;
00162         }
00163         // 3) Apply operator A.
00164         problem_->applyOp( *wR_, *wL_ );
00165         // 4) Apply left preconditioner, if we have one.
00166         if (LP_ != Teuchos::null) {
00167           problem_->applyLeftPrec( *wL_, *v_next );
00168         }
00169 
00170         // Compute A*v_curr - v_prev*H(1:i,i)
00171         Teuchos::SerialDenseMatrix<int,ScalarType> h(Teuchos::View,*H_,i+1,1,0,i);
00172         MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
00173 
00174         // Scale by H(i+1,i)
00175         MVT::MvScale( *v_next, SCT::one()/(*H_)(i+1,i) );  
00176       }
00177 
00178       // Compute output y = V*y_./r0_
00179       if (RP_ != Teuchos::null) {
00180         MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *wR_ );
00181         problem_->applyRightPrec( *wR_, *y_view );
00182       } 
00183       else {
00184         MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *y_view );
00185       }
00186     } // (int j=0; j<n; ++j)
00187   } // end Apply()
00188 
00190   //
00191   // Implementation of the Belos::OperatorTraits for Belos::GmresPolyOp 
00192   //
00194 
00198   template <class ScalarType, class MV, class OP> 
00199   class OperatorTraits < ScalarType, MV, GmresPolyOp<ScalarType,MV,OP> > 
00200   {
00201   public:
00202     
00204     static void Apply ( const GmresPolyOp<ScalarType,MV,OP>& Op, 
00205       const MV& x, MV& y,
00206       ETrans trans=NOTRANS )
00207     { Op.Apply( x, y, trans ); }
00208     
00209   };
00210   
00211 } // end Belos namespace
00212 
00213 #endif
00214 
00215 // end of file BelosGmresPolyOp.hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:06 2011 for Belos by  doxygen 1.6.3