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