Anasazi Version of the Day
AnasaziBlockKrylovSchur.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_KRYLOV_SCHUR_HPP
00034 #define ANASAZI_BLOCK_KRYLOV_SCHUR_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 #include "AnasaziHelperTraits.hpp"
00043 
00044 #include "AnasaziOrthoManager.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 
00066 namespace Anasazi {
00067 
00069 
00070 
00075   template <class ScalarType, class MulVec>
00076   struct BlockKrylovSchurState {
00081     int curDim;
00083     Teuchos::RCP<const MulVec> V;
00089     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > H;
00091     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > S;
00093     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > Q;
00094     BlockKrylovSchurState() : curDim(0), V(Teuchos::null),
00095                               H(Teuchos::null), S(Teuchos::null),
00096                               Q(Teuchos::null) {}
00097   };
00098 
00100 
00102 
00103 
00116   class BlockKrylovSchurInitFailure : public AnasaziError {public:
00117     BlockKrylovSchurInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
00118     {}};
00119 
00126   class BlockKrylovSchurOrthoFailure : public AnasaziError {public:
00127     BlockKrylovSchurOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
00128     {}};
00129   
00131 
00132 
00133   template <class ScalarType, class MV, class OP>
00134   class BlockKrylovSchur : public Eigensolver<ScalarType,MV,OP> { 
00135   public:
00137 
00138     
00150     BlockKrylovSchur( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 
00151                       const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00152                       const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00153                       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00154                       const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
00155                       Teuchos::ParameterList &params 
00156                     );
00157     
00159     virtual ~BlockKrylovSchur() {};
00161 
00162 
00164 
00165     
00187     void iterate();
00188 
00210     void initialize(BlockKrylovSchurState<ScalarType,MV> state);
00211 
00215     void initialize();
00216 
00225     bool isInitialized() const { return initialized_; }
00226 
00234     BlockKrylovSchurState<ScalarType,MV> getState() const {
00235       BlockKrylovSchurState<ScalarType,MV> state;
00236       state.curDim = curDim_;
00237       state.V = V_;
00238       state.H = H_;
00239       state.Q = Q_;
00240       state.S = schurH_;
00241       return state;
00242     }
00243     
00245 
00246 
00248 
00249 
00251     int getNumIters() const { return(iter_); }
00252 
00254     void resetNumIters() { iter_=0; }
00255 
00263     Teuchos::RCP<const MV> getRitzVectors() { return ritzVectors_; }
00264 
00272     std::vector<Value<ScalarType> > getRitzValues() { 
00273       std::vector<Value<ScalarType> > ret = ritzValues_;
00274       ret.resize(ritzIndex_.size());
00275       return ret;
00276     }
00277 
00283     std::vector<int> getRitzIndex() { return ritzIndex_; }
00284 
00289     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms() {
00290       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
00291       return ret;
00292     }
00293 
00298     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms() {
00299       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
00300       return ret;
00301     }
00302 
00307     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms() {
00308       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
00309       ret.resize(ritzIndex_.size());
00310       return ret;
00311     }
00312 
00314 
00316 
00317 
00319     void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
00320 
00322     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
00323 
00325     const Eigenproblem<ScalarType,MV,OP>& getProblem() const { return(*problem_); };
00326 
00333     void setSize(int blockSize, int numBlocks);
00334 
00336     void setBlockSize(int blockSize);
00337 
00339     void setStepSize(int stepSize);
00340 
00342     void setNumRitzVectors(int numRitzVecs);
00343 
00345     int getStepSize() const { return(stepSize_); }
00346 
00348     int getBlockSize() const { return(blockSize_); }
00349 
00351     int getNumRitzVectors() const { return(numRitzVecs_); }
00352 
00358     int getCurSubspaceDim() const {
00359       if (!initialized_) return 0;
00360       return curDim_;
00361     }
00362 
00364     int getMaxSubspaceDim() const { return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
00365 
00366 
00379     void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
00380 
00382     Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const {return auxVecs_;}
00383 
00385 
00387 
00388     
00390     void currentStatus(std::ostream &os);
00391 
00393 
00395 
00396     
00398     bool isRitzVecsCurrent() const { return ritzVecsCurrent_; }
00399 
00401     bool isRitzValsCurrent() const { return ritzValsCurrent_; }
00402     
00404     bool isSchurCurrent() const { return schurCurrent_; }
00405     
00407 
00409 
00410     
00412     void computeRitzVectors();
00413 
00415     void computeRitzValues();
00416     
00418     void computeSchurForm( const bool sort = true );
00419 
00421 
00422   private:
00423     //
00424     // Convenience typedefs
00425     //
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     typedef typename std::vector<ScalarType>::iterator STiter;
00431     typedef typename std::vector<MagnitudeType>::iterator MTiter;
00432     const MagnitudeType MT_ONE;  
00433     const MagnitudeType MT_ZERO; 
00434     const MagnitudeType NANVAL;
00435     const ScalarType ST_ONE;
00436     const ScalarType ST_ZERO;
00437     //
00438     // Internal structs
00439     //
00440     struct CheckList {
00441       bool checkV;
00442       bool checkArn;
00443       bool checkAux;
00444       CheckList() : checkV(false), checkArn(false), checkAux(false) {};
00445     };
00446     //
00447     // Internal methods
00448     //
00449     std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
00450     void sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
00451                         Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
00452                         std::vector<int>& order );
00453     //
00454     // Classes inputed through constructor that define the eigenproblem to be solved.
00455     //
00456     const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> >     problem_;
00457     const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
00458     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00459     Teuchos::RCP<StatusTest<ScalarType,MV,OP> >             tester_;
00460     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        orthman_;
00461     //
00462     // Information obtained from the eigenproblem
00463     //
00464     Teuchos::RCP<const OP> Op_;
00465     //
00466     // Internal timers
00467     //
00468     Teuchos::RCP<Teuchos::Time> timerOp_, timerSortRitzVal_,
00469                                         timerCompSF_, timerSortSF_,
00470                                         timerCompRitzVec_, timerOrtho_;
00471     //
00472     // Counters
00473     //
00474     int count_ApplyOp_;
00475 
00476     //
00477     // Algorithmic parameters.
00478     //
00479     // blockSize_ is the solver block size; it controls the number of eigenvectors that 
00480     // we compute, the number of residual vectors that we compute, and therefore the number
00481     // of vectors added to the basis on each iteration.
00482     int blockSize_;
00483     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00484     int numBlocks_; 
00485     // stepSize_ dictates how many iterations are performed before eigenvectors and eigenvalues
00486     // are computed again
00487     int stepSize_;
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: for Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_
00499     //   for non-Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_ + 1
00500     // this also tells us how many of the values in _theta are valid Ritz values
00501     int curDim_;
00502     //
00503     // State Multivecs
00504     Teuchos::RCP<MV> ritzVectors_, V_;
00505     int numRitzVecs_;
00506     //
00507     // Projected matrices
00508     // H_ : Projected matrix from the Krylov-Schur factorization AV = VH + FB^T
00509     //
00510     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00511     // 
00512     // Schur form of Projected matrices (these are only updated when the Ritz values/vectors are updated).
00513     // schurH_: Schur form reduction of H
00514     // Q_: Schur vectors of H
00515     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > schurH_;
00516     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Q_;
00517     // 
00518     // Auxiliary vectors
00519     Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
00520     int numAuxVecs_;
00521     //
00522     // Number of iterations that have been performed.
00523     int iter_;
00524     //
00525     // State flags
00526     bool ritzVecsCurrent_, ritzValsCurrent_, schurCurrent_;
00527     // 
00528     // Current eigenvalues, residual norms
00529     std::vector<Value<ScalarType> > ritzValues_;
00530     std::vector<MagnitudeType> ritzResiduals_;
00531     // 
00532     // Current index vector for Ritz values and vectors
00533     std::vector<int> ritzIndex_;  // computed by BKS
00534     std::vector<int> ritzOrder_;  // returned from sort manager
00535     //
00536     // Number of Ritz pairs to be printed upon output, if possible
00537     int numRitzPrint_;
00538   };
00539 
00540 
00542   // Constructor
00543   template <class ScalarType, class MV, class OP>
00544   BlockKrylovSchur<ScalarType,MV,OP>::BlockKrylovSchur(
00545         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 
00546         const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00547         const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00548         const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00549         const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
00550         Teuchos::ParameterList &params
00551         ) :
00552     MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
00553     MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00554     NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
00555     ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
00556     ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
00557     // problem, tools
00558     problem_(problem), 
00559     sm_(sorter),
00560     om_(printer),
00561     tester_(tester),
00562     orthman_(ortho),
00563     // timers, counters
00564 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00565     timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Operation Op*x")),
00566     timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Ritz values")),
00567     timerCompSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Schur form")),
00568     timerSortSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Schur form")),
00569     timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Ritz vectors")),
00570     timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Orthogonalization")),
00571 #endif
00572     count_ApplyOp_(0),
00573     // internal data
00574     blockSize_(0),
00575     numBlocks_(0),
00576     stepSize_(0),
00577     initialized_(false),
00578     curDim_(0),
00579     numRitzVecs_(0),
00580     auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ), 
00581     numAuxVecs_(0),
00582     iter_(0),
00583     ritzVecsCurrent_(false),
00584     ritzValsCurrent_(false),
00585     schurCurrent_(false),
00586     numRitzPrint_(0)
00587   {     
00588     TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
00589                        "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
00590     TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
00591                        "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
00592     TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
00593                        "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
00594     TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
00595                        "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
00596     TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
00597                        "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
00598     TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
00599                        "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
00600     TEUCHOS_TEST_FOR_EXCEPTION(sorter == Teuchos::null,std::invalid_argument,
00601                        "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
00602     TEUCHOS_TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
00603                        "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
00604     TEUCHOS_TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
00605                        "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
00606     TEUCHOS_TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
00607                        "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
00608 
00609     // Get problem operator
00610     Op_ = problem_->getOperator();
00611 
00612     // get the step size
00613     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Step Size"), std::invalid_argument,
00614                        "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
00615     int ss = params.get("Step Size",numBlocks_);
00616     setStepSize(ss);
00617 
00618     // set the block size and allocate data
00619     int bs = params.get("Block Size", 1);
00620     int nb = params.get("Num Blocks", 3*problem_->getNEV());
00621     setSize(bs,nb);
00622 
00623     // get the number of Ritz vectors to compute and allocate data.
00624     // --> if this parameter is not specified in the parameter list, then it's assumed that no Ritz vectors will be computed.
00625     int numRitzVecs = params.get("Number of Ritz Vectors", 0);
00626     setNumRitzVectors( numRitzVecs );
00627 
00628     // get the number of Ritz values to print out when currentStatus is called.
00629     numRitzPrint_ = params.get("Print Number of Ritz Values", bs);
00630   }
00631 
00632 
00634   // Set the block size
00635   // This simply calls setSize(), modifying the block size while retaining the number of blocks.
00636   template <class ScalarType, class MV, class OP>
00637   void BlockKrylovSchur<ScalarType,MV,OP>::setBlockSize (int blockSize) 
00638   {
00639     setSize(blockSize,numBlocks_);
00640   }
00641 
00642 
00644   // Set the step size.
00645   template <class ScalarType, class MV, class OP>
00646   void BlockKrylovSchur<ScalarType,MV,OP>::setStepSize (int stepSize)
00647   {
00648     TEUCHOS_TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
00649     stepSize_ = stepSize;
00650   }
00651 
00653   // Set the number of Ritz vectors to compute.
00654   template <class ScalarType, class MV, class OP>
00655   void BlockKrylovSchur<ScalarType,MV,OP>::setNumRitzVectors (int numRitzVecs)
00656   {
00657     // This routine only allocates space; it doesn't not perform any computation
00658     // any change in size will invalidate the state of the solver.
00659 
00660     TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
00661 
00662     // Check to see if the number of requested Ritz vectors has changed.
00663     if (numRitzVecs != numRitzVecs_) {
00664       if (numRitzVecs) {
00665         ritzVectors_ = Teuchos::null;
00666         ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
00667       } else {
00668         ritzVectors_ = Teuchos::null;
00669       }
00670       numRitzVecs_ = numRitzVecs;
00671       ritzVecsCurrent_ = false;
00672     }      
00673   }
00674 
00676   // Set the block size and make necessary adjustments.
00677   template <class ScalarType, class MV, class OP>
00678   void BlockKrylovSchur<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks) 
00679   {
00680     // This routine only allocates space; it doesn't not perform any computation
00681     // any change in size will invalidate the state of the solver.
00682 
00683     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
00684     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
00685     if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00686       // do nothing
00687       return;
00688     }
00689 
00690     blockSize_ = blockSize;
00691     numBlocks_ = numBlocks;
00692 
00693     Teuchos::RCP<const MV> tmp;
00694     // grab some Multivector to Clone
00695     // in practice, getInitVec() should always provide this, but it is possible to use a 
00696     // Eigenproblem with nothing in getInitVec() by manually initializing with initialize(); 
00697     // in case of that strange scenario, we will try to Clone from V_; first resort to getInitVec(), 
00698     // because we would like to clear the storage associated with V_ so we have room for the new V_
00699     if (problem_->getInitVec() != Teuchos::null) {
00700       tmp = problem_->getInitVec();
00701     }
00702     else {
00703       tmp = V_;
00704       TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00705           "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
00706     }
00707 
00708 
00710     // blockSize*numBlocks dependent
00711     //
00712     int newsd;
00713     if (problem_->isHermitian()) {
00714       newsd = blockSize_*numBlocks_;
00715     } else {
00716       newsd = blockSize_*numBlocks_+1;
00717     }
00718     // check that new size is valid
00719     TEUCHOS_TEST_FOR_EXCEPTION(newsd > MVT::GetVecLength(*tmp),std::invalid_argument,
00720         "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
00721 
00722     ritzValues_.resize(newsd);
00723     ritzResiduals_.resize(newsd,MT_ONE);
00724     ritzOrder_.resize(newsd);
00725     V_ = Teuchos::null;
00726     V_ = MVT::Clone(*tmp,newsd+blockSize_);
00727     H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd+blockSize_,newsd) );
00728     Q_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
00729 
00730     initialized_ = false;
00731     curDim_ = 0;
00732   }
00733 
00734 
00736   // Set the auxiliary vectors
00737   template <class ScalarType, class MV, class OP>
00738   void BlockKrylovSchur<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00739     typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
00740 
00741     // set new auxiliary vectors
00742     auxVecs_ = auxvecs;
00743     
00744     if (om_->isVerbosity( Debug ) ) {
00745       // Check almost everything here
00746       CheckList chk;
00747       chk.checkAux = true;
00748       om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
00749     }
00750 
00751     numAuxVecs_ = 0;
00752     for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
00753       numAuxVecs_ += MVT::GetNumberVecs(**i);
00754     }
00755     
00756     // If the solver has been initialized, X and P are not necessarily orthogonal to new auxiliary vectors
00757     if (numAuxVecs_ > 0 && initialized_) {
00758       initialized_ = false;
00759     }
00760   }
00761 
00763   /* Initialize the state of the solver
00764    * 
00765    * POST-CONDITIONS:
00766    *
00767    * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
00768    *
00769    */
00770 
00771   template <class ScalarType, class MV, class OP>
00772   void BlockKrylovSchur<ScalarType,MV,OP>::initialize(BlockKrylovSchurState<ScalarType,MV> newstate)
00773   {
00774     // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
00775 
00776     std::vector<int> bsind(blockSize_);
00777     for (int i=0; i<blockSize_; i++) bsind[i] = i;
00778 
00779     // in BlockKrylovSchur, V and H are required.  
00780     // if either doesn't exist, then we will start with the initial vector.
00781     //
00782     // inconsistent multivectors widths and lengths will not be tolerated, and
00783     // will be treated with exceptions.
00784     //
00785     std::string errstr("Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
00786 
00787     // set up V,H: if the user doesn't specify both of these these, 
00788     // we will start over with the initial vector.
00789     if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
00790 
00791       // initialize V_,H_, and curDim_
00792 
00793       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_),
00794                           std::invalid_argument, errstr );
00795       if (newstate.V != V_) {
00796         TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
00797             std::invalid_argument, errstr );
00798         TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) > getMaxSubspaceDim(),
00799             std::invalid_argument, errstr );
00800       }
00801       TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > getMaxSubspaceDim(),
00802                           std::invalid_argument, errstr );
00803 
00804       curDim_ = newstate.curDim;
00805       int lclDim = MVT::GetNumberVecs(*newstate.V);
00806 
00807       // check size of H
00808       TEUCHOS_TEST_FOR_EXCEPTION(newstate.H->numRows() < curDim_ || newstate.H->numCols() < curDim_, std::invalid_argument, errstr);
00809       
00810       if (curDim_ == 0 && lclDim > blockSize_) {
00811         om_->stream(Warnings) << "Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
00812                                   << "The block size however is only " << blockSize_ << std::endl
00813                                   << "The last " << lclDim - blockSize_ << " vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
00814       }
00815 
00816 
00817       // copy basis vectors from newstate into V
00818       if (newstate.V != V_) {
00819         std::vector<int> nevind(lclDim);
00820         for (int i=0; i<lclDim; i++) nevind[i] = i;
00821         MVT::SetBlock(*newstate.V,nevind,*V_);
00822       }
00823 
00824       // put data into H_, make sure old information is not still hanging around.
00825       if (newstate.H != H_) {
00826         H_->putScalar( ST_ZERO );
00827         Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H,curDim_+blockSize_,curDim_);
00828         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
00829         lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_,curDim_+blockSize_,curDim_) );
00830         lclH->assign(newH);
00831 
00832         // done with local pointers
00833         lclH = Teuchos::null;
00834       }
00835 
00836     }
00837     else {
00838       // user did not specify a basis V
00839       // get vectors from problem or generate something, projectAndNormalize, call initialize() recursively
00840       Teuchos::RCP<const MV> ivec = problem_->getInitVec();
00841       TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
00842                          "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
00843 
00844       int lclDim = MVT::GetNumberVecs(*ivec);
00845       bool userand = false;
00846       if (lclDim < blockSize_) {
00847         // we need at least blockSize_ vectors
00848         // use a random multivec
00849         userand = true;
00850       }
00851               
00852       if (userand) {
00853         // make an index
00854         std::vector<int> dimind2(lclDim);
00855         for (int i=0; i<lclDim; i++) { dimind2[i] = i; }
00856 
00857         // alloc newV as a view of the first block of V
00858         Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,dimind2);
00859 
00860         // copy the initial vectors into the first lclDim vectors of V
00861         MVT::SetBlock(*ivec,dimind2,*newV1);
00862 
00863         // resize / reinitialize the index vector        
00864         dimind2.resize(blockSize_-lclDim);
00865         for (int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
00866 
00867         // initialize the rest of the vectors with random vectors
00868         Teuchos::RCP<MV> newV2 = MVT::CloneViewNonConst(*V_,dimind2);
00869         MVT::MvRandom(*newV2);
00870       }
00871       else {
00872         // alloc newV as a view of the first block of V
00873         Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,bsind);
00874        
00875         // get a view of the first block of initial vectors
00876         Teuchos::RCP<const MV> ivecV = MVT::CloneView(*ivec,bsind);
00877  
00878         // assign ivec to first part of newV
00879         MVT::SetBlock(*ivecV,bsind,*newV1);
00880       }
00881 
00882       // get pointer into first block of V
00883       Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*V_,bsind);
00884 
00885       // remove auxVecs from newV and normalize newV
00886       if (auxVecs_.size() > 0) {
00887 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00888         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
00889 #endif
00890         
00891         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
00892         int rank = orthman_->projectAndNormalize(*newV,auxVecs_);
00893         TEUCHOS_TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
00894                             "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
00895       }
00896       else {
00897 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00898         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
00899 #endif
00900 
00901         int rank = orthman_->normalize(*newV);
00902         TEUCHOS_TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
00903                             "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
00904       }
00905 
00906       // set curDim
00907       curDim_ = 0;
00908 
00909       // clear pointer
00910       newV = Teuchos::null;
00911     }
00912 
00913     // The Ritz vectors/values and Schur form are no longer current.
00914     ritzVecsCurrent_ = false;
00915     ritzValsCurrent_ = false;
00916     schurCurrent_ = false;
00917 
00918     // the solver is initialized
00919     initialized_ = true;
00920 
00921     if (om_->isVerbosity( Debug ) ) {
00922       // Check almost everything here
00923       CheckList chk;
00924       chk.checkV = true;
00925       chk.checkArn = true;
00926       chk.checkAux = true;
00927       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00928     }
00929 
00930     // Print information on current status
00931     if (om_->isVerbosity(Debug)) {
00932       currentStatus( om_->stream(Debug) );
00933     }
00934     else if (om_->isVerbosity(IterationDetails)) {
00935       currentStatus( om_->stream(IterationDetails) );
00936     }
00937   }
00938 
00939 
00941   // initialize the solver with default state
00942   template <class ScalarType, class MV, class OP>
00943   void BlockKrylovSchur<ScalarType,MV,OP>::initialize()
00944   {
00945     BlockKrylovSchurState<ScalarType,MV> empty;
00946     initialize(empty);
00947   }
00948 
00949 
00951   // Perform BlockKrylovSchur iterations until the StatusTest tells us to stop.
00952   template <class ScalarType, class MV, class OP>
00953   void BlockKrylovSchur<ScalarType,MV,OP>::iterate()
00954   {
00955     //
00956     // Allocate/initialize data structures
00957     //
00958     if (initialized_ == false) {
00959       initialize();
00960     }
00961 
00962     // Compute the current search dimension. 
00963     // If the problem is non-Hermitian and the blocksize is one, let the solver use the extra vector.
00964     int searchDim = blockSize_*numBlocks_;
00965     if (problem_->isHermitian() == false) {
00966       searchDim++;
00967     } 
00968 
00970     // iterate until the status test tells us to stop.
00971     //
00972     // also break if our basis is full
00973     //
00974     while (tester_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
00975 
00976       iter_++;
00977 
00978       // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
00979       int lclDim = curDim_ + blockSize_; 
00980 
00981       // Get the current part of the basis.
00982       std::vector<int> curind(blockSize_);
00983       for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
00984       Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
00985 
00986       // Get a view of the previous vectors
00987       // this is used for orthogonalization and for computing V^H K H
00988       for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
00989       Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
00990       // om_->stream(Debug) << "Vprev: " << std::endl;
00991       // MVT::MvPrint(*Vprev,om_->stream(Debug));
00992 
00993       // Compute the next vector in the Krylov basis:  Vnext = Op*Vprev
00994       {
00995 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00996         Teuchos::TimeMonitor lcltimer( *timerOp_ );
00997 #endif
00998         OPT::Apply(*Op_,*Vprev,*Vnext);
00999         count_ApplyOp_ += blockSize_;
01000       }
01001       // om_->stream(Debug) << "Vnext: " << std::endl;
01002       // MVT::MvPrint(*Vnext,om_->stream(Debug));
01003       Vprev = Teuchos::null;
01004       
01005       // Remove all previous Krylov-Schur basis vectors and auxVecs from Vnext
01006       {
01007 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
01008         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01009 #endif
01010         
01011         // Get a view of all the previous vectors
01012         std::vector<int> prevind(lclDim);
01013         for (int i=0; i<lclDim; i++) { prevind[i] = i; }
01014         Vprev = MVT::CloneView(*V_,prevind);
01015         Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
01016         
01017         // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
01018         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
01019           subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
01020                                ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
01021         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
01022         AsubH.append( subH );
01023         
01024         // Add the auxiliary vectors to the current basis vectors if any exist
01025         if (auxVecs_.size() > 0) {
01026           for (Array_size_type i=0; i<auxVecs_.size(); i++) {
01027             AVprev.append( auxVecs_[i] );
01028             AsubH.append( Teuchos::null );
01029           }
01030         }
01031         
01032         // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
01033         // om_->stream(Debug) << "V before ortho: " << std::endl; 
01034         // MVT::MvPrint(*Vprev,om_->stream(Debug));
01035         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
01036           subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
01037                                ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
01038         int rank = orthman_->projectAndNormalize(*Vnext,AVprev,AsubH,subR);
01039         // om_->stream(Debug) << "Vnext after ortho: " << std::endl;
01040         // MVT::MvPrint(*Vnext,om_->stream(Debug));
01041         // om_->stream(Debug) << "subH: " << std::endl << *subH << std::endl;
01042         // om_->stream(Debug) << "subR: " << std::endl << *subR << std::endl;
01043         // om_->stream(Debug) << "H:    " << std::endl << *H_ << std::endl;
01044         TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockKrylovSchurOrthoFailure,
01045                            "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
01046       }
01047       //
01048       // V has been extended, and H has been extended. 
01049       //
01050       // Update basis dim and release all pointers.
01051       Vnext = Teuchos::null;
01052       curDim_ += blockSize_;
01053       // The Ritz vectors/values and Schur form are no longer current.
01054       ritzVecsCurrent_ = false;
01055       ritzValsCurrent_ = false;
01056       schurCurrent_ = false;
01057       //
01058       // Update Ritz values and residuals if needed
01059       if (!(iter_%stepSize_)) {
01060         computeRitzValues();
01061       }
01062       
01063       // When required, monitor some orthogonalities
01064       if (om_->isVerbosity( Debug ) ) {
01065         // Check almost everything here
01066         CheckList chk;
01067         chk.checkV = true;
01068         chk.checkArn = true;
01069         om_->print( Debug, accuracyCheck(chk, ": after local update") );
01070       }
01071       else if (om_->isVerbosity( OrthoDetails ) ) {
01072         CheckList chk;
01073         chk.checkV = true;
01074         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
01075       }
01076       
01077       // Print information on current iteration
01078       if (om_->isVerbosity(Debug)) {
01079         currentStatus( om_->stream(Debug) );
01080       }
01081       else if (om_->isVerbosity(IterationDetails)) {
01082         currentStatus( om_->stream(IterationDetails) );
01083       }
01084       
01085     } // end while (statusTest == false)
01086     
01087   } // end of iterate()
01088 
01089 
01091   // Check accuracy, orthogonality, and other debugging stuff
01092   // 
01093   // bools specify which tests we want to run (instead of running more than we actually care about)
01094   //
01095   // checkV : V orthonormal
01096   //          orthogonal to auxvecs
01097   // checkAux: check that auxiliary vectors are actually orthonormal
01098   //
01099   // checkArn: check the Arnoldi factorization
01100   //
01101   // NOTE:  This method needs to check the current dimension of the subspace, since it is possible to
01102   //        call this method when curDim_ = 0 (after initialization).
01103   template <class ScalarType, class MV, class OP>
01104   std::string BlockKrylovSchur<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const 
01105   {
01106     std::stringstream os;
01107     os.precision(2);
01108     os.setf(std::ios::scientific, std::ios::floatfield);
01109     MagnitudeType tmp;
01110 
01111     os << " Debugging checks: iteration " << iter_ << where << std::endl;
01112 
01113     // index vectors for V and F
01114     std::vector<int> lclind(curDim_);
01115     for (int i=0; i<curDim_; i++) lclind[i] = i;
01116     std::vector<int> bsind(blockSize_);
01117     for (int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
01118     
01119     Teuchos::RCP<const MV> lclV,lclF;
01120     Teuchos::RCP<MV> lclAV;
01121     if (curDim_)
01122       lclV = MVT::CloneView(*V_,lclind);
01123     lclF = MVT::CloneView(*V_,bsind);
01124 
01125     if (chk.checkV) {
01126       if (curDim_) {
01127         tmp = orthman_->orthonormError(*lclV);
01128         os << " >> Error in V^H M V == I  : " << tmp << std::endl;
01129       }
01130       tmp = orthman_->orthonormError(*lclF);
01131       os << " >> Error in F^H M F == I  : " << tmp << std::endl;
01132       if (curDim_) {
01133         tmp = orthman_->orthogError(*lclV,*lclF);
01134         os << " >> Error in V^H M F == 0  : " << tmp << std::endl;
01135       }
01136       for (Array_size_type i=0; i<auxVecs_.size(); i++) {
01137         if (curDim_) {
01138           tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
01139           os << " >> Error in V^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01140         }
01141         tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
01142         os << " >> Error in F^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01143       }
01144     }
01145     
01146     if (chk.checkArn) {
01147 
01148       if (curDim_) {
01149         // Compute AV      
01150         lclAV = MVT::Clone(*V_,curDim_);
01151         {
01152 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
01153           Teuchos::TimeMonitor lcltimer( *timerOp_ );
01154 #endif
01155           OPT::Apply(*Op_,*lclV,*lclAV);
01156         }
01157         
01158         // Compute AV - VH
01159         Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
01160         MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
01161         
01162         // Compute FB_k^T - (AV-VH)
01163         Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
01164                                                         blockSize_,curDim_, curDim_ );
01165         MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
01166         
01167         // Compute || FE_k^T - (AV-VH) ||
01168         std::vector<MagnitudeType> arnNorms( curDim_ );
01169         orthman_->norm( *lclAV, arnNorms );
01170         
01171         for (int i=0; i<curDim_; i++) {        
01172         os << " >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
01173         }
01174       }
01175     }
01176 
01177     if (chk.checkAux) {
01178       for (Array_size_type i=0; i<auxVecs_.size(); i++) {
01179         tmp = orthman_->orthonormError(*auxVecs_[i]);
01180         os << " >> Error in Aux[" << i << "]^H M Aux[" << i << "] == I : " << tmp << std::endl;
01181         for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
01182           tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01183           os << " >> Error in Aux[" << i << "]^H M Aux[" << j << "] == 0 : " << tmp << std::endl;
01184         }
01185       }
01186     }
01187 
01188     os << std::endl;
01189 
01190     return os.str();
01191   }
01192 
01194   /* Get the current approximate eigenvalues, i.e. Ritz values.
01195    * 
01196    * POST-CONDITIONS:
01197    *
01198    * ritzValues_ contains Ritz w.r.t. V, H
01199    * Q_ contains the Schur vectors w.r.t. H
01200    * schurH_ contains the Schur matrix w.r.t. H
01201    * ritzOrder_ contains the current ordering from sort manager
01202    */
01203 
01204   template <class ScalarType, class MV, class OP>  
01205   void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzValues()
01206   {
01207     // Can only call this if the solver is initialized
01208     if (initialized_) {
01209 
01210       // This just updates the Ritz values and residuals.
01211       // --> ritzValsCurrent_ will be set to 'true' by this method.
01212       if (!ritzValsCurrent_) {
01213         // Compute the current Ritz values, through computing the Schur form
01214         //   without updating the current projection matrix or sorting the Schur form.
01215         computeSchurForm( false );
01216       }
01217     }
01218   }
01219 
01221   /* Get the current approximate eigenvectors, i.e. Ritz vectors.
01222    * 
01223    * POST-CONDITIONS:
01224    *
01225    * ritzValues_ contains Ritz w.r.t. V, H
01226    * ritzVectors_ is first blockSize_ Ritz vectors w.r.t. V, H
01227    * Q_ contains the Schur vectors w.r.t. H
01228    * schurH_ contains the Schur matrix w.r.t. H
01229    * ritzOrder_ contains the current ordering from sort manager
01230    */
01231 
01232   template <class ScalarType, class MV, class OP>  
01233   void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzVectors()
01234   {
01235 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
01236     Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
01237 #endif
01238 
01239     TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
01240                        "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
01241 
01242     TEUCHOS_TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
01243                        "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
01244 
01245 
01246     // Check to see if the current subspace dimension is non-trivial and the solver is initialized
01247     if (curDim_ && initialized_) {
01248 
01249       // Check to see if the Ritz vectors are current.
01250       if (!ritzVecsCurrent_) {
01251         
01252         // Check to see if the Schur factorization of H (schurH_, Q) is current and sorted.
01253         if (!schurCurrent_) {
01254           // Compute the Schur factorization of the current H, which will not directly change H,
01255           // the factorization will be sorted and placed in (schurH_, Q)
01256           computeSchurForm( true );
01257         }
01258         
01259         // After the Schur form is computed, then the Ritz values are current.
01260         // Thus, I can check the Ritz index vector to see if I have enough space for the Ritz vectors requested.
01261         TEUCHOS_TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
01262                            "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
01263 
01264         Teuchos::LAPACK<int,ScalarType> lapack;
01265         Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01266 
01267         // Compute the Ritz vectors.
01268         // --> For a Hermitian problem this is simply the current basis times the first numRitzVecs_ Schur vectors
01269         //     
01270         // --> For a non-Hermitian problem, this involves solving the projected eigenproblem, then
01271         //     placing the product of the current basis times the first numRitzVecs_ Schur vectors times the
01272         //     eigenvectors of interest into the Ritz vectors.
01273 
01274         // Get a view of the current Krylov-Schur basis vectors and Schur vectors
01275         std::vector<int> curind( curDim_ );
01276         for (int i=0; i<curDim_; i++) { curind[i] = i; }
01277         Teuchos::RCP<const MV> Vtemp = MVT::CloneView( *V_, curind );
01278         if (problem_->isHermitian()) {
01279           // Get a view into the current Schur vectors
01280           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
01281 
01282           // Compute the current Ritz vectors      
01283           MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
01284           
01285         } else {
01286 
01287           // Get a view into the current Schur vectors.
01288           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01289           
01290           // Get a set of work vectors to hold the current Ritz vectors.
01291           Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
01292 
01293           // Compute the current Krylov-Schur vectors.
01294           MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );          
01295 
01296           //  Now compute the eigenvectors of the Schur form
01297           //  Reset the dense matrix and compute the eigenvalues of the Schur form.
01298           //
01299           // Allocate the work space. This space will be used below for calls to:
01300           // * TREVC (requires 3*N for real, 2*N for complex) 
01301 
01302           int lwork = 3*curDim_;
01303           std::vector<ScalarType> work( lwork );
01304           std::vector<MagnitudeType> rwork( curDim_ );
01305           char side = 'R';
01306           int mm, info = 0; 
01307           const int ldvl = 1;
01308           ScalarType vl[ ldvl ];
01309           Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
01310           lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01311                         copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01312           TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01313                              "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
01314 
01315           // Get a view into the eigenvectors of the Schur form
01316           Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
01317           
01318           // Convert back to Ritz vectors of the operator.
01319           curind.resize( numRitzVecs_ );  // This is already initialized above
01320           Teuchos::RCP<MV> view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, curind );
01321           MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
01322 
01323           // Compute the norm of the new Ritz vectors
01324           std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
01325           MVT::MvNorm( *view_ritzVectors, ritzNrm );
01326 
01327           // Release memory used to compute Ritz vectors before scaling the current vectors.
01328           tmpritzVectors_ = Teuchos::null;
01329           view_ritzVectors = Teuchos::null;
01330           
01331           // Scale the Ritz vectors to have Euclidean norm.
01332           ScalarType ritzScale = ST_ONE;
01333           for (int i=0; i<numRitzVecs_; i++) {
01334             
01335             // If this is a conjugate pair then normalize by the real and imaginary parts.
01336             if (ritzIndex_[i] == 1 ) {
01337               ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
01338               std::vector<int> newind(2);
01339               newind[0] = i; newind[1] = i+1;
01340               tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01341               view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
01342               MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01343 
01344               // Increment counter for imaginary part
01345               i++;
01346             } else {
01347 
01348               // This is a real Ritz value, normalize the vector
01349               std::vector<int> newind(1);
01350               newind[0] = i;
01351               tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01352               view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
01353               MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01354             }              
01355           }
01356           
01357         } // if (problem_->isHermitian()) 
01358         
01359         // The current Ritz vectors have been computed.
01360         ritzVecsCurrent_ = true;
01361         
01362       } // if (!ritzVecsCurrent_)      
01363     } // if (curDim_)    
01364   } // computeRitzVectors()
01365 
01366   
01368   // Set a new StatusTest for the solver.
01369   template <class ScalarType, class MV, class OP>
01370   void BlockKrylovSchur<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
01371     TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
01372         "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
01373     tester_ = test;
01374   }
01375 
01377   // Get the current StatusTest used by the solver.
01378   template <class ScalarType, class MV, class OP>
01379   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockKrylovSchur<ScalarType,MV,OP>::getStatusTest() const {
01380     return tester_;
01381   }
01382   
01383 
01385   /* Get the current approximate eigenvalues, i.e. Ritz values.
01386    * 
01387    * POST-CONDITIONS:
01388    *
01389    * ritzValues_ contains Ritz w.r.t. V, H
01390    * Q_ contains the Schur vectors w.r.t. H
01391    * schurH_ contains the Schur matrix w.r.t. H
01392    * ritzOrder_ contains the current ordering from sort manager
01393    * schurCurrent_ = true if sort = true; i.e. the Schur form is sorted according to the index
01394    *  vector returned by the sort manager.
01395    */
01396   template <class ScalarType, class MV, class OP>
01397   void BlockKrylovSchur<ScalarType,MV,OP>::computeSchurForm( const bool sort )
01398   {
01399     // local timer
01400 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
01401     Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
01402 #endif
01403 
01404     // Check to see if the dimension of the factorization is greater than zero.
01405     if (curDim_) {
01406 
01407       // Check to see if the Schur factorization is current.
01408       if (!schurCurrent_) {
01409         
01410         // Check to see if the Ritz values are current
01411         // --> If they are then the Schur factorization is current but not sorted.
01412         if (!ritzValsCurrent_) {
01413           Teuchos::LAPACK<int,ScalarType> lapack; 
01414           Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01415           Teuchos::BLAS<int,ScalarType> blas;
01416           Teuchos::BLAS<int,MagnitudeType> blas_mag;
01417           
01418           // Get a view into Q, the storage for H's Schur vectors.
01419           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01420           
01421           // Get a copy of H to compute/sort the Schur form.
01422           schurH_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
01423           //
01424           //---------------------------------------------------
01425           // Compute the Schur factorization of subH
01426           // ---> Use driver GEES to first reduce to upper Hessenberg 
01427           //         form and then compute Schur form, outputting Ritz values
01428           //---------------------------------------------------
01429           //
01430           // Allocate the work space. This space will be used below for calls to:
01431           // * GEES  (requires 3*N for real, 2*N for complex)
01432           // * TREVC (requires 3*N for real, 2*N for complex) 
01433           // * TREXC (requires N for real, none for complex)
01434           // Furthermore, GEES requires a real array of length curDim_ (for complex datatypes)
01435           //
01436           int lwork = 3*curDim_;
01437           std::vector<ScalarType> work( lwork );
01438           std::vector<MagnitudeType> rwork( curDim_ );
01439           std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
01440           std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
01441           std::vector<int> bwork( curDim_ );
01442           int info = 0, sdim = 0; 
01443           char jobvs = 'V';
01444           lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
01445                        &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork, 
01446                        &rwork[0], &bwork[0], &info );
01447           
01448           TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01449                              "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ << ") returned info " << info << " != 0.");
01450           //
01451           //---------------------------------------------------
01452           // Use the Krylov-Schur factorization to compute the current Ritz residuals 
01453           // for ALL the eigenvalues estimates (Ritz values)
01454           //           || Ax - x\theta || = || U_m+1*B_m+1^H*Q*s || 
01455           //                              = || B_m+1^H*Q*s ||
01456           //
01457           // where U_m+1 is the current Krylov-Schur basis, Q are the Schur vectors, and x = U_m+1*Q*s
01458           // NOTE: This means that s = e_i if the problem is hermitian, else the eigenvectors
01459           //       of the Schur form need to be computed.
01460           //
01461           // First compute H_{m+1,m}*B_m^T, then determine what 's' is.
01462           //---------------------------------------------------
01463           //
01464           // Get current B_m+1
01465           Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View, *H_,
01466                                                           blockSize_, curDim_, curDim_ );
01467           //
01468           // Compute B_m+1^H*Q
01469           Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
01470           blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE, 
01471                      curB.values(), curB.stride(), subQ.values(), subQ.stride(), 
01472                      ST_ZERO, subB.values(), subB.stride() );
01473           //
01474           // Determine what 's' is and compute Ritz residuals.
01475           //
01476           ScalarType* b_ptr = subB.values();
01477           if (problem_->isHermitian()) {
01478             //
01479             // 's' is the i-th canonical basis vector.
01480             //
01481             for (int i=0; i<curDim_ ; i++) {
01482               ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
01483             }   
01484           } else {
01485             //
01486             //  Compute S: the eigenvectors of the block upper triangular, Schur matrix.
01487             //
01488             char side = 'R';
01489             int mm;
01490             const int ldvl = 1;
01491             ScalarType vl[ ldvl ];
01492             Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
01493             lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01494                           S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01495             
01496             TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01497                                "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
01498             //
01499             // Scale the eigenvectors so that their Euclidean norms are all one.
01500             //
01501             HelperTraits<ScalarType>::scaleRitzVectors( tmp_iRitzValues, &S );
01502             //
01503             // Compute ritzRes = *B_m+1^H*Q*S where the i-th column of S is 's' for the i-th Ritz-value
01504             //
01505             Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
01506             blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE, 
01507                        subB.values(), subB.stride(), S.values(), S.stride(), 
01508                        ST_ZERO, ritzRes.values(), ritzRes.stride() );
01509 
01510             /* TO DO:  There's be an incorrect assumption made in the computation of the Ritz residuals.
01511                        This assumption is that the next vector in the Krylov subspace is Euclidean orthonormal.
01512                        It may not be normalized using Euclidean norm.
01513             Teuchos::RCP<MV> ritzResVecs = MVT::Clone( *V_, curDim_ );
01514             std::vector<int> curind(blockSize_);
01515             for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
01516             Teuchos::RCP<MV> Vtemp = MVT::CloneView(*V_,curind);     
01517             
01518             MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, ritzRes, ST_ZERO, *ritzResVecs );
01519             std::vector<MagnitudeType> ritzResNrms(curDim_);
01520             MVT::MvNorm( *ritzResVecs, ritzResNrms );
01521             i = 0;
01522             while( i < curDim_ ) {
01523               if ( tmp_ritzValues[curDim_+i] != MT_ZERO ) {
01524                 ritzResiduals_[i] = lapack_mag.LAPY2( ritzResNrms[i], ritzResNrms[i+1] );
01525                 ritzResiduals_[i+1] = ritzResiduals_[i];
01526                 i = i+2;
01527               } else {
01528                 ritzResiduals_[i] = ritzResNrms[i];
01529                 i++;
01530               }
01531             }
01532             */
01533             //
01534             // Compute the Ritz residuals for each Ritz value.
01535             // 
01536             HelperTraits<ScalarType>::computeRitzResiduals( tmp_iRitzValues, ritzRes, &ritzResiduals_ );
01537           }
01538           //
01539           // Sort the Ritz values.
01540           //
01541           {
01542 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
01543             Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
01544 #endif
01545             int i=0;
01546             if (problem_->isHermitian()) {
01547               //
01548               // Sort using just the real part of the Ritz values.
01549               sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_); // don't catch exception
01550               ritzIndex_.clear();
01551               while ( i < curDim_ ) {
01552                 // The Ritz value is not complex.
01553                 ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
01554                 ritzIndex_.push_back(0);
01555                 i++;
01556               }
01557             }
01558             else {
01559               //
01560               // Sort using both the real and imaginary parts of the Ritz values.
01561               sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
01562               HelperTraits<ScalarType>::sortRitzValues( tmp_rRitzValues, tmp_iRitzValues, &ritzValues_, &ritzOrder_, &ritzIndex_ );
01563             }
01564             //
01565             // Sort the ritzResiduals_ based on the ordering from the Sort Manager.
01566             std::vector<MagnitudeType> ritz2( curDim_ );
01567             for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
01568             blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
01569             
01570             // The Ritz values have now been updated.
01571             ritzValsCurrent_ = true;
01572           }
01573 
01574         } // if (!ritzValsCurrent_) ...
01575         // 
01576         //---------------------------------------------------
01577         // * The Ritz values and residuals have been updated at this point.
01578         // 
01579         // * The Schur factorization of the projected matrix has been computed,
01580         //   and is stored in (schurH_, Q_).
01581         //
01582         // Now the Schur factorization needs to be sorted.
01583         //---------------------------------------------------
01584         //
01585         // Sort the Schur form using the ordering from the Sort Manager.
01586         if (sort) {
01587           sortSchurForm( *schurH_, *Q_, ritzOrder_ );    
01588           //
01589           // Indicate the Schur form in (schurH_, Q_) is current and sorted
01590           schurCurrent_ = true;
01591         }
01592       } // if (!schurCurrent_) ...
01593   
01594     } // if (curDim_) ...
01595   
01596   } // computeSchurForm( ... )
01597   
01598 
01600   // Sort the Schur form of H stored in (H,Q) using the ordering vector.
01601   template <class ScalarType, class MV, class OP>
01602   void BlockKrylovSchur<ScalarType,MV,OP>::sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
01603                                                           Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
01604                                                           std::vector<int>& order ) 
01605   {
01606     // local timer
01607 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
01608     Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
01609 #endif
01610     //
01611     //---------------------------------------------------
01612     // Reorder real Schur factorization, remember to add one to the indices for the
01613     // fortran call and determine offset.  The offset is necessary since the TREXC
01614     // method reorders in a nonsymmetric fashion, thus we use the reordering in
01615     // a stack-like fashion.  Also take into account conjugate pairs, which may mess
01616     // up the reordering, since the pair is moved if one of the pair is moved.
01617     //---------------------------------------------------
01618     //
01619     int i = 0, nevtemp = 0;
01620     char compq = 'V';
01621     std::vector<int> offset2( curDim_ );
01622     std::vector<int> order2( curDim_ );
01623 
01624     // LAPACK objects.
01625     Teuchos::LAPACK<int,ScalarType> lapack; 
01626     int lwork = 3*curDim_;
01627     std::vector<ScalarType> work( lwork );
01628 
01629     while (i < curDim_) {
01630       if ( ritzIndex_[i] != 0 ) { // This is the first value of a complex conjugate pair
01631         offset2[nevtemp] = 0;
01632         for (int j=i; j<curDim_; j++) {
01633           if (order[j] > order[i]) { offset2[nevtemp]++; }
01634         }
01635         order2[nevtemp] = order[i];
01636         i = i+2;
01637       } else {
01638         offset2[nevtemp] = 0;
01639         for (int j=i; j<curDim_; j++) {
01640           if (order[j] > order[i]) { offset2[nevtemp]++; }
01641         }
01642         order2[nevtemp] = order[i];
01643         i++;
01644       }
01645       nevtemp++;
01646     }
01647     ScalarType *ptr_h = H.values();
01648     ScalarType *ptr_q = Q.values();
01649     int ldh = H.stride(), ldq = Q.stride();
01650     int info = 0;
01651     for (i=nevtemp-1; i>=0; i--) {
01652       lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, order2[i]+1+offset2[i], 
01653                     1, &work[0], &info );
01654       TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01655                          "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ << ") returned info " << info << " != 0.");
01656     }
01657   }
01658 
01660   // Print the current status of the solver
01661   template <class ScalarType, class MV, class OP>
01662   void BlockKrylovSchur<ScalarType,MV,OP>::currentStatus(std::ostream &os) 
01663   {
01664     using std::endl;
01665 
01666     os.setf(std::ios::scientific, std::ios::floatfield);
01667     os.precision(6);
01668     os <<"================================================================================" << endl;
01669     os << endl;
01670     os <<"                         BlockKrylovSchur Solver Status" << endl;
01671     os << endl;
01672     os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01673     os <<"The number of iterations performed is " <<iter_<<endl;
01674     os <<"The block size is         " << blockSize_<<endl;
01675     os <<"The number of blocks is   " << numBlocks_<<endl;
01676     os <<"The current basis size is " << curDim_<<endl;
01677     os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
01678     os <<"The number of operations Op*x   is "<<count_ApplyOp_<<endl;
01679 
01680     os.setf(std::ios_base::right, std::ios_base::adjustfield);
01681 
01682     os << endl;
01683     if (initialized_) {
01684       os <<"CURRENT RITZ VALUES             "<<endl;
01685       if (ritzIndex_.size() != 0) {
01686         int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
01687         if (problem_->isHermitian()) {
01688           os << std::setw(20) << "Ritz Value" 
01689              << std::setw(20) << "Ritz Residual"
01690              << endl;
01691           os <<"--------------------------------------------------------------------------------"<<endl;
01692           for (int i=0; i<numPrint; i++) {
01693             os << std::setw(20) << ritzValues_[i].realpart 
01694                << std::setw(20) << ritzResiduals_[i] 
01695                << endl;
01696           }
01697         } else {
01698           os << std::setw(24) << "Ritz Value" 
01699              << std::setw(30) << "Ritz Residual"
01700              << endl;
01701           os <<"--------------------------------------------------------------------------------"<<endl;
01702           for (int i=0; i<numPrint; i++) {
01703             // Print out the real eigenvalue.
01704             os << std::setw(15) << ritzValues_[i].realpart;
01705             if (ritzValues_[i].imagpart < MT_ZERO) {
01706               os << " - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
01707             } else {
01708               os << " + i" << std::setw(15) << ritzValues_[i].imagpart;
01709             }              
01710             os << std::setw(20) << ritzResiduals_[i] << endl;
01711           }
01712         }
01713       } else {
01714         os << std::setw(20) << "[ NONE COMPUTED ]" << endl;
01715       }
01716     }
01717     os << endl;
01718     os <<"================================================================================" << endl;
01719     os << endl;
01720   }
01721   
01722 } // End of namespace Anasazi
01723 
01724 #endif
01725 
01726 // End of file AnasaziBlockKrylovSchur.hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends