Belos Package Browser (Single Doxygen Collection) Development
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_GMRESPOLYOP_HPP
00043 #define BELOS_GMRESPOLYOP_HPP
00044 
00049 #include "BelosOperatorTraits.hpp"
00050 #include "BelosMultiVecTraits.hpp"
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosConfigDefs.hpp"
00053 #include "Teuchos_RCP.hpp"
00054 #include "Teuchos_SerialDenseMatrix.hpp"
00055 #include "Teuchos_SerialDenseVector.hpp"
00056 
00068 namespace Belos {
00069   
00070   template <class ScalarType, class MV, class OP>
00071   class GmresPolyOp {
00072   public:
00073     
00075 
00076     
00078     GmresPolyOp() {}
00079     
00081     GmresPolyOp( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem_in, 
00082                  const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> >& hess,
00083                  const Teuchos::RCP<Teuchos::SerialDenseMatrix<int, ScalarType> >& comb,
00084                  const Teuchos::RCP<Teuchos::SerialDenseVector<int, ScalarType> >& scal
00085                ) 
00086       : problem_(problem_in), LP_(problem_in->getLeftPrec()), RP_(problem_in->getRightPrec()), H_(hess), y_(comb), r0_(scal)
00087     {
00088       dim_ = y_->numRows();
00089     }
00090 
00092     virtual ~GmresPolyOp() {};
00094     
00096 
00097     
00103     void Apply ( const MV& x, MV& y, ETrans trans=NOTRANS );
00104 
00105     private:
00106 
00107     typedef MultiVecTraits<ScalarType,MV> MVT;
00108     typedef Teuchos::ScalarTraits<ScalarType> SCT ;
00109 
00110     int dim_;
00111     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00112     Teuchos::RCP<const OP> LP_, RP_;
00113     Teuchos::RCP<MV> V_, wL_, wR_;
00114     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_, y_;
00115     Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > r0_;
00116   };
00117   
00118   template <class ScalarType, class MV, class OP>
00119   void GmresPolyOp<ScalarType, MV, OP>::Apply( const MV& x, MV& y, ETrans trans ) 
00120   {
00121     // Initialize vector storage.
00122     if (V_ == Teuchos::null) {
00123       V_ = MVT::Clone( x, dim_ );
00124       if (LP_ != Teuchos::null) {
00125         wL_ = MVT::Clone( y, 1 );
00126       }
00127       if (RP_ != Teuchos::null) {
00128         wR_ = MVT::Clone( y, 1 );
00129       }
00130     }
00131     //
00132     // Apply polynomial to x.
00133     // 
00134     int n = MVT::GetNumberVecs( x );
00135     std::vector<int> idxi(1), idxi2, idxj(1);
00136 
00137     // Select vector x[j].
00138     for (int j=0; j<n; ++j) {
00139 
00140       idxi[0] = 0;
00141       idxj[0] = j;
00142       Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
00143       Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
00144       if (LP_ != Teuchos::null) {
00145         Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
00146         problem_->applyLeftPrec( *x_view, *v_curr ); // Left precondition x into the first vector of V
00147       } else {
00148         MVT::SetBlock( *x_view, idxi, *V_ );  // Set x as the first vector of V
00149       }
00150 
00151       for (int i=0; i<dim_-1; ++i) {
00152 
00153         // Get views into the current and next vectors
00154         idxi2.resize(i+1);
00155         for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
00156         Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
00157         // the tricks below with wR_ and wL_ (potentially set to v_curr and v_next) unfortunately imply that 
00158         // v_curr and v_next must be non-const views.
00159         Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
00160         idxi[0] = i+1;
00161         Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
00162 
00163         //---------------------------------------------
00164         // Apply operator to next vector
00165         //---------------------------------------------
00166         // 1) Apply right preconditioner, if we have one.
00167         if (RP_ != Teuchos::null) {
00168           problem_->applyRightPrec( *v_curr, *wR_ );
00169         } else {
00170           wR_ = v_curr;
00171         }
00172         // 2) Check for right preconditioner, if none exists, point at the next vector.
00173         if (LP_ == Teuchos::null) {
00174           wL_ = v_next;
00175         }
00176         // 3) Apply operator A.
00177         problem_->applyOp( *wR_, *wL_ );
00178         // 4) Apply left preconditioner, if we have one.
00179         if (LP_ != Teuchos::null) {
00180           problem_->applyLeftPrec( *wL_, *v_next );
00181         }
00182 
00183         // Compute A*v_curr - v_prev*H(1:i,i)
00184         Teuchos::SerialDenseMatrix<int,ScalarType> h(Teuchos::View,*H_,i+1,1,0,i);
00185         MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
00186 
00187         // Scale by H(i+1,i)
00188         MVT::MvScale( *v_next, SCT::one()/(*H_)(i+1,i) );  
00189       }
00190 
00191       // Compute output y = V*y_./r0_
00192       if (RP_ != Teuchos::null) {
00193         MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *wR_ );
00194         problem_->applyRightPrec( *wR_, *y_view );
00195       } 
00196       else {
00197         MVT::MvTimesMatAddMv( SCT::one()/(*r0_)(0), *V_, *y_, SCT::zero(), *y_view );
00198       }
00199     } // (int j=0; j<n; ++j)
00200   } // end Apply()
00201 
00203   //
00204   // Implementation of the Belos::OperatorTraits for Belos::GmresPolyOp 
00205   //
00207 
00211   template <class ScalarType, class MV, class OP> 
00212   class OperatorTraits < ScalarType, MV, GmresPolyOp<ScalarType,MV,OP> > 
00213   {
00214   public:
00215     
00217     static void Apply ( const GmresPolyOp<ScalarType,MV,OP>& Op, 
00218       const MV& x, MV& y,
00219       ETrans trans=NOTRANS )
00220     { Op.Apply( x, y, trans ); }
00221     
00222   };
00223   
00224 } // end Belos namespace
00225 
00226 #endif
00227 
00228 // end of file BelosGmresPolyOp.hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines