Belos Package Browser (Single Doxygen Collection) Development
BelosCaGmres.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_CaGmres_hpp
00043 #define __Belos_CaGmres_hpp
00044 
00048 
00049 #include "BelosGmresBase.hpp"
00050 
00051 namespace Belos {
00052 
00054 
00055   
00063   class CaGmresUpdateFailure : public BelosError {
00064   public:
00065     CaGmresUpdateFailure(const std::string& what_arg) : 
00066       BelosError(what_arg) {}
00067   };
00068 
00074   class CaGmresOrthoFailure : public BelosError {
00075   public:
00076     CaGmresOrthoFailure (const int theCurOuterIter,
00077        const int theNumBasisVectors,
00078        const int theCandidateBasisLength,
00079        const int theRank) :
00080       curOuterIter_ (theCurOuterIter),
00081       numBasisVectors_ (theNumBasisVectors)
00082       candidateBasisLength_ (theCandidateBasisLength),
00083       rank_ (theRank)
00084     {
00085       // This may throw an exception if there isn't enough memory with
00086       // which to allocate a string.  std::exception subclass
00087       // constructors really should not throw exceptions.  Howver, if
00088       // there isn't enough memory left to allocate a short string,
00089       // then the program is likely in deep trouble anyway.
00090       std::ostringstream os;
00091       os << "CA-GMRES: After constructing " << numBasisVectors_ 
00092    << " successfully, outer iteration" << (curOuterIter_+1) 
00093    << " failed: the candidate basis of length " 
00094    << candidateBasisLength_ << " only has rank " 
00095    << rank_ << ".";
00096       what_ = os.str();
00097     }
00098 
00100     int curOuterIter() const { return curOuterIter_; }
00101 
00103     int numBasisVectors() const { return numBasisVectors_; }
00104 
00106     int candidateBasisLength() const { return candidateBasisLength_; }
00107 
00109     int rank() const { return rank_; }
00110 
00112     const char* const what() { 
00113       return what_.c_str();
00114     }
00115 
00116   private:
00117     int curOuterIter_, numBasisVectors_, candidateBasisLength_, rank_;
00118     std::string what_;
00119   };
00120   
00122 
00145   template<class Scalar, class MV, class OP>
00146   class CaGmres : public GmresBase<Scalar, MV, OP> {
00147   public:
00148     typedef Scalar scalar_type;
00149     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00150     typedef MV multivector_type;
00151     typedef OP operator_type;
00152 
00153   private:
00156     typedef GmresBase<Scalar, MV, OP> base_type;
00157     typedef Teuchos::ScalarTraits<Scalar> STS;
00158     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00159     typedef MultiVecTraits<Scalar, MV> MVT;
00160     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00161 
00162   public:
00212     CaGmres (const Teuchos::RCP<LinearProblem<Scalar,MV,OP> >& lp,
00213        const Teuchos::RCP<const OrthoManager<Scalar, MV> >& ortho,
00214        const Teuchos::RCP<Akx<Scalar, MV> >& akx,
00215        const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00216        const int maxIterCount,
00217        const bool flexible,
00218        const Teuchos::RCP<const Teuchos::ParameterList>& params)
00219       : GmresBase (lp, ortho, outMan, maxIterCount, flexible),
00220   params_ (validateParameters (params)),
00221   akxParams_ (akxParameters (params_)),
00222   akx_ (initAkx (akx, akxParams_)),
00223   candidateBasisLength_ (initialCandidateBasisLength (akx_, akxParams_))
00224     {}
00225 
00227     virtual ~CaGmres() {}
00228 
00238     static Teuchos::RCP<const Teuchos::ParameterList>
00239     akxParameters (AkxFactory<Scalar, MV, OP>& akxFactory,
00240        const Teuchos::RCP<const Teuchos::ParameterList>& params)
00241     {
00242       using Teuchos::ParameterList;
00243       using Teuchos::parameterList;
00244       using Teuchos::RCP;
00245       using Teuchos::rcp;
00246 
00247       if (params.is_null())
00248   return akxFactory.getDefaultParameters ();
00249       else
00250   {
00251     RCP<ParameterList> akxParams;
00252     bool gotAkxParams = false;
00253     const std::string& listName = akxFactory.parameterListName ();
00254     //
00255     // Users may have specified the parameter list either as a
00256     // sublist, or as an RCP<const ParameterList>.  Try both.
00257     // Regardless, make a deep copy of the result, for safety.
00258     //
00259     try { // Could it be a sublist?
00260       const ParameterList& result = params->sublist (listName);
00261       akxParams = parameterList (result);
00262       gotAkxParams = true;
00263     } catch (Teuchos::Exceptions::InvalidParameter&) {
00264       // No luck; gotAkxParams stays false.
00265     }
00266     try { // Could it be an RCP<const ParameterList>?
00267       RCP<const ParameterList> result = 
00268         params->get<RCP<const ParameterList> (listName);
00269       akxParams = parameterList (*result);
00270       gotAkxParams = true;
00271     } catch (Teuchos::Exceptions::InvalidParameter&) {
00272       // No luck; gotAkxParams stays false.
00273     }
00274     if (gotAkxParams)
00275       return akxParams;
00276     else
00277       return akxFactory.getDefaultParameters ();
00278   }
00279     }
00280 
00288     static Teuchos::RCP<Akx<Scalar, MV> >
00289     initAkx (const Teuchos::RCP<Akx<Scalar, MV> >& akx, 
00290        const Teuchos::RCP<const Teuchos::ParameterList>& akxParams)
00291     {
00292       using Teuchos::ParameterList;
00293       using Teuchos::RCP;
00294 
00295       if (! akx.is_null())
00296   return akx;
00297       else
00298   {
00299     RCP<const OP> A = lp->getOperator();
00300     RCP<const OP> M_left = lp->getLeftPrec();
00301     RCP<const OP> M_right = lp->getRightPrec();
00302     return akxFactory.makeOpAkx (A, M_left, M_right, akxParams);
00303   }
00304     }
00305 
00306     int 
00307     initialCandidateBasisLength (const Teuchos::RCP<Akx<Scalar, MV> >& akx, 
00308          const int maxIterCount,
00309          const Teuchos::RCP<const Teuchos::ParameterList>& akxParams) const
00310     {
00311       TEUCHOS_TEST_FOR_EXCEPTION(akx.is_null(), std::logic_error,
00312        "The matrix powers kernel implementation must be "
00313        "initialized before asking for the initial candidate "
00314        "basis length, since the former gives a strict upper "
00315        "bound on the latter.");
00316       TEUCHOS_TEST_FOR_EXCEPTION(akxParams.is_null(), std::logic_error,
00317        "The list of parameters must be initialized before "
00318        "asking for the initial candidate basis length.");
00319       // An initial basis length of 5 is a reasonable first guess.
00320       // Also make sure you don't exceed the max iteration count.  We
00321       // may validly call maxNumIters() in the CaGmres constructor, as
00322       // long as the GmresBase constructor has been invoked first.
00323       int initLength = std::min (5, this->maxNumIters ());
00324 
00325       // Look for the basis length parameter.  ("s" is a shorthand
00326       // notation from my PhD dissertation.)
00327       const char* candidateNames[] = {"Basis Length", "s"};
00328       const int numCandidateNames = 2;
00329       bool gotLength = false;
00330       for (int k = 0; k < numCandidateNames && ! gotLength; ++k)
00331   {
00332     try {
00333       initLength = akxParams->get<int> (candidateNames[k]);
00334       gotLength = true;
00335     } catch (Teuchos::Exceptions::InvalidParameter&) {
00336       // No luck; gotLength stays false.
00337     }
00338   }
00339       // The Akx implementation may only bound the candidate basis
00340       // length loosely from above, but it may also have insights into
00341       // the numerical properties of the matrix and the largest
00342       // reasonable basis length.
00343       if (gotLength)
00344   return std::min (initLength, akx_->maxCandidateBasisLength());
00345       else
00346   return std::min (akx_->recommendedCandidateBasisLength(),
00347        akx_->maxCandidateBasisLength());
00348     }
00349 
00350     virtual bool canExtendBasis() const { 
00351       return this->getNumIters() < this->maxNumIters();
00352     }
00353 
00354     int 
00355     newCandidateBasisLength () const
00356     {
00357       const int k = this->getNumIters(); 
00358       const int m = this->maxNumIters();
00359 
00360       // V[0 : k] is already occupied, and V[0 : m] is the max
00361       // (inclusive) range.  V[k+1 : k+s] will be used to store the
00362       // candidate basis.
00363       //
00364       // Reduce candidate basis length as necessary, if we are running
00365       // out of iterations to perform.
00366       return (k + candidateBasisLength_ > m) ? (m - k) : candidateBasisLength_;
00367     }
00368 
00369     void
00370     extendBasis (Teuchos::RCP<MV>& V_cur, 
00371      Teuchos::RCP<MV>& Z_cur)
00372     {
00373       using Teuchos::null;
00374       using Teuchos::Range1D;
00375       using Teuchos::RCP;
00376       using std::endl;
00377       const bool verboseDebug = false;
00378       //
00379       // mfh 28 Feb 2011: The use of "this->..." here and elsewhere to
00380       // refer to GmresBase member data is obligatory, since we are
00381       // inheriting from a templated class.  See the C++ FAQ:
00382       //
00383       // http://www.parashift.com/c++-faq-lite/templates.html#faq-35.19
00384       // 
00385       std::ostream& dbg = this->outMan_->stream(Debug);
00386       const int k = this->getNumIters(); 
00387       const int m = this->maxNumIters();
00388       TEUCHOS_TEST_FOR_EXCEPTION(k >= m, GmresCantExtendBasis,
00389        "Belos::CaGmres::extendBasis: "
00390        "Maximum number of iterations " << getNumIters() 
00391        << " reached; cannot extend basis further.");
00392       RCP<LinearProblem<Scalar, MV, OP> > lp = this->lp_;
00393       RCP<MV> V = this->V_;
00394       RCP<MV> Z = this->Z_;
00395       const int s = newCandidateBasisLength();
00396 
00397       // The most recently orthogonalized and accepted basis vector in
00398       // the Q / V basis.
00399       RCP<const MV> q_prv = MVT::CloneView (*V, Range1D (k, k));
00400       V_cur = MVT::CloneViewNonConst (*V, Range1D (k+1, k+s));
00401       if (this->flexible_)
00402   {
00403     Z_cur = MVT::CloneViewNonConst (*Z, Range1D (k, k+s-1));
00404     akx_->computeFlexibleBasis (q_prv, Z_cur, V_cur, s);
00405   }
00406       else
00407   {
00408     if (! Z_cur.is_null())
00409       Z_cur = null;
00410     akx_->computeBasis (q_prv, V_cur, s);
00411   }
00412       // In case we need to retry with a shorter candidate basis
00413       // length, remember what basis length we used this time.
00414       candidateBasisLength_ = s;
00415     }
00416 
00417     bool
00418     acceptedCandidateBasis () const 
00419     {
00420       return candidateBasisLength_ > 0 && 
00421   candidateBasisRank_ == candidateBasisLength_;
00422     }
00423 
00424     virtual void
00425     acceptCandidateBasis (const int newNumVectors)
00426     {
00427       (void) newNumVectors; // Ignore the input
00428 
00429       const int k = getNumIters(); 
00430       const int s = candidateBasisLength_;
00431       // In iteration 0, we want H(0:1, 0:0), which is 2 x 1.
00432       mat_type H_view (Teuchos::View, *(this->H_), k+2, k+1);
00433 
00434       // FIXME (mfh 28 Feb 2011) We should only do this a few times;
00435       // don't need to do this every outer iteration.
00436       //
00437       // FIXME (mfh 28 Feb 2011) In order to avoid problems with
00438       // heterogeneous nodes, we should really compute all these on a
00439       // single node and broadcast the result.
00440       updateEigenInfo (H_view);
00441       akx_->updateBasis (eigRealParts_, eigImagParts_, 
00442        eigMults_, eigNumUnique_);
00443 
00444       curNumIters_ += s;
00445       numInnerIters_.push_back (s);
00446       numOuterIters_++;
00447     }
00448 
00449     void
00450     updateEigenInfo ();
00451 
00452     virtual void
00453     rejectCandidateBasis (const int newNumVectors)
00454     {
00455       // Whether we should invoke advance() recursively.  The
00456       // alternative is to throw an exception signalling the inability
00457       // of CA-GMRES to make progress, given the linear system to
00458       // solve and the matrix powers kernel implementation.
00459       const bool shouldRecurse = 
00460   candidateBasisRank_ >= 2 && candidateBasisLength_ > 2;
00461 
00462       // It's not worth running CA-GMRES with a candidate basis length
00463       // of 1, so give up if that would happen on recursion.
00464       if (! shouldRecurse)
00465   {
00466     const std::ostringstream os;
00467     os << "You should either use standard GMRES for solving this "
00468       "problem, or try setting up the matrix powers kernel differently.";
00469     // Rank of the current candidate basis is < 2.
00470     TEUCHOS_TEST_FOR_EXCEPTION(lastCandidateBasisRank_ < 2, 
00471            GmresRejectsCandidateBasis,
00472            "CA-GMRES cannot produce a candidate basis of rank at "
00473            "least 2.  " << os.str());
00474     // Candidate basis length (not the "last candidate basis
00475     // length," which may be different) is <= 2.
00476     TEUCHOS_TEST_FOR_EXCEPTION(candidateBasisLength_ <= 2, 
00477            GmresRejectsCandidateBasis,
00478            "CA-GMRES: the current candidate basis length " 
00479            << candidateBasisLength_ << " is <= 2.  This means"
00480            " we cannot attempt recovery by recursing with a "
00481            "shorter candidate basis length." << os.str());
00482     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00483   }
00484 
00485       // Decrease the candidate basis length and try again.  It could
00486       // have been that lastCandidateBasisLength_ was quite a bit less
00487       // than candidateBasisLength_, if we were at the end of a
00488       // restart cycle and candidateBasisLength_ does not divide the
00489       // restart length evenly.  However, since the last candidate
00490       // basis was still not full rank, we should still use
00491       // lastCandidateBasisLength_ as a guide.
00492       //
00493       // The max is to ensure the second term of the min is at least 2,
00494       // for example if candidateBasisLength == 3.
00495       const int nextCandidateBasisLength = 
00496   std::min (lastCandidateBasisRank_, 
00497       std::max (2, candidateBasisLength_ / 2));
00498       // Throw an exception if a programming bug would otherwise have
00499       // caused an infinite loop.
00500       TEUCHOS_TEST_FOR_EXCEPTION(nextCandidateBasisLength < 2, std::logic_error, 
00501        "CA-GMRES should never get here!");
00502       candidateBasisLength_ = nextCandidateBasisLength;
00503 
00504       // Invoke advance() recursively.  The cost of recursion should
00505       // be minimal, since typical candidate basis lengths are short,
00506       // and advance() (by itself, not counting subclass
00507       // implementations of hooks, etc.) doesn't allocate any dense
00508       // vector or sparse matrix data.
00509       advance();
00510     }
00511 
00512     virtual void 
00513     updateUpperHessenbergMatrix (const Teuchos::RCP<mat_type>& C_Q,
00514          const Teuchos::RCP<mat_type>& B_Q,
00515          const Teuchos::RCP<mat_type>& C_Z,
00516          const Teuchos::RCP<mat_type>& B_Z)
00517     {
00518       using Teuchos::Copy, Teuchos::View, Teuchos::NO_TRANS, Teuchos::RCP;
00519       const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00520       const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00521 
00522       // Sanity checks
00523       TEUCHOS_TEST_FOR_EXCEPTION(C_Q.is_null(), std::logic_error, "C_Q is null");
00524       TEUCHOS_TEST_FOR_EXCEPTION(B_Q.is_null(), std::logic_error, "B_Q is null");
00525       TEUCHOS_TEST_FOR_EXCEPTION(flexible_ && C_Z.is_null(), std::logic_error, 
00526        "C_Z is null; not allowed when Flexible GMRES "
00527        "is being used");
00528       TEUCHOS_TEST_FOR_EXCEPTION(flexible_ && B_Z.is_null(), std::logic_error, 
00529        "B_Z is null; not allowed when Flexible GMRES "
00530        "is being used");
00531       // Even on the first outer iteration, C_Q should be non-null; it
00532       // should contain exactly one row.
00533       const int m = C_Q->numRows();
00534       const int s = C_Q->numCols();
00535       if (flexible_)
00536   {
00537     const int m_z = C_Z->numRows();
00538     const int s_z = C_Q->numCols();
00539     TEUCHOS_TEST_FOR_EXCEPTION(m_z != m || s_z != s-1, std::logic_error, 
00540            "C_Z has the wrong dimensions: C_Q is " << m 
00541            << " x " << s << " and so C_Z should be " << m 
00542            << " x " << (s-1) << ", but is " << m_z << " x " 
00543            << s_z << " instead.");
00544   }
00545       // Helper matrices make the algebra easier, though they take up
00546       // a modest amount of extra space.  Cleverness could make them
00547       // go away, but we are more interested in clarity at this point.
00548       mat_type F_Q (m-1, s+1);
00549       {
00550   if (m > 1)
00551     {
00552       mat_type F_Q_left (View, F_Q, m-1, 1, 0, 0);
00553       F_Q_left.putScalar (zero);
00554       mat_type F_Q_right (View, F_Q, m-1, s, 0, 1);
00555       F_Q_right = mat_type (View, *C_Q, m-1, s);
00556     }
00557       }
00558       mat_type F_Z (m-1, s);
00559       {
00560   if (m > 1)
00561     {
00562       mat_type F_Z_left (View, F_Z, m-1, 1, 0, 0);
00563       F_Z_left.putScalar (zero);
00564       mat_type F_Z_right (View, F_Z, m-1, s-1, 0, 1);
00565       if (flexible_)
00566         F_Z_right = mat_type (View, *C_Z, m-1, s-1);
00567       else
00568         F_Z_right = mat_type (View, *C_Q, m-1, s-1);
00569     }
00570       }
00571       mat_type R_Q (s+1, s+1);
00572       {
00573   R_Q(0,0) = one;
00574   mat_type R_Q_top (View, R_Q, 1, s, 0, 1);
00575   R_Q_top = mat_type (View, *C_Q, 1, s, m-1, 1);
00576   mat_type R_Q_bot (View, R_Q, s, s, 1, 1);
00577   R_Q_bot = *B_Q;
00578       }
00579       mat_type R_Z (s, s);
00580       {
00581   R_Z(0,0) = one;
00582   mat_type R_Z_top (View, R_Z, 1, s-1, 0, 1);
00583   if (flexible_)
00584     R_Z_top = mat_type (View, *C_Z, 1, s-1, m-1, 1);
00585   else
00586     R_Z_top = mat_type (View, *C_Q, 1, s-1, m-1, 1);
00587   mat_type R_Z_bot (View, R_Z, s-1, s-1, 1, 1);
00588   R_Z_bot = *B_Z;
00589       }
00590       Teuchos::BLAS<int,Scalar> blas;
00591       
00592       // E := B * inv(R_Z), an s+1 by s matrix, is a handy
00593       // intermediate quantity.  In the previous sentence, B is the
00594       // s+1 by s change-of-basis matrix (yes, naming it B is
00595       // confusing).
00596       mat_type E (s+1, s);
00597       {
00598   using Teuchos::RIGHT_SIDE, Teuchos::UPPER_TRI, Teuchos::NON_UNIT_DIAG;
00599   // The change-of-basis matrix is s+1 by s.  We confusingly also
00600   // name it B.  It is the same whether or not we are doing
00601   // Flexible GMRES.
00602   RCP<const mat_type> B = akx_->changeOfBasisMatrix (s);
00603   // BLAS' TRSM overwrites its second matrix argument, so copy
00604   // the right-hand side into E before solving.
00605   E = *B; 
00606   // E := B * inv(R_Z)
00607   blas.TRSM (RIGHT_SIDE, UPPER_TRI, NON_TRANS, NON_UNIT_DIAG, s+1, s,
00608        one, R_Z.values(), R_Z.stride(), E.values(), E.stride());
00609       }
00610       //
00611       // Compute H_kk, the lower right s+1 by s part of the upper
00612       // Hessenberg matrix.  Note that m-1 >= 0, even on the first
00613       // outer iteration (where m==1, since we still have to project
00614       // against the initial basis vector).
00615       //
00616       mat_type H_kk (View, *H_, s+1, s, m-1, m-1);
00617       // H_kk := R_Q * E (result is s+1 by s)
00618       blas.GEMM (NO_TRANS, NO_TRANS, s+1, s, s+1,
00619      one, R_Q.values(), R.stride(), 
00620      E.values(), E.stride(),
00621      zero, H_kk.values(), H_kk.stride());
00622       // If this is the first outer iteration, then we're done.
00623       if (getNumIters() == 0)
00624   return;
00625       // H_kk(1, 2:s) -= beta_km1 F_Z(m-1, 2:s)
00626       const Scalar beta_km1 = (*H_)(m-1, m-2);
00627       {
00628   mat_type H_kk_top (View, H_kk, 1, s-1, 0, 1);
00629   const mat_type F_Z_bot (View, F_Z, 1, s-1, m-2, 1);
00630   H_kk_top = F_Z_bot;
00631   H_kk_top *= -beta_km1; // elementwise multiplication
00632       }
00633       //
00634       // Compute H_km1_k
00635       //
00636       mat_type H_km1_k = mat_type (View, *H_, m-1, s, 0, m-1);
00637       // H_km1_k := -H_km1 * F_Z (result is m-1 by s)
00638       const mat_type H_km1 = mat_type (View, *H_, m-1, m-1);
00639       blas.GEMM (NO_TRANS, NO_TRANS, m-1, s, m-1, 
00640      -one, H_km1.values(), H_km1.stride(), 
00641      F_Z.values(), F_Z.stride(),
00642      zero, H_km1_k.values(), H_km1_k.stride());
00643       // H_km1_k += F_Q * E (result is m-1 by s)
00644       blas.GEMM (NO_TRANS, NO_TRANS, m-1, s, s+1,
00645      one, F_Q.values(), F.stride(), 
00646      E.values(), E.stride(),
00647      one, H_km1_k.values(), H_km1_k.stride());
00648     } 
00649 
00650 
00651   private:
00653     RCP<ParameterList> params_;
00654 
00656     RCP<ParameterList> akxParams_;
00657 
00680     Teuchos::RCP<Akx<Scalar, MV> > akx_;
00681 
00689     int candidateBasisLength_;
00690 
00692     int candidateBasisRank_;
00693 
00700     int numOuterIters_;
00701 
00703     int maxNumOuterIters_;
00704 
00714     std::vector<int> numInnerIters_;
00715   };
00716 } // namespace Belos
00717 
00718 #endif //  __Belos_CaGmres_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines