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