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