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

Generated on Thu Dec 17 11:02:20 2009 for Anasazi by  doxygen 1.5.9