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