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       if (RP_ != Teuchos::null)
00114         wR_ = MVT::Clone( y, 1 );
00115     }
00116     //
00117     // Apply polynomial to x.
00118     // 
00119     int n = MVT::GetNumberVecs( x );
00120     std::vector<int> idxi(1), idxi2, idxj(1);
00121     Teuchos::RCP<MV> v_curr, v_next, v_prev;
00122 
00123     // Select vector x[j].
00124     for (int j=0; j<n; ++j) {
00125 
00126       idxi[0] = 0;
00127       idxj[0] = j;
00128       Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
00129       Teuchos::RCP<MV> y_view = MVT::CloneView( y, idxj );
00130       if (LP_ != Teuchos::null) {
00131   v_curr = MVT::CloneView( *V_, idxi );
00132   problem_->applyLeftPrec( *x_view, *v_curr ); // Left precondition x into the first vector of V
00133       } else {
00134   MVT::SetBlock( *x_view, idxi, *V_ );  // Set x as the first vector of V
00135       }
00136       
00137       for (int i=0; i<dim_-1; ++i) {
00138   
00139   // Get views into the current and next vectors
00140   idxi2.resize(i+1);
00141   for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
00142   v_prev = MVT::CloneView( *V_, idxi2 );
00143   v_curr = MVT::CloneView( *V_, idxi );
00144   idxi[0] = i+1;
00145   v_next = MVT::CloneView( *V_, idxi );
00146   
00147   //---------------------------------------------
00148   // Apply operator to next vector
00149   //---------------------------------------------
00150   // 1) Apply right preconditioner, if we have one.
00151   if (RP_ != Teuchos::null) {
00152     problem_->applyRightPrec( *v_curr, *wR_ );
00153   } else {
00154     wR_ = v_curr;
00155   }
00156   // 2) Check for right preconditioner, if none exists, point at the next vector.
00157   if (LP_ == Teuchos::null) {
00158     wL_ = v_next;
00159   }
00160   // 3) Apply operator A.
00161   problem_->applyOp( *wR_, *wL_ );
00162   // 4) Apply left preconditioner, if we have one.
00163   if (LP_ != Teuchos::null) {
00164     problem_->applyLeftPrec( *wL_, *v_next );
00165   }
00166   
00167   // Compute A*v_curr - v_prev*H(1:i,i)
00168   Teuchos::SerialDenseMatrix<int,ScalarType> h(Teuchos::View,*H_,i+1,1,0,i);
00169   MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
00170   
00171   // Scale by H(i+1,i)
00172   MVT::MvScale( *v_next, SCT::one()/(*H_)(i+1,i) );  
00173       }
00174       
00175       // Compute output y = V*y_./r0_
00176       if (RP_ != Teuchos::null) {
00177   MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *wR_ );
00178   problem_->applyRightPrec( *wR_, *y_view );
00179       } else {
00180   MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *y_view );
00181       }
00182     } // (int j=0; j<n; ++j)
00183   } // end Apply()
00184     
00186   //
00187   // Implementation of the Belos::OperatorTraits for Belos::GmresPolyOp 
00188   //
00190   
00194   template <class ScalarType, class MV, class OP> 
00195   class OperatorTraits < ScalarType, MV, GmresPolyOp<ScalarType,MV,OP> > 
00196   {
00197   public:
00198     
00200     static void Apply ( const GmresPolyOp<ScalarType,MV,OP>& Op, 
00201       const MV& x, MV& y,
00202       ETrans trans=NOTRANS )
00203     { Op.Apply( x, y, trans ); }
00204     
00205   };
00206   
00207 } // end Belos namespace
00208 
00209 #endif
00210 
00211 // end of file BelosGmresPolyOp.hpp

Generated on Tue Jul 13 09:27:02 2010 for Belos by  doxygen 1.4.7