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