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 
00053 #include <Teuchos_BLAS.hpp>
00054 #include <Teuchos_LAPACK.hpp>
00055 #include <Teuchos_ScalarTraits.hpp>
00056 #include <Teuchos_SerialDenseVector.hpp>
00057 #include <Teuchos_TypeNameTraits.hpp>
00058 
00059 #include <algorithm>
00060 #include <iterator>
00061 #include <utility> // std::pair
00062 
00063 // Anonymous namespace is a more portable way than "static" to define
00064 // functions that won't escape the current file's scope.
00065 namespace {
00066   template<class Ordinal, class Scalar>
00067   std::ostream& 
00068   operator<< (std::ostream& out,
00069         const Teuchos::SerialDenseMatrix<Ordinal, Scalar>& A)
00070   {
00071     const Ordinal m = A.numRows();
00072     const Ordinal n = A.numCols();
00073     using std::endl;
00074 
00075     // Summarize the matrix if it's too big to print comfortably on
00076     // the terminal.
00077     //
00078     // FIXME (mfh 25 Feb 2011) Obviously this is no substitute for a
00079     // "curses"-type terminal interface that can query the terminal
00080     // for its dimensions.
00081     if (m > 10 || n > 10)
00082       {
00083   out << "[<" << m << " x " << n << " SerialDenseMatrix<int," 
00084       << Teuchos::TypeNameTraits<Scalar>::name() << ">]";
00085   return out;
00086       }
00087     else if (m == 0 || n == 0)
00088       { // Empty matrix: either zero rows or zeo columns
00089   out << "[]";
00090   return out;
00091       }
00092     out << "[";
00093     if (m > 1)
00094       out << endl;
00095     for (Ordinal i = 0; i < m; ++i)
00096       { // Indent if printing on more than one row.
00097   if (m > 1)
00098     out << "  ";
00099   // Print row i of the matrix: comma-delimited, like Matlab.
00100   for (Ordinal j = 0; j < n; ++j)
00101     {
00102       out << A(i,j);
00103       if (j < n-1)
00104         out << ", ";
00105     }
00106   // End of the current row, which we know is not empty.
00107   if (i < m-1 && m > 1)
00108     out << ";" << endl;
00109       }
00110     if (m > 1)
00111       out << endl;
00112     out << "]";
00113     return out;
00114   }
00115 
00116 } // (anonymous namespace)
00117 
00118 namespace Belos {
00119 
00152   class GmresCantExtendBasis : public BelosError {
00153   public:
00154     GmresCantExtendBasis(const std::string& what_arg) : BelosError(what_arg) {}
00155   };
00156 
00190   class GmresRejectsCandidateBasis : public BelosError {
00191   public:
00192     GmresRejectsCandidateBasis(const std::string& what_arg) : BelosError(what_arg) {}
00193   };
00194 
00195 
00218   template<class Scalar, class MV, class OP>
00219   class GmresBase {
00220   public:
00221     typedef Scalar scalar_type;
00222     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00223     typedef MV multivector_type;
00224     typedef OP operator_type;
00225 
00226   private:
00227     typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00228     typedef Teuchos::SerialDenseVector<int,Scalar> vec_type;
00229     typedef Teuchos::BLAS<int, scalar_type> blas_type;
00230     typedef Teuchos::LAPACK<int, scalar_type> lapack_type;
00231     typedef MultiVecTraits<scalar_type, MV> MVT;
00232     typedef Teuchos::ScalarTraits<Scalar> STS;
00233     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00234 
00235   public:
00272     GmresBase (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
00273          const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00274          const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00275          const int maxIterCount,
00276          const bool flexible);
00277 
00279     bool canAdvance () const { return canExtendBasis(); }
00280 
00294     void advance();
00295 
00301     virtual void 
00302     acceptCandidateBasis (const int newNumVectors = 1) 
00303     {
00304       TEST_FOR_EXCEPTION(newNumVectors <= 0, std::logic_error,
00305        "The number of candidate basis vectors to accept, " 
00306        << newNumVectors << ", is nonpositive, which is not "
00307        "allowed.  This is an std:logic_error because it "
00308        "relates to the implementation of the GmresBase "
00309        "subclass.");
00310       curNumIters_ += newNumVectors;
00311     }
00312 
00317     virtual void 
00318     rejectCandidateBasis (const int newNumVectors = 1) 
00319     {
00320       std::ostringstream os;
00321       os << (flexible_ ? "Flexible GMRES" : "GMRES") 
00322    << " rejects the computed candidate basis vector"
00323    << (newNumVectors != 1 ? "s" : "");
00324       throw GmresRejectsCandidateBasis (os.str());
00325     }
00326 
00330     int getNumIters() const { return curNumIters_; }
00331 
00335     int maxNumIters() const { return maxNumIters_; }
00336 
00360     void restart (const int maxIterCount);
00361 
00363     void restart () { restart (maxNumIters_); }
00364 
00373     void backOut (const int numIters);
00374 
00379     Teuchos::RCP<MV> getCurrentUpdate ();
00380 
00388     Teuchos::RCP<MV> currentNativeResidualVector ();
00389 
00396     typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00397     currentNativeResidualNorm ();
00398 
00416     typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00417     arnoldiRelationError () const;
00418 
00425     void
00426     updateSolution ()
00427     {
00428       (void) lp_->updateSolution (xUpdate_, true, STS::one());
00429     }
00430 
00431   protected:
00432 
00437     //{@
00438 
00446     virtual bool canExtendBasis() const = 0;
00447 
00473     virtual void 
00474     extendBasis (Teuchos::RCP<MV>& V_cur, 
00475      Teuchos::RCP<MV>& Z_cur) = 0;
00476 
00522     virtual void
00523     orthogonalize (const Teuchos::RCP<MV>& V_cur,
00524        const Teuchos::RCP<const MV>& V_prv,
00525        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_V,
00526        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_V,
00527        const Teuchos::RCP<MV>& Z_cur,
00528        const Teuchos::RCP<const MV>& Z_prv,
00529        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_Z,
00530        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_Z) = 0;
00531     
00559     virtual void 
00560     updateUpperHessenbergMatrix (const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_V,
00561          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_V,
00562          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& C_Z,
00563          const Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> >& B_Z) = 0;
00564 
00566     virtual bool acceptedCandidateBasis() const = 0;
00567 
00580     virtual bool shouldRetryAdvance() { return false; }
00581 
00583 
00613 
00614 
00626     virtual void preRestartHook (const int maxIterCount) {}
00627 
00638     virtual void postRestartHook() {}
00639 
00649     virtual void 
00650     preOrthogonalizeHook (const Teuchos::RCP<MV>& V_cur, 
00651         const Teuchos::RCP<MV>& Z_cur) {}
00652 
00662     void
00663     postOrthogonalizeHook (const Teuchos::RCP<MV>& V_cur, 
00664          const Teuchos::RCP<mat_type>& C_V,
00665          const Teuchos::RCP<mat_type>& B_V,
00666          const Teuchos::RCP<MV>& Z_cur,
00667          const Teuchos::RCP<mat_type>& C_Z,
00668          const Teuchos::RCP<mat_type>& B_Z) {}
00670 
00671   private:
00672 
00674 
00675 
00696     void 
00697     setInitialResidual(const int maxIterCount,
00698            const bool reallocateBasis = false);
00699 
00713     typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00714     updateProjectedLeastSquaresProblem ();
00715 
00728     void 
00729     previousVectors (Teuchos::RCP<const MV>& V_prv,
00730          Teuchos::RCP<const MV>& Z_prv) const;
00731 
00745     Teuchos::RCP<const MV> initResVec() const {
00746       return initResVec_;
00747     }
00748 
00750 
00751   protected:
00753 
00754 
00756     Teuchos::RCP<LinearProblem<Scalar, MV, OP> > lp_;
00757 
00759     Teuchos::RCP<const OrthoManager<Scalar, MV> > ortho_;
00760 
00762     Teuchos::RCP<OutputManager<Scalar> > outMan_;
00763 
00785     Teuchos::RCP<MV> nativeResVec_;
00786 
00792     Teuchos::RCP<MV> initResVec_;
00793     
00801     Teuchos::RCP<MV> xUpdate_;
00802 
00809     Teuchos::RCP<MV> V_;
00810 
00816     Teuchos::RCP<MV> Z_;
00817 
00828     Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > H_;
00829 
00831     Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > R_;
00832 
00834     Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > y_;
00835 
00843     Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > z_;
00844 
00845     Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> theCosines_;
00846     Teuchos::Array<Scalar> theSines_;
00847 
00856     typename Teuchos::ScalarTraits<Scalar>::magnitudeType initialResidualNorm_;
00857 
00860     // This is the absolute residual error from solving the projected
00861     // least-squares problem, if we've finished at least one iteration
00862     // of (F)GMRES.  Otherwise, it's the same as the initial residual
00863     // norm.
00864     typename Teuchos::ScalarTraits<Scalar>::magnitudeType nativeResidualNorm_;
00865 
00870     int lastUpdatedCol_;
00871 
00877     int curNumIters_;
00878 
00887     int maxNumIters_;
00888 
00919     bool flexible_;
00920 
00922   };
00923 
00926 
00927   template<class Scalar, class MV, class OP>
00928   void
00929   GmresBase<Scalar, MV, OP>::
00930   previousVectors (Teuchos::RCP<const MV>& V_prv,
00931        Teuchos::RCP<const MV>& Z_prv) const
00932   {
00933     using Teuchos::Range1D;
00934     const bool verboseDebug = false;
00935 
00936     // The number of iterations (returned by getNumIters()) does not
00937     // include the initial basis vector, which is the first column of
00938     // V_.  Hence, the range of previous basis vectors (in the V
00939     // basis) is [0,k] inclusive.
00940     const int k = getNumIters(); 
00941 
00942     if (verboseDebug)
00943       {
00944   using std::endl;
00945   std::ostream& dbg = outMan_->stream(Debug);
00946   dbg << "---- previousVectors: V_prv = V[0:" << k << "]";
00947   if (flexible_)
00948     dbg << ", Z_prv = Z[0:" << k << "]" << endl;
00949   else 
00950     dbg << endl;
00951       }
00952     
00953     V_prv = MVT::CloneView (*V_, Range1D(0,k));
00954     if (flexible_)
00955       // FIXME (mfh 24 Feb 2011) Is this right?  When should we apply
00956       // the right preconditioner in Flexible GMRES?
00957       Z_prv = MVT::CloneView (*Z_, Range1D(0,k));
00958     else
00959       Z_prv = Teuchos::null;
00960   }
00961 
00962   template<class Scalar, class MV, class OP>
00963   void
00964   GmresBase<Scalar, MV, OP>::advance () 
00965   {
00966     using Teuchos::is_null;
00967     using Teuchos::null;
00968     using Teuchos::Range1D;
00969     using Teuchos::RCP;
00970     using Teuchos::rcp;
00971     using Teuchos::tuple;
00972     using std::endl;
00973     typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00974     const bool verboseDebug = false;
00975 
00976     std::ostream& dbg = outMan_->stream(Debug);
00977     dbg << "Belos::GmresBase::advance():" << endl
00978   << "--" << endl
00979   << "-- Iteration " << getNumIters() << endl
00980   << "--" << endl;
00981 
00982     TEST_FOR_EXCEPTION( !canExtendBasis(), GmresCantExtendBasis,
00983       "GMRES (iteration " << getNumIters() << ") cannot "
00984       "add any more basis vectors." );
00985 
00986     RCP<const MV> V_prv, Z_prv;
00987     RCP<MV> V_cur, Z_cur;
00988     // Fetch the previous basis vector(s), against which to project.
00989     previousVectors (V_prv, Z_prv);
00990     if (verboseDebug)
00991     dbg << "---- Number of previous vectors: " 
00992       << MVT::GetNumberVecs(*V_prv) << endl;
00993     // Compute new basis vector(s).
00994     extendBasis (V_cur, Z_cur);
00995 
00996     // Do whatever the subclass wants to do before orthogonalizing the
00997     // new basis vectors.  This might include projecting them against
00998     // a fixed basis (e.g., for Krylov subspace recycling).
00999     preOrthogonalizeHook (V_cur, Z_cur);
01000 
01001     // Flexible CA-GMRES must orthogonalize the Z basis along with the
01002     // V basis, and it must keep both sets of orthogonalization
01003     // coefficients.  Flexible standard GMRES only needs to
01004     // orthogonalize the new V basis vector; it doesn't need to
01005     // orthogonalize the Z basis vector.  This is why we let the
01006     // subclass implement this functionality, rather than calling the
01007     // OrthoManager's projectAndNormalize() method directly; the
01008     // subclass has to decide whether it needs to orthogonalize Z_cur
01009     // as well as V_cur.
01010     //
01011     // Subclasses may want to save the rank (if applicable) from
01012     // ortho_->projectAndNormalize() somewhere.  In the case of
01013     // Flexible CA-GMRES, there are _two_ ranks -- one for the V
01014     // basis, and the other for the Z basis.
01015     //
01016     // C_V and B_V are the "C" (projection) resp. "B" (normalization)
01017     // block coefficients from ortho_->projectAndNormalize() on the
01018     // new V basis vectors.  If the new Z basis vectors need to be
01019     // orthogonalized as well, then C_Z and B_Z are the "C" resp. "B"
01020     // block coefficients from ortho_->projectAndNormalize() on the
01021     // new Z basis vectors.  Subclasses' implementations of
01022     // orthogonalize() are responsible for allocating space for
01023     // C_V and B_V, and C_Z and B_Z if necessary (otherwise they are
01024     // set to Teuchos::null, as the RCPs are passed by reference).
01025     RCP<mat_type> C_V, B_V, C_Z, B_Z;
01026     orthogonalize (V_cur, V_prv, C_V, B_V, Z_cur, Z_prv, C_Z, B_Z);
01027     if (false && verboseDebug)
01028       {
01029   dbg << "---- Projection coefficient(s): ";
01030   dbg << *C_V << endl;
01031   dbg << "---- Normalization coefficient(s): ";
01032   dbg << *B_V << endl;
01033       }
01034     // Do whatever the subclass wants to do after projecting the new
01035     // basis vectors against the previous ones and normalizing them.
01036     postOrthogonalizeHook (V_cur, C_V, B_V, Z_cur, C_Z, B_Z);
01037 
01038     // Implemented by subclasses.  C_Z and B_Z need not be used,
01039     // regardless of whether we are computing Flexible GMRES.  (Only
01040     // Flexible CA-GMRES needs them, as far as we know.)  If they are
01041     // not used, orthogonalize() may set them to null (the RCPs are
01042     // passed by reference).  Some subclasses may choose to update the
01043     // upper Hessenberg matrix "in place" in their implementation of
01044     // the orthogonalize() method (see above).
01045     updateUpperHessenbergMatrix (C_V, B_V, C_Z, B_Z);
01046 
01047     // The subclass decides whether or not to accept the candidate
01048     // basis.
01049     //
01050     // We update the upper Hessenberg matrix before deciding
01051     // acceptance, because the acceptance critera might depend on the
01052     // entire upper Hessenberg matrix (e.g., (an estimate of) its
01053     // condition number), not just the condition number or rank of the
01054     // candidate basis vector(s).
01055     //
01056     // Subclasses' implementations of orthogonalize() and/or
01057     // updateUpperHessenbergMatrix() might wish to set state to
01058     // indicate acceptance or rejection of the candidate basis
01059     // vector(s).  They may also call advance() recursively (e.g.,
01060     // setting the "last rejected rank" so that the next call of
01061     // advance() computes fewer candidate basis vectors).
01062     if (! acceptedCandidateBasis())
01063       rejectCandidateBasis (MVT::GetNumberVecs(*V_cur));
01064     else 
01065       // This increments the current iteration count by the number of
01066       // accepted candidate basis vectors.
01067       acceptCandidateBasis (MVT::GetNumberVecs(*V_cur));
01068 
01069     // We don't update the projected least-squares problem here.
01070     // This is done lazily, whenever the "native" residual norm
01071     // or the current solution update are requested.
01072   }
01073 
01074 
01075   template<class Scalar, class MV, class OP>
01076   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
01077   GmresBase<Scalar, MV, OP>::arnoldiRelationError () const
01078   {
01079     using Teuchos::Range1D;
01080     using Teuchos::RCP;
01081     const Scalar one = STS::one();
01082     const int m = getNumIters();
01083 
01084     if (m == 0) // No iterations completed means no error (yet)
01085       return magnitude_type(0);
01086     RCP<const MV> V_mp1 = MVT::CloneView (*V_, Range1D(0,m));
01087     const mat_type H_m (Teuchos::View, *H_, m+1, m);
01088     std::vector<magnitude_type> norms (m);
01089     if (flexible_)
01090       {
01091   RCP<const MV> Z_view = MVT::CloneView (*Z_, Range1D(0,m-1));
01092   // A*Z is in the same space as V, so we clone from V
01093   RCP<MV> AZ = MVT::Clone (*V_, m);
01094   // Compute AZ := A*Z_m
01095   lp_->applyOp (Z_view, AZ);
01096   // Compute AZ := AZ - V_{m+1} \underline{H}_m
01097   MVT::MvTimesMatAddMv (-one, V_mp1, H_m, one, AZ);
01098   MVT::MvNorm (AZ, norms);
01099       }
01100     else
01101       {
01102   RCP<const MV> V_m = MVT::CloneView (*V_, Range1D(0,m-1));
01103   RCP<MV> AV = MVT::Clone (*V_, m);
01104   // Compute AV := A*V_m.  Here, "A" means the preconditioned
01105   // operator.
01106   lp_->apply (V_m, AV);
01107   // Compute AV := A V_m - V_{m+1} \underline{H}_m
01108   MVT::MvTimesMatAddMv (-one, V_mp1, H_m, one, AV);
01109   MVT::MvNorm (AV, norms);
01110       }
01111     // Compute and return the Frobenius norm of the above (either AZ
01112     // or AV, depending on whether we are performing Flexible GMRES).
01113     //
01114     // FIXME (mfh 16 Feb 2011) This computation might overflow for
01115     // unjustifiable reasons.  I should really include intermediate
01116     // scaling.
01117     magnitude_type result (0);
01118     for (int k = 0; k < m; ++k)
01119       result += norms[k];
01120     return STM::squareroot(result);
01121   }
01122 
01123 
01124   template<class Scalar, class MV, class OP>
01125   GmresBase<Scalar, MV, OP>::
01126   GmresBase (const Teuchos::RCP<LinearProblem<Scalar, MV, OP> >& problem,
01127        const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
01128        const Teuchos::RCP<OutputManager<Scalar> >& outMan,
01129        const int maxIterCount,
01130        const bool flexible) :
01131     lp_ (problem), 
01132     ortho_ (ortho),
01133     outMan_ (outMan),
01134     lastUpdatedCol_ (-1), // column updates start with zero
01135     curNumIters_ (0),
01136     maxNumIters_ (maxIterCount),
01137     flexible_ (flexible && ! lp_.is_null() && ! lp_->getRightPrec().is_null())
01138   {
01139     using Teuchos::is_null;
01140     using Teuchos::null;
01141     using Teuchos::RCP;
01142     using Teuchos::rcp;
01143     using std::endl;
01144 
01145     std::ostream& dbg = outMan_->stream(Debug);
01146     dbg << "Belos::GmresBase:" << endl;
01147 
01148     // Fragments of error messages for use in sanity checks.
01149     const char prefix[] = "Belos::GmresBase constructor: ";
01150     const char postfix[] = "  Please call setProblem() with non-null data "
01151       "before passing the LinearProblem to the GmresBase constructor.";
01152     const char postfix2[] = "  That may mean that the LinearProblem's "
01153       "setLSIndex() method has not yet been called, even though setProblem() "
01154       "has been called.  After calling setProblem(), you should call "
01155       "setLSIndex() with the ind{ex,ices} of the right-hand side(s) to solve.";
01156 
01157     // Sanity checks on the linear problem.
01158     //
01159     // First, make sure A, X, and B are all non-null.
01160     TEST_FOR_EXCEPTION(lp_.is_null(), std::invalid_argument,
01161            prefix << "The given LinearProblem is null.");
01162     TEST_FOR_EXCEPTION(! lp_->isProblemSet(), std::invalid_argument,
01163                prefix << "The given LinearProblem has not yet been set."
01164            << postfix);
01165     TEST_FOR_EXCEPTION(lp_->getLHS().is_null(), std::invalid_argument,
01166                prefix << "The given LinearProblem has null initial guess"
01167            " (getLHS()) for all right-hand sides." << postfix);
01168     TEST_FOR_EXCEPTION(lp_->getRHS().is_null(), std::invalid_argument,
01169                prefix << "The given LinearProblem's right-hand side(s) "
01170            "are all null (getRHS().is_null() == true)." << postfix);
01171     TEST_FOR_EXCEPTION(lp_->getOperator().is_null(), std::invalid_argument,
01172                prefix << "The given LinearProblem's operator (the "
01173            "matrix A) is null." << postfix);
01174     // Next, make sure that setLSIndex() has been called on the linear
01175     // problem instance, so that the "current" right-hand side and
01176     // "current" initial guess have been set.
01177     TEST_FOR_EXCEPTION(lp_->getCurrLHSVec().is_null(), 
01178            std::invalid_argument,
01179                prefix << "Although the given LinearProblem has non-null "
01180            "initial guess (getLHS()) for all right-hand sides, the "
01181            "current initial guess (getCurrLHSVec()) is null."
01182            << postfix2);
01183     TEST_FOR_EXCEPTION(lp_->getCurrRHSVec().is_null(), 
01184            std::invalid_argument,
01185                prefix << "Although the given LinearProblem has non-null "
01186            "initial guess (getLHS()) for all right-hand sides, the "
01187            "current right-hand side to solve (getCurrRHSVec()) is "
01188            "null." << postfix2);
01189     // Get the initial guess x0, and allocate the solution update
01190     // (xUpdate_) from the same vector space as x0.
01191     RCP<const MV> x0 = lp_->getCurrLHSVec();
01192     {
01193       const int numLHS = MVT::GetNumberVecs (*x0);
01194       TEST_FOR_EXCEPTION(numLHS != 1, std::invalid_argument,
01195        "Our GMRES implementation only works for single-"
01196        "vector problems, but the supplied initial guess has "
01197        << numLHS << " columns.");
01198       const int numRHS = MVT::GetNumberVecs (*(lp_->getCurrRHSVec()));
01199       TEST_FOR_EXCEPTION(numRHS != 1, std::invalid_argument,
01200        "Our GMRES implementation only works for single-"
01201        "vector problems, but the current right-hand side has "
01202        << numRHS << " columns.");
01203     }
01204     dbg << "-- Initializing xUpdate_ and filling with zeros" << endl;
01205     xUpdate_ = MVT::Clone (*x0, 1);
01206     MVT::MvInit (*xUpdate_, STS::zero());
01207 
01208     // Get the (left preconditioned, if applicable) residual vector,
01209     // and make a deep copy of it in initResVec_.  Allocate the V_
01210     // basis vectors, compute the initial residual norm, and compute
01211     // the first basis vector (first column of V_).
01212     setInitialResidual (maxIterCount);
01213 
01214     // The Z_ vectors, if we need them, are always in the same vector
01215     // space as the initial solution guess (x0).  Even if the caller
01216     // asks for the flexible option, we don't allocate space for Z_
01217     // unless the caller specifies a right preconditioner.
01218     Z_ = (flexible && ! lp_->getRightPrec().is_null()) ? 
01219       MVT::Clone(*x0, maxIterCount+1) : null;
01220 
01221     // These (small dense) matrices and vectors encode the projected
01222     // least-squares problem.  z_ was already allocated and
01223     // initialized in setInitialResidual() above.
01224     H_ = rcp (new mat_type (maxIterCount+1, maxIterCount));
01225     R_ = rcp (new mat_type (maxIterCount+1, maxIterCount));
01226     // If we choose to use _GELS to solve the projected least-squares
01227     // problem, we need to leave an extra entry in y so that we can
01228     // copy the right-hand side z into it.  This is because _GELS
01229     // overwrites the input right-hand side with the output
01230     // least-squares solution, so our typical procedure is to copy the
01231     // right-hand side into the solution vector on input to _GELS.
01232     y_ = rcp (new mat_type (maxIterCount+1, 1));
01233 
01234     // These cosines and sines encode the Q factor in the QR
01235     // factorization of the upper Hessenberg matrix H_.  We compute
01236     // this by computing H_ into R_ and operating on R_ in place; H_
01237     // itself is left alone.
01238     theCosines_.resize (maxIterCount);
01239     theSines_.resize (maxIterCount);
01240   }
01241 
01242   template<class Scalar, class MV, class OP>
01243   Teuchos::RCP<MV> 
01244   GmresBase<Scalar, MV, OP>::currentNativeResidualVector ()
01245   {
01246     const Scalar zero = STS::zero();
01247     const Scalar one = STS::one();
01248     using Teuchos::is_null;
01249     using Teuchos::Range1D;
01250     using Teuchos::RCP;
01251     using std::endl;
01252 
01253     const bool verboseDebug = true;
01254     const char prefix[] = "Belos::GmresBase::currentNativeResidualVector: ";
01255     const int m = getNumIters();
01256 
01257     std::ostream& dbg = outMan_->stream(Debug);
01258     if (verboseDebug)
01259       dbg << "-- currentNativeResidualVector: getNumIters() = " << m << endl;
01260 
01261     TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*initResVec_) != 1,
01262            std::logic_error,
01263            prefix << "Initial residual vector (initResVec_) has " 
01264            << MVT::GetNumberVecs(*initResVec_)
01265            << " columns, but should only have 1 column.  "
01266            "This should never happen.");
01267     // If left preconditioning is used, then the "native" residual
01268     // vector is in the same vector space as the left-preconditioned
01269     // initial residual vector.  Otherwise, it is in the same space as
01270     // the right-hand side b and the unpreconditioned initial residual
01271     // vector.  Regardless, it's always in the same space as
01272     // *initResVec_, from which we clone it.  Allocate only if
01273     // necessary.
01274     if (nativeResVec_.is_null() || MVT::GetNumberVecs (*nativeResVec_) != 1)
01275       nativeResVec_ = MVT::Clone (*initResVec_, 1);
01276 
01277     // The native residual vector is the same as the initial
01278     // residual vector when the number of iterations is 0.
01279     // Otherwise, we will update *nativeResVec_ below.
01280     MVT::Assign (*initResVec_, *nativeResVec_);
01281     if (m > 0)
01282       {
01283   TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*V_) < m+1,
01284          std::logic_error,
01285          "Only " << MVT::GetNumberVecs(*V_) << " basis vectors "
01286          "were given, but getNumIters()+1=" << (m+1) 
01287          << " of them are required.  "
01288          "This likely indicates a bug in Belos.");
01289   TEST_FOR_EXCEPTION(H_->numRows() < m+1, 
01290          std::logic_error,
01291          "H only has " << H_->numRows() << " rows, but " 
01292          << (m+1) << " rows are required.  "
01293          "This likely indicates a bug in Belos.");
01294   TEST_FOR_EXCEPTION(H_->numCols() < m,
01295          std::logic_error,
01296          "H only has " << H_->numCols() << " columns, but "
01297          << (m+1) << " columns are required.  "
01298          "This likely indicates a bug in Belos.");
01299   TEST_FOR_EXCEPTION(y_->numRows() < m, 
01300          std::logic_error,
01301          "y only has " << y_->numRows() << " entries, but "
01302          << m << " entries are required.  "
01303          "This likely indicates a bug in Belos.");
01304   // Update the QR factorization of the upper Hessenberg matrix,
01305   // and let y[0:m-1] be the projected least-squares problem's
01306   // solution.  This won't do any work if the work has already
01307   // been done.
01308   updateProjectedLeastSquaresProblem();
01309   RCP<const MV> V_view = MVT::CloneView (*V_, Range1D(0, m));
01310   const mat_type H_view (Teuchos::View, *H_, m+1, m);
01311   const mat_type y_view (Teuchos::View, *y_, m, 1);
01312   mat_type H_times_y (m+1, 1);
01313   {
01314     const int err = 
01315       H_times_y.multiply (Teuchos::NO_TRANS, Teuchos::NO_TRANS, 
01316         one, H_view, y_view, zero);
01317     TEST_FOR_EXCEPTION(err != 0, std::logic_error,
01318            "In (F)GMRES, when computing the current native "
01319            "residual vector via the Arnoldi relation, H*y "
01320            "failed due to incompatible dimensions.  "
01321            "This likely indicates a bug in Belos.");
01322   }
01323   // nativeResVec_ := -V * (H * y) + r0 = r0 - V (H y)
01324   // (where we stored r0 in nativeResVec_).
01325   MVT::MvTimesMatAddMv (-one, *V_view, H_times_y, one, *nativeResVec_);
01326       }
01327     // Recport the "native" residual vector's norm, scaled by the
01328     // initial residual norm (if the latter is not zero).
01329     if (verboseDebug)
01330       {
01331   std::vector<magnitude_type> theNorm (1);
01332   MVT::MvNorm (*nativeResVec_, theNorm);
01333   dbg << "---- ||r_0||_2 = " << initialResidualNorm_ << endl;
01334   if (initialResidualNorm_ == STM::zero())
01335     {
01336       dbg << "---- ||Native residual vector||_2 = " << theNorm[0] << endl;
01337     }
01338   else
01339     {
01340       const magnitude_type relResNorm = theNorm[0] / initialResidualNorm_;
01341       dbg << "---- ||Native residual vector||_2 / ||r_0||_2 = " 
01342     << relResNorm << endl;
01343     }
01344   RCP<MV> xUpdate = getCurrentUpdate();
01345   // FIXME (mfh 28 Feb 2011) Is this correct for Flexible GMRES
01346   // too?
01347   RCP<MV> r = MVT::CloneCopy (*initResVec_);
01348   RCP<MV> AxUpdate = MVT::CloneCopy (*initResVec_);
01349   lp_->apply (*xUpdate, *AxUpdate);
01350   // r := initResVec_ - A * xUpdate
01351   MVT::MvAddMv (STS::one(), *initResVec_, -STS::one(), *AxUpdate, *r);
01352   MVT::MvNorm (*r, theNorm);
01353 
01354   dbg << "---- ||b-Ax||_2 = " << theNorm[0] << endl;
01355   if (initialResidualNorm_ != STM::zero())
01356     {
01357       const magnitude_type relResNorm = theNorm[0] / initialResidualNorm_;
01358       dbg << "---- ||b-Ax||_2 / ||r_0|| = " << relResNorm << endl;
01359     }
01360       }
01361     return nativeResVec_;
01362   }
01363 
01364   template<class Scalar, class MV, class OP>
01365   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
01366   GmresBase<Scalar, MV, OP>::currentNativeResidualNorm ()
01367   {
01368     // Update the least-squares problem if necessary.  This computes
01369     // the current native residual norm for "free."
01370     const magnitude_type theNorm = updateProjectedLeastSquaresProblem();
01371     return theNorm;
01372   }
01373 
01374   template<class Scalar, class MV, class OP>
01375   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
01376   GmresBase<Scalar, MV, OP>::updateProjectedLeastSquaresProblem()
01377   {
01378     using std::endl;
01379     const bool verboseDebug = false;
01380 
01381     const int startCol = (lastUpdatedCol_ < 0) ? 0 : lastUpdatedCol_ + 1;
01382     const int endCol = getNumIters() - 1;
01383     const int m = getNumIters();
01384     blas_type blas;
01385     lapack_type lapack;
01386 
01387     const char prefix[] = "Belos::GmresBase::updateProjectedLeastSquaresProblem: ";
01388     std::ostream& dbg = outMan_->stream(Debug);
01389     dbg << "---- updateProjectedLeastSquaresProblem: ";
01390 
01391     TEST_FOR_EXCEPTION(startCol > endCol+1, std::logic_error,
01392            prefix << "Somehow, GmresBase updated the QR "
01393            "factorization of the upper Hessenberg matrix past the "
01394            "last updated column, which was " << lastUpdatedCol_ 
01395            << ". The rightmost column to update is " << endCol 
01396            << ".  This is likely a bug in GmresBase.");
01397     if (m == 0)
01398       {
01399   dbg << "Nothing to do: # iters = 0" << endl;
01400   // Assign here, just in case we haven't done this yet.
01401   nativeResidualNorm_ = initialResidualNorm_; 
01402   return initialResidualNorm_;
01403       }
01404     else if (startCol == endCol+1) // Nothing to do
01405       {
01406   dbg << "Nothing to do: No columns to update (startCol == endCol+1 "
01407     "== " << startCol << ")" << endl;
01408   return nativeResidualNorm_;
01409       }
01410     else
01411       dbg << endl;
01412 
01413     const int numCols = endCol - startCol + 1;
01414     dbg << "------ # cols to update: " << numCols << endl;
01415 
01416     // Define some references to ease notation.
01417     mat_type& H = *H_;
01418     mat_type& R = *R_;
01419     mat_type& y = *y_;
01420     mat_type& z = *z_;
01421 
01422     // FIXME (mfh 24 Feb 2011) 
01423     //
01424     // The "if 0" part attempts to update the QR factorization of the
01425     // upper Hessenberg matrix using givens rotations.  I haven't
01426     // quite gotten that to work yet, so I've provided as an
01427     // alternative a least-squares solve using LAPACK's _GELS routine.
01428     // This will be slightly more expensive, but it should only matter
01429     // when solving very small problems or when the iteration count is
01430     // large.
01431 #if 0    
01432     {
01433       // Copy columns [startCol, endCol] (inclusive) of the current
01434       // upper Hessenberg matrix H, into the upper triangular matrix
01435       // R.  (We avoid overwriting H because having the original upper
01436       // Hessenberg matrix around can save us some work when computing
01437       // the "native" residual.)  Columns [0, startCol-1] of R are
01438       // already upper triangular; we don't have to do any more work
01439       // there.
01440       const int numRows = endCol + 2;
01441       const mat_type H_view (Teuchos::View, *H_, numRows, numCols, 0, startCol);
01442       mat_type R_view (Teuchos::View, *R_, numRows, numCols, 0, startCol);
01443       R_view.assign (H_view);
01444 
01445       if (verboseDebug)
01446   dbg << "------ R factor of upper Hessenberg matrix before rotations:" << endl
01447       << mat_type(Teuchos::View, *R_, m+1, m) << endl;
01448     }
01449 
01450     const scalar_type zero = STS::zero();
01451     const scalar_type one = STS::one();
01452 
01453     // Update columns [startCol, endCol] of the QR factorization of
01454     // the upper Hessenberg matrix.  Use Givens rotations to maintain
01455     // the factorization.  We use a left-looking update procedure,
01456     // since Arnoldi is always adding new column(s) on the right of
01457     // the upper Hessenberg matrix.
01458     const int LDR = R.stride();
01459     if (startCol > 0)
01460       dbg << "-- Applying " << startCol << " previous Givens rotation"
01461     << (startCol != 1 ? "s" : "") << " to new columns of H (in R)"
01462     << endl;
01463     for (int j = 0; j < startCol; ++j)
01464       {
01465   // Apply all previous Givens rotations to the new columns of
01466   // the upper Hessenberg matrix (as stored in R).
01467   //
01468   // FIXME (mfh 25 Dec 2010) Teuchos::BLAS wants nonconst
01469   // pointers for the sine and cosine arguments to ROT().  They
01470   // should be const, since ROT doesn't change them (???).  Fix
01471   // this in Teuchos.
01472   magnitude_type theCosine = theCosines_[j];
01473   scalar_type theSine = theSines_[j];
01474   blas.ROT (numCols, &R(j,j), LDR, &R(j+1,j), LDR, 
01475       &theCosine, &theSine);
01476       }
01477     for (int j = startCol; j < endCol; ++j)
01478       {
01479   // Calculate new Givens rotation for [R(j,j); R(j+1,j)]
01480   magnitude_type theCosine;
01481   scalar_type theSine;
01482   blas.ROTG (&R(j,j), &R(j+1,j), &theCosine, &theSine);
01483   theCosines_[j] = theCosine;
01484   theSines_[j] = theSine;
01485   // Clear the subdiagonal element R(j+1,j)
01486   R(j+1,j) = zero;
01487   // Update the "trailing matrix."  The "if" check is not
01488   // strictly necessary, but ensures that the references to
01489   // column j+1 of R always makes sense (even though ROT
01490   // shouldn't access column j+1 if endCol-1-j==0).
01491   if (j < endCol)
01492     blas.ROT (endCol-1-j, &R(j,j+1), LDR, &R(j+1,j+1), LDR, 
01493         &theCosine, &theSine);
01494   // Update the right-hand side of the least-squares problem
01495   // with the new Givens rotation.
01496   blas.ROT (1, &z(j,0), 1, &z(j+1,0), 1, &theCosine, &theSine);
01497       }
01498 
01499     // The absolute value of the last element of z gives the current
01500     // "native" residual norm.
01501     nativeResidualNorm_ = STS::magnitude ((*z_)(m,0));
01502 
01503     // Now that we have the updated R factor of H, and the updated
01504     // right-hand side z, solve the least-squares problem by solving
01505     // the linear system Ry=z.
01506     // 
01507     // TRSM() overwrites the right-hand side with the solution, so
01508     // copy z into y.
01509     {
01510       mat_type y_view (Teuchos::View, *y_, m, 1);
01511       mat_type z_view (Teuchos::View, *z_, m, 1);
01512       y_view.assign (z_view);
01513     }
01514     // Solve Ry = z for y.
01515     {
01516       dbg << "-- R factor of upper Hessenberg matrix:" << endl
01517     << mat_type(Teuchos::View, *R_, m, m) << endl
01518     << "-- Right-hand side of Ry=z:" << endl
01519     << mat_type(Teuchos::View, *y_, m, 1) << endl;
01520     }
01521     blas.TRSM(Teuchos::RIGHT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
01522         Teuchos::NON_UNIT_DIAG, m, 1,
01523         one, R_->values(), R_->stride(), y_->values(), y_->stride());
01524 #else // not 0
01525     // Let's solve the projected least-squares problem the easy way:
01526     // Rather than messing around with Givens rotations, let's just
01527     // invoke the appropriate LAPACK routine.
01528     {
01529       mat_type H_view (Teuchos::View, H, m+1, m);
01530       mat_type R_view (Teuchos::View, R, m+1, m);
01531       // y gets an extra entry, since it needs to hold a copy of the
01532       // whole right-hand side (see below).
01533       mat_type y_view (Teuchos::View, y, m+1, 1);
01534       mat_type z_view (Teuchos::View, z, m+1, 1);
01535 
01536       // Copy _all_ of H into R.  Remember that _GELS overwrites the
01537       // entire input matrix with the implicitly stored Q factor.  If
01538       // we call _GELS each time, we then have to recopy all of H into
01539       // R each time.
01540       R_view.assign (H_view);
01541       // The least-squares solver overwrites the right-hand side with
01542       // the solution, so first copy z into y.
01543       y_view.assign (z_view);
01544 
01545       dbg << "------ Current upper Hessenberg matrix (" << m+1 << " by " 
01546     << m << "): " << R_view << endl;
01547       dbg << "------ Current projected right-hand side (" << m+1 << " by " 
01548     << 1 << "): " << y_view << endl;
01549 
01550       // Workspace query.
01551       int info = 0;
01552       Scalar lwork_scalar = STS::zero();
01553       if (verboseDebug)
01554   dbg << "------ _GELS workspace query: " << endl;
01555       // Recall that _GELS overwrites the right-hand side with the
01556       // solution.  We've copied z_view into y_view in preparation for
01557       // this.  Note that this means the stride of y_view must be no
01558       // less than the number of rows in R_view.
01559       TEST_FOR_EXCEPTION(y_view.stride() < m+1, std::logic_error,
01560        "Before calling LAPACK _GELS: y_view has stride "
01561        << y_view.stride() << ", which is less than the "
01562        "number of rows " << m+1 << " in R_view.  This is "
01563        "an internal Belos bug.");
01564       lapack.GELS ('N', m+1, m, 1, 
01565        R_view.values(), R_view.stride(),
01566        y_view.values(), y_view.stride(),
01567        &lwork_scalar, -1, &info);
01568       TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01569        "LAPACK _GELS workspace query failed with INFO = " 
01570        << info << ", for an " << (m+1) << " x " << m 
01571        << " matrix with 1 right-hand side.");
01572       // Cast workspace from Scalar (which may be complex, hence the
01573       // request for the real part -- don't call magnitude() since it
01574       // may overflow due to squaring and square root) to int.
01575       // Hopefully LAPACK doesn't ever overflow int this way.
01576       const int lwork = static_cast<int> (STS::real (lwork_scalar));
01577       dbg << "------ LAPACK _GELS workspace size: " << lwork << endl;
01578       TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
01579        "LAPACK _GELS workspace query returned LWORK = " 
01580        << lwork << " < 0.  This should never happen.");
01581       std::vector<Scalar> work (lwork, STS::zero()); // Workspace
01582       // Solve the least-squares problem.  The ?: operator prevents
01583       // accessing the first element of the work array, if it has
01584       // length zero.
01585       lapack.GELS ('N', m+1, m, 1, 
01586        R_view.values(), R_view.stride(),
01587        y_view.values(), y_view.stride(),
01588        (lwork > 0 ? &work[0] : (Scalar*) NULL), 
01589        lwork, &info);
01590       TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01591        "Solving projected least-squares problem with LAPACK "
01592        "_GELS failed with INFO = " << info << ", for an " 
01593        << (m+1) << " x " << m << " matrix with 1 right-hand "
01594        "side.");
01595       // Extract the projected least-squares problem's residual error.
01596       // It's the magnitude of the last entry of y_view on output from
01597       // LAPACK's least-squares solver.  That is the absolute residual
01598       // norm, so scale by the initial residual norm.  Be sure not to
01599       // divide by zero (though if the initial residual norm is zero,
01600       // GMRES shouldn't progress anyway).
01601       nativeResidualNorm_ = (initialResidualNorm_ == STM::zero()) ? 
01602   STS::magnitude (y_view(m, 0)) : 
01603   STS::magnitude (y_view(m, 0)) / initialResidualNorm_;
01604       dbg << "Native residual norm (via _GELS): " << nativeResidualNorm_ << endl;
01605       dbg << "------ Current projected least-squares solution (" << m << " by " 
01606     << 1 << "): " << mat_type(Teuchos::View, y_view, m, 1) << endl;
01607     }
01608 #endif // 0
01609 
01610     // Remember the rightmost updated column.
01611     lastUpdatedCol_ = endCol;
01612 
01613     // Return the "native" residual norm (which is the absolute
01614     // residual error from solving the projected least-squares
01615     // problem).
01616     return nativeResidualNorm_;
01617   }
01618 
01619 
01620   template<class Scalar, class MV, class OP>
01621 #if 0
01622   std::pair< Teuchos::RCP<MV>, typename Teuchos::ScalarTraits<Scalar>::magnitudeType >
01623 #else // not 0
01624   Teuchos::RCP<MV>
01625 #endif // 0
01626   GmresBase<Scalar, MV, OP>::getCurrentUpdate ()
01627   {
01628     using Teuchos::is_null;
01629     using Teuchos::Range1D;
01630     using Teuchos::RCP;
01631     const Scalar one = STS::one();
01632     const Scalar zero = STS::zero();
01633 
01634     if (is_null (xUpdate_))
01635       {
01636   // xUpdate_ comes from the Z_ space and not the V_ space in
01637   // the flexible variant of GMRES.  This makes a difference
01638   // only if there are multiple vector spaces involved, that
01639   // is if the preconditioner maps V-space vectors to Z-space
01640   // vectors, and the matrix maps Z-space vectors to V-space
01641   // vectors (where V and Z are different spaces).  This is
01642   // legitimate as long as the initial guess for the solution
01643   // (x0) comes from the Z space.
01644   if (flexible_)
01645     xUpdate_ = MVT::Clone (*Z_, 1);
01646   else
01647     xUpdate_ = MVT::Clone (*V_, 1);
01648       }
01649     if (getNumIters() == 0)
01650       {
01651   MVT::MvInit (*xUpdate_, zero);
01652   return xUpdate_;
01653       }
01654 
01655     // When one iteration of GMRES is completed, the upper
01656     // Hessenberg matrix H is 2 by 1.  Thus, we want to update
01657     // columns up to and including curNumIters_ - 1: actually,
01658     // columns [lastUpdatedCol_+1, curNumIters_-1] (inclusive).  The
01659     // solution update gives you the current "native" residual norm
01660     // for free.  We could save and return that as well...
01661     //const magnitude_type nativeResNorm = updateProjectedLeastSquaresProblem();
01662     (void) updateProjectedLeastSquaresProblem();
01663 
01664     // Current iteration count does not include the initial basis vector.
01665     const int m = getNumIters();
01666     // Basis for solution update coefficients; Flexible GMRES uses the
01667     // preconditioned basis here.
01668     RCP<const MV> Z_view = flexible_ ? 
01669       MVT::CloneView (*Z_, Range1D(0, m-1)) : 
01670       MVT::CloneView (*V_, Range1D(0, m-1));
01671     const mat_type y_view (Teuchos::View, *y_, m, 1);
01672     // xUpdate_ := Z_view * y_view
01673     MVT::MvTimesMatAddMv (one, *Z_view, y_view, zero, *xUpdate_);
01674     return xUpdate_;
01675   }
01676 
01677   template<class Scalar, class MV, class OP>
01678   void
01679   GmresBase<Scalar, MV, OP>::
01680   setInitialResidual(const int maxIterCount,
01681          const bool reallocateBasis)
01682   {
01683     using Teuchos::Range1D;
01684     using Teuchos::RCP;
01685     using std::endl;
01686 
01687     // Message fragment for error messages.
01688     const char prefix[] = "Belos::GmresBase::setInitialResidual: ";
01689     std::ostream& dbg = outMan_->stream(Debug);
01690 
01691     TEST_FOR_EXCEPTION(maxIterCount < 0, std::invalid_argument,
01692            prefix << "maxIterCount = " << maxIterCount 
01693            << " < 0.");
01694     // Do we have a left preconditioner?
01695     const bool haveLeftPrec = ! lp_->getLeftPrec().is_null();
01696     if (haveLeftPrec)
01697       dbg << "-- Have a left preconditioner" << endl;
01698     else
01699       dbg << "-- No left preconditioner" << endl;
01700 
01701     // mfh 22 Feb 2011: LinearProblem has the annoying quirk that
01702     // getInit(Prec)ResVec() return the residual vector for the whole
01703     // linear problem (X_ and B_), not for the "current" linear
01704     // problem (curX_ and curB_, returned by getCurrLHSVec()
01705     // resp. getCurrRHSVec()).  In order to get the residual vector(s)
01706     // for the _current_ linear problem, we have to use getLSIndex()
01707     // to extract the column ind(ex/ices) and take the corresponding
01708     // column of getInit(Prec)ResVec().  It should be exactly one
01709     // column, since this GMRES solver only knows how to solve for one
01710     // right-hand side at a time.
01711     RCP<const MV> R_full = 
01712       haveLeftPrec ? lp_->getInitPrecResVec() : lp_->getInitResVec();
01713     std::vector<int> inputIndices = lp_->getLSIndex();
01714     TEST_FOR_EXCEPTION(inputIndices.size() == 0, 
01715            std::invalid_argument,
01716            "The LinearProblem claims that there are zero linear "
01717            "systems to solve: getLSIndex() returns an index "
01718            "vector of length zero.");
01719     dbg << "-- Input indices: [";
01720     for (std::vector<int>::size_type k = 0; k < inputIndices.size(); ++k)
01721       {
01722   dbg << inputIndices[k];
01723   if (k < inputIndices.size() - 1)
01724     dbg << ", ";
01725       }
01726     dbg << "]" << endl;
01727 
01728     // Some of the indices returned by getLSIndex() might be -1.
01729     // These are for column(s) of getCurr{L,R}HSVec() that don't
01730     // correspond to any actual linear system(s) that we want to
01731     // solve.  They may be filled in by block iterative methods in
01732     // order to make the blocks full rank or otherwise improve
01733     // convergence.  We shouldn't have any of those columns here, but
01734     // we check just to make sure.
01735     std::vector<int> outputIndices;
01736     std::remove_copy_if (inputIndices.begin(), inputIndices.end(), 
01737        std::back_inserter (outputIndices), 
01738        std::bind2nd (std::equal_to<int>(), -1));
01739     dbg << "-- Output indices: [";
01740     for (std::vector<int>::size_type k = 0; k < outputIndices.size(); ++k)
01741       {
01742   dbg << outputIndices[k];
01743   if (k < outputIndices.size() - 1)
01744     dbg << ", ";
01745       }
01746     dbg << "]" << endl;
01747 
01748     // outputIndices should have exactly one entry, corresponding to
01749     // the one linear system to solve.  Our GMRES implementation can
01750     // only solve for one right-hand side at a time.
01751     if (outputIndices.size() != 1 || outputIndices[0] == -1)
01752       {
01753   std::string inputIndicesAsString;
01754   {
01755     std::ostringstream os;
01756     os << "[";
01757     std::copy (inputIndices.begin(), inputIndices.end(),
01758          std::ostream_iterator<int>(os, " "));
01759     os << "]";
01760     inputIndicesAsString = os.str();
01761   }
01762   std::string outputIndicesAsString;
01763   {
01764     std::ostringstream os;
01765     os << "[";
01766     std::copy (outputIndices.begin(), outputIndices.end(),
01767          std::ostream_iterator<int>(os, " "));
01768     os << "]";
01769     outputIndicesAsString = os.str();
01770   }
01771   std::ostringstream os;
01772   os << prefix << "The LinearProblem instance's getLSIndex() method "
01773     "returns indices " << inputIndicesAsString << ", of which the "
01774     "following are not -1: " << outputIndicesAsString << ".";
01775   TEST_FOR_EXCEPTION(outputIndices.size() != 1,
01776          std::invalid_argument, 
01777          os.str() << "  The latter list should have length "
01778          "exactly 1, since our GMRES implementation only "
01779          "knows how to solve for one right-hand side at a "
01780          "time.");
01781   TEST_FOR_EXCEPTION(outputIndices[0] == -1,
01782          std::invalid_argument,
01783          os.str() << "  the latter list contains no "
01784          "nonnegative entries, meaning that there is no "
01785          "valid linear system to solve.");
01786       }
01787     // View of the "current" linear system's residual vector.
01788     RCP<const MV> r0 = MVT::CloneView (*R_full, outputIndices);
01789     // Sanity check that the residual vector has exactly 1 column.
01790     // We've already checked this a different way above.
01791     const int numResidVecs = MVT::GetNumberVecs (*r0);
01792     TEST_FOR_EXCEPTION(numResidVecs != 1, std::logic_error,
01793            prefix << "Residual vector has " << numResidVecs 
01794            << " columns, but should only have one.  "
01795            "Should never get here.");
01796     // Save a deep copy of the initial residual vector.
01797     initResVec_ = MVT::CloneCopy (*r0);
01798 
01799     // Allocate space for the V_ basis vectors if necessary.
01800     //
01801     // If left preconditioning is used, then the V_ vectors are in the
01802     // same vector space as the (left-)preconditioned initial residual
01803     // vector.  Otherwise they are in the same space as the
01804     // unpreconditioned initial residual vector.  That's r0 in either
01805     // case.
01806     //
01807     // FIXME (mfh 22 Feb 2011) We should really check to make sure
01808     // that the columns of V_ (if previously allocated) are in the
01809     // same vector space as r0.  Belos::MultiVecTraits doesn't give us
01810     // a way to do that.  Just checking if the number of rows is the
01811     // same is not enough.  The third Boolean argument
01812     // (reallocateBasis) lets us add in the check externally if
01813     // MultiVecTraits gains that ability sometime.
01814     if (V_.is_null() || 
01815   MVT::GetNumberVecs(*V_) != maxIterCount+1 || 
01816   reallocateBasis)
01817       {
01818   dbg << "-- Reallocating V_ basis with " << (maxIterCount+1) 
01819       << " columns" << endl;
01820   V_ = MVT::Clone (*initResVec_, maxIterCount+1);
01821       }
01822     else
01823       dbg << "-- V_ basis has alread been allocated with " 
01824     << (maxIterCount+1) << " columns." << endl;
01825 
01826     // Initial residual norm is computed with respect to the inner
01827     // product defined by the OrthoManager.  Since the basis
01828     // vectors are orthonormal (resp. unitary) with respect to that
01829     // inner product, this ensures that the least-squares minimization
01830     // is also computed with respect to that inner product.
01831     std::vector<magnitude_type> result (1);
01832     ortho_->norm (*r0, result);
01833     initialResidualNorm_ = result[0];
01834     dbg << "-- Initial residual norm (as computed by the OrthoManager): " 
01835   << initialResidualNorm_ << endl;
01836 
01837     // z_ is the right-hand side of the projected least-squares
01838     // problem.  If we reallocated z_, we have to reset it by setting
01839     // its first entry to the initial residual norm.
01840     if (z_.is_null() || z_->numRows() != maxIterCount+1)
01841       z_ = rcp (new mat_type (maxIterCount+1, 1));
01842     else
01843       (void) z_->putScalar (STS::zero());
01844     // The residual norm is a magnitude_type; promote to Scalar.
01845     (*z_)(0,0) = Scalar (initialResidualNorm_);
01846 
01847     // Copy the initial residual vector into the first column of V_,
01848     // and scale the vector by its norm.
01849     RCP<MV> v1 = MVT::CloneViewNonConst (*V_, Range1D(0,0));
01850     MVT::Assign (*r0, *v1);
01851     MVT::MvScale (*v1, STS::one() / initialResidualNorm_);
01852   }
01853 
01854   template<class Scalar, class MV, class OP>
01855   void 
01856   GmresBase<Scalar, MV, OP>::backOut (const int numIters)
01857   {
01858     const Scalar zero = STS::zero();
01859 
01860     TEST_FOR_EXCEPTION(numIters < 0, std::invalid_argument,
01861            "The GMRES iteration count cannot be less than "
01862            "zero, but you specified numIters = " << numIters);
01863     TEST_FOR_EXCEPTION(numIters > getNumIters(), std::invalid_argument,
01864            "The GMRES iteration count cannot be reset to more than "
01865            "its current iteration count " << getNumIters()
01866            << ", but you specified numIters = " << numIters << ".");
01867     // Reset the R factor in the QR factorization of H_, and the
01868     // right-hand side z_ of the projected least-squares problem.  We
01869     // will "replay" the first numIters steps of the QR factorization
01870     // below.
01871     //
01872     // The integer flag return value of SerialDenseMatrix's
01873     // putScalar() is not informative; it's always zero.
01874     (void) R_->putScalar(zero); // zero out below-diagonal entries
01875     *R_ = *H_; // deep copy
01876     (void) y_->putScalar(zero);
01877     (void) z_->putScalar(zero);
01878     (*z_)(0,0) = Scalar (initialResidualNorm_);
01879 
01880     // Replay the first numIters-1 Givens rotations.
01881     lastUpdatedCol_ = -1;
01882     curNumIters_ = numIters;
01883     (void) updateProjectedLeastSquaresProblem ();
01884   }
01885 
01886   template<class Scalar, class MV, class OP>
01887   void 
01888   GmresBase<Scalar, MV, OP>::
01889   restart (const int maxIterCount)
01890   {
01891     using Teuchos::RCP;
01892     using Teuchos::null;
01893     const Scalar zero = STS::zero();
01894 
01895     // This would be the place where subclasses may implement things
01896     // like recycling basis vectors for the next restart cycle.  Note
01897     // that the linear problem hasn't been updated yet; this happens
01898     // below.
01899     preRestartHook (maxIterCount);
01900 
01901     // Make sure that the LinearProblem object gets the current
01902     // approximate solution (which becomes the initial guess for the
01903     // upcoming new restart cycle), so that the "exact" residuals are
01904     // correct.
01905     (void) getCurrentUpdate (); // results stored in xUpdate_
01906     updateSolution ();
01907 
01908     // Check whether there is a right preconditioner, since (at least
01909     // in theory) the caller could have added a right preconditioner
01910     // after a restart cycle.  (Is that legitimate?  Adding or
01911     // changing a right preconditioner between restarts wouldn't
01912     // change the initial residual vector, and therefore wouldn't
01913     // change the V_ basis.  The V_ basis is all that matters for the
01914     // restart.)
01915     const bool needToAllocateZ = flexible_ && ! lp_->getRightPrec().is_null();
01916 
01917     // This (re)allocates the V_ basis if necessary.  It also
01918     // (re)allocates the projected least-squares problem's right-hand
01919     // side z_ if necessary, and initializes it regardless.
01920     setInitialResidual (maxIterCount);
01921 
01922     // If the new max iteration count has changed, reallocate space
01923     // for the Z_ basis vectors and the projected least-squares
01924     // problem.  setInitialResidual() already does that for the V_
01925     // basis.
01926     if (maxNumIters_ != maxIterCount)
01927       {
01928   maxNumIters_ = maxIterCount;
01929   // Initial guess.
01930   RCP<MV> x0 = lp_->getCurrLHSVec();
01931   Z_ = needToAllocateZ ? MVT::Clone (*x0, maxIterCount+1) : null;
01932   H_ = rcp (new mat_type (maxIterCount+1, maxIterCount));
01933   R_ = rcp (new mat_type (maxIterCount+1, maxIterCount));
01934   y_ = rcp (new mat_type (maxIterCount+1, 1));
01935       }
01936     // Even if the max iteration count hasn't changed, we still want
01937     // to fill in the projected least-squares problem data with zeros,
01938     // to ensure (for example) that R_ is upper triangular and that H_
01939     // is upper Hessenberg.
01940     (void) H_->putScalar(zero);
01941     (void) R_->putScalar(zero);
01942     (void) y_->putScalar(zero);
01943 
01944     lastUpdatedCol_ = -1; // column updates start with zero
01945     curNumIters_ = 0;
01946     flexible_ = needToAllocateZ;
01947 
01948     // This would be the place where subclasses may implement things
01949     // like orthogonalizing the first basis vector (after restart)
01950     // against the recycled subspace.
01951     postRestartHook();
01952   }
01953 
01954 } // namespace Belos
01955 
01956 #endif // __Belos_GmresBase_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines