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