Belos Package Browser (Single Doxygen Collection) Development
BelosStandardGmres.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_StandardGmres_hpp
00043 #define __Belos_StandardGmres_hpp
00044 
00048 
00049 #include <BelosGmresBase.hpp>
00050 
00051 namespace Belos {
00052 
00073   template<class Scalar, class MV, class OP>
00074   class StandardGmres : public GmresBase<Scalar,MV,OP> {
00075   public:
00076     typedef Scalar scalar_type;
00077     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00078     typedef MV multivector_type;
00079     typedef OP operator_type;
00080 
00081   private:
00084     typedef GmresBase<Scalar, MV, OP> base_type;
00085     typedef Teuchos::ScalarTraits<Scalar> STS;
00086     typedef MultiVecTraits<Scalar, MV> MVT;
00087     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00088     
00089   public:
00108     StandardGmres (const Teuchos::RCP<LinearProblem<Scalar,MV,OP> >& lp,
00109        const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00110        const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00111        const int maxIterCount,
00112        const bool flexible) :
00113       GmresBase<Scalar, MV, OP> (lp, ortho, outMan, maxIterCount, flexible) {}
00114 
00116     virtual ~StandardGmres() {}
00117 
00118     virtual bool canExtendBasis() const {
00119       return this->getNumIters() < this->maxNumIters();
00120     }
00121 
00122     virtual void
00123     extendBasis (Teuchos::RCP<MV>& V_cur, 
00124      Teuchos::RCP<MV>& Z_cur)
00125     { 
00126       using Teuchos::RCP;
00127       using Teuchos::Range1D;
00128       using std::endl;
00129       const bool verboseDebug = false;
00130       //
00131       // mfh 16 Feb 2011: The use of "this->..." here and elsewhere to
00132       // refer to GmresBase member data is obligatory, since we are
00133       // inheriting from a templated class.  See the C++ FAQ:
00134       //
00135       // http://www.parashift.com/c++-faq-lite/templates.html#faq-35.19
00136       // 
00137       RCP<LinearProblem<Scalar, MV, OP> > lp = this->lp_;
00138       RCP<MV> V = this->V_;
00139       RCP<MV> Z = this->Z_;
00140       // This does not count the initial basis vector.
00141       const int k = this->getNumIters(); 
00142       const int m = this->maxNumIters();
00143       TEUCHOS_TEST_FOR_EXCEPTION(k >= m, GmresCantExtendBasis,
00144        "Belos::StandardGmres::extendBasis: "
00145        "Maximum number of iterations " << m << "reached; "
00146        "cannot extend basis further.");
00147       std::ostream& dbg = this->outMan_->stream(Debug);
00148       if (verboseDebug)
00149   {
00150     dbg << "---- StandardGmres::extendBasis: "
00151       "V_prv = V[" << k << "," << k << "], "
00152       "V_cur = V[" << k+1 << "," << k+1 << "]";
00153     if (this->flexible_)
00154       dbg << ", Z_cur = Z[" << k << "," << k << "]" << endl;
00155     else
00156       dbg << endl;
00157   }
00158       RCP<const MV> V_prv = MVT::CloneView(*V, Range1D(k, k));
00159       V_cur = MVT::CloneViewNonConst(*V, Range1D(k+1, k+1));
00160       if (this->flexible_)
00161   {
00162     Z_cur = MVT::CloneViewNonConst(*Z, Range1D(k, k));
00163     lp->applyOp (*V_prv, *Z_cur);
00164     lp->applyRightPrec (*Z_cur, *V_cur);
00165   }
00166       else
00167   {
00168     lp->apply (*V_prv, *V_cur);
00169     if (! Z_cur.is_null())
00170       Z_cur = Teuchos::null;
00171   }
00172     }
00173 
00174     virtual void
00175     orthogonalize (const Teuchos::RCP<MV>& V_cur,
00176        const Teuchos::RCP<const MV>& V_prv,
00177        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_V,
00178        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_V,
00179        const Teuchos::RCP<MV>& Z_cur,
00180        const Teuchos::RCP<const MV>& Z_prv,
00181        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_Z,
00182        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_Z)
00183     {
00184       using Teuchos::null;
00185       using Teuchos::rcp;
00186       using Teuchos::tuple;
00187       const bool verboseDebug = true;
00188 
00189       // Standard GMRES doesn't generate projection and normalization
00190       // coefficients for the Z basis, regardless of whether we are
00191       // performing Flexible GMRES.
00192       if (! C_Z.is_null())
00193   C_Z = null;
00194       if (! B_Z.is_null())
00195   B_Z = null;
00196       // Help C++'s type inference by referring to the "this" pointer.
00197       mat_type& H = this->projectedProblem_->H;
00198       mat_type& z = this->projectedProblem_->z;
00199       const int k = this->getNumIters();
00200 
00201       // Flexible standard GMRES only needs to orthogonalize the new V
00202       // basis vector; it doesn't need to orthogonalize the Z basis
00203       // vector.
00204       //
00205       // We don't need to do anything with the "rank" output of
00206       // projectAndNormalize(), since the single value in B will tell
00207       // us (or rather, the caller) what we (rather, they) need to
00208       // know.
00209       C_V = rcp (new mat_type (Teuchos::View, H, k+1, 1, 0, k));
00210       B_V = rcp (new mat_type (Teuchos::View, H, 1, 1, k+1, k));
00211       (void) this->ortho_->projectAndNormalize (*V_cur, tuple(C_V), 
00212             B_V, tuple(V_prv));
00213       if (verboseDebug)
00214   {
00215     using std::endl;
00216     std::ostream& dbg = this->outMan_->stream(Debug);
00217     dbg << "---- Current upper Hessenberg matrix (" 
00218         << k+2 << " x " << k+1 << "): "
00219         << mat_type(Teuchos::View, H, k+2, k+1) << endl;
00220     dbg << "---- Current projected right-hand side ("
00221         << k+2 << " x " << 1 << "): "
00222         << mat_type(Teuchos::View, z, k+2, 1) << endl;
00223   }
00224     }
00225 
00226     virtual bool 
00227     acceptedCandidateBasis () const 
00228     {
00229       // We haven't accepted the candidate basis yet, so the number of
00230       // iterations hasn't yet been incremented.
00231       const int k = this->getNumIters();
00232       const Scalar H_kp1k = this->projectedProblem_->H(k+1, k);
00233       
00234       // NOTE (mfh {15,16} Jan 2011) This test should perhaps be more
00235       // sophisticated.  However, perhaps the right place for such a
00236       // test is the status check, rather than here.  Certainly if
00237       // H_(k+1,k) is zero, the basis cannot be extended, even in
00238       // exact arithmetic.
00239       return STS::magnitude (H_kp1k) > STS::zero();
00240     }
00241 
00242     virtual void 
00243     updateUpperHessenbergMatrix (const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_V,
00244          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_V,
00245          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_Z,
00246          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_Z)
00247     {
00248       // Standard GMRES just writes to the upper Hessenberg matrix in
00249       // place in its implementation of orthogonalize(), so we don't
00250       // need to do anything here.  Just include the usual boilerplate
00251       // to avoid compiler warnings for unused inputs.
00252       (void) C_V;
00253       (void) B_V;
00254       (void) C_Z;
00255       (void) B_Z;
00256     }
00257 
00258   };
00259 
00260 } // namespace Belos
00261 
00262 #endif // __Belos_StandardGmres_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines