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