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