Anasazi Version of the Day
AnasaziRTRBase.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00033 #ifndef ANASAZI_RTRBASE_HPP
00034 #define ANASAZI_RTRBASE_HPP
00035 
00036 #include "AnasaziTypes.hpp"
00037 
00038 #include "AnasaziEigensolver.hpp"
00039 #include "AnasaziMultiVecTraits.hpp"
00040 #include "AnasaziOperatorTraits.hpp"
00041 #include "Teuchos_ScalarTraits.hpp"
00042 
00043 #include "AnasaziGenOrthoManager.hpp"
00044 #include "AnasaziSolverUtils.hpp"
00045 
00046 #include "Teuchos_LAPACK.hpp"
00047 #include "Teuchos_BLAS.hpp"
00048 #include "Teuchos_SerialDenseMatrix.hpp"
00049 #include "Teuchos_ParameterList.hpp"
00050 #include "Teuchos_TimeMonitor.hpp"
00051 
00101 namespace Anasazi {
00102 
00104 
00105 
00110   template <class ScalarType, class MV>
00111   struct RTRState {
00113     Teuchos::RCP<const MV> X; 
00115     Teuchos::RCP<const MV> AX; 
00117     Teuchos::RCP<const MV> BX;
00119     Teuchos::RCP<const MV> R;
00121     Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
00125     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType rho;
00126     RTRState() : X(Teuchos::null),AX(Teuchos::null),BX(Teuchos::null),
00127                  R(Teuchos::null),T(Teuchos::null),rho(0) {};
00128   };
00129 
00131 
00133 
00134 
00148   class RTRRitzFailure : public AnasaziError {public:
00149     RTRRitzFailure(const std::string& what_arg) : AnasaziError(what_arg)
00150     {}};
00151 
00160   class RTRInitFailure : public AnasaziError {public:
00161     RTRInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
00162     {}};
00163 
00180   class RTROrthoFailure : public AnasaziError {public:
00181     RTROrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
00182     {}};
00183 
00184 
00186 
00187 
00188   template <class ScalarType, class MV, class OP>
00189   class RTRBase : public Eigensolver<ScalarType,MV,OP> {
00190   public:
00191 
00193 
00194 
00200     RTRBase(const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 
00201             const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00202             const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00203             const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00204             const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
00205                   Teuchos::ParameterList &params,
00206             const std::string &solverLabel, bool skinnySolver
00207            );
00208 
00210     virtual ~RTRBase() {};
00211 
00213 
00215 
00216 
00238     virtual void iterate() = 0;
00239 
00264     void initialize(RTRState<ScalarType,MV> newstate);
00265 
00269     void initialize();
00270 
00283     bool isInitialized() const;
00284 
00292     RTRState<ScalarType,MV> getState() const;
00293 
00295 
00297 
00298 
00300     int getNumIters() const;
00301 
00303     void resetNumIters();
00304 
00312     Teuchos::RCP<const MV> getRitzVectors();
00313 
00319     std::vector<Value<ScalarType> > getRitzValues();
00320 
00328     std::vector<int> getRitzIndex();
00329 
00335     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
00336 
00337 
00343     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
00344 
00345 
00350     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
00351 
00352 
00361     int getCurSubspaceDim() const;
00362 
00365     int getMaxSubspaceDim() const;
00366 
00368 
00370 
00371 
00373     void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
00374 
00376     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
00377 
00379     const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
00380 
00381 
00390     void setBlockSize(int blockSize);
00391 
00392 
00394     int getBlockSize() const;
00395 
00396 
00417     void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
00418 
00420     Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
00421 
00423 
00425 
00426 
00428     virtual void currentStatus(std::ostream &os);
00429 
00431 
00432   protected:
00433     //
00434     // Convenience typedefs
00435     //
00436     typedef SolverUtils<ScalarType,MV,OP> Utils;
00437     typedef MultiVecTraits<ScalarType,MV> MVT;
00438     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00439     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00440     typedef typename SCT::magnitudeType MagnitudeType;
00441     typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
00442     const MagnitudeType ONE;  
00443     const MagnitudeType ZERO; 
00444     const MagnitudeType NANVAL;
00445     typedef typename std::vector<MagnitudeType>::iterator vecMTiter;
00446     typedef typename std::vector<ScalarType>::iterator    vecSTiter;
00447     //
00448     // Internal structs
00449     //
00450     struct CheckList {
00451       bool checkX, checkAX, checkBX;
00452       bool checkEta, checkAEta, checkBEta;
00453       bool checkR, checkQ, checkBR;
00454       bool checkZ, checkPBX;
00455       CheckList() : checkX(false),checkAX(false),checkBX(false),
00456                     checkEta(false),checkAEta(false),checkBEta(false),
00457                     checkR(false),checkQ(false),checkBR(false),
00458                     checkZ(false), checkPBX(false) {};
00459     };
00460     //
00461     // Internal methods
00462     //
00463     std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
00464     // solves the model minimization
00465     virtual void solveTRSubproblem() = 0;
00466     // Riemannian metric
00467     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi) const;
00468     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi, const MV &zeta) const;
00469     void ginnersep(const MV &xi, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
00470     void ginnersep(const MV &xi, const MV &zeta, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
00471     //
00472     // Classes input through constructor that define the eigenproblem to be solved.
00473     //
00474     const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> >     problem_;
00475     const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> >           sm_;
00476     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00477     Teuchos::RCP<StatusTest<ScalarType,MV,OP> >             tester_;
00478     const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> >  orthman_;
00479     //
00480     // Information obtained from the eigenproblem
00481     //
00482     Teuchos::RCP<const OP> AOp_;
00483     Teuchos::RCP<const OP> BOp_;
00484     Teuchos::RCP<const OP> Prec_;
00485     bool hasBOp_, hasPrec_, olsenPrec_;
00486     //
00487     // Internal timers
00488     //
00489     Teuchos::RCP<Teuchos::Time> timerAOp_, timerBOp_, timerPrec_,
00490                                 timerSort_, 
00491                                 timerLocalProj_, timerDS_,
00492                                 timerLocalUpdate_, timerCompRes_,
00493                                 timerOrtho_, timerInit_;
00494     //
00495     // Counters
00496     //
00497     // Number of operator applications
00498     int counterAOp_, counterBOp_, counterPrec_;
00499 
00500     //
00501     // Algorithmic parameters.
00502     //
00503     // blockSize_ is the solver block size
00504     int blockSize_;
00505     //
00506     // Current solver state
00507     //
00508     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00509     // is capable of running; _initialize is controlled  by the initialize() member method
00510     // For the implications of the state of initialized_, please see documentation for initialize()
00511     bool initialized_;
00512     //
00513     // nevLocal_ reflects how much of the current basis is valid (0 <= nevLocal_ <= blockSize_)
00514     // this tells us how many of the values in theta_ are valid Ritz values
00515     int nevLocal_;
00516     // 
00517     // are we implementing a skinny solver? (SkinnyIRTR)
00518     bool isSkinny_;
00519     // 
00520     // are we computing leftmost or rightmost eigenvalue?
00521     bool leftMost_;
00522     //
00523     // State Multivecs
00524     //
00525     // if we are implementing a skinny solver (SkinnyIRTR), 
00526     // then some of these will never be allocated
00527     // 
00528     // In order to handle auxiliary vectors, we need to handle the projector
00529     //   P_{[BQ BX],[BQ BX]}
00530     // Using an orthomanager with B-inner product, this requires calling with multivectors
00531     // [BQ,BX] and [Q,X].
00532     // These multivectors must be combined because <[BQ,BX],[Q,X]>_B != I
00533     // Hence, we will create two multivectors V and BV, which store
00534     //   V = [Q,X]
00535     //  BV = [BQ,BX]
00536     // 
00537     //  In the context of preconditioning, we may need to apply the projector
00538     //   P_{prec*[BQ,BX],[BQ,BX]}
00539     //  Because [BQ,BX] do not change during the supproblem solver, we will apply 
00540     //  the preconditioner to [BQ,BX] only once. This result is stored in PBV.
00541     // 
00542     // X,BX are views into V,BV
00543     // We don't need views for Q,BQ
00544     // Inside the subproblem solver, X,BX are constant, so we can allow these
00545     // views to exist alongside the full view of V,BV without violating
00546     // view semantics.
00547     // 
00548     // Skinny solver allocates 6/7/8 multivectors:
00549     //    V_, BV_ (if hasB)
00550     //    PBV_ (if hasPrec and olsenPrec)
00551     //    R_, Z_  (regardless of hasPrec)
00552     //    eta_, delta_, Hdelta_
00553     //
00554     // Hefty solver allocates 10/11/12/13 multivectors:
00555     //    V_, AX_, BV_ (if hasB)
00556     //    PBV_ (if hasPrec and olsenPrec)
00557     //    R_, Z_ (if hasPrec)
00558     //    eta_, Aeta_, Beta_
00559     //    delta_, Adelta_, Bdelta_, Hdelta_
00560     //
00561     Teuchos::RCP<MV> V_, BV_, PBV_,                     // V = [Q,X]; B*V; Prec*B*V
00562                      AX_, R_,                           // A*X_; residual, gradient, and residual of model minimization
00563                      eta_, Aeta_, Beta_,                // update vector from model minimization
00564                      delta_, Adelta_, Bdelta_, Hdelta_, // search direction in model minimization
00565                      Z_;                                // preconditioned residual
00566     Teuchos::RCP<const MV> X_, BX_;
00567     // 
00568     // auxiliary vectors
00569     Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
00570     int numAuxVecs_;
00571     //
00572     // Number of iterations that have been performed.
00573     int iter_;
00574     // 
00575     // Current eigenvalues, residual norms
00576     std::vector<MagnitudeType> theta_, Rnorms_, R2norms_, ritz2norms_;
00577     // 
00578     // are the residual norms current with the residual?
00579     bool Rnorms_current_, R2norms_current_;
00580     // 
00581     // parameters solver and inner solver
00582     MagnitudeType conv_kappa_, conv_theta_;
00583     MagnitudeType rho_;
00584     // 
00585     // current objective function value
00586     MagnitudeType fx_;
00587   };
00588 
00589 
00590 
00591 
00593   // Constructor
00594   template <class ScalarType, class MV, class OP>
00595   RTRBase<ScalarType,MV,OP>::RTRBase(
00596         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> >    &problem, 
00597         const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00598         const Teuchos::RCP<OutputManager<ScalarType> >         &printer,
00599         const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >      &tester,
00600         const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
00601         Teuchos::ParameterList                                 &params,
00602         const std::string &solverLabel, bool skinnySolver
00603         ) :
00604     ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
00605     ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00606     NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
00607     // problem, tools
00608     problem_(problem), 
00609     sm_(sorter),
00610     om_(printer),
00611     tester_(tester),
00612     orthman_(ortho),
00613     // timers, counters
00614     timerAOp_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation A*x")),
00615     timerBOp_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation B*x")),
00616     timerPrec_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation Prec*x")),
00617     timerSort_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Sorting eigenvalues")),
00618     timerLocalProj_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Local projection")),
00619     timerDS_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Direct solve")),
00620     timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Local update")),
00621     timerCompRes_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Computing residuals")),
00622     timerOrtho_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Orthogonalization")),
00623     timerInit_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Initialization")),
00624     counterAOp_(0),
00625     counterBOp_(0),
00626     counterPrec_(0),
00627     // internal data
00628     blockSize_(0),
00629     initialized_(false),
00630     nevLocal_(0),
00631     isSkinny_(skinnySolver),
00632     leftMost_(true),
00633     auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ), 
00634     numAuxVecs_(0),
00635     iter_(0),
00636     Rnorms_current_(false),
00637     R2norms_current_(false),
00638     conv_kappa_(.1), 
00639     conv_theta_(1),
00640     rho_( MAT::nan() ),
00641     fx_( MAT::nan() )
00642   {
00643     TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
00644         "Anasazi::RTRBase::constructor: user passed null problem pointer.");
00645     TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
00646         "Anasazi::RTRBase::constructor: user passed null sort manager pointer.");
00647     TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
00648         "Anasazi::RTRBase::constructor: user passed null output manager pointer.");
00649     TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
00650         "Anasazi::RTRBase::constructor: user passed null status test pointer.");
00651     TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
00652         "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer.");
00653     TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
00654         "Anasazi::RTRBase::constructor: problem is not set.");
00655     TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
00656         "Anasazi::RTRBase::constructor: problem is not Hermitian.");
00657 
00658     // get the problem operators
00659     AOp_   = problem_->getOperator();
00660     TEST_FOR_EXCEPTION(AOp_ == Teuchos::null, std::invalid_argument,
00661                        "Anasazi::RTRBase::constructor: problem provides no A matrix.");
00662     BOp_  = problem_->getM();
00663     Prec_ = problem_->getPrec();
00664     hasBOp_ = (BOp_ != Teuchos::null);
00665     hasPrec_ = (Prec_ != Teuchos::null);
00666     olsenPrec_ = params.get<bool>("Olsen Prec", true);
00667 
00668     TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument,
00669         "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product.");
00670 
00671     // set the block size and allocate data
00672     int bs = params.get("Block Size", problem_->getNEV());
00673     setBlockSize(bs);
00674 
00675     // leftmost or rightmost?
00676     leftMost_ = params.get("Leftmost",leftMost_);
00677 
00678     conv_kappa_ = params.get("Kappa Convergence",conv_kappa_);
00679     TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument,
00680                        "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
00681     conv_theta_ = params.get("Theta Convergence",conv_theta_);
00682     TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument,
00683                        "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
00684   }
00685 
00686 
00688   // Set the block size and make necessary adjustments.
00689   template <class ScalarType, class MV, class OP>
00690   void RTRBase<ScalarType,MV,OP>::setBlockSize (int blockSize) 
00691   {
00692     // time spent here counts towards timerInit_
00693     Teuchos::TimeMonitor lcltimer( *timerInit_ );
00694 
00695     // This routine only allocates space; it doesn't not perform any computation
00696     // if solver is initialized and size is to be decreased, take the first blockSize vectors of all to preserve state
00697     // otherwise, shrink/grow/allocate space and set solver to unitialized
00698 
00699     Teuchos::RCP<const MV> tmp;
00700     // grab some Multivector to Clone
00701     // in practice, getInitVec() should always provide this, but it is possible to use a 
00702     // Eigenproblem with nothing in getInitVec() by manually initializing with initialize(); 
00703     // in case of that strange scenario, we will try to Clone from R_
00704     // we like R_ for this, because it has minimal size (blockSize_), as opposed to V_ (blockSize_+numAuxVecs_)
00705     if (blockSize_ > 0) {
00706       tmp = R_;
00707     }
00708     else {
00709       tmp = problem_->getInitVec();
00710       TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
00711           "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
00712     }
00713 
00714     TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetVecLength(*tmp), std::invalid_argument, 
00715         "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
00716 
00717     // last chance to quit before causing side-effects
00718     if (blockSize == blockSize_) {
00719       // do nothing
00720       return;
00721     }
00722 
00723     // clear views
00724     X_  = Teuchos::null;
00725     BX_ = Teuchos::null;
00726 
00727     // regardless of whether we preserve any data, any content in V, BV and PBV corresponding to the
00728     // auxiliary vectors must be retained
00729     // go ahead and do these first
00730     // 
00731     // two cases here: 
00732     // * we are growing (possibly, from empty) 
00733     //   any aux data must be copied over, nothing from X need copying
00734     // * we are shrinking
00735     //   any aux data must be copied over, go ahead and copy X material (initialized or not)
00736     //
00737     if (blockSize > blockSize_)
00738     {
00739       // GROWING 
00740       // get a pointer for Q-related material, and an index for extracting/setting it
00741       Teuchos::RCP<const MV> Q;
00742       std::vector<int> indQ(numAuxVecs_);
00743       for (int i=0; i<numAuxVecs_; i++) indQ[i] = i;
00744       // if numAuxVecs_ > 0, then necessarily blockSize_ > 0 (we have already been allocated once)
00745       TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error,
00746           "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team.");
00747       // V
00748       if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ);
00749       V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00750       if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_);
00751       // BV
00752       if (hasBOp_) {
00753         if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ);
00754         BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00755         if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_);
00756       }
00757       else {
00758         BV_ = V_;
00759       }
00760       // PBV
00761       if (hasPrec_ && olsenPrec_) {
00762         if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ);
00763         PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00764         if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_);
00765       }
00766     }
00767     else 
00768     {
00769       // SHRINKING
00770       std::vector<int> indV(numAuxVecs_+blockSize);
00771       for (int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
00772       // V
00773       V_ = MVT::CloneCopy(*V_,indV);
00774       // BV
00775       if (hasBOp_) {
00776         BV_ = MVT::CloneCopy(*BV_,indV);
00777       }
00778       else {
00779         BV_ = V_;
00780       }
00781       // PBV
00782       if (hasPrec_ && olsenPrec_) {
00783         PBV_ = MVT::CloneCopy(*PBV_,indV);
00784       }
00785     }
00786 
00787     if (blockSize < blockSize_) {
00788       // shrink vectors
00789       blockSize_ = blockSize;
00790 
00791       theta_.resize(blockSize_);
00792       ritz2norms_.resize(blockSize_);
00793       Rnorms_.resize(blockSize_);
00794       R2norms_.resize(blockSize_);
00795 
00796       if (initialized_) {
00797         // shrink multivectors with copy
00798         std::vector<int> ind(blockSize_);
00799         for (int i=0; i<blockSize_; i++) ind[i] = i;
00800 
00801         // Z can be deleted, no important info there
00802         Z_ = Teuchos::null;
00803         
00804         // we will not use tmp for cloning, clear it and free some space
00805         tmp = Teuchos::null;
00806 
00807         R_      = MVT::CloneCopy(*R_     ,ind);
00808         eta_    = MVT::CloneCopy(*eta_   ,ind);
00809         delta_  = MVT::CloneCopy(*delta_ ,ind);
00810         Hdelta_ = MVT::CloneCopy(*Hdelta_,ind);
00811         if (!isSkinny_) {
00812           AX_     = MVT::CloneCopy(*AX_    ,ind);
00813           Aeta_   = MVT::CloneCopy(*Aeta_  ,ind);
00814           Adelta_ = MVT::CloneCopy(*Adelta_,ind);
00815         }
00816         else {
00817           AX_     = Teuchos::null;
00818           Aeta_   = Teuchos::null;
00819           Adelta_ = Teuchos::null;
00820         }
00821 
00822         if (hasBOp_) {
00823           if (!isSkinny_) {
00824             Beta_   = MVT::CloneCopy(*Beta_,ind);
00825             Bdelta_ = MVT::CloneCopy(*Bdelta_,ind);
00826           }
00827           else {
00828             Beta_   = Teuchos::null;
00829             Bdelta_ = Teuchos::null;
00830           }
00831         }
00832         else {
00833           Beta_   = eta_;
00834           Bdelta_ = delta_;
00835         }
00836         
00837         // we need Z if we have a preconditioner
00838         // we also use Z for temp storage in the skinny solvers
00839         if (hasPrec_ || isSkinny_) {
00840           Z_ = MVT::Clone(*V_,blockSize_);
00841         }
00842         else {
00843           Z_ = R_;
00844         }
00845 
00846       }
00847       else {
00848         // release previous multivectors
00849         // shrink multivectors without copying
00850         AX_     = Teuchos::null;
00851         R_      = Teuchos::null;
00852         eta_    = Teuchos::null;
00853         Aeta_   = Teuchos::null;
00854         delta_  = Teuchos::null;
00855         Adelta_ = Teuchos::null;
00856         Hdelta_ = Teuchos::null;
00857         Beta_   = Teuchos::null;
00858         Bdelta_ = Teuchos::null;
00859         Z_      = Teuchos::null;
00860 
00861         R_      = MVT::Clone(*tmp,blockSize_);
00862         eta_    = MVT::Clone(*tmp,blockSize_);
00863         delta_  = MVT::Clone(*tmp,blockSize_);
00864         Hdelta_ = MVT::Clone(*tmp,blockSize_);
00865         if (!isSkinny_) {
00866           AX_     = MVT::Clone(*tmp,blockSize_);
00867           Aeta_   = MVT::Clone(*tmp,blockSize_);
00868           Adelta_ = MVT::Clone(*tmp,blockSize_);
00869         }
00870 
00871         if (hasBOp_) {
00872           if (!isSkinny_) {
00873             Beta_   = MVT::Clone(*tmp,blockSize_);
00874             Bdelta_ = MVT::Clone(*tmp,blockSize_);
00875           }
00876         }
00877         else {
00878           Beta_   = eta_;
00879           Bdelta_ = delta_;
00880         }
00881 
00882         // we need Z if we have a preconditioner
00883         // we also use Z for temp storage in the skinny solvers
00884         if (hasPrec_ || isSkinny_) {
00885           Z_ = MVT::Clone(*tmp,blockSize_);
00886         }
00887         else {
00888           Z_ = R_;
00889         }
00890       } // if initialized_
00891     } // if blockSize is shrinking
00892     else {  // if blockSize > blockSize_
00893       // this is also the scenario for our initial call to setBlockSize(), in the constructor
00894       initialized_ = false;
00895 
00896       // grow/allocate vectors
00897       theta_.resize(blockSize,NANVAL);
00898       ritz2norms_.resize(blockSize,NANVAL);
00899       Rnorms_.resize(blockSize,NANVAL);
00900       R2norms_.resize(blockSize,NANVAL);
00901 
00902       // deallocate old multivectors
00903       AX_     = Teuchos::null;
00904       R_      = Teuchos::null;
00905       eta_    = Teuchos::null;
00906       Aeta_   = Teuchos::null;
00907       delta_  = Teuchos::null;
00908       Adelta_ = Teuchos::null;
00909       Hdelta_ = Teuchos::null;
00910       Beta_   = Teuchos::null;
00911       Bdelta_ = Teuchos::null;
00912       Z_      = Teuchos::null;
00913 
00914       // clone multivectors off of tmp
00915       R_      = MVT::Clone(*tmp,blockSize);
00916       eta_    = MVT::Clone(*tmp,blockSize);
00917       delta_  = MVT::Clone(*tmp,blockSize);
00918       Hdelta_ = MVT::Clone(*tmp,blockSize);
00919       if (!isSkinny_) {
00920         AX_     = MVT::Clone(*tmp,blockSize);
00921         Aeta_   = MVT::Clone(*tmp,blockSize);
00922         Adelta_ = MVT::Clone(*tmp,blockSize);
00923       }
00924 
00925       if (hasBOp_) {
00926         if (!isSkinny_) {
00927           Beta_   = MVT::Clone(*tmp,blockSize);
00928           Bdelta_ = MVT::Clone(*tmp,blockSize);
00929         }
00930       }
00931       else {
00932         Beta_   = eta_;
00933         Bdelta_ = delta_;
00934       }
00935       if (hasPrec_ || isSkinny_) {
00936         Z_ = MVT::Clone(*tmp,blockSize);
00937       }
00938       else {
00939         Z_ = R_;
00940       }
00941       blockSize_ = blockSize;
00942     }
00943 
00944     // get view of X from V, BX from BV
00945     // these are located after the first numAuxVecs columns
00946     {
00947       std::vector<int> indX(blockSize_);
00948       for (int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
00949       X_ = MVT::CloneView(*V_,indX);
00950       if (hasBOp_) {
00951         BX_ = MVT::CloneView(*BV_,indX);
00952       }
00953       else {
00954         BX_ = X_;
00955       }
00956     }
00957   }
00958 
00959 
00961   // Set a new StatusTest for the solver.
00962   template <class ScalarType, class MV, class OP>
00963   void RTRBase<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
00964     TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
00965         "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest.");
00966     tester_ = test;
00967   }
00968 
00969 
00971   // Get the current StatusTest used by the solver.
00972   template <class ScalarType, class MV, class OP>
00973   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > RTRBase<ScalarType,MV,OP>::getStatusTest() const {
00974     return tester_;
00975   }
00976   
00977 
00979   // Set the auxiliary vectors
00980   template <class ScalarType, class MV, class OP>
00981   void RTRBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00982     typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv;
00983 
00984     // set new auxiliary vectors
00985     auxVecs_.resize(0);
00986     auxVecs_.reserve(auxvecs.size());
00987 
00988     numAuxVecs_ = 0;
00989     for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
00990       numAuxVecs_ += MVT::GetNumberVecs(**v);
00991     }
00992 
00993     // If the solver has been initialized, X is not necessarily orthogonal to new auxiliary vectors
00994     if (numAuxVecs_ > 0 && initialized_) {
00995       initialized_ = false;
00996     }
00997 
00998     // clear X,BX views
00999     X_   = Teuchos::null;
01000     BX_  = Teuchos::null;
01001     // deallocate, we'll clone off R_ below
01002     V_   = Teuchos::null;
01003     BV_  = Teuchos::null;
01004     PBV_ = Teuchos::null;
01005 
01006     // put auxvecs into V, update BV and PBV if necessary
01007     if (numAuxVecs_ > 0) {
01008       V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01009       int numsofar = 0;
01010       for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
01011         std::vector<int> ind(MVT::GetNumberVecs(**v));
01012         for (size_t j=0; j<ind.size(); j++) ind[j] = numsofar++;
01013         MVT::SetBlock(**v,ind,*V_);
01014         auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind));
01015       }
01016       TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error,
01017           "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team.");
01018       // compute B*V, Prec*B*V
01019       if (hasBOp_) {
01020         BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01021         OPT::Apply(*BOp_,*V_,*BV_);
01022       }
01023       else {
01024         BV_ = V_;
01025       }
01026       if (hasPrec_ && olsenPrec_) {
01027         PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01028         OPT::Apply(*Prec_,*BV_,*V_);
01029       }
01030     }
01031     // 
01032 
01033     if (om_->isVerbosity( Debug ) ) {
01034       // Check almost everything here
01035       CheckList chk;
01036       chk.checkQ = true;
01037       om_->print( Debug, accuracyCheck(chk, "in setAuxVecs()") );
01038     }
01039   }
01040 
01041 
01043   /* Initialize the state of the solver
01044    * 
01045    * POST-CONDITIONS:
01046    *
01047    * initialized_ == true
01048    * X is orthonormal, orthogonal to auxVecs_
01049    * AX = A*X if not skinnny
01050    * BX = B*X if hasBOp_
01051    * theta_ contains Ritz values of X
01052    * R = AX - BX*diag(theta_)
01053    */
01054   template <class ScalarType, class MV, class OP>
01055   void RTRBase<ScalarType,MV,OP>::initialize(RTRState<ScalarType,MV> newstate)
01056   {
01057     // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
01058     // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
01059 
01060     // clear const views to X,BX in V,BV
01061     // work with temporary non-const views
01062     X_  = Teuchos::null;
01063     BX_ = Teuchos::null;
01064     Teuchos::RCP<MV> X, BX;
01065     {
01066       std::vector<int> ind(blockSize_);
01067       for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
01068       X = MVT::CloneViewNonConst(*V_,ind);
01069       if (hasBOp_) {
01070         BX = MVT::CloneViewNonConst(*BV_,ind);
01071       }
01072       else {
01073         BX = X;
01074       }
01075     }
01076 
01077     Teuchos::TimeMonitor inittimer( *timerInit_ );
01078 
01079     std::vector<int> bsind(blockSize_);
01080     for (int i=0; i<blockSize_; i++) bsind[i] = i;
01081 
01082     // in RTR, X (the subspace iterate) is primary
01083     // the order of dependence follows like so.
01084     // --init->                 X
01085     //    --op apply->          AX,BX
01086     //       --ritz analysis->  theta
01087     // 
01088     // if the user specifies all data for a level, we will accept it.
01089     // otherwise, we will generate the whole level, and all subsequent levels.
01090     //
01091     // the data members are ordered based on dependence, and the levels are
01092     // partitioned according to the amount of work required to produce the
01093     // items in a level.
01094     //
01095     // inconsitent multivectors widths and lengths will not be tolerated, and
01096     // will be treated with exceptions.
01097 
01098     // set up X, AX, BX: get them from "state" if user specified them
01099     if (newstate.X != Teuchos::null) {
01100       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.X) != MVT::GetVecLength(*X),
01101                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
01102       // newstate.X must have blockSize_ vectors; any more will be ignored
01103       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
01104                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
01105 
01106       // put data in X
01107       MVT::SetBlock(*newstate.X,bsind,*X);
01108 
01109       // put data in AX
01110       // if we are implementing a skinny solver, then we don't have memory allocated for AX
01111       // in this case, point AX at Z (skinny solvers allocate Z) and use it for temporary storage
01112       // we will clear the pointer later
01113       if (isSkinny_) {
01114         AX_ = Z_;
01115       }
01116       if (newstate.AX != Teuchos::null) {
01117         TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.AX) != MVT::GetVecLength(*X),
01118                             std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
01119         // newstate.AX must have blockSize_ vectors; any more will be ignored
01120         TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.AX) < blockSize_,
01121                             std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
01122         MVT::SetBlock(*newstate.AX,bsind,*AX_);
01123       }
01124       else {
01125         {
01126           Teuchos::TimeMonitor lcltimer( *timerAOp_ );
01127           OPT::Apply(*AOp_,*X,*AX_);
01128           counterAOp_ += blockSize_;
01129         }
01130         // we generated AX; we will generate R as well
01131         newstate.R = Teuchos::null;
01132       }
01133 
01134       // put data in BX
01135       // skinny solvers always allocate BX if hasB, so this is unconditionally appropriate
01136       if (hasBOp_) {
01137         if (newstate.BX != Teuchos::null) {
01138           TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.BX) != MVT::GetVecLength(*X),
01139                               std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
01140           // newstate.BX must have blockSize_ vectors; any more will be ignored
01141           TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.BX) < blockSize_,
01142                               std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
01143           MVT::SetBlock(*newstate.BX,bsind,*BX);
01144         }
01145         else {
01146           {
01147             Teuchos::TimeMonitor lcltimer( *timerBOp_ );
01148             OPT::Apply(*BOp_,*X,*BX);
01149             counterBOp_ += blockSize_;
01150           }
01151           // we generated BX; we will generate R as well
01152           newstate.R = Teuchos::null;
01153         }
01154       }
01155       else {
01156         // the assignment BX_==X_ would be redundant; take advantage of this opportunity to debug a little
01157         TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
01158       }
01159 
01160     }
01161     else {
01162       // user did not specify X
01163 
01164       // clear state so we won't use any data from it below
01165       newstate.R  = Teuchos::null;
01166       newstate.T  = Teuchos::null;
01167 
01168       // generate something and projectAndNormalize
01169       Teuchos::RCP<const MV> ivec = problem_->getInitVec();
01170       TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
01171                          "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from.");
01172 
01173       int initSize = MVT::GetNumberVecs(*ivec);
01174       if (initSize > blockSize_) {
01175         // we need only the first blockSize_ vectors from ivec; get a view of them
01176         initSize = blockSize_;
01177         std::vector<int> ind(blockSize_);
01178         for (int i=0; i<blockSize_; i++) ind[i] = i;
01179         ivec = MVT::CloneView(*ivec,ind);
01180       }
01181 
01182       // assign ivec to first part of X
01183       if (initSize > 0) {
01184         std::vector<int> ind(initSize);
01185         for (int i=0; i<initSize; i++) ind[i] = i;
01186         MVT::SetBlock(*ivec,ind,*X);
01187       }
01188       // fill the rest of X with random
01189       if (blockSize_ > initSize) {
01190         std::vector<int> ind(blockSize_ - initSize);
01191         for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
01192         Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X,ind);
01193         MVT::MvRandom(*rX);
01194         rX = Teuchos::null;
01195       }
01196 
01197       // put data in BX
01198       if (hasBOp_) {
01199         Teuchos::TimeMonitor lcltimer( *timerBOp_ );
01200         OPT::Apply(*BOp_,*X,*BX);
01201         counterBOp_ += blockSize_;
01202       }
01203       else {
01204         // the assignment BX==X would be redundant; take advantage of this opportunity to debug a little
01205         TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
01206       }
01207   
01208       // remove auxVecs from X and normalize it
01209       if (numAuxVecs_ > 0) {
01210         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01211         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01212         int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX);
01213         TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
01214                            "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
01215       }
01216       else {
01217         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01218         int rank = orthman_->normalizeMat(*X,Teuchos::null,BX);
01219         TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
01220                            "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
01221       }
01222 
01223       // put data in AX
01224       if (isSkinny_) {
01225         AX_ = Z_;
01226       }
01227       {
01228         Teuchos::TimeMonitor lcltimer( *timerAOp_ );
01229         OPT::Apply(*AOp_,*X,*AX_);
01230         counterAOp_ += blockSize_;
01231       }
01232 
01233     } // end if (newstate.X != Teuchos::null)
01234 
01235 
01236     // set up Ritz values
01237     if (newstate.T != Teuchos::null) {
01238       TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_,
01239                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values.");
01240       for (int i=0; i<blockSize_; i++) {
01241         theta_[i] = (*newstate.T)[i];
01242       }
01243     }
01244     else {
01245       // get ritz vecs/vals
01246       Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_),
01247                                                  BB(blockSize_,blockSize_),
01248                                                   S(blockSize_,blockSize_);
01249       // project A
01250       {
01251         Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
01252         MVT::MvTransMv(ONE,*X,*AX_,AA);
01253         if (hasBOp_) {
01254           MVT::MvTransMv(ONE,*X,*BX,BB);
01255         }
01256       }
01257       nevLocal_ = blockSize_;
01258 
01259       // solve the projected problem
01260       // project B
01261       int ret;
01262       if (hasBOp_) {
01263         Teuchos::TimeMonitor lcltimer( *timerDS_ );
01264         ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
01265       }
01266       else {
01267         Teuchos::TimeMonitor lcltimer( *timerDS_ );
01268         ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10);
01269       }
01270       TEST_FOR_EXCEPTION(ret != 0,RTRInitFailure,
01271           "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret);
01272       TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,RTRInitFailure,"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis.");
01273 
01274       // We only have blockSize_ ritz pairs, ergo we do not need to select.
01275       // However, we still require them to be ordered correctly
01276       {
01277         Teuchos::TimeMonitor lcltimer( *timerSort_ );
01278         
01279         std::vector<int> order(blockSize_);
01280         // 
01281         // sort the first blockSize_ values in theta_
01282         sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);   // don't catch exception
01283         //
01284         // apply the same ordering to the primitive ritz vectors
01285         Utils::permuteVectors(order,S);
01286       }
01287 
01288       // compute ritz residual norms
01289       {
01290         Teuchos::BLAS<int,ScalarType> blas;
01291         Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_);
01292         // RR = AA*S - BB*S*diag(theta)
01293         int info;
01294         if (hasBOp_) {
01295           info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO);
01296           TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
01297         }
01298         else {
01299           RR.assign(S);
01300         }
01301         for (int i=0; i<blockSize_; i++) {
01302           blas.SCAL(blockSize_,theta_[i],RR[i],1);
01303         }
01304         info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE);
01305         TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
01306         for (int i=0; i<blockSize_; i++) {
01307           ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1);
01308         }
01309       }
01310 
01311       // update the solution
01312       // use R as local work space
01313       // Z may already be in use as work space (holding AX if isSkinny)
01314       {
01315         Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
01316         // X <- X*S
01317         MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ );        
01318         MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X );
01319         // AX <- AX*S
01320         MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ );        
01321         MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ );
01322         if (hasBOp_) {
01323           // BX <- BX*S
01324           MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ );        
01325           MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX );
01326         }
01327       }
01328     }
01329     
01330     // done modifying X,BX; clear temp views and setup const views
01331     X  = Teuchos::null;
01332     BX = Teuchos::null;
01333     {
01334       std::vector<int> ind(blockSize_);
01335       for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
01336       this->X_ = MVT::CloneView(static_cast<const MV&>(*V_),ind);
01337       if (hasBOp_) {
01338         this->BX_ = MVT::CloneView(static_cast<const MV&>(*BV_),ind);
01339       }
01340       else {
01341         this->BX_ = this->X_;
01342       }
01343     }
01344 
01345 
01346     // get objective function value
01347     fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
01348 
01349     // set up R
01350     if (newstate.R != Teuchos::null) {
01351       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
01352                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
01353       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
01354                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
01355       MVT::SetBlock(*newstate.R,bsind,*R_);
01356     }
01357     else {
01358       Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
01359       // form R <- AX - BX*T
01360       MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
01361       Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
01362       T.putScalar(ZERO);
01363       for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
01364       MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_);
01365     }
01366 
01367     // R has been updated; mark the norms as out-of-date
01368     Rnorms_current_ = false;
01369     R2norms_current_ = false;
01370 
01371     // if isSkinny, then AX currently points to Z for temp storage
01372     // set AX back to null
01373     if (isSkinny_) {
01374       AX_ = Teuchos::null;
01375     }
01376 
01377     // finally, we are initialized
01378     initialized_ = true;
01379 
01380     if (om_->isVerbosity( Debug ) ) {
01381       // Check almost everything here
01382       CheckList chk;
01383       chk.checkX = true;
01384       chk.checkAX = true;
01385       chk.checkBX = true;
01386       chk.checkR = true;
01387       chk.checkQ = true;
01388       om_->print( Debug, accuracyCheck(chk, "after initialize()") );
01389     }
01390   }
01391 
01392   template <class ScalarType, class MV, class OP>
01393   void RTRBase<ScalarType,MV,OP>::initialize()
01394   {
01395     RTRState<ScalarType,MV> empty;
01396     initialize(empty);
01397   }
01398 
01399 
01400 
01401 
01403   // compute/return residual B-norms
01404   template <class ScalarType, class MV, class OP>
01405   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01406   RTRBase<ScalarType,MV,OP>::getResNorms() {
01407     if (Rnorms_current_ == false) {
01408       // Update the residual norms
01409       orthman_->norm(*R_,Rnorms_);
01410       Rnorms_current_ = true;
01411     }
01412     return Rnorms_;
01413   }
01414 
01415   
01417   // compute/return residual 2-norms
01418   template <class ScalarType, class MV, class OP>
01419   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01420   RTRBase<ScalarType,MV,OP>::getRes2Norms() {
01421     if (R2norms_current_ == false) {
01422       // Update the residual 2-norms 
01423       MVT::MvNorm(*R_,R2norms_);
01424       R2norms_current_ = true;
01425     }
01426     return R2norms_;
01427   }
01428 
01429 
01430 
01431 
01433   // Check accuracy, orthogonality, and other debugging stuff
01434   // 
01435   // bools specify which tests we want to run (instead of running more than we actually care about)
01436   //
01437   // we don't bother checking the following because they are computed explicitly:
01438   //   AH == A*H
01439   //
01440   // 
01441   // checkX    : X orthonormal
01442   //             orthogonal to auxvecs
01443   // checkAX   : check AX == A*X
01444   // checkBX   : check BX == B*X
01445   // checkEta  : check that Eta'*B*X == 0
01446   //             orthogonal to auxvecs
01447   // checkAEta : check that AEta = A*Eta
01448   // checkBEta : check that BEta = B*Eta
01449   // checkR    : check R orthogonal to X
01450   // checkBR   : check R B-orthogonal to X, valid in and immediately after solveTRSubproblem
01451   // checkQ    : check that auxiliary vectors are actually orthonormal
01452   //
01453   // TODO: 
01454   //  add checkTheta 
01455   //
01456   template <class ScalarType, class MV, class OP>
01457   std::string RTRBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const 
01458   {
01459     using std::setprecision;
01460     using std::scientific;
01461     using std::setw;
01462     using std::endl;
01463     std::stringstream os;
01464     MagnitudeType tmp;
01465 
01466     os << " Debugging checks: " << where << endl;
01467 
01468     // X and friends
01469     if (chk.checkX && initialized_) {
01470       tmp = orthman_->orthonormError(*X_);
01471       os << " >> Error in X^H B X == I :    " << scientific << setprecision(10) << tmp << endl;
01472       for (Array_size_type i=0; i<auxVecs_.size(); i++) {
01473         tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
01474         os << " >> Error in X^H B Q[" << i << "] == 0 : " << scientific << setprecision(10) << tmp << endl;
01475       }
01476     }
01477     if (chk.checkAX && !isSkinny_ && initialized_) {
01478       tmp = Utils::errorEquality(*X_, *AX_, AOp_);
01479       os << " >> Error in AX == A*X    :    " << scientific << setprecision(10) << tmp << endl;
01480     }
01481     if (chk.checkBX && hasBOp_ && initialized_) {
01482       Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_);
01483       OPT::Apply(*BOp_,*this->X_,*BX);
01484       tmp = Utils::errorEquality(*BX, *BX_);
01485       os << " >> Error in BX == B*X    :    " << scientific << setprecision(10) << tmp << endl;
01486     }
01487 
01488     // Eta and friends
01489     if (chk.checkEta && initialized_) {
01490       tmp = orthman_->orthogError(*X_,*eta_);
01491       os << " >> Error in X^H B Eta == 0 :  " << scientific << setprecision(10) << tmp << endl;
01492       for (Array_size_type i=0; i<auxVecs_.size(); i++) {
01493         tmp = orthman_->orthogError(*eta_,*auxVecs_[i]);
01494         os << " >> Error in Eta^H B Q[" << i << "]==0 : " << scientific << setprecision(10) << tmp << endl;
01495       }
01496     }
01497     if (chk.checkAEta && !isSkinny_ && initialized_) {
01498       tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_);
01499       os << " >> Error in AEta == A*Eta    :    " << scientific << setprecision(10) << tmp << endl;
01500     }
01501     if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
01502       tmp = Utils::errorEquality(*eta_, *Beta_, BOp_);
01503       os << " >> Error in BEta == B*Eta    :    " << scientific << setprecision(10) << tmp << endl;
01504     }
01505 
01506     // R: this is not B-orthogonality, but standard euclidean orthogonality
01507     if (chk.checkR && initialized_) {
01508       Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01509       MVT::MvTransMv(ONE,*X_,*R_,xTx);
01510       tmp = xTx.normFrobenius();
01511       os << " >> Error in X^H R == 0   :    " << scientific << setprecision(10) << tmp << endl;
01512     }
01513 
01514     // BR: this is B-orthogonality: this is only valid inside and immediately after solveTRSubproblem 
01515     // only check if B != I, otherwise, it is equivalent to the above test
01516     if (chk.checkBR && hasBOp_ && initialized_) {
01517       Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01518       MVT::MvTransMv(ONE,*BX_,*R_,xTx);
01519       tmp = xTx.normFrobenius();
01520       os << " >> Error in X^H B R == 0 :    " << scientific << setprecision(10) << tmp << endl;
01521     }
01522 
01523     // Z: Z is preconditioned R, should be on tangent plane
01524     if (chk.checkZ && initialized_) {
01525       Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01526       MVT::MvTransMv(ONE,*BX_,*Z_,xTx);
01527       tmp = xTx.normFrobenius();
01528       os << " >> Error in X^H B Z == 0 :    " << scientific << setprecision(10) << tmp << endl;
01529     }
01530 
01531     // Q
01532     if (chk.checkQ) {
01533       for (Array_size_type i=0; i<auxVecs_.size(); i++) {
01534         tmp = orthman_->orthonormError(*auxVecs_[i]);
01535         os << " >> Error in Q[" << i << "]^H B Q[" << i << "]==I: " << scientific << setprecision(10) << tmp << endl;
01536         for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
01537           tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01538           os << " >> Error in Q[" << i << "]^H B Q[" << j << "]==0: " << scientific << setprecision(10) << tmp << endl;
01539         }
01540       }
01541     }
01542     os << endl;
01543     return os.str();
01544   }
01545 
01546 
01548   // Print the current status of the solver
01549   template <class ScalarType, class MV, class OP>
01550   void 
01551   RTRBase<ScalarType,MV,OP>::currentStatus(std::ostream &os) 
01552   {
01553     using std::setprecision;
01554     using std::scientific;
01555     using std::setw;
01556     using std::endl;
01557 
01558     os <<endl;
01559     os <<"================================================================================" << endl;
01560     os << endl;
01561     os <<"                              RTRBase Solver Status" << endl;
01562     os << endl;
01563     os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01564     os <<"The number of iterations performed is " << iter_       << endl;
01565     os <<"The current block size is             " << blockSize_  << endl;
01566     os <<"The number of auxiliary vectors is    " << numAuxVecs_ << endl;
01567     os <<"The number of operations A*x    is " << counterAOp_   << endl;
01568     os <<"The number of operations B*x    is " << counterBOp_    << endl;
01569     os <<"The number of operations Prec*x is " << counterPrec_ << endl;
01570     os <<"The most recent rho was " << scientific << setprecision(10) << rho_ << endl;
01571     os <<"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl;
01572 
01573     if (initialized_) {
01574       os << endl;
01575       os <<"CURRENT EIGENVALUE ESTIMATES             "<<endl;
01576       os << setw(20) << "Eigenvalue" 
01577          << setw(20) << "Residual(B)"
01578          << setw(20) << "Residual(2)"
01579          << endl;
01580       os <<"--------------------------------------------------------------------------------"<<endl;
01581       for (int i=0; i<blockSize_; i++) {
01582         os << scientific << setprecision(10) << setw(20) << theta_[i];
01583         if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i];
01584         else os << scientific << setprecision(10) << setw(20) << "not current";
01585         if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i];
01586         else os << scientific << setprecision(10) << setw(20) << "not current";
01587         os << endl;
01588       }
01589     }
01590     os <<"================================================================================" << endl;
01591     os << endl;
01592   }
01593 
01594 
01596   // Inner product 1
01597   template <class ScalarType, class MV, class OP>
01598   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
01599   RTRBase<ScalarType,MV,OP>::ginner(const MV &xi) const 
01600   {
01601     std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi));
01602     MVT::MvNorm(xi,d);
01603     MagnitudeType ret = 0;
01604     for (vecMTiter i=d.begin(); i != d.end(); ++i) {
01605       ret += (*i)*(*i);
01606     }
01607     return ret;
01608   }
01609 
01610 
01612   // Inner product 2
01613   template <class ScalarType, class MV, class OP>
01614   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
01615   RTRBase<ScalarType,MV,OP>::ginner(const MV &xi, const MV &zeta) const 
01616   {
01617     std::vector<ScalarType> d(MVT::GetNumberVecs(xi));
01618     MVT::MvDot(xi,zeta,d);
01619     return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
01620   }
01621 
01622 
01624   // Inner product 1 without trace accumulation
01625   template <class ScalarType, class MV, class OP>
01626   void RTRBase<ScalarType,MV,OP>::ginnersep(
01627       const MV &xi, 
01628       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const 
01629   {
01630     MVT::MvNorm(xi,d);
01631     for (vecMTiter i=d.begin(); i != d.end(); ++i) {
01632       (*i) = (*i)*(*i);
01633     }
01634   }
01635 
01636 
01638   // Inner product 2 without trace accumulation
01639   template <class ScalarType, class MV, class OP>
01640   void RTRBase<ScalarType,MV,OP>::ginnersep(
01641       const MV &xi, const MV &zeta, 
01642       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const 
01643   {
01644     std::vector<ScalarType> dC(MVT::GetNumberVecs(xi));
01645     MVT::MvDot(xi,zeta,dC);
01646     vecMTiter iR=d.begin(); 
01647     vecSTiter iS=dC.begin();
01648     for (; iR != d.end(); iR++, iS++) {
01649       (*iR) = SCT::real(*iS);
01650     }
01651   }
01652 
01653   template <class ScalarType, class MV, class OP>
01654   Teuchos::Array<Teuchos::RCP<const MV> > RTRBase<ScalarType,MV,OP>::getAuxVecs() const {
01655     return auxVecs_;
01656   }
01657 
01658   template <class ScalarType, class MV, class OP>
01659   int RTRBase<ScalarType,MV,OP>::getBlockSize() const { 
01660     return(blockSize_); 
01661   }
01662 
01663   template <class ScalarType, class MV, class OP>
01664   const Eigenproblem<ScalarType,MV,OP>& RTRBase<ScalarType,MV,OP>::getProblem() const { 
01665     return(*problem_); 
01666   }
01667 
01668   template <class ScalarType, class MV, class OP>
01669   int RTRBase<ScalarType,MV,OP>::getMaxSubspaceDim() const {
01670     return blockSize_;
01671   }
01672 
01673   template <class ScalarType, class MV, class OP>
01674   int RTRBase<ScalarType,MV,OP>::getCurSubspaceDim() const 
01675   {
01676     if (!initialized_) return 0;
01677     return nevLocal_;
01678   }
01679 
01680   template <class ScalarType, class MV, class OP>
01681   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01682   RTRBase<ScalarType,MV,OP>::getRitzRes2Norms() 
01683   {
01684     std::vector<MagnitudeType> ret = ritz2norms_;
01685     ret.resize(nevLocal_);
01686     return ret;
01687   }
01688 
01689   template <class ScalarType, class MV, class OP>
01690   std::vector<Value<ScalarType> > 
01691   RTRBase<ScalarType,MV,OP>::getRitzValues() 
01692   {
01693     std::vector<Value<ScalarType> > ret(nevLocal_);
01694     for (int i=0; i<nevLocal_; i++) {
01695       ret[i].realpart = theta_[i];
01696       ret[i].imagpart = ZERO;
01697     }
01698     return ret;
01699   }
01700 
01701   template <class ScalarType, class MV, class OP>
01702   Teuchos::RCP<const MV> 
01703   RTRBase<ScalarType,MV,OP>::getRitzVectors() 
01704   {
01705     return X_;
01706   }
01707 
01708   template <class ScalarType, class MV, class OP>
01709   void RTRBase<ScalarType,MV,OP>::resetNumIters() 
01710   { 
01711     iter_=0; 
01712   }
01713 
01714   template <class ScalarType, class MV, class OP>
01715   int RTRBase<ScalarType,MV,OP>::getNumIters() const 
01716   { 
01717     return iter_; 
01718   }
01719 
01720   template <class ScalarType, class MV, class OP>
01721   RTRState<ScalarType,MV> RTRBase<ScalarType,MV,OP>::getState() const 
01722   {
01723     RTRState<ScalarType,MV> state;
01724     state.X = X_;
01725     state.AX = AX_;
01726     if (hasBOp_) {
01727       state.BX = BX_;
01728     }
01729     else {
01730       state.BX = Teuchos::null;
01731     }
01732     state.rho = rho_;
01733     state.R = R_;
01734     state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_));
01735     return state;
01736   }
01737 
01738   template <class ScalarType, class MV, class OP>
01739   bool RTRBase<ScalarType,MV,OP>::isInitialized() const 
01740   { 
01741     return initialized_; 
01742   }
01743 
01744   template <class ScalarType, class MV, class OP>
01745   std::vector<int> RTRBase<ScalarType,MV,OP>::getRitzIndex() 
01746   {
01747     std::vector<int> ret(nevLocal_,0);
01748     return ret;
01749   }
01750 
01751 
01752 } // end Anasazi namespace
01753 
01754 #endif // ANASAZI_RTRBASE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends