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::CloneView(*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::CloneView(*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::CloneView(*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::CloneView(*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::CloneView(*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<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 (unsigned int 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<MV> lclV,lclF,lclAV;
01110     if (curDim_)
01111       lclV = MVT::CloneView(*V_,lclind);
01112     lclF = MVT::CloneView(*V_,bsind);
01113 
01114     if (chk.checkV) {
01115       if (curDim_) {
01116         tmp = orthman_->orthonormError(*lclV);
01117         os << " >> Error in V^H M V == I  : " << tmp << std::endl;
01118       }
01119       tmp = orthman_->orthonormError(*lclF);
01120       os << " >> Error in F^H M F == I  : " << tmp << std::endl;
01121       if (curDim_) {
01122         tmp = orthman_->orthogError(*lclV,*lclF);
01123         os << " >> Error in V^H M F == 0  : " << tmp << std::endl;
01124       }
01125       for (unsigned int i=0; i<auxVecs_.size(); i++) {
01126         if (curDim_) {
01127           tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
01128           os << " >> Error in V^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01129         }
01130         tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
01131         os << " >> Error in F^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01132       }
01133     }
01134     
01135     if (chk.checkArn) {
01136 
01137       if (curDim_) {
01138         // Compute AV      
01139         lclAV = MVT::Clone(*V_,curDim_);
01140         {
01141           Teuchos::TimeMonitor lcltimer( *timerOp_ );
01142           OPT::Apply(*Op_,*lclV,*lclAV);
01143         }
01144         
01145         // Compute AV - VH
01146         Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
01147         MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
01148         
01149         // Compute FB_k^T - (AV-VH)
01150         Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
01151                                                         blockSize_,curDim_, curDim_ );
01152         MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
01153         
01154         // Compute || FE_k^T - (AV-VH) ||
01155         std::vector<MagnitudeType> arnNorms( curDim_ );
01156         orthman_->norm( *lclAV, arnNorms );
01157         
01158         for (int i=0; i<curDim_; i++) {        
01159         os << " >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
01160         }
01161       }
01162     }
01163 
01164     if (chk.checkAux) {
01165       for (unsigned int i=0; i<auxVecs_.size(); i++) {
01166         tmp = orthman_->orthonormError(*auxVecs_[i]);
01167         os << " >> Error in Aux[" << i << "]^H M Aux[" << i << "] == I : " << tmp << std::endl;
01168         for (unsigned int j=i+1; j<auxVecs_.size(); j++) {
01169           tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01170           os << " >> Error in Aux[" << i << "]^H M Aux[" << j << "] == 0 : " << tmp << std::endl;
01171         }
01172       }
01173     }
01174 
01175     os << std::endl;
01176 
01177     return os.str();
01178   }
01179 
01181   /* Get the current approximate eigenvalues, i.e. Ritz values.
01182    * 
01183    * POST-CONDITIONS:
01184    *
01185    * ritzValues_ contains Ritz w.r.t. V, H
01186    * Q_ contains the Schur vectors w.r.t. H
01187    * schurH_ contains the Schur matrix w.r.t. H
01188    * ritzOrder_ contains the current ordering from sort manager
01189    */
01190 
01191   template <class ScalarType, class MV, class OP>  
01192   void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzValues()
01193   {
01194     // Can only call this if the solver is initialized
01195     if (initialized_) {
01196 
01197       // This just updates the Ritz values and residuals.
01198       // --> ritzValsCurrent_ will be set to 'true' by this method.
01199       if (!ritzValsCurrent_) {
01200         // Compute the current Ritz values, through computing the Schur form
01201         //   without updating the current projection matrix or sorting the Schur form.
01202         computeSchurForm( false );
01203       }
01204     }
01205   }
01206 
01208   /* Get the current approximate eigenvectors, i.e. Ritz vectors.
01209    * 
01210    * POST-CONDITIONS:
01211    *
01212    * ritzValues_ contains Ritz w.r.t. V, H
01213    * ritzVectors_ is first blockSize_ Ritz vectors w.r.t. V, H
01214    * Q_ contains the Schur vectors w.r.t. H
01215    * schurH_ contains the Schur matrix w.r.t. H
01216    * ritzOrder_ contains the current ordering from sort manager
01217    */
01218 
01219   template <class ScalarType, class MV, class OP>  
01220   void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzVectors()
01221   {
01222     Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
01223 
01224     TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
01225                        "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
01226 
01227     TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
01228                        "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
01229 
01230 
01231     // Check to see if the current subspace dimension is non-trivial and the solver is initialized
01232     if (curDim_ && initialized_) {
01233 
01234       // Check to see if the Ritz vectors are current.
01235       if (!ritzVecsCurrent_) {
01236         
01237         // Check to see if the Schur factorization of H (schurH_, Q) is current and sorted.
01238         if (!schurCurrent_) {
01239           // Compute the Schur factorization of the current H, which will not directly change H,
01240           // the factorization will be sorted and placed in (schurH_, Q)
01241           computeSchurForm( true );
01242         }
01243         
01244         // After the Schur form is computed, then the Ritz values are current.
01245         // Thus, I can check the Ritz index vector to see if I have enough space for the Ritz vectors requested.
01246         TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
01247                            "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
01248 
01249         Teuchos::LAPACK<int,ScalarType> lapack;
01250         Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01251 
01252         // Compute the Ritz vectors.
01253         // --> For a Hermitian problem this is simply the current basis times the first numRitzVecs_ Schur vectors
01254         //     
01255         // --> For a non-Hermitian problem, this involves solving the projected eigenproblem, then
01256         //     placing the product of the current basis times the first numRitzVecs_ Schur vectors times the
01257         //     eigenvectors of interest into the Ritz vectors.
01258 
01259         // Get a view of the current Krylov-Schur basis vectors and Schur vectors
01260         std::vector<int> curind( curDim_ );
01261         for (int i=0; i<curDim_; i++) { curind[i] = i; }
01262         Teuchos::RCP<MV> Vtemp = MVT::CloneView( *V_, curind );
01263         if (problem_->isHermitian()) {
01264           // Get a view into the current Schur vectors
01265           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
01266 
01267           // Compute the current Ritz vectors      
01268           MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
01269           
01270         } else {
01271 
01272           // Get a view into the current Schur vectors.
01273           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01274           
01275           // Get a set of work vectors to hold the current Ritz vectors.
01276           Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
01277 
01278           // Compute the current Krylov-Schur vectors.
01279           MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );          
01280 
01281           //  Now compute the eigenvectors of the Schur form
01282           //  Reset the dense matrix and compute the eigenvalues of the Schur form.
01283           //
01284           // Allocate the work space. This space will be used below for calls to:
01285           // * TREVC (requires 3*N for real, 2*N for complex) 
01286 
01287           int lwork = 3*curDim_;
01288           std::vector<ScalarType> work( lwork );
01289           std::vector<MagnitudeType> rwork( curDim_ );
01290           char side = 'R';
01291           int mm, info = 0; 
01292           const int ldvl = 1;
01293           ScalarType vl[ ldvl ];
01294           Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
01295           lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01296                         copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01297           TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01298                              "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
01299 
01300           // Get a view into the eigenvectors of the Schur form
01301           Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
01302           
01303           // Convert back to Ritz vectors of the operator.
01304           curind.resize( numRitzVecs_ );  // This is already initialized above
01305           Teuchos::RCP<MV> view_ritzVectors = MVT::CloneView( *ritzVectors_, curind );
01306           MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
01307 
01308           // Compute the norm of the new Ritz vectors
01309           std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
01310           MVT::MvNorm( *view_ritzVectors, ritzNrm );
01311 
01312           // Release memory used to compute Ritz vectors before scaling the current vectors.
01313           tmpritzVectors_ = Teuchos::null;
01314           view_ritzVectors = Teuchos::null;
01315           
01316           // Scale the Ritz vectors to have Euclidean norm.
01317           ScalarType ritzScale = ST_ONE;
01318           for (int i=0; i<numRitzVecs_; i++) {
01319             
01320             // If this is a conjugate pair then normalize by the real and imaginary parts.
01321             if (ritzIndex_[i] == 1 ) {
01322               ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
01323               std::vector<int> newind(2);
01324               newind[0] = i; newind[1] = i+1;
01325               tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01326               view_ritzVectors = MVT::CloneView( *ritzVectors_, newind );
01327               MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01328 
01329               // Increment counter for imaginary part
01330               i++;
01331             } else {
01332 
01333               // This is a real Ritz value, normalize the vector
01334               std::vector<int> newind(1);
01335               newind[0] = i;
01336               tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01337               view_ritzVectors = MVT::CloneView( *ritzVectors_, newind );
01338               MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01339             }              
01340           }
01341           
01342         } // if (problem_->isHermitian()) 
01343         
01344         // The current Ritz vectors have been computed.
01345         ritzVecsCurrent_ = true;
01346         
01347       } // if (!ritzVecsCurrent_)      
01348     } // if (curDim_)    
01349   } // computeRitzVectors()
01350 
01351   
01353   // Set a new StatusTest for the solver.
01354   template <class ScalarType, class MV, class OP>
01355   void BlockKrylovSchur<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
01356     TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
01357         "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
01358     tester_ = test;
01359   }
01360 
01362   // Get the current StatusTest used by the solver.
01363   template <class ScalarType, class MV, class OP>
01364   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockKrylovSchur<ScalarType,MV,OP>::getStatusTest() const {
01365     return tester_;
01366   }
01367   
01368 
01370   /* Get the current approximate eigenvalues, i.e. Ritz values.
01371    * 
01372    * POST-CONDITIONS:
01373    *
01374    * ritzValues_ contains Ritz w.r.t. V, H
01375    * Q_ contains the Schur vectors w.r.t. H
01376    * schurH_ contains the Schur matrix w.r.t. H
01377    * ritzOrder_ contains the current ordering from sort manager
01378    * schurCurrent_ = true if sort = true; i.e. the Schur form is sorted according to the index
01379    *  vector returned by the sort manager.
01380    */
01381   template <class ScalarType, class MV, class OP>
01382   void BlockKrylovSchur<ScalarType,MV,OP>::computeSchurForm( const bool sort )
01383   {
01384     // local timer
01385     Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
01386 
01387     // Check to see if the dimension of the factorization is greater than zero.
01388     if (curDim_) {
01389 
01390       // Check to see if the Schur factorization is current.
01391       if (!schurCurrent_) {
01392         
01393         // Check to see if the Ritz values are current
01394         // --> If they are then the Schur factorization is current but not sorted.
01395         if (!ritzValsCurrent_) {
01396           Teuchos::LAPACK<int,ScalarType> lapack; 
01397           Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01398           Teuchos::BLAS<int,ScalarType> blas;
01399           Teuchos::BLAS<int,MagnitudeType> blas_mag;
01400           
01401           // Get a view into Q, the storage for H's Schur vectors.
01402           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01403           
01404           // Get a copy of H to compute/sort the Schur form.
01405           schurH_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
01406           //
01407           //---------------------------------------------------
01408           // Compute the Schur factorization of subH
01409           // ---> Use driver GEES to first reduce to upper Hessenberg 
01410           //         form and then compute Schur form, outputting Ritz values
01411           //---------------------------------------------------
01412           //
01413           // Allocate the work space. This space will be used below for calls to:
01414           // * GEES  (requires 3*N for real, 2*N for complex)
01415           // * TREVC (requires 3*N for real, 2*N for complex) 
01416           // * TREXC (requires N for real, none for complex)
01417           // Furthermore, GEES requires a real array of length curDim_ (for complex datatypes)
01418           //
01419           int lwork = 3*curDim_;
01420           std::vector<ScalarType> work( lwork );
01421           std::vector<MagnitudeType> rwork( curDim_ );
01422           std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
01423           std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
01424           std::vector<int> bwork( curDim_ );
01425           int info = 0, sdim = 0; 
01426           char jobvs = 'V';
01427           lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
01428                        &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork, 
01429                        &rwork[0], &bwork[0], &info );
01430           
01431           TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01432                              "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ << ") returned info " << info << " != 0.");
01433           //
01434           //---------------------------------------------------
01435           // Use the Krylov-Schur factorization to compute the current Ritz residuals 
01436           // for ALL the eigenvalues estimates (Ritz values)
01437           //           || Ax - x\theta || = || U_m+1*B_m+1^H*Q*s || 
01438           //                              = || B_m+1^H*Q*s ||
01439           //
01440           // where U_m+1 is the current Krylov-Schur basis, Q are the Schur vectors, and x = U_m+1*Q*s
01441           // NOTE: This means that s = e_i if the problem is hermitian, else the eigenvectors
01442           //       of the Schur form need to be computed.
01443           //
01444           // First compute H_{m+1,m}*B_m^T, then determine what 's' is.
01445           //---------------------------------------------------
01446           //
01447           // Get current B_m+1
01448           Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View, *H_,
01449                                                           blockSize_, curDim_, curDim_ );
01450           //
01451           // Compute B_m+1^H*Q
01452           Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
01453           blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE, 
01454                      curB.values(), curB.stride(), subQ.values(), subQ.stride(), 
01455                      ST_ZERO, subB.values(), subB.stride() );
01456           //
01457           // Determine what 's' is and compute Ritz residuals.
01458           //
01459           ScalarType* b_ptr = subB.values();
01460           if (problem_->isHermitian()) {
01461             //
01462             // 's' is the i-th canonical basis vector.
01463             //
01464             for (int i=0; i<curDim_ ; i++) {
01465               ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
01466             }   
01467           } else {
01468             //
01469             //  Compute S: the eigenvectors of the block upper triangular, Schur matrix.
01470             //
01471             char side = 'R';
01472             int mm;
01473             const int ldvl = 1;
01474             ScalarType vl[ ldvl ];
01475             Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
01476             lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01477                           S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01478             
01479             TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01480                                "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
01481             //
01482             // Scale the eigenvectors so that their Euclidean norms are all one.
01483             //
01484             HelperTraits<ScalarType>::scaleRitzVectors( tmp_iRitzValues, &S );
01485             //
01486             // Compute ritzRes = *B_m+1^H*Q*S where the i-th column of S is 's' for the i-th Ritz-value
01487             //
01488             Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
01489             blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE, 
01490                        subB.values(), subB.stride(), S.values(), S.stride(), 
01491                        ST_ZERO, ritzRes.values(), ritzRes.stride() );
01492 
01493             /* TO DO:  There's be an incorrect assumption made in the computation of the Ritz residuals.
01494                        This assumption is that the next vector in the Krylov subspace is Euclidean orthonormal.
01495                        It may not be normalized using Euclidean norm.
01496             Teuchos::RCP<MV> ritzResVecs = MVT::Clone( *V_, curDim_ );
01497             std::vector<int> curind(blockSize_);
01498             for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
01499             Teuchos::RCP<MV> Vtemp = MVT::CloneView(*V_,curind);     
01500             
01501             MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, ritzRes, ST_ZERO, *ritzResVecs );
01502             std::vector<MagnitudeType> ritzResNrms(curDim_);
01503             MVT::MvNorm( *ritzResVecs, ritzResNrms );
01504             i = 0;
01505             while( i < curDim_ ) {
01506               if ( tmp_ritzValues[curDim_+i] != MT_ZERO ) {
01507                 ritzResiduals_[i] = lapack_mag.LAPY2( ritzResNrms[i], ritzResNrms[i+1] );
01508                 ritzResiduals_[i+1] = ritzResiduals_[i];
01509                 i = i+2;
01510               } else {
01511                 ritzResiduals_[i] = ritzResNrms[i];
01512                 i++;
01513               }
01514             }
01515             */
01516             //
01517             // Compute the Ritz residuals for each Ritz value.
01518             // 
01519             HelperTraits<ScalarType>::computeRitzResiduals( tmp_iRitzValues, ritzRes, &ritzResiduals_ );
01520           }
01521           //
01522           // Sort the Ritz values.
01523           //
01524           {
01525             Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
01526             int i=0;
01527             if (problem_->isHermitian()) {
01528               //
01529               // Sort using just the real part of the Ritz values.
01530               sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_); // don't catch exception
01531               ritzIndex_.clear();
01532               while ( i < curDim_ ) {
01533                 // The Ritz value is not complex.
01534                 ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
01535                 ritzIndex_.push_back(0);
01536                 i++;
01537               }
01538             }
01539             else {
01540               //
01541               // Sort using both the real and imaginary parts of the Ritz values.
01542               sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
01543               HelperTraits<ScalarType>::sortRitzValues( tmp_rRitzValues, tmp_iRitzValues, &ritzValues_, &ritzOrder_, &ritzIndex_ );
01544             }
01545             //
01546             // Sort the ritzResiduals_ based on the ordering from the Sort Manager.
01547             std::vector<MagnitudeType> ritz2( curDim_ );
01548             for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
01549             blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
01550             
01551             // The Ritz values have now been updated.
01552             ritzValsCurrent_ = true;
01553           }
01554 
01555         } // if (!ritzValsCurrent_) ...
01556         // 
01557         //---------------------------------------------------
01558         // * The Ritz values and residuals have been updated at this point.
01559         // 
01560         // * The Schur factorization of the projected matrix has been computed,
01561         //   and is stored in (schurH_, Q_).
01562         //
01563         // Now the Schur factorization needs to be sorted.
01564         //---------------------------------------------------
01565         //
01566         // Sort the Schur form using the ordering from the Sort Manager.
01567         if (sort) {
01568           sortSchurForm( *schurH_, *Q_, ritzOrder_ );    
01569           //
01570           // Indicate the Schur form in (schurH_, Q_) is current and sorted
01571           schurCurrent_ = true;
01572         }
01573       } // if (!schurCurrent_) ...
01574   
01575     } // if (curDim_) ...
01576   
01577   } // computeSchurForm( ... )
01578   
01579 
01581   // Sort the Schur form of H stored in (H,Q) using the ordering vector.
01582   template <class ScalarType, class MV, class OP>
01583   void BlockKrylovSchur<ScalarType,MV,OP>::sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
01584                                                           Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
01585                                                           std::vector<int>& order ) 
01586   {
01587     // local timer
01588     Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
01589     //
01590     //---------------------------------------------------
01591     // Reorder real Schur factorization, remember to add one to the indices for the
01592     // fortran call and determine offset.  The offset is necessary since the TREXC
01593     // method reorders in a nonsymmetric fashion, thus we use the reordering in
01594     // a stack-like fashion.  Also take into account conjugate pairs, which may mess
01595     // up the reordering, since the pair is moved if one of the pair is moved.
01596     //---------------------------------------------------
01597     //
01598     int i = 0, nevtemp = 0;
01599     char compq = 'V';
01600     std::vector<int> offset2( curDim_ );
01601     std::vector<int> order2( curDim_ );
01602 
01603     // LAPACK objects.
01604     Teuchos::LAPACK<int,ScalarType> lapack; 
01605     int lwork = 3*curDim_;
01606     std::vector<ScalarType> work( lwork );
01607 
01608     while (i < curDim_) {
01609       if ( ritzIndex_[i] != 0 ) { // This is the first value of a complex conjugate pair
01610         offset2[nevtemp] = 0;
01611         for (int j=i; j<curDim_; j++) {
01612           if (order[j] > order[i]) { offset2[nevtemp]++; }
01613         }
01614         order2[nevtemp] = order[i];
01615         i = i+2;
01616       } else {
01617         offset2[nevtemp] = 0;
01618         for (int j=i; j<curDim_; j++) {
01619           if (order[j] > order[i]) { offset2[nevtemp]++; }
01620         }
01621         order2[nevtemp] = order[i];
01622         i++;
01623       }
01624       nevtemp++;
01625     }
01626     ScalarType *ptr_h = H.values();
01627     ScalarType *ptr_q = Q.values();
01628     int ldh = H.stride(), ldq = Q.stride();
01629     int info = 0;
01630     for (i=nevtemp-1; i>=0; i--) {
01631       lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, order2[i]+1+offset2[i], 
01632                     1, &work[0], &info );
01633       TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01634                          "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ << ") returned info " << info << " != 0.");
01635     }
01636   }
01637 
01639   // Print the current status of the solver
01640   template <class ScalarType, class MV, class OP>
01641   void BlockKrylovSchur<ScalarType,MV,OP>::currentStatus(std::ostream &os) 
01642   {
01643     using std::endl;
01644 
01645     os.setf(std::ios::scientific, std::ios::floatfield);
01646     os.precision(6);
01647     os <<"================================================================================" << endl;
01648     os << endl;
01649     os <<"                         BlockKrylovSchur Solver Status" << endl;
01650     os << endl;
01651     os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01652     os <<"The number of iterations performed is " <<iter_<<endl;
01653     os <<"The block size is         " << blockSize_<<endl;
01654     os <<"The number of blocks is   " << numBlocks_<<endl;
01655     os <<"The current basis size is " << curDim_<<endl;
01656     os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
01657     os <<"The number of operations Op*x   is "<<count_ApplyOp_<<endl;
01658 
01659     os.setf(std::ios_base::right, std::ios_base::adjustfield);
01660 
01661     os << endl;
01662     if (initialized_) {
01663       os <<"CURRENT RITZ VALUES             "<<endl;
01664       if (ritzIndex_.size() != 0) {
01665         int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
01666         if (problem_->isHermitian()) {
01667           os << std::setw(20) << "Ritz Value" 
01668              << std::setw(20) << "Ritz Residual"
01669              << endl;
01670           os <<"--------------------------------------------------------------------------------"<<endl;
01671           for (int i=0; i<numPrint; i++) {
01672             os << std::setw(20) << ritzValues_[i].realpart 
01673                << std::setw(20) << ritzResiduals_[i] 
01674                << endl;
01675           }
01676         } else {
01677           os << std::setw(24) << "Ritz Value" 
01678              << std::setw(30) << "Ritz Residual"
01679              << endl;
01680           os <<"--------------------------------------------------------------------------------"<<endl;
01681           for (int i=0; i<numPrint; i++) {
01682             // Print out the real eigenvalue.
01683             os << std::setw(15) << ritzValues_[i].realpart;
01684             if (ritzValues_[i].imagpart < MT_ZERO) {
01685               os << " - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
01686             } else {
01687               os << " + i" << std::setw(15) << ritzValues_[i].imagpart;
01688             }              
01689             os << std::setw(20) << ritzResiduals_[i] << endl;
01690           }
01691         }
01692       } else {
01693         os << std::setw(20) << "[ NONE COMPUTED ]" << endl;
01694       }
01695     }
01696     os << endl;
01697     os <<"================================================================================" << endl;
01698     os << endl;
01699   }
01700   
01701 } // End of namespace Anasazi
01702 
01703 #endif
01704 
01705 // End of file AnasaziBlockKrylovSchur.hpp

Generated on Wed May 12 21:24:34 2010 for Anasazi by  doxygen 1.4.7