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