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_, olsenPrec_;
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 may 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 constant, 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 and olsenPrec)
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 and olsenPrec)
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     olsenPrec_ = params.get<bool>("Olsen Prec", true);
00669 
00670     TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument,
00671         "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product.");
00672 
00673     // set the block size and allocate data
00674     int bs = params.get("Block Size", problem_->getNEV());
00675     setBlockSize(bs);
00676 
00677     // leftmost or rightmost?
00678     leftMost_ = params.get("Leftmost",leftMost_);
00679 
00680     conv_kappa_ = params.get("Kappa Convergence",conv_kappa_);
00681     TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument,
00682                        "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
00683     conv_theta_ = params.get("Theta Convergence",conv_theta_);
00684     TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument,
00685                        "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
00686   }
00687 
00688 
00690   // Set the block size and make necessary adjustments.
00691   template <class ScalarType, class MV, class OP>
00692   void RTRBase<ScalarType,MV,OP>::setBlockSize (int blockSize) 
00693   {
00694     // time spent here counts towards timerInit_
00695     Teuchos::TimeMonitor lcltimer( *timerInit_ );
00696 
00697     // This routine only allocates space; it doesn't not perform any computation
00698     // if solver is initialized and size is to be decreased, take the first blockSize vectors of all to preserve state
00699     // otherwise, shrink/grow/allocate space and set solver to unitialized
00700 
00701     Teuchos::RCP<const MV> tmp;
00702     // grab some Multivector to Clone
00703     // in practice, getInitVec() should always provide this, but it is possible to use a 
00704     // Eigenproblem with nothing in getInitVec() by manually initializing with initialize(); 
00705     // in case of that strange scenario, we will try to Clone from R_
00706     // we like R_ for this, because it has minimal size (blockSize_), as opposed to V_ (blockSize_+numAuxVecs_)
00707     if (blockSize_ > 0) {
00708       tmp = R_;
00709     }
00710     else {
00711       tmp = problem_->getInitVec();
00712       TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
00713           "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
00714     }
00715 
00716     TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetVecLength(*tmp), std::invalid_argument, 
00717         "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
00718 
00719     // last chance to quit before causing side-effects
00720     if (blockSize == blockSize_) {
00721       // do nothing
00722       return;
00723     }
00724 
00725     // clear views
00726     X_  = Teuchos::null;
00727     BX_ = Teuchos::null;
00728 
00729     // regardless of whether we preserve any data, any content in V, BV and PBV corresponding to the
00730     // auxiliary vectors must be retained
00731     // go ahead and do these first
00732     // 
00733     // two cases here: 
00734     // * we are growing (possibly, from empty) 
00735     //   any aux data must be copied over, nothing from X need copying
00736     // * we are shrinking
00737     //   any aux data must be copied over, go ahead and copy X material (initialized or not)
00738     //
00739     if (blockSize > blockSize_)
00740     {
00741       // GROWING 
00742       // get a pointer for Q-related material, and an index for extracting/setting it
00743       Teuchos::RCP<const MV> Q;
00744       std::vector<int> indQ(numAuxVecs_);
00745       for (int i=0; i<numAuxVecs_; i++) indQ[i] = i;
00746       // if numAuxVecs_ > 0, then necessarily blockSize_ > 0 (we have already been allocated once)
00747       TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error,
00748           "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team.");
00749       // V
00750       if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ);
00751       V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00752       if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_);
00753       // BV
00754       if (hasBOp_) {
00755         if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ);
00756         BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00757         if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_);
00758       }
00759       else {
00760         BV_ = V_;
00761       }
00762       // PBV
00763       if (hasPrec_ && olsenPrec_) {
00764         if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ);
00765         PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00766         if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_);
00767       }
00768     }
00769     else 
00770     {
00771       // SHRINKING
00772       std::vector<int> indV(numAuxVecs_+blockSize);
00773       for (int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
00774       // V
00775       V_ = MVT::CloneCopy(*V_,indV);
00776       // BV
00777       if (hasBOp_) {
00778         BV_ = MVT::CloneCopy(*BV_,indV);
00779       }
00780       else {
00781         BV_ = V_;
00782       }
00783       // PBV
00784       if (hasPrec_ && olsenPrec_) {
00785         PBV_ = MVT::CloneCopy(*PBV_,indV);
00786       }
00787     }
00788 
00789     if (blockSize < blockSize_) {
00790       // shrink vectors
00791       blockSize_ = blockSize;
00792 
00793       theta_.resize(blockSize_);
00794       ritz2norms_.resize(blockSize_);
00795       Rnorms_.resize(blockSize_);
00796       R2norms_.resize(blockSize_);
00797 
00798       if (initialized_) {
00799         // shrink multivectors with copy
00800         std::vector<int> ind(blockSize_);
00801         for (int i=0; i<blockSize_; i++) ind[i] = i;
00802 
00803         // Z can be deleted, no important info there
00804         Z_ = Teuchos::null;
00805         
00806         // we will not use tmp for cloning, clear it and free some space
00807         tmp = Teuchos::null;
00808 
00809         R_      = MVT::CloneCopy(*R_     ,ind);
00810         eta_    = MVT::CloneCopy(*eta_   ,ind);
00811         delta_  = MVT::CloneCopy(*delta_ ,ind);
00812         Hdelta_ = MVT::CloneCopy(*Hdelta_,ind);
00813         if (!isSkinny_) {
00814           AX_     = MVT::CloneCopy(*AX_    ,ind);
00815           Aeta_   = MVT::CloneCopy(*Aeta_  ,ind);
00816           Adelta_ = MVT::CloneCopy(*Adelta_,ind);
00817         }
00818         else {
00819           AX_     = Teuchos::null;
00820           Aeta_   = Teuchos::null;
00821           Adelta_ = Teuchos::null;
00822         }
00823 
00824         if (hasBOp_) {
00825           if (!isSkinny_) {
00826             Beta_   = MVT::CloneCopy(*Beta_,ind);
00827             Bdelta_ = MVT::CloneCopy(*Bdelta_,ind);
00828           }
00829           else {
00830             Beta_   = Teuchos::null;
00831             Bdelta_ = Teuchos::null;
00832           }
00833         }
00834         else {
00835           Beta_   = eta_;
00836           Bdelta_ = delta_;
00837         }
00838         
00839         // we need Z if we have a preconditioner
00840         // we also use Z for temp storage in the skinny solvers
00841         if (hasPrec_ || isSkinny_) {
00842           Z_ = MVT::Clone(*V_,blockSize_);
00843         }
00844         else {
00845           Z_ = R_;
00846         }
00847 
00848       }
00849       else {
00850         // release previous multivectors
00851         // shrink multivectors without copying
00852         AX_     = Teuchos::null;
00853         R_      = Teuchos::null;
00854         eta_    = Teuchos::null;
00855         Aeta_   = Teuchos::null;
00856         delta_  = Teuchos::null;
00857         Adelta_ = Teuchos::null;
00858         Hdelta_ = Teuchos::null;
00859         Beta_   = Teuchos::null;
00860         Bdelta_ = Teuchos::null;
00861         Z_      = Teuchos::null;
00862 
00863         R_      = MVT::Clone(*tmp,blockSize_);
00864         eta_    = MVT::Clone(*tmp,blockSize_);
00865         delta_  = MVT::Clone(*tmp,blockSize_);
00866         Hdelta_ = MVT::Clone(*tmp,blockSize_);
00867         if (!isSkinny_) {
00868           AX_     = MVT::Clone(*tmp,blockSize_);
00869           Aeta_   = MVT::Clone(*tmp,blockSize_);
00870           Adelta_ = MVT::Clone(*tmp,blockSize_);
00871         }
00872 
00873         if (hasBOp_) {
00874           if (!isSkinny_) {
00875             Beta_   = MVT::Clone(*tmp,blockSize_);
00876             Bdelta_ = MVT::Clone(*tmp,blockSize_);
00877           }
00878         }
00879         else {
00880           Beta_   = eta_;
00881           Bdelta_ = delta_;
00882         }
00883 
00884         // we need Z if we have a preconditioner
00885         // we also use Z for temp storage in the skinny solvers
00886         if (hasPrec_ || isSkinny_) {
00887           Z_ = MVT::Clone(*tmp,blockSize_);
00888         }
00889         else {
00890           Z_ = R_;
00891         }
00892       } // if initialized_
00893     } // if blockSize is shrinking
00894     else {  // if blockSize > blockSize_
00895       // this is also the scenario for our initial call to setBlockSize(), in the constructor
00896       initialized_ = false;
00897 
00898       // grow/allocate vectors
00899       theta_.resize(blockSize,NANVAL);
00900       ritz2norms_.resize(blockSize,NANVAL);
00901       Rnorms_.resize(blockSize,NANVAL);
00902       R2norms_.resize(blockSize,NANVAL);
00903 
00904       // deallocate old multivectors
00905       AX_     = Teuchos::null;
00906       R_      = Teuchos::null;
00907       eta_    = Teuchos::null;
00908       Aeta_   = Teuchos::null;
00909       delta_  = Teuchos::null;
00910       Adelta_ = Teuchos::null;
00911       Hdelta_ = Teuchos::null;
00912       Beta_   = Teuchos::null;
00913       Bdelta_ = Teuchos::null;
00914       Z_      = Teuchos::null;
00915 
00916       // clone multivectors off of tmp
00917       R_      = MVT::Clone(*tmp,blockSize);
00918       eta_    = MVT::Clone(*tmp,blockSize);
00919       delta_  = MVT::Clone(*tmp,blockSize);
00920       Hdelta_ = MVT::Clone(*tmp,blockSize);
00921       if (!isSkinny_) {
00922         AX_     = MVT::Clone(*tmp,blockSize);
00923         Aeta_   = MVT::Clone(*tmp,blockSize);
00924         Adelta_ = MVT::Clone(*tmp,blockSize);
00925       }
00926 
00927       if (hasBOp_) {
00928         if (!isSkinny_) {
00929           Beta_   = MVT::Clone(*tmp,blockSize);
00930           Bdelta_ = MVT::Clone(*tmp,blockSize);
00931         }
00932       }
00933       else {
00934         Beta_   = eta_;
00935         Bdelta_ = delta_;
00936       }
00937       if (hasPrec_ || isSkinny_) {
00938         Z_ = MVT::Clone(*tmp,blockSize);
00939       }
00940       else {
00941         Z_ = R_;
00942       }
00943       blockSize_ = blockSize;
00944     }
00945 
00946     // get view of X from V, BX from BV
00947     // these are located after the first numAuxVecs columns
00948     {
00949       std::vector<int> indX(blockSize_);
00950       for (int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
00951       X_ = MVT::CloneView(*V_,indX);
00952       if (hasBOp_) {
00953         BX_ = MVT::CloneView(*BV_,indX);
00954       }
00955       else {
00956         BX_ = X_;
00957       }
00958     }
00959   }
00960 
00961 
00963   // Set a new StatusTest for the solver.
00964   template <class ScalarType, class MV, class OP>
00965   void RTRBase<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
00966     TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
00967         "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest.");
00968     tester_ = test;
00969   }
00970 
00971 
00973   // Get the current StatusTest used by the solver.
00974   template <class ScalarType, class MV, class OP>
00975   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > RTRBase<ScalarType,MV,OP>::getStatusTest() const {
00976     return tester_;
00977   }
00978   
00979 
00981   // Set the auxiliary vectors
00982   template <class ScalarType, class MV, class OP>
00983   void RTRBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00984     typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv;
00985 
00986     // set new auxiliary vectors
00987     auxVecs_.resize(0);
00988     auxVecs_.reserve(auxvecs.size());
00989 
00990     numAuxVecs_ = 0;
00991     for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
00992       numAuxVecs_ += MVT::GetNumberVecs(**v);
00993     }
00994 
00995     // If the solver has been initialized, X is not necessarily orthogonal to new auxiliary vectors
00996     if (numAuxVecs_ > 0 && initialized_) {
00997       initialized_ = false;
00998     }
00999 
01000     // clear X,BX views
01001     X_   = Teuchos::null;
01002     BX_  = Teuchos::null;
01003     // deallocate, we'll clone off R_ below
01004     V_   = Teuchos::null;
01005     BV_  = Teuchos::null;
01006     PBV_ = Teuchos::null;
01007 
01008     // put auxvecs into V, update BV and PBV if necessary
01009     if (numAuxVecs_ > 0) {
01010       V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01011       int numsofar = 0;
01012       for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
01013         std::vector<int> ind(MVT::GetNumberVecs(**v));
01014         for (unsigned int j=0; j<ind.size(); j++) ind[j] = numsofar++;
01015         MVT::SetBlock(**v,ind,*V_);
01016         auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind));
01017       }
01018       TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error,
01019           "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team.");
01020       // compute B*V, Prec*B*V
01021       if (hasBOp_) {
01022         BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01023         OPT::Apply(*BOp_,*V_,*BV_);
01024       }
01025       else {
01026         BV_ = V_;
01027       }
01028       if (hasPrec_ && olsenPrec_) {
01029         PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01030         OPT::Apply(*Prec_,*BV_,*V_);
01031       }
01032     }
01033     // 
01034 
01035     if (om_->isVerbosity( Debug ) ) {
01036       // Check almost everything here
01037       CheckList chk;
01038       chk.checkQ = true;
01039       om_->print( Debug, accuracyCheck(chk, "in setAuxVecs()") );
01040     }
01041   }
01042 
01043 
01045   /* Initialize the state of the solver
01046    * 
01047    * POST-CONDITIONS:
01048    *
01049    * initialized_ == true
01050    * X is orthonormal, orthogonal to auxVecs_
01051    * AX = A*X if not skinnny
01052    * BX = B*X if hasBOp_
01053    * theta_ contains Ritz values of X
01054    * R = AX - BX*diag(theta_)
01055    */
01056   template <class ScalarType, class MV, class OP>
01057   void RTRBase<ScalarType,MV,OP>::initialize(RTRState<ScalarType,MV> newstate)
01058   {
01059     // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
01060     // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
01061 
01062     // clear const views to X,BX in V,BV
01063     // work with temporary non-const views
01064     X_  = Teuchos::null;
01065     BX_ = Teuchos::null;
01066     Teuchos::RCP<MV> X, BX;
01067     {
01068       std::vector<int> ind(blockSize_);
01069       for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
01070       X = MVT::CloneView(*V_,ind);
01071       if (hasBOp_) {
01072         BX = MVT::CloneView(*BV_,ind);
01073       }
01074       else {
01075         BX = X;
01076       }
01077     }
01078 
01079     Teuchos::TimeMonitor inittimer( *timerInit_ );
01080 
01081     std::vector<int> bsind(blockSize_);
01082     for (int i=0; i<blockSize_; i++) bsind[i] = i;
01083 
01084     // in RTR, X (the subspace iterate) is primary
01085     // the order of dependence follows like so.
01086     // --init->                 X
01087     //    --op apply->          AX,BX
01088     //       --ritz analysis->  theta
01089     // 
01090     // if the user specifies all data for a level, we will accept it.
01091     // otherwise, we will generate the whole level, and all subsequent levels.
01092     //
01093     // the data members are ordered based on dependence, and the levels are
01094     // partitioned according to the amount of work required to produce the
01095     // items in a level.
01096     //
01097     // inconsitent multivectors widths and lengths will not be tolerated, and
01098     // will be treated with exceptions.
01099 
01100     // set up X, AX, BX: get them from "state" if user specified them
01101     if (newstate.X != Teuchos::null) {
01102       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.X) != MVT::GetVecLength(*X),
01103                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
01104       // newstate.X must have blockSize_ vectors; any more will be ignored
01105       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
01106                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
01107 
01108       // put data in X
01109       MVT::SetBlock(*newstate.X,bsind,*X);
01110 
01111       // put data in AX
01112       // if we are implementing a skinny solver, then we don't have memory allocated for AX
01113       // in this case, point AX at Z (skinny solvers allocate Z) and use it for temporary storage
01114       // we will clear the pointer later
01115       if (isSkinny_) {
01116         AX_ = Z_;
01117       }
01118       if (newstate.AX != Teuchos::null) {
01119         TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.AX) != MVT::GetVecLength(*X),
01120                             std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
01121         // newstate.AX must have blockSize_ vectors; any more will be ignored
01122         TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.AX) < blockSize_,
01123                             std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
01124         MVT::SetBlock(*newstate.AX,bsind,*AX_);
01125       }
01126       else {
01127         {
01128           Teuchos::TimeMonitor lcltimer( *timerAOp_ );
01129           OPT::Apply(*AOp_,*X,*AX_);
01130           counterAOp_ += blockSize_;
01131         }
01132         // we generated AX; we will generate R as well
01133         newstate.R = Teuchos::null;
01134       }
01135 
01136       // put data in BX
01137       // skinny solvers always allocate BX if hasB, so this is unconditionally appropriate
01138       if (hasBOp_) {
01139         if (newstate.BX != Teuchos::null) {
01140           TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.BX) != MVT::GetVecLength(*X),
01141                               std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
01142           // newstate.BX must have blockSize_ vectors; any more will be ignored
01143           TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.BX) < blockSize_,
01144                               std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
01145           MVT::SetBlock(*newstate.BX,bsind,*BX);
01146         }
01147         else {
01148           {
01149             Teuchos::TimeMonitor lcltimer( *timerBOp_ );
01150             OPT::Apply(*BOp_,*X,*BX);
01151             counterBOp_ += blockSize_;
01152           }
01153           // we generated BX; we will generate R as well
01154           newstate.R = Teuchos::null;
01155         }
01156       }
01157       else {
01158         // the assignment BX_==X_ would be redundant; take advantage of this opportunity to debug a little
01159         TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
01160       }
01161 
01162     }
01163     else {
01164       // user did not specify X
01165 
01166       // clear state so we won't use any data from it below
01167       newstate.R  = Teuchos::null;
01168       newstate.T  = Teuchos::null;
01169 
01170       // generate something and projectAndNormalize
01171       Teuchos::RCP<const MV> ivec = problem_->getInitVec();
01172       TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
01173                          "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from.");
01174 
01175       int initSize = MVT::GetNumberVecs(*ivec);
01176       if (initSize > blockSize_) {
01177         // we need only the first blockSize_ vectors from ivec; get a view of them
01178         initSize = blockSize_;
01179         std::vector<int> ind(blockSize_);
01180         for (int i=0; i<blockSize_; i++) ind[i] = i;
01181         ivec = MVT::CloneView(*ivec,ind);
01182       }
01183 
01184       // assign ivec to first part of X
01185       if (initSize > 0) {
01186         std::vector<int> ind(initSize);
01187         for (int i=0; i<initSize; i++) ind[i] = i;
01188         MVT::SetBlock(*ivec,ind,*X);
01189       }
01190       // fill the rest of X with random
01191       if (blockSize_ > initSize) {
01192         std::vector<int> ind(blockSize_ - initSize);
01193         for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
01194         Teuchos::RCP<MV> rX = MVT::CloneView(*X,ind);
01195         MVT::MvRandom(*rX);
01196         rX = Teuchos::null;
01197       }
01198 
01199       // put data in BX
01200       if (hasBOp_) {
01201         Teuchos::TimeMonitor lcltimer( *timerBOp_ );
01202         OPT::Apply(*BOp_,*X,*BX);
01203         counterBOp_ += blockSize_;
01204       }
01205       else {
01206         // the assignment BX==X would be redundant; take advantage of this opportunity to debug a little
01207         TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
01208       }
01209   
01210       // remove auxVecs from X and normalize it
01211       if (numAuxVecs_ > 0) {
01212         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01213         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01214         int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX);
01215         TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
01216                            "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
01217       }
01218       else {
01219         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01220         int rank = orthman_->normalizeMat(*X,Teuchos::null,BX);
01221         TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
01222                            "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
01223       }
01224 
01225       // put data in AX
01226       if (isSkinny_) {
01227         AX_ = Z_;
01228       }
01229       {
01230         Teuchos::TimeMonitor lcltimer( *timerAOp_ );
01231         OPT::Apply(*AOp_,*X,*AX_);
01232         counterAOp_ += blockSize_;
01233       }
01234 
01235     } // end if (newstate.X != Teuchos::null)
01236 
01237 
01238     // set up Ritz values
01239     if (newstate.T != Teuchos::null) {
01240       TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_,
01241                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values.");
01242       for (int i=0; i<blockSize_; i++) {
01243         theta_[i] = (*newstate.T)[i];
01244       }
01245     }
01246     else {
01247       // get ritz vecs/vals
01248       Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_),
01249                                                  BB(blockSize_,blockSize_),
01250                                                   S(blockSize_,blockSize_);
01251       // project A
01252       {
01253         Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
01254         MVT::MvTransMv(ONE,*X,*AX_,AA);
01255         if (hasBOp_) {
01256           MVT::MvTransMv(ONE,*X,*BX,BB);
01257         }
01258       }
01259       nevLocal_ = blockSize_;
01260 
01261       // solve the projected problem
01262       // project B
01263       int ret;
01264       if (hasBOp_) {
01265         Teuchos::TimeMonitor lcltimer( *timerDS_ );
01266         ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
01267       }
01268       else {
01269         Teuchos::TimeMonitor lcltimer( *timerDS_ );
01270         ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10);
01271       }
01272       TEST_FOR_EXCEPTION(ret != 0,RTRInitFailure,
01273           "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret);
01274       TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,RTRInitFailure,"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis.");
01275 
01276       // We only have blockSize_ ritz pairs, ergo we do not need to select.
01277       // However, we still require them to be ordered correctly
01278       {
01279         Teuchos::TimeMonitor lcltimer( *timerSort_ );
01280         
01281         std::vector<int> order(blockSize_);
01282         // 
01283         // sort the first blockSize_ values in theta_
01284         sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);   // don't catch exception
01285         //
01286         // apply the same ordering to the primitive ritz vectors
01287         Utils::permuteVectors(order,S);
01288       }
01289 
01290       // compute ritz residual norms
01291       {
01292         Teuchos::BLAS<int,ScalarType> blas;
01293         Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_);
01294         // RR = AA*S - BB*S*diag(theta)
01295         int info;
01296         if (hasBOp_) {
01297           info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO);
01298           TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
01299         }
01300         else {
01301           RR.assign(S);
01302         }
01303         for (int i=0; i<blockSize_; i++) {
01304           blas.SCAL(blockSize_,theta_[i],RR[i],1);
01305         }
01306         info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE);
01307         TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
01308         for (int i=0; i<blockSize_; i++) {
01309           ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1);
01310         }
01311       }
01312 
01313       // update the solution
01314       // use R as local work space
01315       // Z may already be in use as work space (holding AX if isSkinny)
01316       {
01317         Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
01318         // X <- X*S
01319         MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ );        
01320         MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X );
01321         // AX <- AX*S
01322         MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ );        
01323         MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ );
01324         if (hasBOp_) {
01325           // BX <- BX*S
01326           MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ );        
01327           MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX );
01328         }
01329       }
01330     }
01331     
01332     // done modifying X,BX; clear temp views and setup const views
01333     X  = Teuchos::null;
01334     BX = Teuchos::null;
01335     {
01336       std::vector<int> ind(blockSize_);
01337       for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
01338       this->X_ = MVT::CloneView(static_cast<const MV&>(*V_),ind);
01339       if (hasBOp_) {
01340         this->BX_ = MVT::CloneView(static_cast<const MV&>(*BV_),ind);
01341       }
01342       else {
01343         this->BX_ = this->X_;
01344       }
01345     }
01346 
01347 
01348     // get objective function value
01349     fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
01350 
01351     // set up R
01352     if (newstate.R != Teuchos::null) {
01353       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
01354                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
01355       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
01356                           std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
01357       MVT::SetBlock(*newstate.R,bsind,*R_);
01358     }
01359     else {
01360       Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
01361       // form R <- AX - BX*T
01362       MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
01363       Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
01364       T.putScalar(ZERO);
01365       for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
01366       MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_);
01367     }
01368 
01369     // R has been updated; mark the norms as out-of-date
01370     Rnorms_current_ = false;
01371     R2norms_current_ = false;
01372 
01373     // if isSkinny, then AX currently points to Z for temp storage
01374     // set AX back to null
01375     if (isSkinny_) {
01376       AX_ = Teuchos::null;
01377     }
01378 
01379     // finally, we are initialized
01380     initialized_ = true;
01381 
01382     if (om_->isVerbosity( Debug ) ) {
01383       // Check almost everything here
01384       CheckList chk;
01385       chk.checkX = true;
01386       chk.checkAX = true;
01387       chk.checkBX = true;
01388       chk.checkR = true;
01389       chk.checkQ = true;
01390       om_->print( Debug, accuracyCheck(chk, "after initialize()") );
01391     }
01392   }
01393 
01394   template <class ScalarType, class MV, class OP>
01395   void RTRBase<ScalarType,MV,OP>::initialize()
01396   {
01397     RTRState<ScalarType,MV> empty;
01398     initialize(empty);
01399   }
01400 
01401 
01402 
01403 
01405   // compute/return residual B-norms
01406   template <class ScalarType, class MV, class OP>
01407   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01408   RTRBase<ScalarType,MV,OP>::getResNorms() {
01409     if (Rnorms_current_ == false) {
01410       // Update the residual norms
01411       orthman_->norm(*R_,Rnorms_);
01412       Rnorms_current_ = true;
01413     }
01414     return Rnorms_;
01415   }
01416 
01417   
01419   // compute/return residual 2-norms
01420   template <class ScalarType, class MV, class OP>
01421   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01422   RTRBase<ScalarType,MV,OP>::getRes2Norms() {
01423     if (R2norms_current_ == false) {
01424       // Update the residual 2-norms 
01425       MVT::MvNorm(*R_,R2norms_);
01426       R2norms_current_ = true;
01427     }
01428     return R2norms_;
01429   }
01430 
01431 
01432 
01433 
01435   // Check accuracy, orthogonality, and other debugging stuff
01436   // 
01437   // bools specify which tests we want to run (instead of running more than we actually care about)
01438   //
01439   // we don't bother checking the following because they are computed explicitly:
01440   //   AH == A*H
01441   //
01442   // 
01443   // checkX    : X orthonormal
01444   //             orthogonal to auxvecs
01445   // checkAX   : check AX == A*X
01446   // checkBX   : check BX == B*X
01447   // checkEta  : check that Eta'*B*X == 0
01448   //             orthogonal to auxvecs
01449   // checkAEta : check that AEta = A*Eta
01450   // checkBEta : check that BEta = B*Eta
01451   // checkR    : check R orthogonal to X
01452   // checkBR   : check R B-orthogonal to X, valid in and immediately after solveTRSubproblem
01453   // checkQ    : check that auxiliary vectors are actually orthonormal
01454   //
01455   // TODO: 
01456   //  add checkTheta 
01457   //
01458   template <class ScalarType, class MV, class OP>
01459   std::string RTRBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const 
01460   {
01461     using std::setprecision;
01462     using std::scientific;
01463     using std::setw;
01464     using std::endl;
01465     std::stringstream os;
01466     MagnitudeType tmp;
01467 
01468     os << " Debugging checks: " << where << endl;
01469 
01470     // X and friends
01471     if (chk.checkX && initialized_) {
01472       tmp = orthman_->orthonormError(*X_);
01473       os << " >> Error in X^H B X == I :    " << scientific << setprecision(10) << tmp << endl;
01474       for (unsigned int i=0; i<auxVecs_.size(); i++) {
01475         tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
01476         os << " >> Error in X^H B Q[" << i << "] == 0 : " << scientific << setprecision(10) << tmp << endl;
01477       }
01478     }
01479     if (chk.checkAX && !isSkinny_ && initialized_) {
01480       tmp = Utils::errorEquality(*X_, *AX_, AOp_);
01481       os << " >> Error in AX == A*X    :    " << scientific << setprecision(10) << tmp << endl;
01482     }
01483     if (chk.checkBX && hasBOp_ && initialized_) {
01484       Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_);
01485       OPT::Apply(*BOp_,*this->X_,*BX);
01486       tmp = Utils::errorEquality(*BX, *BX_);
01487       os << " >> Error in BX == B*X    :    " << scientific << setprecision(10) << tmp << endl;
01488     }
01489 
01490     // Eta and friends
01491     if (chk.checkEta && initialized_) {
01492       tmp = orthman_->orthogError(*X_,*eta_);
01493       os << " >> Error in X^H B Eta == 0 :  " << scientific << setprecision(10) << tmp << endl;
01494       for (unsigned int i=0; i<auxVecs_.size(); i++) {
01495         tmp = orthman_->orthogError(*eta_,*auxVecs_[i]);
01496         os << " >> Error in Eta^H B Q[" << i << "]==0 : " << scientific << setprecision(10) << tmp << endl;
01497       }
01498     }
01499     if (chk.checkAEta && !isSkinny_ && initialized_) {
01500       tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_);
01501       os << " >> Error in AEta == A*Eta    :    " << scientific << setprecision(10) << tmp << endl;
01502     }
01503     if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
01504       tmp = Utils::errorEquality(*eta_, *Beta_, BOp_);
01505       os << " >> Error in BEta == B*Eta    :    " << scientific << setprecision(10) << tmp << endl;
01506     }
01507 
01508     // R: this is not B-orthogonality, but standard euclidean orthogonality
01509     if (chk.checkR && initialized_) {
01510       Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01511       MVT::MvTransMv(ONE,*X_,*R_,xTx);
01512       tmp = xTx.normFrobenius();
01513       os << " >> Error in X^H R == 0   :    " << scientific << setprecision(10) << tmp << endl;
01514     }
01515 
01516     // BR: this is B-orthogonality: this is only valid inside and immediately after solveTRSubproblem 
01517     // only check if B != I, otherwise, it is equivalent to the above test
01518     if (chk.checkBR && hasBOp_ && initialized_) {
01519       Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01520       MVT::MvTransMv(ONE,*BX_,*R_,xTx);
01521       tmp = xTx.normFrobenius();
01522       os << " >> Error in X^H B R == 0 :    " << scientific << setprecision(10) << tmp << endl;
01523     }
01524 
01525     // Z: Z is preconditioned R, should be on tangent plane
01526     if (chk.checkZ && initialized_) {
01527       Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01528       MVT::MvTransMv(ONE,*BX_,*Z_,xTx);
01529       tmp = xTx.normFrobenius();
01530       os << " >> Error in X^H B Z == 0 :    " << scientific << setprecision(10) << tmp << endl;
01531     }
01532 
01533     // Q
01534     if (chk.checkQ) {
01535       for (unsigned int i=0; i<auxVecs_.size(); i++) {
01536         tmp = orthman_->orthonormError(*auxVecs_[i]);
01537         os << " >> Error in Q[" << i << "]^H B Q[" << i << "]==I: " << scientific << setprecision(10) << tmp << endl;
01538         for (unsigned int j=i+1; j<auxVecs_.size(); j++) {
01539           tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01540           os << " >> Error in Q[" << i << "]^H B Q[" << j << "]==0: " << scientific << setprecision(10) << tmp << endl;
01541         }
01542       }
01543     }
01544     os << endl;
01545     return os.str();
01546   }
01547 
01548 
01550   // Print the current status of the solver
01551   template <class ScalarType, class MV, class OP>
01552   void 
01553   RTRBase<ScalarType,MV,OP>::currentStatus(std::ostream &os) 
01554   {
01555     using std::setprecision;
01556     using std::scientific;
01557     using std::setw;
01558     using std::endl;
01559 
01560     os <<endl;
01561     os <<"================================================================================" << endl;
01562     os << endl;
01563     os <<"                              RTRBase Solver Status" << endl;
01564     os << endl;
01565     os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01566     os <<"The number of iterations performed is " << iter_       << endl;
01567     os <<"The current block size is             " << blockSize_  << endl;
01568     os <<"The number of auxiliary vectors is    " << numAuxVecs_ << endl;
01569     os <<"The number of operations A*x    is " << counterAOp_   << endl;
01570     os <<"The number of operations B*x    is " << counterBOp_    << endl;
01571     os <<"The number of operations Prec*x is " << counterPrec_ << endl;
01572     os <<"The most recent rho was " << scientific << setprecision(10) << rho_ << endl;
01573     os <<"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl;
01574 
01575     if (initialized_) {
01576       os << endl;
01577       os <<"CURRENT EIGENVALUE ESTIMATES             "<<endl;
01578       os << setw(20) << "Eigenvalue" 
01579          << setw(20) << "Residual(B)"
01580          << setw(20) << "Residual(2)"
01581          << endl;
01582       os <<"--------------------------------------------------------------------------------"<<endl;
01583       for (int i=0; i<blockSize_; i++) {
01584         os << scientific << setprecision(10) << setw(20) << theta_[i];
01585         if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i];
01586         else os << scientific << setprecision(10) << setw(20) << "not current";
01587         if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i];
01588         else os << scientific << setprecision(10) << setw(20) << "not current";
01589         os << endl;
01590       }
01591     }
01592     os <<"================================================================================" << endl;
01593     os << endl;
01594   }
01595 
01596 
01598   // Inner product 1
01599   template <class ScalarType, class MV, class OP>
01600   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
01601   RTRBase<ScalarType,MV,OP>::ginner(const MV &xi) const 
01602   {
01603     std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi));
01604     MVT::MvNorm(xi,d);
01605     MagnitudeType ret = 0;
01606     for (vecMTiter i=d.begin(); i != d.end(); ++i) {
01607       ret += (*i)*(*i);
01608     }
01609     return ret;
01610   }
01611 
01612 
01614   // Inner product 2
01615   template <class ScalarType, class MV, class OP>
01616   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
01617   RTRBase<ScalarType,MV,OP>::ginner(const MV &xi, const MV &zeta) const 
01618   {
01619     std::vector<ScalarType> d(MVT::GetNumberVecs(xi));
01620     MVT::MvDot(xi,zeta,d);
01621     return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
01622   }
01623 
01624 
01626   // Inner product 1 without trace accumulation
01627   template <class ScalarType, class MV, class OP>
01628   void RTRBase<ScalarType,MV,OP>::ginnersep(
01629       const MV &xi, 
01630       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const 
01631   {
01632     MVT::MvNorm(xi,d);
01633     for (vecMTiter i=d.begin(); i != d.end(); ++i) {
01634       (*i) = (*i)*(*i);
01635     }
01636   }
01637 
01638 
01640   // Inner product 2 without trace accumulation
01641   template <class ScalarType, class MV, class OP>
01642   void RTRBase<ScalarType,MV,OP>::ginnersep(
01643       const MV &xi, const MV &zeta, 
01644       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const 
01645   {
01646     std::vector<ScalarType> dC(MVT::GetNumberVecs(xi));
01647     MVT::MvDot(xi,zeta,dC);
01648     vecMTiter iR=d.begin(); 
01649     vecSTiter iS=dC.begin();
01650     for (; iR != d.end(); iR++, iS++) {
01651       (*iR) = SCT::real(*iS);
01652     }
01653   }
01654 
01655   template <class ScalarType, class MV, class OP>
01656   Teuchos::Array<Teuchos::RCP<const MV> > RTRBase<ScalarType,MV,OP>::getAuxVecs() const {
01657     return auxVecs_;
01658   }
01659 
01660   template <class ScalarType, class MV, class OP>
01661   int RTRBase<ScalarType,MV,OP>::getBlockSize() const { 
01662     return(blockSize_); 
01663   }
01664 
01665   template <class ScalarType, class MV, class OP>
01666   const Eigenproblem<ScalarType,MV,OP>& RTRBase<ScalarType,MV,OP>::getProblem() const { 
01667     return(*problem_); 
01668   }
01669 
01670   template <class ScalarType, class MV, class OP>
01671   int RTRBase<ScalarType,MV,OP>::getMaxSubspaceDim() const {
01672     return blockSize_;
01673   }
01674 
01675   template <class ScalarType, class MV, class OP>
01676   int RTRBase<ScalarType,MV,OP>::getCurSubspaceDim() const 
01677   {
01678     if (!initialized_) return 0;
01679     return nevLocal_;
01680   }
01681 
01682   template <class ScalarType, class MV, class OP>
01683   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01684   RTRBase<ScalarType,MV,OP>::getRitzRes2Norms() 
01685   {
01686     std::vector<MagnitudeType> ret = ritz2norms_;
01687     ret.resize(nevLocal_);
01688     return ret;
01689   }
01690 
01691   template <class ScalarType, class MV, class OP>
01692   std::vector<Value<ScalarType> > 
01693   RTRBase<ScalarType,MV,OP>::getRitzValues() 
01694   {
01695     std::vector<Value<ScalarType> > ret(nevLocal_);
01696     for (int i=0; i<nevLocal_; i++) {
01697       ret[i].realpart = theta_[i];
01698       ret[i].imagpart = ZERO;
01699     }
01700     return ret;
01701   }
01702 
01703   template <class ScalarType, class MV, class OP>
01704   Teuchos::RCP<const MV> 
01705   RTRBase<ScalarType,MV,OP>::getRitzVectors() 
01706   {
01707     return X_;
01708   }
01709 
01710   template <class ScalarType, class MV, class OP>
01711   void RTRBase<ScalarType,MV,OP>::resetNumIters() 
01712   { 
01713     iter_=0; 
01714   }
01715 
01716   template <class ScalarType, class MV, class OP>
01717   int RTRBase<ScalarType,MV,OP>::getNumIters() const 
01718   { 
01719     return iter_; 
01720   }
01721 
01722   template <class ScalarType, class MV, class OP>
01723   RTRState<ScalarType,MV> RTRBase<ScalarType,MV,OP>::getState() const 
01724   {
01725     RTRState<ScalarType,MV> state;
01726     state.X = X_;
01727     state.AX = AX_;
01728     if (hasBOp_) {
01729       state.BX = BX_;
01730     }
01731     else {
01732       state.BX = Teuchos::null;
01733     }
01734     state.rho = rho_;
01735     state.R = R_;
01736     state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_));
01737     return state;
01738   }
01739 
01740   template <class ScalarType, class MV, class OP>
01741   bool RTRBase<ScalarType,MV,OP>::isInitialized() const 
01742   { 
01743     return initialized_; 
01744   }
01745 
01746   template <class ScalarType, class MV, class OP>
01747   std::vector<int> RTRBase<ScalarType,MV,OP>::getRitzIndex() 
01748   {
01749     std::vector<int> ret(nevLocal_,0);
01750     return ret;
01751   }
01752 
01753 
01754 } // end Anasazi namespace
01755 
01756 #endif // ANASAZI_RTRBASE_HPP

Generated on Wed May 12 21:24:34 2010 for Anasazi by  doxygen 1.4.7