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

Generated on Thu Dec 17 11:02:21 2009 for Anasazi by  doxygen 1.5.9