Anasazi Version of the Day
AnasaziBlockDavidson.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_BLOCK_DAVIDSON_HPP
00034 #define ANASAZI_BLOCK_DAVIDSON_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 "AnasaziMatOrthoManager.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 
00067 namespace Anasazi {
00068 
00070 
00071 
00076   template <class ScalarType, class MV>
00077   struct BlockDavidsonState {
00082     int curDim;
00087     Teuchos::RCP<const MV> V;
00089     Teuchos::RCP<const MV> X; 
00091     Teuchos::RCP<const MV> KX; 
00093     Teuchos::RCP<const MV> MX;
00095     Teuchos::RCP<const MV> R;
00100     Teuchos::RCP<const MV> H;
00102     Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
00108     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > KK;
00109     BlockDavidsonState() : curDim(0), V(Teuchos::null),
00110                            X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null),
00111                            R(Teuchos::null), H(Teuchos::null),
00112                            T(Teuchos::null), KK(Teuchos::null) {}
00113   };
00114 
00116 
00118 
00119 
00132   class BlockDavidsonInitFailure : public AnasaziError {public:
00133     BlockDavidsonInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
00134     {}};
00135 
00143   class BlockDavidsonOrthoFailure : public AnasaziError {public:
00144     BlockDavidsonOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
00145     {}};
00146   
00148 
00149 
00150   template <class ScalarType, class MV, class OP>
00151   class BlockDavidson : public Eigensolver<ScalarType,MV,OP> { 
00152   public:
00154 
00155     
00163     BlockDavidson( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> >    &problem, 
00164                    const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00165                    const Teuchos::RCP<OutputManager<ScalarType> >         &printer,
00166                    const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >      &tester,
00167                    const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00168                    Teuchos::ParameterList &params 
00169                  );
00170     
00172     virtual ~BlockDavidson();
00174 
00175 
00177 
00178 
00202     void iterate();
00203 
00237     void initialize(BlockDavidsonState<ScalarType,MV> newstate);
00238 
00242     void initialize();
00243 
00259     bool isInitialized() const;
00260 
00273     BlockDavidsonState<ScalarType,MV> getState() const;
00274     
00276 
00277 
00279 
00280 
00282     int getNumIters() const;
00283 
00285     void resetNumIters();
00286 
00294     Teuchos::RCP<const MV> getRitzVectors();
00295 
00301     std::vector<Value<ScalarType> > getRitzValues();
00302 
00303 
00311     std::vector<int> getRitzIndex();
00312 
00313 
00319     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
00320 
00321 
00327     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
00328 
00329 
00337     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
00338 
00344     int getCurSubspaceDim() const;
00345 
00347     int getMaxSubspaceDim() const;
00348 
00350 
00351 
00353 
00354 
00356     void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
00357 
00359     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
00360 
00362     const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
00363 
00373     void setBlockSize(int blockSize);
00374 
00376     int getBlockSize() const;
00377 
00390     void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
00391 
00393     Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
00394 
00396 
00398 
00399 
00409     void setSize(int blockSize, int numBlocks);
00410 
00412 
00414 
00415 
00417     void currentStatus(std::ostream &os);
00418 
00420 
00421   private:
00422     //
00423     // Convenience typedefs
00424     //
00425     typedef SolverUtils<ScalarType,MV,OP> Utils;
00426     typedef MultiVecTraits<ScalarType,MV> MVT;
00427     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00428     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00429     typedef typename SCT::magnitudeType MagnitudeType;
00430     const MagnitudeType ONE;  
00431     const MagnitudeType ZERO; 
00432     const MagnitudeType NANVAL;
00433     //
00434     // Internal structs
00435     //
00436     struct CheckList {
00437       bool checkV;
00438       bool checkX, checkMX, checkKX;
00439       bool checkH, checkMH, checkKH;
00440       bool checkR, checkQ;
00441       bool checkKK;
00442       CheckList() : checkV(false),
00443                     checkX(false),checkMX(false),checkKX(false),
00444                     checkH(false),checkMH(false),checkKH(false),
00445                     checkR(false),checkQ(false),checkKK(false) {};
00446     };
00447     //
00448     // Internal methods
00449     //
00450     std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
00451     //
00452     // Classes inputed through constructor that define the eigenproblem to be solved.
00453     //
00454     const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> >     problem_;
00455     const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
00456     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00457     Teuchos::RCP<StatusTest<ScalarType,MV,OP> >             tester_;
00458     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >  orthman_;
00459     //
00460     // Information obtained from the eigenproblem
00461     //
00462     Teuchos::RCP<const OP> Op_;
00463     Teuchos::RCP<const OP> MOp_;
00464     Teuchos::RCP<const OP> Prec_;
00465     bool hasM_;
00466     //
00467     // Internal timers
00468     //
00469     Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
00470                                         timerSortEval_, timerDS_,
00471                                         timerLocal_, timerCompRes_, 
00472                                         timerOrtho_, timerInit_;
00473     //
00474     // Counters
00475     //
00476     int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
00477 
00478     //
00479     // Algorithmic parameters.
00480     //
00481     // blockSize_ is the solver block size; it controls the number of eigenvectors that 
00482     // we compute, the number of residual vectors that we compute, and therefore the number
00483     // of vectors added to the basis on each iteration.
00484     int blockSize_;
00485     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00486     int numBlocks_; 
00487     
00488     // 
00489     // Current solver state
00490     //
00491     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00492     // is capable of running; _initialize is controlled  by the initialize() member method
00493     // For the implications of the state of initialized_, please see documentation for initialize()
00494     bool initialized_;
00495     //
00496     // curDim_ reflects how much of the current basis is valid 
00497     // NOTE: 0 <= curDim_ <= blockSize_*numBlocks_
00498     // this also tells us how many of the values in theta_ are valid Ritz values
00499     int curDim_;
00500     //
00501     // State Multivecs
00502     // H_,KH_,MH_ will not own any storage
00503     // H_ will occasionally point at the current block of vectors in the basis V_
00504     // MH_,KH_ will occasionally point at MX_,KX_ when they are used as temporary storage
00505     Teuchos::RCP<MV> X_, KX_, MX_, R_,
00506                      H_, KH_, MH_,
00507                      V_;
00508     //
00509     // Projected matrices
00510     //
00511     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_;
00512     // 
00513     // auxiliary vectors
00514     Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
00515     int numAuxVecs_;
00516     //
00517     // Number of iterations that have been performed.
00518     int iter_;
00519     // 
00520     // Current eigenvalues, residual norms
00521     std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
00522     // 
00523     // are the residual norms current with the residual?
00524     bool Rnorms_current_, R2norms_current_;
00525 
00526   };
00527 
00530   //
00531   // Implementations
00532   //
00535 
00536 
00538   // Constructor
00539   template <class ScalarType, class MV, class OP>
00540   BlockDavidson<ScalarType,MV,OP>::BlockDavidson(
00541         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> >    &problem, 
00542         const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00543         const Teuchos::RCP<OutputManager<ScalarType> >         &printer,
00544         const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >      &tester,
00545         const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00546         Teuchos::ParameterList &params
00547         ) :
00548     ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
00549     ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00550     NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
00551     // problem, tools
00552     problem_(problem), 
00553     sm_(sorter),
00554     om_(printer),
00555     tester_(tester),
00556     orthman_(ortho),
00557     // timers, counters
00558     timerOp_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Operation Op*x")),
00559     timerMOp_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Operation M*x")),
00560     timerPrec_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Operation Prec*x")),
00561     timerSortEval_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Sorting eigenvalues")),
00562     timerDS_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Direct solve")),
00563     timerLocal_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Local update")),
00564     timerCompRes_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Computing residuals")),
00565     timerOrtho_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Orthogonalization")),
00566     timerInit_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Initialization")),
00567     count_ApplyOp_(0),
00568     count_ApplyM_(0),
00569     count_ApplyPrec_(0),
00570     // internal data
00571     blockSize_(0),
00572     numBlocks_(0),
00573     initialized_(false),
00574     curDim_(0),
00575     auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ), 
00576     numAuxVecs_(0),
00577     iter_(0),
00578     Rnorms_current_(false),
00579     R2norms_current_(false)
00580   {     
00581     TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
00582                        "Anasazi::BlockDavidson::constructor: user passed null problem pointer.");
00583     TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
00584                        "Anasazi::BlockDavidson::constructor: user passed null sort manager pointer.");
00585     TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
00586                        "Anasazi::BlockDavidson::constructor: user passed null output manager pointer.");
00587     TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
00588                        "Anasazi::BlockDavidson::constructor: user passed null status test pointer.");
00589     TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
00590                        "Anasazi::BlockDavidson::constructor: user passed null orthogonalization manager pointer.");
00591     TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
00592                        "Anasazi::BlockDavidson::constructor: problem is not set.");
00593     TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
00594                        "Anasazi::BlockDavidson::constructor: problem is not hermitian.");
00595 
00596     // get the problem operators
00597     Op_   = problem_->getOperator();
00598     TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
00599                        "Anasazi::BlockDavidson::constructor: problem provides no operator.");
00600     MOp_  = problem_->getM();
00601     Prec_ = problem_->getPrec();
00602     hasM_ = (MOp_ != Teuchos::null);
00603 
00604     // set the block size and allocate data
00605     int bs = params.get("Block Size", problem_->getNEV());
00606     int nb = params.get("Num Blocks", 2);
00607     setSize(bs,nb);
00608   }
00609 
00610 
00612   // Destructor
00613   template <class ScalarType, class MV, class OP>
00614   BlockDavidson<ScalarType,MV,OP>::~BlockDavidson() {}
00615 
00616 
00618   // Set the block size
00619   // This simply calls setSize(), modifying the block size while retaining the number of blocks.
00620   template <class ScalarType, class MV, class OP>
00621   void BlockDavidson<ScalarType,MV,OP>::setBlockSize (int blockSize) 
00622   {
00623     setSize(blockSize,numBlocks_);
00624   }
00625 
00626 
00628   // Return the current auxiliary vectors
00629   template <class ScalarType, class MV, class OP>
00630   Teuchos::Array<Teuchos::RCP<const MV> > BlockDavidson<ScalarType,MV,OP>::getAuxVecs() const {
00631     return auxVecs_;
00632   }
00633 
00634 
00635 
00637   // return the current block size
00638   template <class ScalarType, class MV, class OP>
00639   int BlockDavidson<ScalarType,MV,OP>::getBlockSize() const {
00640     return(blockSize_); 
00641   }
00642 
00643 
00645   // return eigenproblem
00646   template <class ScalarType, class MV, class OP>
00647   const Eigenproblem<ScalarType,MV,OP>& BlockDavidson<ScalarType,MV,OP>::getProblem() const { 
00648     return(*problem_); 
00649   }
00650 
00651 
00653   // return max subspace dim
00654   template <class ScalarType, class MV, class OP>
00655   int BlockDavidson<ScalarType,MV,OP>::getMaxSubspaceDim() const {
00656     return blockSize_*numBlocks_;
00657   }
00658 
00660   // return current subspace dim
00661   template <class ScalarType, class MV, class OP>
00662   int BlockDavidson<ScalarType,MV,OP>::getCurSubspaceDim() const {
00663     if (!initialized_) return 0;
00664     return curDim_;
00665   }
00666 
00667 
00669   // return ritz residual 2-norms
00670   template <class ScalarType, class MV, class OP>
00671   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
00672   BlockDavidson<ScalarType,MV,OP>::getRitzRes2Norms() {
00673     return this->getRes2Norms();
00674   }
00675 
00676 
00678   // return ritz index
00679   template <class ScalarType, class MV, class OP>
00680   std::vector<int> BlockDavidson<ScalarType,MV,OP>::getRitzIndex() {
00681     std::vector<int> ret(curDim_,0);
00682     return ret;
00683   }
00684 
00685 
00687   // return ritz values
00688   template <class ScalarType, class MV, class OP>
00689   std::vector<Value<ScalarType> > BlockDavidson<ScalarType,MV,OP>::getRitzValues() { 
00690     std::vector<Value<ScalarType> > ret(curDim_);
00691     for (int i=0; i<curDim_; ++i) {
00692       ret[i].realpart = theta_[i];
00693       ret[i].imagpart = ZERO;
00694     }
00695     return ret;
00696   }
00697 
00698 
00700   // return pointer to ritz vectors
00701   template <class ScalarType, class MV, class OP>
00702   Teuchos::RCP<const MV> BlockDavidson<ScalarType,MV,OP>::getRitzVectors() {
00703     return X_;
00704   }
00705 
00706 
00708   // reset number of iterations
00709   template <class ScalarType, class MV, class OP>
00710   void BlockDavidson<ScalarType,MV,OP>::resetNumIters() { 
00711     iter_=0; 
00712   }
00713 
00714 
00716   // return number of iterations
00717   template <class ScalarType, class MV, class OP>
00718   int BlockDavidson<ScalarType,MV,OP>::getNumIters() const { 
00719     return(iter_); 
00720   }
00721 
00722 
00724   // return state pointers
00725   template <class ScalarType, class MV, class OP>
00726   BlockDavidsonState<ScalarType,MV> BlockDavidson<ScalarType,MV,OP>::getState() const {
00727     BlockDavidsonState<ScalarType,MV> state;
00728     state.curDim = curDim_;
00729     state.V = V_;
00730     state.X = X_;
00731     state.KX = KX_;
00732     if (hasM_) {
00733       state.MX = MX_;
00734     }
00735     else {
00736       state.MX = Teuchos::null;
00737     }
00738     state.R = R_;
00739     state.H = H_;
00740     state.KK = KK_;
00741     if (curDim_ > 0) {
00742       state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
00743     } else {
00744       state.T = Teuchos::rcp(new std::vector<MagnitudeType>(0));
00745     }
00746     return state;
00747   }
00748 
00749 
00751   // Return initialized state
00752   template <class ScalarType, class MV, class OP>
00753   bool BlockDavidson<ScalarType,MV,OP>::isInitialized() const { return initialized_; }
00754 
00755 
00757   // Set the block size and make necessary adjustments.
00758   template <class ScalarType, class MV, class OP>
00759   void BlockDavidson<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks) 
00760   {
00761     // time spent here counts towards timerInit_
00762     Teuchos::TimeMonitor initimer( *timerInit_ );
00763 
00764     // This routine only allocates space; it doesn't not perform any computation
00765     // any change in size will invalidate the state of the solver.
00766 
00767     TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): blocksize must be strictly positive.");
00768     TEST_FOR_EXCEPTION(numBlocks < 2, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): numblocks must be greater than one.");
00769     if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00770       // do nothing
00771       return;
00772     }
00773 
00774     blockSize_ = blockSize;
00775     numBlocks_ = numBlocks;
00776 
00777     Teuchos::RCP<const MV> tmp;
00778     // grab some Multivector to Clone
00779     // in practice, getInitVec() should always provide this, but it is possible to use a 
00780     // Eigenproblem with nothing in getInitVec() by manually initializing with initialize(); 
00781     // in case of that strange scenario, we will try to Clone from X_ first, then resort to getInitVec()
00782     if (X_ != Teuchos::null) { // this is equivalent to blockSize_ > 0
00783       tmp = X_;
00784     }
00785     else {
00786       tmp = problem_->getInitVec();
00787       TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00788                          "Anasazi::BlockDavidson::setSize(): eigenproblem did not specify initial vectors to clone from.");
00789     }
00790 
00791     TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*numBlocks > MVT::GetVecLength(*tmp),std::invalid_argument,
00792                        "Anasazi::BlockDavidson::setSize(): max subspace dimension and auxilliary subspace too large.");
00793 
00794 
00796     // blockSize dependent
00797     //
00798     // grow/allocate vectors
00799     Rnorms_.resize(blockSize_,NANVAL);
00800     R2norms_.resize(blockSize_,NANVAL);
00801     //
00802     // clone multivectors off of tmp
00803     //
00804     // free current allocation first, to make room for new allocation
00805     X_ = Teuchos::null;
00806     KX_ = Teuchos::null;
00807     MX_ = Teuchos::null;
00808     R_ = Teuchos::null;
00809     V_ = Teuchos::null;
00810 
00811     om_->print(Debug," >> Allocating X_\n");
00812     X_ = MVT::Clone(*tmp,blockSize_);
00813     om_->print(Debug," >> Allocating KX_\n");
00814     KX_ = MVT::Clone(*tmp,blockSize_);
00815     if (hasM_) {
00816       om_->print(Debug," >> Allocating MX_\n");
00817       MX_ = MVT::Clone(*tmp,blockSize_);
00818     }
00819     else {
00820       MX_ = X_;
00821     }
00822     om_->print(Debug," >> Allocating R_\n");
00823     R_ = MVT::Clone(*tmp,blockSize_);
00824 
00826     // blockSize*numBlocks dependent
00827     //
00828     int newsd = blockSize_*numBlocks_;
00829     theta_.resize(blockSize_*numBlocks_,NANVAL);
00830     om_->print(Debug," >> Allocating V_\n");
00831     V_ = MVT::Clone(*tmp,newsd);
00832     KK_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
00833 
00834     om_->print(Debug," >> done allocating.\n");
00835 
00836     initialized_ = false;
00837     curDim_ = 0;
00838   }
00839 
00840 
00842   // Set the auxiliary vectors
00843   template <class ScalarType, class MV, class OP>
00844   void BlockDavidson<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00845     typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
00846 
00847     // set new auxiliary vectors
00848     auxVecs_ = auxvecs;
00849     numAuxVecs_ = 0;
00850     for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
00851       numAuxVecs_ += MVT::GetNumberVecs(**i);
00852     }
00853 
00854     // If the solver has been initialized, V is not necessarily orthogonal to new auxiliary vectors
00855     if (numAuxVecs_ > 0 && initialized_) {
00856       initialized_ = false;
00857     }
00858 
00859     if (om_->isVerbosity( Debug ) ) {
00860       CheckList chk;
00861       chk.checkQ = true;
00862       om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
00863     }
00864   }
00865 
00866 
00868   /* Initialize the state of the solver
00869    * 
00870    * POST-CONDITIONS:
00871    *
00872    * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
00873    * theta_ contains Ritz w.r.t. V_(1:curDim_)
00874    * X is Ritz vectors w.r.t. V_(1:curDim_)
00875    * KX = Op*X
00876    * MX = M*X if hasM_
00877    * R = KX - MX*diag(theta_)
00878    *
00879    */
00880   template <class ScalarType, class MV, class OP>
00881   void BlockDavidson<ScalarType,MV,OP>::initialize(BlockDavidsonState<ScalarType,MV> newstate)
00882   {
00883     // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
00884     // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
00885 
00886     Teuchos::TimeMonitor inittimer( *timerInit_ );
00887 
00888     std::vector<int> bsind(blockSize_);
00889     for (int i=0; i<blockSize_; ++i) bsind[i] = i;
00890 
00891     Teuchos::BLAS<int,ScalarType> blas;
00892 
00893     // in BlockDavidson, V is primary
00894     // the order of dependence follows like so.
00895     // --init->               V,KK
00896     //    --ritz analysis->   theta,X  
00897     //       --op apply->     KX,MX  
00898     //          --compute->   R
00899     // 
00900     // if the user specifies all data for a level, we will accept it.
00901     // otherwise, we will generate the whole level, and all subsequent levels.
00902     //
00903     // the data members are ordered based on dependence, and the levels are
00904     // partitioned according to the amount of work required to produce the
00905     // items in a level.
00906     //
00907     // inconsistent multivectors widths and lengths will not be tolerated, and
00908     // will be treated with exceptions.
00909     //
00910     // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to 
00911     // multivectors in the solver, the copy will not be affected.
00912 
00913     // set up V and KK: get them from newstate if user specified them
00914     // otherwise, set them manually
00915     Teuchos::RCP<MV> lclV;
00916     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK;
00917 
00918     if (newstate.V != Teuchos::null && newstate.KK != Teuchos::null) {
00919       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_), std::invalid_argument, 
00920                           "Anasazi::BlockDavidson::initialize(newstate): Vector length of V not correct." );
00921       TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument, 
00922                           "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be at least blockSize().");
00923       TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument, 
00924                           "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
00925       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument, 
00926                           "Anasazi::BlockDavidson::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
00927 
00928       curDim_ = newstate.curDim;
00929       // pick an integral amount
00930       curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
00931 
00932       TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument, 
00933                           "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
00934 
00935       // check size of KK
00936       TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument, 
00937                           "Anasazi::BlockDavidson::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
00938 
00939       // put data in V
00940       std::vector<int> nevind(curDim_);
00941       for (int i=0; i<curDim_; ++i) nevind[i] = i;
00942       if (newstate.V != V_) {
00943         if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
00944           newstate.V = MVT::CloneView(*newstate.V,nevind);
00945         }
00946         MVT::SetBlock(*newstate.V,nevind,*V_);
00947       }
00948       lclV = MVT::CloneViewNonConst(*V_,nevind);
00949 
00950       // put data into KK_
00951       lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
00952       if (newstate.KK != KK_) {
00953         if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
00954           newstate.KK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
00955         }
00956         lclKK->assign(*newstate.KK);
00957       }
00958       //
00959       // make lclKK Hermitian in memory (copy the upper half to the lower half)
00960       for (int j=0; j<curDim_-1; ++j) {
00961         for (int i=j+1; i<curDim_; ++i) {
00962           (*lclKK)(i,j) = SCT::conjugate((*lclKK)(j,i));
00963         }
00964       }
00965     }
00966     else {
00967       // user did not specify one of V or KK
00968       // get vectors from problem or generate something, projectAndNormalize
00969       Teuchos::RCP<const MV> ivec = problem_->getInitVec();
00970       TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
00971                          "Anasazi::BlockDavdison::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
00972       // clear newstate so we won't use any data from it below
00973       newstate.X      = Teuchos::null;
00974       newstate.MX     = Teuchos::null;
00975       newstate.KX     = Teuchos::null;
00976       newstate.R      = Teuchos::null;
00977       newstate.H      = Teuchos::null;
00978       newstate.T      = Teuchos::null;
00979       newstate.KK     = Teuchos::null;
00980       newstate.V      = Teuchos::null;
00981       newstate.curDim = 0;
00982 
00983       curDim_ = MVT::GetNumberVecs(*ivec);
00984       // pick the largest multiple of blockSize_
00985       curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
00986       if (curDim_ > blockSize_*numBlocks_) {
00987         // user specified too many vectors... truncate
00988         // this produces a full subspace, but that is okay
00989         curDim_ = blockSize_*numBlocks_;
00990       }
00991       bool userand = false;
00992       if (curDim_ == 0) {
00993         // we need at least blockSize_ vectors
00994         // use a random multivec: ignore everything from InitVec
00995         userand = true;
00996         curDim_ = blockSize_;
00997       }
00998 
00999       // get pointers into V,KV,MV
01000       // tmpVecs will be used below for M*V and K*V (not simultaneously)
01001       // lclV has curDim vectors
01002       // if there is space for lclV and tmpVecs in V_, point tmpVecs into V_
01003       // otherwise, we must allocate space for these products
01004       //
01005       // get pointer to first curDim vector in V_
01006       std::vector<int> dimind(curDim_);
01007       for (int i=0; i<curDim_; ++i) dimind[i] = i;
01008       lclV = MVT::CloneViewNonConst(*V_,dimind);
01009       if (userand) {
01010         // generate random vector data
01011         MVT::MvRandom(*lclV);
01012       }
01013       else {
01014         if (MVT::GetNumberVecs(*ivec) > curDim_) {
01015           ivec = MVT::CloneView(*ivec,dimind);
01016         }
01017         // assign ivec to first part of lclV
01018         MVT::SetBlock(*ivec,dimind,*lclV);
01019       }
01020       Teuchos::RCP<MV> tmpVecs; 
01021       if (curDim_*2 <= blockSize_*numBlocks_) {
01022         // partition V_ = [lclV tmpVecs _leftover_]
01023         std::vector<int> block2(curDim_);
01024         for (int i=0; i<curDim_; ++i) block2[i] = curDim_+i;
01025         tmpVecs = MVT::CloneViewNonConst(*V_,block2);
01026       }
01027       else {
01028         // allocate space for tmpVecs
01029         tmpVecs = MVT::Clone(*V_,curDim_);
01030       }
01031 
01032       // compute M*lclV if hasM_
01033       if (hasM_) {
01034         Teuchos::TimeMonitor lcltimer( *timerMOp_ );
01035         OPT::Apply(*MOp_,*lclV,*tmpVecs);
01036         count_ApplyM_ += curDim_;
01037       }
01038 
01039       // remove auxVecs from lclV and normalize it
01040       if (auxVecs_.size() > 0) {
01041         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01042 
01043         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01044         int rank = orthman_->projectAndNormalizeMat(*lclV,auxVecs_,dummyC,Teuchos::null,tmpVecs);
01045         TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
01046                            "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
01047       }
01048       else {
01049         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01050 
01051         int rank = orthman_->normalizeMat(*lclV,Teuchos::null,tmpVecs);
01052         TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
01053                            "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
01054       }
01055 
01056       // compute K*lclV: we are re-using tmpVecs to store the result
01057       {
01058         Teuchos::TimeMonitor lcltimer( *timerOp_ );
01059         OPT::Apply(*Op_,*lclV,*tmpVecs);
01060         count_ApplyOp_ += curDim_;
01061       }
01062 
01063       // generate KK
01064       lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
01065       MVT::MvTransMv(ONE,*lclV,*tmpVecs,*lclKK);
01066 
01067       // clear tmpVecs
01068       tmpVecs = Teuchos::null;
01069     }
01070 
01071     // X,theta require Ritz analysis; if we have to generate one of these, we might as well generate both
01072     if (newstate.X != Teuchos::null && newstate.T != Teuchos::null) {
01073       TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetVecLength(*newstate.X) != MVT::GetVecLength(*X_),
01074                           std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of X must be consistent with block size and length of V.");
01075       TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
01076                           std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of T must be consistent with dimension of V.");
01077 
01078       if (newstate.X != X_) {
01079         MVT::SetBlock(*newstate.X,bsind,*X_);
01080       }
01081 
01082       std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
01083     }
01084     else {
01085       // compute ritz vecs/vals
01086       Teuchos::SerialDenseMatrix<int,ScalarType> S(curDim_,curDim_);
01087       {
01088         Teuchos::TimeMonitor lcltimer( *timerDS_ );
01089         int rank = curDim_;
01090         Utils::directSolver(curDim_, *lclKK, Teuchos::null, S, theta_, rank, 10);
01091         // we want all ritz values back
01092         TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
01093                            "Anasazi::BlockDavidson::initialize(newstate): Not enough Ritz vectors to initialize algorithm.");
01094       }
01095       // sort ritz pairs
01096       {
01097         Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
01098 
01099         std::vector<int> order(curDim_);
01100         //
01101         // sort the first curDim_ values in theta_
01102         sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_);   // don't catch exception
01103         //
01104         // apply the same ordering to the primitive ritz vectors
01105         Utils::permuteVectors(order,S);
01106       }
01107 
01108       // compute eigenvectors
01109       Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,curDim_,blockSize_);
01110       {
01111         Teuchos::TimeMonitor lcltimer( *timerLocal_ );
01112 
01113         // X <- lclV*S
01114         MVT::MvTimesMatAddMv( ONE, *lclV, S1, ZERO, *X_ );
01115       }
01116       // we generated theta,X so we don't want to use the user's KX,MX
01117       newstate.KX = Teuchos::null;
01118       newstate.MX = Teuchos::null;
01119     }
01120 
01121     // done with local pointers
01122     lclV = Teuchos::null;
01123     lclKK = Teuchos::null;
01124 
01125     // set up KX
01126     if ( newstate.KX != Teuchos::null ) {
01127       TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_,
01128                           std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
01129       TEST_FOR_EXCEPTION(MVT::GetVecLength(*newstate.KX) != MVT::GetVecLength(*X_),
01130                           std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
01131       if (newstate.KX != KX_) {
01132         MVT::SetBlock(*newstate.KX,bsind,*KX_);
01133       }
01134     }
01135     else {
01136       // generate KX
01137       {
01138         Teuchos::TimeMonitor lcltimer( *timerOp_ );
01139         OPT::Apply(*Op_,*X_,*KX_);
01140         count_ApplyOp_ += blockSize_;
01141       }
01142       // we generated KX; we will generate R as well
01143       newstate.R = Teuchos::null;
01144     }
01145 
01146     // set up MX
01147     if (hasM_) {
01148       if ( newstate.MX != Teuchos::null ) {
01149         TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_,
01150                             std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
01151         TEST_FOR_EXCEPTION(MVT::GetVecLength(*newstate.MX) != MVT::GetVecLength(*X_),
01152                             std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
01153         if (newstate.MX != MX_) {
01154           MVT::SetBlock(*newstate.MX,bsind,*MX_);
01155         }
01156       }
01157       else {
01158         // generate MX
01159         {
01160           Teuchos::TimeMonitor lcltimer( *timerOp_ );
01161           OPT::Apply(*MOp_,*X_,*MX_);
01162           count_ApplyOp_ += blockSize_;
01163         }
01164         // we generated MX; we will generate R as well
01165         newstate.R = Teuchos::null;
01166       }
01167     }
01168     else {
01169       // the assignment MX_==X_ would be redundant; take advantage of this opportunity to debug a little
01170       TEST_FOR_EXCEPTION(MX_ != X_, std::logic_error, "Anasazi::BlockDavidson::initialize(): solver invariant not satisfied (MX==X).");
01171     }
01172 
01173     // set up R
01174     if (newstate.R != Teuchos::null) {
01175       TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
01176                           std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
01177       TEST_FOR_EXCEPTION(MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*X_),
01178                           std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
01179       if (newstate.R != R_) {
01180         MVT::SetBlock(*newstate.R,bsind,*R_);
01181       }
01182     }
01183     else {
01184       Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
01185       
01186       // form R <- KX - MX*T
01187       MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
01188       Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
01189       T.putScalar(ZERO);
01190       for (int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
01191       MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
01192 
01193     }
01194 
01195     // R has been updated; mark the norms as out-of-date
01196     Rnorms_current_ = false;
01197     R2norms_current_ = false;
01198 
01199     // finally, we are initialized
01200     initialized_ = true;
01201 
01202     if (om_->isVerbosity( Debug ) ) {
01203       // Check almost everything here
01204       CheckList chk;
01205       chk.checkV = true;
01206       chk.checkX = true;
01207       chk.checkKX = true;
01208       chk.checkMX = true;
01209       chk.checkR = true;
01210       chk.checkQ = true;
01211       chk.checkKK = true;
01212       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
01213     }
01214 
01215     // Print information on current status
01216     if (om_->isVerbosity(Debug)) {
01217       currentStatus( om_->stream(Debug) );
01218     }
01219     else if (om_->isVerbosity(IterationDetails)) {
01220       currentStatus( om_->stream(IterationDetails) );
01221     }
01222   }
01223 
01224 
01226   // initialize the solver with default state
01227   template <class ScalarType, class MV, class OP>
01228   void BlockDavidson<ScalarType,MV,OP>::initialize()
01229   {
01230     BlockDavidsonState<ScalarType,MV> empty;
01231     initialize(empty);
01232   }
01233 
01234 
01236   // Perform BlockDavidson iterations until the StatusTest tells us to stop.
01237   template <class ScalarType, class MV, class OP>
01238   void BlockDavidson<ScalarType,MV,OP>::iterate ()
01239   {
01240     //
01241     // Initialize solver state
01242     if (initialized_ == false) {
01243       initialize();
01244     }
01245 
01246     // as a data member, this would be redundant and require synchronization with
01247     // blockSize_ and numBlocks_; we'll use a constant here.
01248     const int searchDim = blockSize_*numBlocks_;
01249 
01250     Teuchos::BLAS<int,ScalarType> blas;
01251 
01252     //
01253     // The projected matrices are part of the state, but the eigenvectors are defined locally.
01254     //    S = Local eigenvectors         (size: searchDim * searchDim
01255     Teuchos::SerialDenseMatrix<int,ScalarType> S( searchDim, searchDim );
01256 
01257 
01259     // iterate until the status test tells us to stop.
01260     // also break if our basis is full
01261     while (tester_->checkStatus(this) != Passed && curDim_ < searchDim) {
01262 
01263       // Print information on current iteration
01264       if (om_->isVerbosity(Debug)) {
01265         currentStatus( om_->stream(Debug) );
01266       }
01267       else if (om_->isVerbosity(IterationDetails)) {
01268         currentStatus( om_->stream(IterationDetails) );
01269       }
01270 
01271       ++iter_;
01272 
01273       // get the current part of the basis
01274       std::vector<int> curind(blockSize_);
01275       for (int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
01276       H_ = MVT::CloneViewNonConst(*V_,curind);
01277       
01278       // Apply the preconditioner on the residuals: H <- Prec*R
01279       // H = Prec*R
01280       if (Prec_ != Teuchos::null) {
01281         Teuchos::TimeMonitor lcltimer( *timerPrec_ );
01282         OPT::Apply( *Prec_, *R_, *H_ );   // don't catch the exception
01283         count_ApplyPrec_ += blockSize_;
01284       }
01285       else {
01286         std::vector<int> bsind(blockSize_);
01287         for (int i=0; i<blockSize_; ++i) { bsind[i] = i; }
01288         MVT::SetBlock(*R_,bsind,*H_);
01289       }
01290 
01291       // Apply the mass matrix on H
01292       if (hasM_) {
01293         // use memory at MX_ for temporary storage
01294         MH_ = MX_;
01295         Teuchos::TimeMonitor lcltimer( *timerMOp_ );
01296         OPT::Apply( *MOp_, *H_, *MH_);    // don't catch the exception
01297         count_ApplyM_ += blockSize_;
01298       }
01299       else  {
01300         MH_ = H_;
01301       }
01302 
01303       // Get a view of the previous vectors
01304       // this is used for orthogonalization and for computing V^H K H
01305       std::vector<int> prevind(curDim_);
01306       for (int i=0; i<curDim_; ++i) prevind[i] = i;
01307       Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,prevind);
01308 
01309       // Orthogonalize H against the previous vectors and the auxiliary vectors, and normalize
01310       {
01311         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01312 
01313         Teuchos::Array<Teuchos::RCP<const MV> > against = auxVecs_;
01314         against.push_back(Vprev);
01315         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01316         int rank = orthman_->projectAndNormalizeMat(*H_,against,
01317                                             dummyC,
01318                                             Teuchos::null,MH_);
01319         TEST_FOR_EXCEPTION(rank != blockSize_,BlockDavidsonOrthoFailure,
01320                            "Anasazi::BlockDavidson::iterate(): unable to compute orthonormal basis for H.");
01321       }
01322 
01323       // Apply the stiffness matrix to H
01324       {
01325         // use memory at KX_ for temporary storage
01326         KH_ = KX_;
01327         Teuchos::TimeMonitor lcltimer( *timerOp_ );
01328         OPT::Apply( *Op_, *H_, *KH_);    // don't catch the exception
01329         count_ApplyOp_ += blockSize_;
01330       }
01331 
01332       if (om_->isVerbosity( Debug ) ) {
01333         CheckList chk;
01334         chk.checkH = true;
01335         chk.checkMH = true;
01336         chk.checkKH = true;
01337         om_->print( Debug, accuracyCheck(chk, ": after ortho H") );
01338       }
01339       else if (om_->isVerbosity( OrthoDetails ) ) {
01340         CheckList chk;
01341         chk.checkH = true;
01342         chk.checkMH = true;
01343         chk.checkKH = true;
01344         om_->print( OrthoDetails, accuracyCheck(chk,": after ortho H") );
01345       }
01346 
01347       // compute next part of the projected matrices
01348       // this this in two parts
01349       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > nextKK;
01350       // Vprev*K*H
01351       nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_) );
01352       MVT::MvTransMv(ONE,*Vprev,*KH_,*nextKK);
01353       // H*K*H
01354       nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
01355       MVT::MvTransMv(ONE,*H_,*KH_,*nextKK);
01356       // 
01357       // make sure that KK_ is Hermitian in memory
01358       nextKK = Teuchos::null;
01359       for (int i=curDim_; i<curDim_+blockSize_; ++i) {
01360         for (int j=0; j<i; ++j) {
01361           (*KK_)(i,j) = SCT::conjugate((*KK_)(j,i));
01362         }
01363       }
01364 
01365       // V has been extended, and KK has been extended. Update basis dim and release all pointers.
01366       curDim_ += blockSize_;
01367       H_ = KH_ = MH_ = Teuchos::null;
01368       Vprev = Teuchos::null;
01369 
01370       if (om_->isVerbosity( Debug ) ) {
01371         CheckList chk;
01372         chk.checkKK = true;
01373         om_->print( Debug, accuracyCheck(chk, ": after expanding KK") );
01374       }
01375 
01376       // Get pointer to complete basis
01377       curind.resize(curDim_);
01378       for (int i=0; i<curDim_; ++i) curind[i] = i;
01379       Teuchos::RCP<const MV> curV = MVT::CloneView(*V_,curind);
01380 
01381       // Perform spectral decomposition
01382       {
01383         Teuchos::TimeMonitor lcltimer(*timerDS_);
01384         int nevlocal = curDim_;
01385         int info = Utils::directSolver(curDim_,*KK_,Teuchos::null,S,theta_,nevlocal,10);
01386         TEST_FOR_EXCEPTION(info != 0,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve returned error code.");
01387         // we did not ask directSolver to perform deflation, so nevLocal better be curDim_
01388         TEST_FOR_EXCEPTION(nevlocal != curDim_,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve did not compute all eigenvectors."); // this should never happen
01389       }
01390 
01391       // Sort ritz pairs
01392       { 
01393         Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
01394 
01395         std::vector<int> order(curDim_);
01396         // 
01397         // sort the first curDim_ values in theta_
01398         sm_->sort(theta_, Teuchos::rcp(&order,false), curDim_);   // don't catch exception
01399         //
01400         // apply the same ordering to the primitive ritz vectors
01401         Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,curDim_,curDim_);
01402         Utils::permuteVectors(order,curS);
01403       }
01404 
01405       // Create a view matrix of the first blockSize_ vectors
01406       Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
01407 
01408       // Compute the new Ritz vectors
01409       {
01410         Teuchos::TimeMonitor lcltimer( *timerLocal_ );
01411         MVT::MvTimesMatAddMv(ONE,*curV,S1,ZERO,*X_);
01412       }
01413 
01414       // Apply the stiffness matrix for the Ritz vectors
01415       {
01416         Teuchos::TimeMonitor lcltimer( *timerOp_ );
01417         OPT::Apply( *Op_, *X_, *KX_);    // don't catch the exception
01418         count_ApplyOp_ += blockSize_;
01419       }
01420       // Apply the mass matrix for the Ritz vectors
01421       if (hasM_) {
01422         Teuchos::TimeMonitor lcltimer( *timerMOp_ );
01423         OPT::Apply(*MOp_,*X_,*MX_);
01424         count_ApplyM_ += blockSize_;
01425       }
01426       else {
01427         MX_ = X_;
01428       }
01429 
01430       // Compute the residual
01431       // R = KX - MX*diag(theta)
01432       {
01433         Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
01434         
01435         MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
01436         Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
01437         for (int i = 0; i < blockSize_; ++i) {
01438           T(i,i) = theta_[i];
01439         }
01440         MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
01441       }
01442 
01443       // R has been updated; mark the norms as out-of-date
01444       Rnorms_current_ = false;
01445       R2norms_current_ = false;
01446 
01447 
01448       // When required, monitor some orthogonalities
01449       if (om_->isVerbosity( Debug ) ) {
01450         // Check almost everything here
01451         CheckList chk;
01452         chk.checkV = true;
01453         chk.checkX = true;
01454         chk.checkKX = true;
01455         chk.checkMX = true;
01456         chk.checkR = true;
01457         om_->print( Debug, accuracyCheck(chk, ": after local update") );
01458       }
01459       else if (om_->isVerbosity( OrthoDetails )) {
01460         CheckList chk;
01461         chk.checkX = true;
01462         chk.checkKX = true;
01463         chk.checkMX = true;
01464         chk.checkR = true;
01465         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
01466       }
01467     } // end while (statusTest == false)
01468 
01469   } // end of iterate()
01470 
01471 
01472 
01474   // compute/return residual M-norms
01475   template <class ScalarType, class MV, class OP>
01476   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01477   BlockDavidson<ScalarType,MV,OP>::getResNorms() {
01478     if (Rnorms_current_ == false) {
01479       // Update the residual norms
01480       orthman_->norm(*R_,Rnorms_);
01481       Rnorms_current_ = true;
01482     }
01483     return Rnorms_;
01484   }
01485 
01486   
01488   // compute/return residual 2-norms
01489   template <class ScalarType, class MV, class OP>
01490   std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
01491   BlockDavidson<ScalarType,MV,OP>::getRes2Norms() {
01492     if (R2norms_current_ == false) {
01493       // Update the residual 2-norms 
01494       MVT::MvNorm(*R_,R2norms_);
01495       R2norms_current_ = true;
01496     }
01497     return R2norms_;
01498   }
01499 
01501   // Set a new StatusTest for the solver.
01502   template <class ScalarType, class MV, class OP>
01503   void BlockDavidson<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
01504     TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
01505         "Anasazi::BlockDavidson::setStatusTest() was passed a null StatusTest.");
01506     tester_ = test;
01507   }
01508 
01510   // Get the current StatusTest used by the solver.
01511   template <class ScalarType, class MV, class OP>
01512   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockDavidson<ScalarType,MV,OP>::getStatusTest() const {
01513     return tester_;
01514   }
01515   
01516 
01518   // Check accuracy, orthogonality, and other debugging stuff
01519   // 
01520   // bools specify which tests we want to run (instead of running more than we actually care about)
01521   //
01522   // we don't bother checking the following because they are computed explicitly:
01523   //    H == Prec*R
01524   //   KH == K*H
01525   //
01526   // 
01527   // checkV : V orthonormal
01528   //          orthogonal to auxvecs
01529   // checkX : X orthonormal
01530   //          orthogonal to auxvecs
01531   // checkMX: check MX == M*X
01532   // checkKX: check KX == K*X
01533   // checkH : H orthonormal 
01534   //          orthogonal to V and H and auxvecs
01535   // checkMH: check MH == M*H
01536   // checkR : check R orthogonal to X
01537   // checkQ : check that auxiliary vectors are actually orthonormal
01538   // checkKK: check that KK is symmetric in memory 
01539   //
01540   // TODO: 
01541   //  add checkTheta 
01542   //
01543   template <class ScalarType, class MV, class OP>
01544   std::string BlockDavidson<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const 
01545   {
01546     using std::endl;
01547 
01548     std::stringstream os;
01549     os.precision(2);
01550     os.setf(std::ios::scientific, std::ios::floatfield);
01551 
01552     os << " Debugging checks: iteration " << iter_ << where << endl;
01553 
01554     // V and friends
01555     std::vector<int> lclind(curDim_);
01556     for (int i=0; i<curDim_; ++i) lclind[i] = i;
01557     Teuchos::RCP<const MV> lclV;
01558     if (initialized_) {
01559       lclV = MVT::CloneView(*V_,lclind);
01560     }
01561     if (chk.checkV && initialized_) {
01562       MagnitudeType err = orthman_->orthonormError(*lclV);
01563       os << " >> Error in V^H M V == I  : " << err << endl;
01564       for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
01565         err = orthman_->orthogError(*lclV,*auxVecs_[i]);
01566         os << " >> Error in V^H M Q[" << i << "] == 0 : " << err << endl;
01567       }
01568       Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
01569       Teuchos::RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
01570       OPT::Apply(*Op_,*lclV,*lclKV);
01571       MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
01572       Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
01573       curKK -= subKK;
01574       // dup the lower tri part
01575       for (int j=0; j<curDim_; ++j) {
01576         for (int i=j+1; i<curDim_; ++i) {
01577           curKK(i,j) = curKK(j,i);
01578         }
01579       }
01580       os << " >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
01581     }
01582 
01583     // X and friends
01584     if (chk.checkX && initialized_) {
01585       MagnitudeType err = orthman_->orthonormError(*X_);
01586       os << " >> Error in X^H M X == I  : " << err << endl;
01587       for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
01588         err = orthman_->orthogError(*X_,*auxVecs_[i]);
01589         os << " >> Error in X^H M Q[" << i << "] == 0 : " << err << endl;
01590       }
01591     }
01592     if (chk.checkMX && hasM_ && initialized_) {
01593       MagnitudeType err = Utils::errorEquality(*X_, *MX_, MOp_);
01594       os << " >> Error in MX == M*X     : " << err << endl;
01595     }
01596     if (chk.checkKX && initialized_) {
01597       MagnitudeType err = Utils::errorEquality(*X_, *KX_, Op_);
01598       os << " >> Error in KX == K*X     : " << err << endl;
01599     }
01600 
01601     // H and friends
01602     if (chk.checkH && initialized_) {
01603       MagnitudeType err = orthman_->orthonormError(*H_);
01604       os << " >> Error in H^H M H == I  : " << err << endl;
01605       err = orthman_->orthogError(*H_,*lclV);
01606       os << " >> Error in H^H M V == 0  : " << err << endl;
01607       err = orthman_->orthogError(*H_,*X_);
01608       os << " >> Error in H^H M X == 0  : " << err << endl;
01609       for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
01610         err = orthman_->orthogError(*H_,*auxVecs_[i]);
01611         os << " >> Error in H^H M Q[" << i << "] == 0 : " << err << endl;
01612       }
01613     }
01614     if (chk.checkKH && initialized_) {
01615       MagnitudeType err = Utils::errorEquality(*H_, *KH_, Op_);
01616       os << " >> Error in KH == K*H     : " << err << endl;
01617     }
01618     if (chk.checkMH && hasM_ && initialized_) {
01619       MagnitudeType err = Utils::errorEquality(*H_, *MH_, MOp_);
01620       os << " >> Error in MH == M*H     : " << err << endl;
01621     }
01622 
01623     // R: this is not M-orthogonality, but standard Euclidean orthogonality
01624     if (chk.checkR && initialized_) {
01625       Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01626       MVT::MvTransMv(ONE,*X_,*R_,xTx);
01627       MagnitudeType err = xTx.normFrobenius();
01628       os << " >> Error in X^H R == 0    : " << err << endl;
01629     }
01630 
01631     // KK
01632     if (chk.checkKK && initialized_) {
01633       Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_), lclKK(Teuchos::View,*KK_,curDim_,curDim_);
01634       for (int j=0; j<curDim_; ++j) {
01635         for (int i=0; i<curDim_; ++i) {
01636           SDMerr(i,j) = lclKK(i,j) - SCT::conjugate(lclKK(j,i));
01637         }
01638       }
01639       os << " >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
01640     }
01641 
01642     // Q
01643     if (chk.checkQ) {
01644       for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
01645         MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
01646         os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << err << endl;
01647         for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
01648           err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01649           os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << err << endl;
01650         }
01651       }
01652     }
01653 
01654     os << endl;
01655 
01656     return os.str();
01657   }
01658 
01659 
01661   // Print the current status of the solver
01662   template <class ScalarType, class MV, class OP>
01663   void 
01664   BlockDavidson<ScalarType,MV,OP>::currentStatus(std::ostream &os) 
01665   {
01666     using std::endl;
01667 
01668     os.setf(std::ios::scientific, std::ios::floatfield);
01669     os.precision(6);
01670     os <<endl;
01671     os <<"================================================================================" << endl;
01672     os << endl;
01673     os <<"                          BlockDavidson Solver Status" << endl;
01674     os << endl;
01675     os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01676     os <<"The number of iterations performed is " <<iter_<<endl;
01677     os <<"The block size is         " << blockSize_<<endl;
01678     os <<"The number of blocks is   " << numBlocks_<<endl;
01679     os <<"The current basis size is " << curDim_<<endl;
01680     os <<"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
01681     os <<"The number of operations Op*x   is "<<count_ApplyOp_<<endl;
01682     os <<"The number of operations M*x    is "<<count_ApplyM_<<endl;
01683     os <<"The number of operations Prec*x is "<<count_ApplyPrec_<<endl;
01684 
01685     os.setf(std::ios_base::right, std::ios_base::adjustfield);
01686 
01687     if (initialized_) {
01688       os << endl;
01689       os <<"CURRENT EIGENVALUE ESTIMATES             "<<endl;
01690       os << std::setw(20) << "Eigenvalue" 
01691          << std::setw(20) << "Residual(M)"
01692          << std::setw(20) << "Residual(2)"
01693          << endl;
01694       os <<"--------------------------------------------------------------------------------"<<endl;
01695       for (int i=0; i<blockSize_; ++i) {
01696         os << std::setw(20) << theta_[i];
01697         if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
01698         else os << std::setw(20) << "not current";
01699         if (R2norms_current_) os << std::setw(20) << R2norms_[i];
01700         else os << std::setw(20) << "not current";
01701         os << endl;
01702       }
01703     }
01704     os <<"================================================================================" << endl;
01705     os << endl;
01706   }
01707   
01708 } // End of namespace Anasazi
01709 
01710 #endif
01711 
01712 // End of file AnasaziBlockDavidson.hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends