Belos Version of the Day
BelosGmresBase.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_GmresBase_hpp
00043 #define __Belos_GmresBase_hpp
00044 
00048 
00049 #include <BelosLinearProblem.hpp>
00050 #include <BelosOrthoManager.hpp>
00051 #include <BelosMultiVecTraits.hpp>
00052 #include <BelosProjectedLeastSquaresSolver.hpp>
00053 
00054 #include <Teuchos_BLAS.hpp>
00055 #include <Teuchos_LAPACK.hpp>
00056 #include <Teuchos_ScalarTraits.hpp>
00057 #include <Teuchos_SerialDenseVector.hpp>
00058 #include <Teuchos_TypeNameTraits.hpp>
00059 
00060 #include <algorithm>
00061 #include <iterator>
00062 #include <utility> // std::pair
00063 
00064 // Anonymous namespace is a more portable way than "static" to define
00065 // functions that won't escape the current file's scope.
00066 namespace {
00067   template<class Ordinal, class Scalar>
00068   std::ostream& 
00069   operator<< (std::ostream& out,
00070         const Teuchos::SerialDenseMatrix<Ordinal, Scalar>& A)
00071   {
00072     const Ordinal m = A.numRows();
00073     const Ordinal n = A.numCols();
00074     using std::endl;
00075 
00076     // Summarize the matrix if it's too big to print comfortably on
00077     // the terminal.
00078     //
00079     // FIXME (mfh 25 Feb 2011) Obviously this is no substitute for a
00080     // "curses"-type terminal interface that can query the terminal
00081     // for its dimensions.
00082     if (m > 10 || n > 10)
00083       {
00084   out << "[<" << m << " x " << n << " SerialDenseMatrix<int," 
00085       << Teuchos::TypeNameTraits<Scalar>::name() << ">]";
00086   return out;
00087       }
00088     else if (m == 0 || n == 0)
00089       { // Empty matrix: either zero rows or zeo columns
00090   out << "[]";
00091   return out;
00092       }
00093     out << "[";
00094     if (m > 1)
00095       out << endl;
00096     for (Ordinal i = 0; i < m; ++i)
00097       { // Indent if printing on more than one row.
00098   if (m > 1)
00099     out << "  ";
00100   // Print row i of the matrix: comma-delimited, like Matlab.
00101   for (Ordinal j = 0; j < n; ++j)
00102     {
00103       out << A(i,j);
00104       if (j < n-1)
00105         out << ", ";
00106     }
00107   // End of the current row, which we know is not empty.
00108   if (i < m-1 && m > 1)
00109     out << ";" << endl;
00110       }
00111     if (m > 1)
00112       out << endl;
00113     out << "]";
00114     return out;
00115   }
00116 
00117 } // (anonymous namespace)
00118 
00119 namespace Belos {
00120 
00142   class GmresCantExtendBasis : public BelosError {
00143   public:
00144     GmresCantExtendBasis (const std::string& what_arg) : 
00145       BelosError(what_arg) {}
00146   };
00147 
00186   class GmresRejectsCandidateBasis : public BelosError {
00187   public:
00188     GmresRejectsCandidateBasis (const std::string& what_arg) : 
00189       BelosError(what_arg) {}
00190   };
00191 
00192 
00225   template<class Scalar, class MV, class OP>
00226   class GmresBase {
00227   public:
00228     typedef Scalar scalar_type;
00229     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00230     typedef MV multivector_type;
00231     typedef OP operator_type;
00232 
00233   private:
00234     typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00235     typedef Teuchos::SerialDenseVector<int,Scalar> vec_type;
00236     typedef Teuchos::BLAS<int, scalar_type> blas_type;
00237     typedef Teuchos::LAPACK<int, scalar_type> lapack_type;
00238     typedef MultiVecTraits<scalar_type, MV> MVT;
00239     typedef Teuchos::ScalarTraits<Scalar> STS;
00240     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00241 
00242   public:
00279     GmresBase (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
00280          const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00281          const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00282          const int maxIterCount,
00283          const bool flexible);
00284 
00286     bool canAdvance () const { return canExtendBasis(); }
00287 
00301     void advance();
00302 
00308     virtual void 
00309     acceptCandidateBasis (const int newNumVectors = 1) 
00310     {
00311       TEUCHOS_TEST_FOR_EXCEPTION(newNumVectors <= 0, std::logic_error,
00312        "The number of candidate basis vectors to accept, " 
00313        << newNumVectors << ", is nonpositive, which is not "
00314        "allowed.  This is an std:logic_error because it "
00315        "relates to the implementation of the GmresBase "
00316        "subclass.");
00317       curNumIters_ += newNumVectors;
00318     }
00319 
00324     virtual void 
00325     rejectCandidateBasis (const int newNumVectors = 1) 
00326     {
00327       std::ostringstream os;
00328       os << (flexible_ ? "Flexible GMRES" : "GMRES") 
00329    << " rejects the computed candidate basis vector"
00330    << (newNumVectors != 1 ? "s" : "");
00331       throw GmresRejectsCandidateBasis (os.str());
00332     }
00333 
00337     int getNumIters() const { return curNumIters_; }
00338 
00342     int maxNumIters() const { return maxNumIters_; }
00343 
00367     void restart (const int maxIterCount);
00368 
00370     void restart () { restart (maxNumIters_); }
00371 
00380     void backOut (const int numIters);
00381 
00386     Teuchos::RCP<MV> getCurrentUpdate ();
00387 
00395     Teuchos::RCP<MV> currentNativeResidualVector ();
00396 
00403     typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00404     currentNativeResidualNorm ();
00405 
00423     typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00424     arnoldiRelationError () const;
00425 
00432     void
00433     updateSolution ()
00434     {
00435       (void) lp_->updateSolution (xUpdate_, true, STS::one());
00436     }
00437 
00438   protected:
00439 
00444     //{@
00445 
00453     virtual bool canExtendBasis() const = 0;
00454 
00480     virtual void 
00481     extendBasis (Teuchos::RCP<MV>& V_cur, 
00482      Teuchos::RCP<MV>& Z_cur) = 0;
00483 
00529     virtual void
00530     orthogonalize (const Teuchos::RCP<MV>& V_cur,
00531        const Teuchos::RCP<const MV>& V_prv,
00532        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_V,
00533        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_V,
00534        const Teuchos::RCP<MV>& Z_cur,
00535        const Teuchos::RCP<const MV>& Z_prv,
00536        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_Z,
00537        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_Z) = 0;
00538     
00566     virtual void 
00567     updateUpperHessenbergMatrix (const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_V,
00568          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_V,
00569          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_Z,
00570          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_Z) = 0;
00571 
00573     virtual bool acceptedCandidateBasis() const = 0;
00574 
00587     virtual bool shouldRetryAdvance() { return false; }
00588 
00590 
00620 
00621 
00633     virtual void preRestartHook (const int maxIterCount) {}
00634 
00645     virtual void postRestartHook() {}
00646 
00657     virtual void 
00658     preOrthogonalizeHook (const Teuchos::RCP<MV>& V_cur, 
00659         const Teuchos::RCP<MV>& Z_cur) {}
00660 
00671     void
00672     postOrthogonalizeHook (const Teuchos::RCP<MV>& V_cur, 
00673          const Teuchos::RCP<mat_type>& C_V,
00674          const Teuchos::RCP<mat_type>& B_V,
00675          const Teuchos::RCP<MV>& Z_cur,
00676          const Teuchos::RCP<mat_type>& C_Z,
00677          const Teuchos::RCP<mat_type>& B_Z) {}
00679 
00680   private:
00681 
00683 
00684 
00705     void 
00706     setInitialResidual(const int maxIterCount,
00707            const bool reallocateBasis = false);
00708 
00722     typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00723     updateProjectedLeastSquaresProblem ();
00724 
00737     void 
00738     previousVectors (Teuchos::RCP<const MV>& V_prv,
00739          Teuchos::RCP<const MV>& Z_prv) const;
00740 
00754     Teuchos::RCP<const MV> initResVec() const {
00755       return initResVec_;
00756     }
00757 
00759 
00760   protected:
00762 
00763 
00765     Teuchos::RCP<LinearProblem<Scalar, MV, OP> > lp_;
00766 
00768     Teuchos::RCP<const OrthoManager<Scalar, MV> > ortho_;
00769 
00771     Teuchos::RCP<OutputManager<Scalar> > outMan_;
00772 
00796     Teuchos::RCP<MV> nativeResVec_;
00797 
00803     Teuchos::RCP<MV> initResVec_;
00804     
00812     Teuchos::RCP<MV> xUpdate_;
00813 
00820     Teuchos::RCP<MV> V_;
00821 
00827     Teuchos::RCP<MV> Z_;
00828 
00835     Teuchos::RCP<details::ProjectedLeastSquaresProblem<Scalar> > projectedProblem_;
00836 
00844     typename Teuchos::ScalarTraits<Scalar>::magnitudeType initialResidualNorm_;
00845 
00848     // This is the absolute residual error from solving the projected
00849     // least-squares problem, if we've finished at least one iteration
00850     // of (F)GMRES.  Otherwise, it's the same as the initial residual
00851     // norm.
00852     typename Teuchos::ScalarTraits<Scalar>::magnitudeType nativeResidualNorm_;
00853 
00858     int lastUpdatedCol_;
00859 
00865     int curNumIters_;
00866 
00875     int maxNumIters_;
00876 
00906     bool flexible_;
00907 
00909   };
00910 
00913 
00914   template<class Scalar, class MV, class OP>
00915   void
00916   GmresBase<Scalar, MV, OP>::
00917   previousVectors (Teuchos::RCP<const MV>& V_prv,
00918        Teuchos::RCP<const MV>& Z_prv) const
00919   {
00920     using Teuchos::Range1D;
00921     const bool verboseDebug = false;
00922 
00923     // The number of iterations (returned by getNumIters()) does not
00924     // include the initial basis vector, which is the first column of
00925     // V_.  Hence, the range of previous basis vectors (in the V
00926     // basis) is [0,k] inclusive.
00927     const int k = getNumIters(); 
00928 
00929     if (verboseDebug)
00930       {
00931   using std::endl;
00932   std::ostream& dbg = outMan_->stream(Debug);
00933   dbg << "---- previousVectors: V_prv = V[0:" << k << "]";
00934   if (flexible_)
00935     dbg << ", Z_prv = Z[0:" << k << "]" << endl;
00936   else 
00937     dbg << endl;
00938       }
00939     
00940     V_prv = MVT::CloneView (*V_, Range1D(0,k));
00941     if (flexible_)
00942       // FIXME (mfh 24 Feb 2011) Is this right?  When should we apply
00943       // the right preconditioner in Flexible GMRES?
00944       Z_prv = MVT::CloneView (*Z_, Range1D(0,k));
00945     else
00946       Z_prv = Teuchos::null;
00947   }
00948 
00949   template<class Scalar, class MV, class OP>
00950   void
00951   GmresBase<Scalar, MV, OP>::advance () 
00952   {
00953     using Teuchos::is_null;
00954     using Teuchos::null;
00955     using Teuchos::Range1D;
00956     using Teuchos::RCP;
00957     using Teuchos::rcp;
00958     using Teuchos::tuple;
00959     using std::endl;
00960     typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00961     const bool verboseDebug = false;
00962 
00963     std::ostream& dbg = outMan_->stream(Debug);
00964     dbg << "Belos::GmresBase::advance():" << endl
00965   << "--" << endl
00966   << "-- Iteration " << getNumIters() << endl
00967   << "--" << endl;
00968 
00969     TEUCHOS_TEST_FOR_EXCEPTION( !canExtendBasis(), GmresCantExtendBasis,
00970       "GMRES (iteration " << getNumIters() << ") cannot "
00971       "add any more basis vectors." );
00972 
00973     RCP<const MV> V_prv, Z_prv;
00974     RCP<MV> V_cur, Z_cur;
00975     // Fetch the previous basis vector(s), against which to project.
00976     previousVectors (V_prv, Z_prv);
00977     if (verboseDebug)
00978     dbg << "---- Number of previous vectors: " 
00979       << MVT::GetNumberVecs(*V_prv) << endl;
00980     // Compute new basis vector(s).
00981     extendBasis (V_cur, Z_cur);
00982 
00983     // Do whatever the subclass wants to do before orthogonalizing the
00984     // new basis vectors.  This might include projecting them against
00985     // a fixed basis (e.g., for Krylov subspace recycling).
00986     preOrthogonalizeHook (V_cur, Z_cur);
00987 
00988     // Flexible CA-GMRES must orthogonalize the Z basis along with the
00989     // V basis, and it must keep both sets of orthogonalization
00990     // coefficients.  Flexible standard GMRES only needs to
00991     // orthogonalize the new V basis vector; it doesn't need to
00992     // orthogonalize the Z basis vector.  This is why we let the
00993     // subclass implement this functionality, rather than calling the
00994     // OrthoManager's projectAndNormalize() method directly; the
00995     // subclass has to decide whether it needs to orthogonalize Z_cur
00996     // as well as V_cur.
00997     //
00998     // Subclasses may want to save the rank (if applicable) from
00999     // ortho_->projectAndNormalize() somewhere.  In the case of
01000     // Flexible CA-GMRES, there are _two_ ranks -- one for the V
01001     // basis, and the other for the Z basis.
01002     //
01003     // C_V and B_V are the "C" (projection) resp. "B" (normalization)
01004     // block coefficients from ortho_->projectAndNormalize() on the
01005     // new V basis vectors.  If the new Z basis vectors need to be
01006     // orthogonalized as well, then C_Z and B_Z are the "C" resp. "B"
01007     // block coefficients from ortho_->projectAndNormalize() on the
01008     // new Z basis vectors.  Subclasses' implementations of
01009     // orthogonalize() are responsible for allocating space for
01010     // C_V and B_V, and C_Z and B_Z if necessary (otherwise they are
01011     // set to Teuchos::null, as the RCPs are passed by reference).
01012     RCP<mat_type> C_V, B_V, C_Z, B_Z;
01013     orthogonalize (V_cur, V_prv, C_V, B_V, Z_cur, Z_prv, C_Z, B_Z);
01014     if (false && verboseDebug)
01015       {
01016   dbg << "---- Projection coefficient(s): ";
01017   dbg << *C_V << endl;
01018   dbg << "---- Normalization coefficient(s): ";
01019   dbg << *B_V << endl;
01020       }
01021     // Do whatever the subclass wants to do after projecting the new
01022     // basis vectors against the previous ones and normalizing them.
01023     postOrthogonalizeHook (V_cur, C_V, B_V, Z_cur, C_Z, B_Z);
01024 
01025     // Implemented by subclasses.  C_Z and B_Z need not be used,
01026     // regardless of whether we are computing Flexible GMRES.  (Only
01027     // Flexible CA-GMRES needs them, as far as we know.)  If they are
01028     // not used, orthogonalize() may set them to null (the RCPs are
01029     // passed by reference).  Some subclasses may choose to update the
01030     // upper Hessenberg matrix "in place" in their implementation of
01031     // the orthogonalize() method (see above).
01032     updateUpperHessenbergMatrix (C_V, B_V, C_Z, B_Z);
01033 
01034     // The subclass decides whether or not to accept the candidate
01035     // basis.
01036     //
01037     // We update the upper Hessenberg matrix before deciding
01038     // acceptance, because the acceptance critera might depend on the
01039     // entire upper Hessenberg matrix (e.g., (an estimate of) its
01040     // condition number), not just the condition number or rank of the
01041     // candidate basis vector(s).
01042     //
01043     // Subclasses' implementations of orthogonalize() and/or
01044     // updateUpperHessenbergMatrix() might wish to set state to
01045     // indicate acceptance or rejection of the candidate basis
01046     // vector(s).  They may also call advance() recursively (e.g.,
01047     // setting the "last rejected rank" so that the next call of
01048     // advance() computes fewer candidate basis vectors).
01049     if (! acceptedCandidateBasis())
01050       rejectCandidateBasis (MVT::GetNumberVecs(*V_cur));
01051     else 
01052       // This increments the current iteration count by the number of
01053       // accepted candidate basis vectors.
01054       acceptCandidateBasis (MVT::GetNumberVecs(*V_cur));
01055 
01056     // We don't update the projected least-squares problem here.
01057     // This is done lazily, whenever the "native" residual norm
01058     // or the current solution update are requested.
01059   }
01060 
01061 
01062   template<class Scalar, class MV, class OP>
01063   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
01064   GmresBase<Scalar, MV, OP>::arnoldiRelationError () const
01065   {
01066     using Teuchos::Range1D;
01067     using Teuchos::RCP;
01068     const Scalar one = STS::one();
01069     const int m = getNumIters();
01070 
01071     if (m == 0) // No iterations completed means no error (yet)
01072       return magnitude_type(0);
01073     RCP<const MV> V_mp1 = MVT::CloneView (*V_, Range1D(0,m));
01074     const mat_type H_m (Teuchos::View, projectedProblem_->H, m+1, m);
01075     std::vector<magnitude_type> norms (m);
01076     if (flexible_)
01077       {
01078   RCP<const MV> Z_view = MVT::CloneView (*Z_, Range1D(0,m-1));
01079   // A*Z is in the same space as V, so we clone from V
01080   RCP<MV> AZ = MVT::Clone (*V_, m);
01081   // Compute AZ := A*Z_m
01082   lp_->applyOp (Z_view, AZ);
01083   // Compute AZ := AZ - V_{m+1} \underline{H}_m
01084   MVT::MvTimesMatAddMv (-one, V_mp1, H_m, one, AZ);
01085   MVT::MvNorm (AZ, norms);
01086       }
01087     else
01088       {
01089   RCP<const MV> V_m = MVT::CloneView (*V_, Range1D(0,m-1));
01090   RCP<MV> AV = MVT::Clone (*V_, m);
01091   // Compute AV := A*V_m.  Here, "A" means the preconditioned
01092   // operator.
01093   lp_->apply (V_m, AV);
01094   // Compute AV := A V_m - V_{m+1} \underline{H}_m
01095   MVT::MvTimesMatAddMv (-one, V_mp1, H_m, one, AV);
01096   MVT::MvNorm (AV, norms);
01097       }
01098     // Compute and return the Frobenius norm of the above (either AZ
01099     // or AV, depending on whether we are performing Flexible GMRES).
01100     //
01101     // FIXME (mfh 16 Feb 2011) This computation might overflow for
01102     // unjustifiable reasons.  I should really include intermediate
01103     // scaling.
01104     magnitude_type result (0);
01105     for (int k = 0; k < m; ++k)
01106       result += norms[k];
01107     return STM::squareroot(result);
01108   }
01109 
01110 
01111   template<class Scalar, class MV, class OP>
01112   GmresBase<Scalar, MV, OP>::
01113   GmresBase (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
01114        const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
01115        const Teuchos::RCP<OutputManager<Scalar> >& outMan,
01116        const int maxIterCount,
01117        const bool flexible) :
01118     lp_ (problem), 
01119     ortho_ (ortho),
01120     outMan_ (outMan),
01121     lastUpdatedCol_ (-1), // column updates start with zero
01122     curNumIters_ (0),
01123     maxNumIters_ (maxIterCount),
01124     flexible_ (flexible && ! lp_.is_null() && ! lp_->getRightPrec().is_null())
01125   {
01126     using Teuchos::is_null;
01127     using Teuchos::null;
01128     using Teuchos::RCP;
01129     using Teuchos::rcp;
01130     using std::endl;
01131 
01132     std::ostream& dbg = outMan_->stream(Debug);
01133     dbg << "Belos::GmresBase:" << endl;
01134 
01135     // Fragments of error messages for use in sanity checks.
01136     const char prefix[] = "Belos::GmresBase constructor: ";
01137     const char postfix[] = "  Please call setProblem() with non-null data "
01138       "before passing the LinearProblem to the GmresBase constructor.";
01139     const char postfix2[] = "  That may mean that the LinearProblem's "
01140       "setLSIndex() method has not yet been called, even though setProblem() "
01141       "has been called.  After calling setProblem(), you should call "
01142       "setLSIndex() with the ind{ex,ices} of the right-hand side(s) to solve.";
01143 
01144     // Sanity checks on the linear problem.
01145     //
01146     // First, make sure A, X, and B are all non-null.
01147     TEUCHOS_TEST_FOR_EXCEPTION(lp_.is_null(), std::invalid_argument,
01148            prefix << "The given LinearProblem is null.");
01149     TEUCHOS_TEST_FOR_EXCEPTION(! lp_->isProblemSet(), std::invalid_argument,
01150                prefix << "The given LinearProblem has not yet been set."
01151            << postfix);
01152     TEUCHOS_TEST_FOR_EXCEPTION(lp_->getLHS().is_null(), std::invalid_argument,
01153                prefix << "The given LinearProblem has null initial guess"
01154            " (getLHS()) for all right-hand sides." << postfix);
01155     TEUCHOS_TEST_FOR_EXCEPTION(lp_->getRHS().is_null(), std::invalid_argument,
01156                prefix << "The given LinearProblem's right-hand side(s) "
01157            "are all null (getRHS().is_null() == true)." << postfix);
01158     TEUCHOS_TEST_FOR_EXCEPTION(lp_->getOperator().is_null(), std::invalid_argument,
01159                prefix << "The given LinearProblem's operator (the "
01160            "matrix A) is null." << postfix);
01161     // Next, make sure that setLSIndex() has been called on the linear
01162     // problem instance, so that the "current" right-hand side and
01163     // "current" initial guess have been set.
01164     TEUCHOS_TEST_FOR_EXCEPTION(lp_->getCurrLHSVec().is_null(), 
01165            std::invalid_argument,
01166                prefix << "Although the given LinearProblem has non-null "
01167            "initial guess (getLHS()) for all right-hand sides, the "
01168            "current initial guess (getCurrLHSVec()) is null."
01169            << postfix2);
01170     TEUCHOS_TEST_FOR_EXCEPTION(lp_->getCurrRHSVec().is_null(), 
01171            std::invalid_argument,
01172                prefix << "Although the given LinearProblem has non-null "
01173            "initial guess (getLHS()) for all right-hand sides, the "
01174            "current right-hand side to solve (getCurrRHSVec()) is "
01175            "null." << postfix2);
01176     // Get the initial guess x0, and allocate the solution update
01177     // (xUpdate_) from the same vector space as x0.
01178     RCP<const MV> x0 = lp_->getCurrLHSVec();
01179     {
01180       const int numLHS = MVT::GetNumberVecs (*x0);
01181       TEUCHOS_TEST_FOR_EXCEPTION(numLHS != 1, std::invalid_argument,
01182        "Our GMRES implementation only works for single-"
01183        "vector problems, but the supplied initial guess has "
01184        << numLHS << " columns.");
01185       const int numRHS = MVT::GetNumberVecs (*(lp_->getCurrRHSVec()));
01186       TEUCHOS_TEST_FOR_EXCEPTION(numRHS != 1, std::invalid_argument,
01187        "Our GMRES implementation only works for single-"
01188        "vector problems, but the current right-hand side has "
01189        << numRHS << " columns.");
01190     }
01191     dbg << "-- Initializing xUpdate_ and filling with zeros" << endl;
01192     xUpdate_ = MVT::Clone (*x0, 1);
01193     MVT::MvInit (*xUpdate_, STS::zero());
01194 
01195     // Initialize the projected least-squares problem.
01196     projectedProblem_ = 
01197       rcp (new details::ProjectedLeastSquaresProblem<Scalar> (maxIterCount));
01198 
01199     // Get the (left preconditioned, if applicable) residual vector,
01200     // and make a deep copy of it in initResVec_.  Allocate the V_
01201     // basis vectors, compute the initial residual norm, and compute
01202     // the first basis vector (first column of V_).
01203     setInitialResidual (maxIterCount);
01204 
01205     // The Z_ vectors, if we need them, are always in the same vector
01206     // space as the initial solution guess (x0).  Even if the caller
01207     // asks for the flexible option, we don't allocate space for Z_
01208     // unless the caller specifies a right preconditioner.
01209     Z_ = (flexible && ! lp_->getRightPrec().is_null()) ? 
01210       MVT::Clone(*x0, maxIterCount+1) : null;
01211   }
01212 
01213   template<class Scalar, class MV, class OP>
01214   Teuchos::RCP<MV> 
01215   GmresBase<Scalar, MV, OP>::currentNativeResidualVector ()
01216   {
01217     const Scalar zero = STS::zero();
01218     const Scalar one = STS::one();
01219     using Teuchos::is_null;
01220     using Teuchos::Range1D;
01221     using Teuchos::RCP;
01222     using std::endl;
01223 
01224     const bool verboseDebug = true;
01225     const char prefix[] = "Belos::GmresBase::currentNativeResidualVector: ";
01226     const int m = getNumIters();
01227 
01228     std::ostream& dbg = outMan_->stream(Debug);
01229     if (verboseDebug)
01230       dbg << "-- currentNativeResidualVector: getNumIters() = " << m << endl;
01231 
01232     TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*initResVec_) != 1,
01233            std::logic_error,
01234            prefix << "Initial residual vector (initResVec_) has " 
01235            << MVT::GetNumberVecs(*initResVec_)
01236            << " columns, but should only have 1 column.  "
01237            "This should never happen.");
01238     // If left preconditioning is used, then the "native" residual
01239     // vector is in the same vector space as the left-preconditioned
01240     // initial residual vector.  Otherwise, it is in the same space as
01241     // the right-hand side b and the unpreconditioned initial residual
01242     // vector.  Regardless, it's always in the same space as
01243     // *initResVec_, from which we clone it.  Allocate only if
01244     // necessary.
01245     if (nativeResVec_.is_null() || MVT::GetNumberVecs (*nativeResVec_) != 1)
01246       nativeResVec_ = MVT::Clone (*initResVec_, 1);
01247 
01248     // The native residual vector is the same as the initial
01249     // residual vector when the number of iterations is 0.
01250     // Otherwise, we will update *nativeResVec_ below.
01251     MVT::Assign (*initResVec_, *nativeResVec_);
01252     if (m > 0)
01253       {
01254   TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*V_) < m+1,
01255          std::logic_error,
01256          "Only " << MVT::GetNumberVecs(*V_) << " basis vectors "
01257          "were given, but getNumIters()+1=" << (m+1) 
01258          << " of them are required.  "
01259          "This likely indicates a bug in Belos.");
01260   // Update the QR factorization of the upper Hessenberg matrix,
01261   // and let y[0:m-1] be the projected least-squares problem's
01262   // solution.  This won't do any work if the work has already
01263   // been done.
01264   updateProjectedLeastSquaresProblem();
01265   RCP<const MV> V_view = MVT::CloneView (*V_, Range1D(0, m));
01266   const mat_type H_view (Teuchos::View, projectedProblem_->H, m+1, m);
01267   const mat_type y_view (Teuchos::View, projectedProblem_->y, m, 1);
01268   mat_type H_times_y (m+1, 1);
01269   {
01270     const int err = 
01271       H_times_y.multiply (Teuchos::NO_TRANS, Teuchos::NO_TRANS, 
01272         one, H_view, y_view, zero);
01273     TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::logic_error,
01274            "In (F)GMRES, when computing the current native "
01275            "residual vector via the Arnoldi relation, H*y "
01276            "failed due to incompatible dimensions.  "
01277            "This likely indicates a bug in Belos.");
01278   }
01279   // nativeResVec_ := -V * (H * y) + r0 = r0 - V (H y)
01280   // (where we stored r0 in nativeResVec_).
01281   MVT::MvTimesMatAddMv (-one, *V_view, H_times_y, one, *nativeResVec_);
01282       }
01283     // Recport the "native" residual vector's norm, scaled by the
01284     // initial residual norm (if the latter is not zero).
01285     if (verboseDebug)
01286       {
01287   std::vector<magnitude_type> theNorm (1);
01288   MVT::MvNorm (*nativeResVec_, theNorm);
01289   dbg << "---- ||r_0||_2 = " << initialResidualNorm_ << endl;
01290   if (initialResidualNorm_ == STM::zero())
01291     {
01292       dbg << "---- ||Native residual vector||_2 = " << theNorm[0] << endl;
01293     }
01294   else
01295     {
01296       const magnitude_type relResNorm = theNorm[0] / initialResidualNorm_;
01297       dbg << "---- ||Native residual vector||_2 / ||r_0||_2 = " 
01298     << relResNorm << endl;
01299     }
01300   RCP<MV> xUpdate = getCurrentUpdate();
01301   // FIXME (mfh 28 Feb 2011) Is this correct for Flexible GMRES
01302   // too?
01303   RCP<MV> r = MVT::CloneCopy (*initResVec_);
01304   RCP<MV> AxUpdate = MVT::CloneCopy (*initResVec_);
01305   lp_->apply (*xUpdate, *AxUpdate);
01306   // r := initResVec_ - A * xUpdate
01307   MVT::MvAddMv (STS::one(), *initResVec_, -STS::one(), *AxUpdate, *r);
01308   MVT::MvNorm (*r, theNorm);
01309 
01310   dbg << "---- ||b-Ax||_2 = " << theNorm[0] << endl;
01311   if (initialResidualNorm_ != STM::zero())
01312     {
01313       const magnitude_type relResNorm = theNorm[0] / initialResidualNorm_;
01314       dbg << "---- ||b-Ax||_2 / ||r_0|| = " << relResNorm << endl;
01315     }
01316       }
01317     return nativeResVec_;
01318   }
01319 
01320   template<class Scalar, class MV, class OP>
01321   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
01322   GmresBase<Scalar, MV, OP>::currentNativeResidualNorm ()
01323   {
01324     // Update the least-squares problem if necessary.  This computes
01325     // the current native residual norm for "free."
01326     const magnitude_type theNorm = updateProjectedLeastSquaresProblem();
01327     return theNorm;
01328   }
01329 
01330   template<class Scalar, class MV, class OP>
01331   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
01332   GmresBase<Scalar, MV, OP>::updateProjectedLeastSquaresProblem()
01333   {
01334     using std::endl;
01335 
01336     // [startCol, endCol] is a zero-based inclusive index range.
01337     const int startCol = (lastUpdatedCol_ < 0) ? 0 : lastUpdatedCol_ + 1;
01338     // getNumIters() returns the number of _completed_ iterations.
01339     const int endCol = getNumIters() - 1;
01340 
01341     const char prefix[] = "Belos::GmresBase::updateProjectedLeastSquaresProblem: ";
01342     std::ostream& dbg = outMan_->stream(Debug);
01343     dbg << "---- updateProjectedLeastSquaresProblem: ";
01344 
01345     TEUCHOS_TEST_FOR_EXCEPTION(startCol > endCol+1, std::logic_error,
01346            prefix << "Somehow, GmresBase updated the QR "
01347            "factorization of the upper Hessenberg matrix past the "
01348            "last updated column, which was " << lastUpdatedCol_ 
01349            << ". The rightmost column to update is " << endCol 
01350            << ".  Please report this bug to the Belos developers.");
01351     if (endCol < 0) {
01352       dbg << "------ Nothing to do: # iters = 0" << endl;
01353       // Assign here, just in case we haven't done this yet.
01354       nativeResidualNorm_ = initialResidualNorm_; 
01355     } else if (startCol == endCol+1) { // Nothing to do
01356       dbg << "------ Nothing to do: No columns to update (startCol == endCol+1 "
01357   "== " << startCol << ")" << endl;
01358     } else {
01359       dbg << endl;
01360 
01361       details::ProjectedLeastSquaresSolver<Scalar> solver (outMan_->stream(Warnings));
01362       const magnitude_type newAbsoluteResidualNorm = 
01363   solver.updateColumns (*projectedProblem_, startCol, endCol);
01364 
01365       // Scale the absolute residual norm by the initial residual norm.
01366       // Be sure not to divide by zero (though if the initial residual
01367       // norm is zero, GMRES shouldn't progress anyway).
01368       if (initialResidualNorm_ == STM::zero()) {
01369   nativeResidualNorm_ = newAbsoluteResidualNorm;
01370       } else {
01371   nativeResidualNorm_ = newAbsoluteResidualNorm / initialResidualNorm_;
01372       }
01373 
01374       // Remember the rightmost updated column.
01375       lastUpdatedCol_ = endCol;
01376     }
01377 
01378     // Return the "native" residual norm (which is the absolute
01379     // residual error from solving the projected least-squares
01380     // problem).
01381     dbg << "------ Native residual norm: " << nativeResidualNorm_ << endl;
01382     return nativeResidualNorm_;
01383   }
01384 
01385 
01386   template<class Scalar, class MV, class OP>
01387 #if 0
01388   std::pair< Teuchos::RCP<MV>, typename Teuchos::ScalarTraits<Scalar>::magnitudeType >
01389 #else // not 0
01390   Teuchos::RCP<MV>
01391 #endif // 0
01392   GmresBase<Scalar, MV, OP>::getCurrentUpdate ()
01393   {
01394     using Teuchos::is_null;
01395     using Teuchos::Range1D;
01396     using Teuchos::RCP;
01397     const Scalar one = STS::one();
01398     const Scalar zero = STS::zero();
01399 
01400     if (is_null (xUpdate_))
01401       {
01402   // xUpdate_ comes from the Z_ space and not the V_ space in
01403   // the flexible variant of GMRES.  This makes a difference
01404   // only if there are multiple vector spaces involved, that
01405   // is if the preconditioner maps V-space vectors to Z-space
01406   // vectors, and the matrix maps Z-space vectors to V-space
01407   // vectors (where V and Z are different spaces).  This is
01408   // legitimate as long as the initial guess for the solution
01409   // (x0) comes from the Z space.
01410   if (flexible_)
01411     xUpdate_ = MVT::Clone (*Z_, 1);
01412   else
01413     xUpdate_ = MVT::Clone (*V_, 1);
01414       }
01415     if (getNumIters() == 0)
01416       {
01417   MVT::MvInit (*xUpdate_, zero);
01418   return xUpdate_;
01419       }
01420 
01421     // When one iteration of GMRES is completed, the upper
01422     // Hessenberg matrix H is 2 by 1.  Thus, we want to update
01423     // columns up to and including curNumIters_ - 1: actually,
01424     // columns [lastUpdatedCol_+1, curNumIters_-1] (inclusive).  The
01425     // solution update gives you the current "native" residual norm
01426     // for free.  We could save and return that as well...
01427     //const magnitude_type nativeResNorm = updateProjectedLeastSquaresProblem();
01428     (void) updateProjectedLeastSquaresProblem();
01429 
01430     // Current iteration count does not include the initial basis vector.
01431     const int m = getNumIters();
01432     // Basis for solution update coefficients; Flexible GMRES uses the
01433     // preconditioned basis here.
01434     RCP<const MV> Z_view = flexible_ ? 
01435       MVT::CloneView (*Z_, Range1D(0, m-1)) : 
01436       MVT::CloneView (*V_, Range1D(0, m-1));
01437     const mat_type y_view (Teuchos::View, projectedProblem_->y, m, 1);
01438     // xUpdate_ := Z_view * y_view
01439     MVT::MvTimesMatAddMv (one, *Z_view, y_view, zero, *xUpdate_);
01440     return xUpdate_;
01441   }
01442 
01443   template<class Scalar, class MV, class OP>
01444   void
01445   GmresBase<Scalar, MV, OP>::
01446   setInitialResidual(const int maxIterCount,
01447          const bool reallocateBasis)
01448   {
01449     using Teuchos::Range1D;
01450     using Teuchos::RCP;
01451     using Teuchos::rcp;
01452     using std::endl;
01453 
01454     // Message fragment for error messages.
01455     const char prefix[] = "Belos::GmresBase::setInitialResidual: ";
01456     std::ostream& dbg = outMan_->stream(Debug);
01457 
01458     TEUCHOS_TEST_FOR_EXCEPTION(maxIterCount < 0, std::invalid_argument,
01459            prefix << "maxIterCount = " << maxIterCount 
01460            << " < 0.");
01461     // Do we have a left preconditioner?
01462     const bool haveLeftPrec = ! lp_->getLeftPrec().is_null();
01463     if (haveLeftPrec)
01464       dbg << "-- Have a left preconditioner" << endl;
01465     else
01466       dbg << "-- No left preconditioner" << endl;
01467 
01468     // mfh 22 Feb 2011: LinearProblem's getInitResVec() and
01469     // getInitPrecResVec() methods both return the residual vector for
01470     // the whole linear problem (X_ and B_), not for the "current"
01471     // linear problem (curX_ and curB_, returned by getCurrLHSVec()
01472     // resp. getCurrRHSVec()).  In order to get the residual vector(s)
01473     // for the _current_ linear problem, we have to use getLSIndex()
01474     // to extract the column ind(ex/ices) and take the corresponding
01475     // column of getInit(Prec)ResVec().  It should be exactly one
01476     // column, since this GMRES solver only knows how to solve for one
01477     // right-hand side at a time.
01478     RCP<const MV> R_full = 
01479       haveLeftPrec ? lp_->getInitPrecResVec() : lp_->getInitResVec();
01480     std::vector<int> inputIndices = lp_->getLSIndex();
01481     TEUCHOS_TEST_FOR_EXCEPTION(inputIndices.size() == 0, 
01482            std::invalid_argument,
01483            "The LinearProblem claims that there are zero linear "
01484            "systems to solve: getLSIndex() returns an index "
01485            "vector of length zero.");
01486     dbg << "-- Input indices: [";
01487     for (std::vector<int>::size_type k = 0; k < inputIndices.size(); ++k)
01488       {
01489   dbg << inputIndices[k];
01490   if (k < inputIndices.size() - 1)
01491     dbg << ", ";
01492       }
01493     dbg << "]" << endl;
01494 
01495     // Some of the indices returned by getLSIndex() might be -1.
01496     // These are for column(s) of getCurr{L,R}HSVec() that don't
01497     // correspond to any actual linear system(s) that we want to
01498     // solve.  They may be filled in by block iterative methods in
01499     // order to make the blocks full rank or otherwise improve
01500     // convergence.  We shouldn't have any of those columns here, but
01501     // we check just to make sure.
01502     std::vector<int> outputIndices;
01503     std::remove_copy_if (inputIndices.begin(), inputIndices.end(), 
01504        std::back_inserter (outputIndices), 
01505        std::bind2nd (std::equal_to<int>(), -1));
01506     dbg << "-- Output indices: [";
01507     for (std::vector<int>::size_type k = 0; k < outputIndices.size(); ++k)
01508       {
01509   dbg << outputIndices[k];
01510   if (k < outputIndices.size() - 1)
01511     dbg << ", ";
01512       }
01513     dbg << "]" << endl;
01514 
01515     // outputIndices should have exactly one entry, corresponding to
01516     // the one linear system to solve.  Our GMRES implementation can
01517     // only solve for one right-hand side at a time.
01518     if (outputIndices.size() != 1 || outputIndices[0] == -1)
01519       {
01520   std::string inputIndicesAsString;
01521   {
01522     std::ostringstream os;
01523     os << "[";
01524     std::copy (inputIndices.begin(), inputIndices.end(),
01525          std::ostream_iterator<int>(os, " "));
01526     os << "]";
01527     inputIndicesAsString = os.str();
01528   }
01529   std::string outputIndicesAsString;
01530   {
01531     std::ostringstream os;
01532     os << "[";
01533     std::copy (outputIndices.begin(), outputIndices.end(),
01534          std::ostream_iterator<int>(os, " "));
01535     os << "]";
01536     outputIndicesAsString = os.str();
01537   }
01538   std::ostringstream os;
01539   os << prefix << "The LinearProblem instance's getLSIndex() method "
01540     "returns indices " << inputIndicesAsString << ", of which the "
01541     "following are not -1: " << outputIndicesAsString << ".";
01542   TEUCHOS_TEST_FOR_EXCEPTION(outputIndices.size() != 1,
01543          std::invalid_argument, 
01544          os.str() << "  The latter list should have length "
01545          "exactly 1, since our GMRES implementation only "
01546          "knows how to solve for one right-hand side at a "
01547          "time.");
01548   TEUCHOS_TEST_FOR_EXCEPTION(outputIndices[0] == -1,
01549          std::invalid_argument,
01550          os.str() << "  the latter list contains no "
01551          "nonnegative entries, meaning that there is no "
01552          "valid linear system to solve.");
01553       }
01554     // View of the "current" linear system's residual vector.
01555     RCP<const MV> r0 = MVT::CloneView (*R_full, outputIndices);
01556     // Sanity check that the residual vector has exactly 1 column.
01557     // We've already checked this a different way above.
01558     const int numResidVecs = MVT::GetNumberVecs (*r0);
01559     TEUCHOS_TEST_FOR_EXCEPTION(numResidVecs != 1, std::logic_error,
01560            prefix << "Residual vector has " << numResidVecs 
01561            << " columns, but should only have one.  "
01562            "Should never get here.");
01563     // Save a deep copy of the initial residual vector.
01564     initResVec_ = MVT::CloneCopy (*r0);
01565 
01566     // Allocate space for the V_ basis vectors if necessary.
01567     //
01568     // If left preconditioning is used, then the V_ vectors are in the
01569     // same vector space as the (left-)preconditioned initial residual
01570     // vector.  Otherwise they are in the same space as the
01571     // unpreconditioned initial residual vector.  That's r0 in either
01572     // case.
01573     //
01574     // FIXME (mfh 22 Feb 2011) We should really check to make sure
01575     // that the columns of V_ (if previously allocated) are in the
01576     // same vector space as r0.  Belos::MultiVecTraits doesn't give us
01577     // a way to do that.  Just checking if the number of rows is the
01578     // same is not enough.  The third Boolean argument
01579     // (reallocateBasis) lets us add in the check externally if
01580     // MultiVecTraits gains that ability sometime.
01581     if (V_.is_null() || 
01582   MVT::GetNumberVecs(*V_) != maxIterCount+1 || 
01583   reallocateBasis)
01584       {
01585   dbg << "-- Reallocating V_ basis with " << (maxIterCount+1) 
01586       << " columns" << endl;
01587   V_ = MVT::Clone (*initResVec_, maxIterCount+1);
01588       }
01589     else
01590       dbg << "-- V_ basis has alread been allocated with " 
01591     << (maxIterCount+1) << " columns." << endl;
01592 
01593     // Initial residual norm is computed with respect to the inner
01594     // product defined by the OrthoManager.  Since the basis
01595     // vectors are orthonormal (resp. unitary) with respect to that
01596     // inner product, this ensures that the least-squares minimization
01597     // is also computed with respect to that inner product.
01598     std::vector<magnitude_type> result (1);
01599     ortho_->norm (*r0, result);
01600     initialResidualNorm_ = result[0];
01601     dbg << "-- Initial residual norm (as computed by the OrthoManager): " 
01602   << initialResidualNorm_ << endl;
01603 
01604     if (projectedProblem_.is_null()) {
01605       projectedProblem_ = 
01606   rcp (new details::ProjectedLeastSquaresProblem<Scalar> (maxIterCount));
01607     }
01608     projectedProblem_->reallocateAndReset (initialResidualNorm_, maxIterCount);
01609 
01610     // Copy the initial residual vector into the first column of V_,
01611     // and scale the vector by its norm.
01612     RCP<MV> v1 = MVT::CloneViewNonConst (*V_, Range1D(0,0));
01613     MVT::Assign (*r0, *v1);
01614     MVT::MvScale (*v1, STS::one() / initialResidualNorm_);
01615   }
01616 
01617   template<class Scalar, class MV, class OP>
01618   void 
01619   GmresBase<Scalar, MV, OP>::backOut (const int numIters)
01620   {
01621     TEUCHOS_TEST_FOR_EXCEPTION(numIters < 0, std::invalid_argument,
01622            "The GMRES iteration count cannot be less than "
01623            "zero, but you specified numIters = " << numIters);
01624     TEUCHOS_TEST_FOR_EXCEPTION(numIters > getNumIters(), std::invalid_argument,
01625            "The GMRES iteration count cannot be reset to more than "
01626            "its current iteration count " << getNumIters()
01627            << ", but you specified numIters = " << numIters << ".");
01628 
01629     // Reset the projected least-squares problem, without reallocating
01630     // its data.  We will "replay" the first numIters steps of the QR
01631     // factorization below.
01632     projectedProblem_->reset (initialResidualNorm_);
01633 
01634     // Replay the first numIters-1 Givens rotations.
01635     lastUpdatedCol_ = -1;
01636     curNumIters_ = numIters;
01637     (void) updateProjectedLeastSquaresProblem ();
01638   }
01639 
01640   template<class Scalar, class MV, class OP>
01641   void 
01642   GmresBase<Scalar, MV, OP>::
01643   restart (const int maxIterCount)
01644   {
01645     using Teuchos::RCP;
01646     using Teuchos::null;
01647 
01648     // This would be the place where subclasses may implement things
01649     // like recycling basis vectors for the next restart cycle.  Note
01650     // that the linear problem hasn't been updated yet; this happens
01651     // below.
01652     preRestartHook (maxIterCount);
01653 
01654     // Make sure that the LinearProblem object gets the current
01655     // approximate solution (which becomes the initial guess for the
01656     // upcoming new restart cycle), so that the "exact" residuals are
01657     // correct.
01658     (void) getCurrentUpdate (); // results stored in xUpdate_
01659     updateSolution ();
01660 
01661     // Check whether there is a right preconditioner, since (at least
01662     // in theory) the caller could have added a right preconditioner
01663     // after a restart cycle.  (Is that legitimate?  Adding or
01664     // changing a right preconditioner between restarts wouldn't
01665     // change the initial residual vector, and therefore wouldn't
01666     // change the V_ basis.  The V_ basis is all that matters for the
01667     // restart.)
01668     const bool needToAllocateZ = flexible_ && ! lp_->getRightPrec().is_null();
01669 
01670     // This (re)allocates the V_ basis if necessary.  It also
01671     // (re)allocates and resets the projected least-squares problem.
01672     setInitialResidual (maxIterCount);
01673 
01674     // If the new max iteration count has changed, reallocate space
01675     // for the Z_ basis vectors.  setInitialResidual() already does
01676     // that for the V_ basis.
01677     if (maxNumIters_ != maxIterCount)
01678       {
01679   maxNumIters_ = maxIterCount;
01680   // Initial guess.
01681   RCP<MV> x0 = lp_->getCurrLHSVec();
01682   Z_ = needToAllocateZ ? MVT::Clone (*x0, maxIterCount+1) : null;
01683       }
01684     lastUpdatedCol_ = -1; // column updates start with zero
01685     curNumIters_ = 0;
01686     flexible_ = needToAllocateZ;
01687 
01688     // This would be the place where subclasses may implement things
01689     // like orthogonalizing the first basis vector (after restart)
01690     // against the recycled subspace.
01691     postRestartHook();
01692   }
01693 
01694 } // namespace Belos
01695 
01696 #endif // __Belos_GmresBase_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines