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