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 
00043 #include "AnasaziOrthoManager.hpp"
00044 
00045 #include "Teuchos_LAPACK.hpp"
00046 #include "Teuchos_BLAS.hpp"
00047 #include "Teuchos_SerialDenseMatrix.hpp"
00048 #include "Teuchos_ParameterList.hpp"
00049 #include "Teuchos_TimeMonitor.hpp"
00050 
00051 #ifdef HAVE_TEUCHOS_COMPLEX
00052 #if defined(HAVE_COMPLEX)
00053 #define ANSZI_CPLX_CLASS std::complex
00054 #elif  defined(HAVE_COMPLEX_H)
00055 #define ANSZI_CPLX_CLASS ::complex
00056 #endif
00057 #endif
00058 
00073 namespace Anasazi {
00074 
00076 
00077 
00082   template <class ScalarType, class MulVec>
00083   struct BlockKrylovSchurState {
00088     int curDim;
00090     Teuchos::RCP<const MulVec> V;
00096     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > H;
00098     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > S;
00100     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > Q;
00101     BlockKrylovSchurState() : curDim(0), V(Teuchos::null),
00102                               H(Teuchos::null), S(Teuchos::null),
00103                               Q(Teuchos::null) {}
00104   };
00105 
00107 
00109 
00110 
00123   class BlockKrylovSchurInitFailure : public AnasaziError {public:
00124     BlockKrylovSchurInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
00125     {}};
00126 
00133   class BlockKrylovSchurOrthoFailure : public AnasaziError {public:
00134     BlockKrylovSchurOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
00135     {}};
00136   
00138 
00139 
00140   template <class ScalarType, class MV, class OP>
00141   class BlockKrylovSchur : public Eigensolver<ScalarType,MV,OP> { 
00142   public:
00144 
00145     
00157     BlockKrylovSchur( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 
00158                    const Teuchos::RCP<SortManager<ScalarType,MV,OP> > &sorter,
00159                    const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00160                    const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00161                    const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
00162                    Teuchos::ParameterList &params 
00163                  );
00164     
00166     virtual ~BlockKrylovSchur() {};
00168 
00169 
00171 
00172     
00194     void iterate();
00195 
00217     void initialize(BlockKrylovSchurState<ScalarType,MV> state);
00218 
00222     void initialize();
00223 
00232     bool isInitialized() const { return initialized_; }
00233 
00241     BlockKrylovSchurState<ScalarType,MV> getState() const {
00242       BlockKrylovSchurState<ScalarType,MV> state;
00243       state.curDim = curDim_;
00244       state.V = V_;
00245       state.H = H_;
00246       state.Q = Q_;
00247       state.S = schurH_;
00248       return state;
00249     }
00250     
00252 
00253 
00255 
00256 
00258     int getNumIters() const { return(iter_); }
00259 
00261     void resetNumIters() { iter_=0; }
00262 
00270     Teuchos::RCP<const MV> getRitzVectors() { return ritzVectors_; }
00271 
00279     std::vector<Value<ScalarType> > getRitzValues() { 
00280       std::vector<Value<ScalarType> > ret = ritzValues_;
00281       ret.resize(ritzIndex_.size());
00282       return ret;
00283     }
00284 
00290     std::vector<int> getRitzIndex() { return ritzIndex_; }
00291 
00296     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms() {
00297       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
00298       return ret;
00299     }
00300 
00305     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms() {
00306       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
00307       return ret;
00308     }
00309 
00314     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms() {
00315       std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
00316       ret.resize(ritzIndex_.size());
00317       return ret;
00318     }
00319 
00321 
00323 
00324 
00326     const Eigenproblem<ScalarType,MV,OP>& getProblem() const { return(*problem_); };
00327 
00334     void setSize(int blockSize, int numBlocks);
00335 
00337     void setBlockSize(int blockSize);
00338 
00340     void setStepSize(int stepSize);
00341 
00343     void setNumRitzVectors(int numRitzVecs);
00344 
00346     int getStepSize() const { return(stepSize_); }
00347 
00349     int getBlockSize() const { return(blockSize_); }
00350 
00352     int getNumRitzVectors() const { return(numRitzVecs_); }
00353 
00359     int getCurSubspaceDim() const {
00360       if (!initialized_) return 0;
00361       return curDim_;
00362     }
00363 
00365     int getMaxSubspaceDim() const { return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
00366 
00367 
00380     void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
00381 
00383     Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const {return auxVecs_;}
00384 
00386 
00388 
00389     
00391     void currentStatus(std::ostream &os);
00392 
00394 
00396 
00397     
00399     bool isRitzVecsCurrent() const { return ritzVecsCurrent_; }
00400 
00402     bool isRitzValsCurrent() const { return ritzValsCurrent_; }
00403     
00405     bool isSchurCurrent() const { return schurCurrent_; }
00406     
00408 
00410 
00411     
00413     void computeRitzVectors();
00414 
00416     void computeRitzValues();
00417     
00419     void computeSchurForm( const bool sort = true );
00420 
00422 
00423   private:
00424     //
00425     // Convenience typedefs
00426     //
00427     typedef MultiVecTraits<ScalarType,MV> MVT;
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<ScalarType,MV,OP> >      sm_;
00459     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00460     const 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   // Helper function for correctly storing the Ritz values when the eigenproblem is non-Hermitian
00544   // This allows us to use template specialization to compute the right index vector and correctly
00545   // handle complex-conjugate pairs.
00546   template<class ScalarType>
00547   void sortRitzValues( const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
00548                        const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00549                        std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI )
00550   {
00551     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00552     MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
00553 
00554     int curDim = (int)rRV.size();
00555     int i = 0;
00556 
00557     // Clear the current index.
00558     RI->clear();
00559 
00560     // Place the Ritz values from rRV and iRV into the RV container.
00561     while( i < curDim ) {
00562       if ( iRV[i] != MT_ZERO ) {
00563         //
00564         // We will have this situation for real-valued, non-Hermitian matrices.
00565         (*RV)[i].set(rRV[i], iRV[i]);
00566         (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
00567         
00568         // Make sure that complex conjugate pairs have their positive imaginary part first.
00569         if ( (*RV)[i].imagpart < MT_ZERO ) {
00570           // The negative imaginary part is first, so swap the order of the ritzValues and ritzOrders.
00571           Anasazi::Value<ScalarType> tmp_ritz( (*RV)[i] );
00572           (*RV)[i] = (*RV)[i+1];
00573           (*RV)[i+1] = tmp_ritz;
00574           
00575           int tmp_order = (*RO)[i];
00576           (*RO)[i] = (*RO)[i+1];
00577           (*RO)[i+1] = tmp_order;
00578           
00579         }
00580         RI->push_back(1); RI->push_back(-1);
00581         i = i+2;
00582       } else {
00583         //
00584         // The Ritz value is not complex.
00585         (*RV)[i].set(rRV[i], MT_ZERO);
00586         RI->push_back(0);
00587         i++;
00588       }
00589     }
00590   }
00591   
00592 #ifdef HAVE_TEUCHOS_COMPLEX
00593   // Template specialization for the complex scalar type.
00594   void sortRitzValues( const std::vector<double>& rRV, 
00595                        const std::vector<double>& iRV,
00596                        std::vector<Value<ANSZI_CPLX_CLASS<double> > >* RV, 
00597                        std::vector<int>* RO, std::vector<int>* RI )
00598   {
00599     int curDim = (int)rRV.size();
00600     int i = 0;
00601 
00602     // Clear the current index.
00603     RI->clear();
00604 
00605     // Place the Ritz values from rRV and iRV into the RV container.
00606     while( i < curDim ) {
00607       (*RV)[i].set(rRV[i], iRV[i]);
00608       RI->push_back(0);
00609       i++;
00610     }    
00611   }
00612 #endif
00613 
00615   // Helper function for correctly scaling the eigenvectors of the projected eigenproblem.
00616   // This allows us to use template specialization to compute the right scaling so the
00617   // Ritz residuals are correct.
00618   template<class ScalarType>
00619   void scaleRitzVectors( const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00620                          Teuchos::SerialDenseMatrix<int, ScalarType>* S )
00621   {
00622     ScalarType ST_ONE = Teuchos::ScalarTraits<ScalarType>::one();
00623     
00624     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00625     MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
00626 
00627     Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
00628     Teuchos::BLAS<int,ScalarType> blas;
00629     
00630     int i = 0, curDim = S->numRows();
00631     ScalarType temp;
00632     ScalarType* s_ptr = S->values();
00633     while( i < curDim ) {
00634       if ( iRV[i] != MT_ZERO ) {
00635         temp = lapack_mag.LAPY2( blas.NRM2( curDim, s_ptr+i*curDim, 1 ), 
00636                                  blas.NRM2( curDim, s_ptr+(i+1)*curDim, 1 ) );
00637         blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
00638         blas.SCAL( curDim, ST_ONE/temp, s_ptr+(i+1)*curDim, 1 );
00639         i = i+2;
00640       } else {
00641         temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
00642         blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
00643         i++;
00644       }
00645     }
00646   }
00647 
00648 #ifdef HAVE_TEUCHOS_COMPLEX
00649   // Template specialization for the complex scalar type.
00650   void scaleRitzVectors( const std::vector<double>& iRV,
00651                          Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<double> >* S )
00652   {
00653     typedef ANSZI_CPLX_CLASS<double> ST;
00654     ST ST_ONE = Teuchos::ScalarTraits<ST>::one();
00655     
00656     Teuchos::BLAS<int,ST> blas;
00657     
00658     int i = 0, curDim = S->numRows();
00659     ST temp;
00660     ST* s_ptr = S->values();
00661     while( i < curDim ) {
00662       temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
00663       blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
00664       i++;
00665     }
00666   }
00667 #endif
00668 
00670   // Helper function for correctly computing the Ritz residuals of the projected eigenproblem.
00671   // This allows us to use template specialization to ensure the Ritz residuals are correct.
00672   template<class ScalarType>
00673   void computeRitzResiduals( const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00674                              const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
00675                              std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR
00676                              )
00677   {
00678     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00679     MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
00680     
00681     Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
00682     Teuchos::BLAS<int,ScalarType> blas;
00683     
00684     int i = 0;
00685     int s_stride = S.stride();
00686     int s_rows = S.numRows();
00687     int s_cols = S.numCols();
00688     ScalarType* s_ptr = S.values();
00689 
00690     while( i < s_cols ) {
00691       if ( iRV[i] != MT_ZERO ) {
00692         (*RR)[i] = lapack_mag.LAPY2( blas.NRM2(s_rows, s_ptr + i*s_stride, 1),
00693                                      blas.NRM2(s_rows, s_ptr + (i+1)*s_stride, 1) );
00694         (*RR)[i+1] = (*RR)[i];
00695         i = i+2;
00696       } else {
00697         (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
00698         i++;
00699       }
00700     }          
00701   }
00702 
00703 #ifdef HAVE_TEUCHOS_COMPLEX
00704   // Template specialization for the complex scalar type.
00705   void computeRitzResiduals( const std::vector<double>& iRV,
00706                              const Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<double> >& S,
00707                              std::vector<double>* RR
00708                              )
00709   {
00710     Teuchos::BLAS<int,ANSZI_CPLX_CLASS<double> > blas;
00711     
00712     int s_stride = S.stride();
00713     int s_rows = S.numRows();
00714     int s_cols = S.numCols();
00715     ANSZI_CPLX_CLASS<double>* s_ptr = S.values();
00716 
00717     for (int i=0; i<s_cols; ++i ) {
00718       (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
00719     }
00720   }          
00721   
00722 #endif
00723 
00725   // Constructor
00726   template <class ScalarType, class MV, class OP>
00727   BlockKrylovSchur<ScalarType,MV,OP>::BlockKrylovSchur(
00728         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 
00729         const Teuchos::RCP<SortManager<ScalarType,MV,OP> > &sorter,
00730         const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00731         const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00732         const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
00733         Teuchos::ParameterList &params
00734         ) :
00735     MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
00736     MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00737     NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
00738     ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
00739     ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
00740     // problem, tools
00741     problem_(problem), 
00742     sm_(sorter),
00743     om_(printer),
00744     tester_(tester),
00745     orthman_(ortho),
00746     // timers, counters
00747     timerOp_(Teuchos::TimeMonitor::getNewTimer("Operation Op*x")),
00748     timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer("Sorting Ritz values")),
00749     timerCompSF_(Teuchos::TimeMonitor::getNewTimer("Computing Schur form")),
00750     timerSortSF_(Teuchos::TimeMonitor::getNewTimer("Sorting Schur form")),
00751     timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer("Computing Ritz vectors")),
00752     timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Orthogonalization")),
00753     count_ApplyOp_(0),
00754     // internal data
00755     blockSize_(0),
00756     numBlocks_(0),
00757     stepSize_(0),
00758     initialized_(false),
00759     curDim_(0),
00760     numRitzVecs_(0),
00761     auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ), 
00762     numAuxVecs_(0),
00763     iter_(0),
00764     ritzVecsCurrent_(false),
00765     ritzValsCurrent_(false),
00766     schurCurrent_(false),
00767     numRitzPrint_(0)
00768   {     
00769     TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
00770                        "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
00771     TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
00772                        "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
00773     TEST_FOR_EXCEPTION(sorter == Teuchos::null,std::invalid_argument,
00774                        "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
00775     TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
00776                        "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
00777     TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
00778                        "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
00779     TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
00780                        "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
00781 
00782     // Get problem operator
00783     Op_ = problem_->getOperator();
00784 
00785     // get the step size
00786     TEST_FOR_EXCEPTION(!params.isParameter("Step Size"), std::invalid_argument,
00787                        "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
00788     int ss = params.get("Step Size",numBlocks_);
00789     setStepSize(ss);
00790 
00791     // set the block size and allocate data
00792     int bs = params.get("Block Size", 1);
00793     int nb = params.get("Num Blocks", 3*problem_->getNEV());
00794     setSize(bs,nb);
00795 
00796     // get the number of Ritz vectors to compute and allocate data.
00797     // --> if this parameter is not specified in the parameter list, then it's assumed that no Ritz vectors will be computed.
00798     int numRitzVecs = params.get("Number of Ritz Vectors", 0);
00799     setNumRitzVectors( numRitzVecs );
00800 
00801     // get the number of Ritz values to print out when currentStatus is called.
00802     numRitzPrint_ = params.get("Print Number of Ritz Values", bs);
00803   }
00804 
00805 
00807   // Set the block size
00808   // This simply calls setSize(), modifying the block size while retaining the number of blocks.
00809   template <class ScalarType, class MV, class OP>
00810   void BlockKrylovSchur<ScalarType,MV,OP>::setBlockSize (int blockSize) 
00811   {
00812     setSize(blockSize,numBlocks_);
00813   }
00814 
00815 
00817   // Set the step size.
00818   template <class ScalarType, class MV, class OP>
00819   void BlockKrylovSchur<ScalarType,MV,OP>::setStepSize (int stepSize)
00820   {
00821     TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
00822     stepSize_ = stepSize;
00823   }
00824 
00826   // Set the number of Ritz vectors to compute.
00827   template <class ScalarType, class MV, class OP>
00828   void BlockKrylovSchur<ScalarType,MV,OP>::setNumRitzVectors (int numRitzVecs)
00829   {
00830     // This routine only allocates space; it doesn't not perform any computation
00831     // any change in size will invalidate the state of the solver.
00832 
00833     TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
00834 
00835     // Check to see if the number of requested Ritz vectors has changed.
00836     if (numRitzVecs != numRitzVecs_) {
00837       if (numRitzVecs) {
00838         ritzVectors_ = Teuchos::null;
00839         ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
00840       } else {
00841         ritzVectors_ = Teuchos::null;
00842       }
00843       numRitzVecs_ = numRitzVecs;
00844       ritzVecsCurrent_ = false;
00845     }      
00846   }
00847 
00849   // Set the block size and make necessary adjustments.
00850   template <class ScalarType, class MV, class OP>
00851   void BlockKrylovSchur<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks) 
00852   {
00853     // This routine only allocates space; it doesn't not perform any computation
00854     // any change in size will invalidate the state of the solver.
00855 
00856     TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
00857     TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
00858     if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00859       // do nothing
00860       return;
00861     }
00862 
00863     blockSize_ = blockSize;
00864     numBlocks_ = numBlocks;
00865 
00866     Teuchos::RCP<const MV> tmp;
00867     // grab some Multivector to Clone
00868     // in practice, getInitVec() should always provide this, but it is possible to use a 
00869     // Eigenproblem with nothing in getInitVec() by manually initializing with initialize(); 
00870     // in case of that strange scenario, we will try to Clone from V_; first resort to getInitVec(), 
00871     // because we would like to clear the storage associated with V_ so we have room for the new V_
00872     if (problem_->getInitVec() != Teuchos::null) {
00873       tmp = problem_->getInitVec();
00874     }
00875     else {
00876       tmp = V_;
00877       TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00878           "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
00879     }
00880 
00881 
00883     // blockSize*numBlocks dependent
00884     //
00885     int newsd;
00886     if (problem_->isHermitian()) {
00887       newsd = blockSize_*numBlocks_;
00888     } else {
00889       newsd = blockSize_*numBlocks_+1;
00890     }
00891     // check that new size is valid
00892     TEST_FOR_EXCEPTION(newsd > MVT::GetVecLength(*tmp),std::invalid_argument,
00893         "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
00894 
00895     ritzValues_.resize(newsd);
00896     ritzResiduals_.resize(newsd,MT_ONE);
00897     ritzOrder_.resize(newsd);
00898     V_ = Teuchos::null;
00899     V_ = MVT::Clone(*tmp,newsd+blockSize_);
00900     H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd+blockSize_,newsd) );
00901     Q_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
00902 
00903     initialized_ = false;
00904     curDim_ = 0;
00905   }
00906 
00907 
00909   // Set the auxiliary vectors
00910   template <class ScalarType, class MV, class OP>
00911   void BlockKrylovSchur<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00912     typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
00913 
00914     // set new auxiliary vectors
00915     auxVecs_ = auxvecs;
00916     
00917     if (om_->isVerbosity( Debug ) ) {
00918       // Check almost everything here
00919       CheckList chk;
00920       chk.checkAux = true;
00921       om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
00922     }
00923 
00924     numAuxVecs_ = 0;
00925     for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
00926       numAuxVecs_ += MVT::GetNumberVecs(**i);
00927     }
00928     
00929     // If the solver has been initialized, X and P are not necessarily orthogonal to new auxiliary vectors
00930     if (numAuxVecs_ > 0 && initialized_) {
00931       initialized_ = false;
00932     }
00933   }
00934 
00936   /* Initialize the state of the solver
00937    * 
00938    * POST-CONDITIONS:
00939    *
00940    * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
00941    *
00942    */
00943 
00944   template <class ScalarType, class MV, class OP>
00945   void BlockKrylovSchur<ScalarType,MV,OP>::initialize(BlockKrylovSchurState<ScalarType,MV> newstate)
00946   {
00947     // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
00948 
00949     std::vector<int> bsind(blockSize_);
00950     for (int i=0; i<blockSize_; i++) bsind[i] = i;
00951 
00952     // in BlockKrylovSchur, V and H are required.  
00953     // if either doesn't exist, then we will start with the initial vector.
00954     //
00955     // inconsistent multivectors widths and lengths will not be tolerated, and
00956     // will be treated with exceptions.
00957     //
00958     std::string errstr("Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
00959 
00960     // set up V,H: if the user doesn't specify both of these these, 
00961     // we will start over with the initial vector.
00962     if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
00963 
00964       // initialize V_,H_, and curDim_
00965 
00966       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_),
00967                           std::invalid_argument, errstr );
00968       if (newstate.V != V_) {
00969         TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
00970             std::invalid_argument, errstr );
00971         TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) > getMaxSubspaceDim(),
00972             std::invalid_argument, errstr );
00973       }
00974       TEST_FOR_EXCEPTION( newstate.curDim > getMaxSubspaceDim(),
00975                           std::invalid_argument, errstr );
00976 
00977       curDim_ = newstate.curDim;
00978       int lclDim = MVT::GetNumberVecs(*newstate.V);
00979 
00980       // check size of H
00981       TEST_FOR_EXCEPTION(newstate.H->numRows() < curDim_ || newstate.H->numCols() < curDim_, std::invalid_argument, errstr);
00982       
00983       if (curDim_ == 0 && lclDim > blockSize_) {
00984         om_->stream(Warnings) << "Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
00985                                   << "The block size however is only " << blockSize_ << std::endl
00986                                   << "The last " << lclDim - blockSize_ << " vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
00987       }
00988 
00989 
00990       // copy basis vectors from newstate into V
00991       if (newstate.V != V_) {
00992         std::vector<int> nevind(lclDim);
00993         for (int i=0; i<lclDim; i++) nevind[i] = i;
00994         MVT::SetBlock(*newstate.V,nevind,*V_);
00995       }
00996 
00997       // put data into H_, make sure old information is not still hanging around.
00998       if (newstate.H != H_) {
00999         H_->putScalar( ST_ZERO );
01000         Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H,curDim_+blockSize_,curDim_);
01001         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
01002         lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_,curDim_+blockSize_,curDim_) );
01003         lclH->assign(newH);
01004 
01005         // done with local pointers
01006         lclH = Teuchos::null;
01007       }
01008 
01009     }
01010     else {
01011       // user did not specify a basis V
01012       // get vectors from problem or generate something, projectAndNormalize, call initialize() recursively
01013       Teuchos::RCP<const MV> ivec = problem_->getInitVec();
01014       TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
01015                          "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
01016 
01017       int lclDim = MVT::GetNumberVecs(*ivec);
01018       bool userand = false;
01019       if (lclDim < blockSize_) {
01020         // we need at least blockSize_ vectors
01021         // use a random multivec
01022         userand = true;
01023       }
01024               
01025       if (userand) {
01026         // make an index
01027         std::vector<int> dimind2(lclDim);
01028         for (int i=0; i<lclDim; i++) { dimind2[i] = i; }
01029 
01030         // alloc newV as a view of the first block of V
01031         Teuchos::RCP<MV> newV1 = MVT::CloneView(*V_,dimind2);
01032 
01033         // copy the initial vectors into the first lclDim vectors of V
01034         MVT::SetBlock(*ivec,dimind2,*newV1);
01035 
01036         // resize / reinitialize the index vector        
01037         dimind2.resize(blockSize_-lclDim);
01038         for (int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
01039 
01040         // initialize the rest of the vectors with random vectors
01041         Teuchos::RCP<MV> newV2 = MVT::CloneView(*V_,dimind2);
01042         MVT::MvRandom(*newV2);
01043       }
01044       else {
01045         // alloc newV as a view of the first block of V
01046         Teuchos::RCP<MV> newV1 = MVT::CloneView(*V_,bsind);
01047        
01048         // get a view of the first block of initial vectors
01049         Teuchos::RCP<const MV> ivecV = MVT::CloneView(*ivec,bsind);
01050  
01051         // assign ivec to first part of newV
01052         MVT::SetBlock(*ivecV,bsind,*newV1);
01053       }
01054 
01055       // get pointer into first block of V
01056       Teuchos::RCP<MV> newV = MVT::CloneView(*V_,bsind);
01057 
01058       // remove auxVecs from newV and normalize newV
01059       if (auxVecs_.size() > 0) {
01060         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01061         
01062         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
01063         int rank = orthman_->projectAndNormalize(*newV,dummy,Teuchos::null,auxVecs_);
01064         TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
01065                             "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
01066       }
01067       else {
01068         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01069 
01070         int rank = orthman_->normalize(*newV,Teuchos::null);
01071         TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
01072                             "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
01073       }
01074 
01075       // set curDim
01076       curDim_ = 0;
01077 
01078       // clear pointer
01079       newV = Teuchos::null;
01080     }
01081 
01082     // The Ritz vectors/values and Schur form are no longer current.
01083     ritzVecsCurrent_ = false;
01084     ritzValsCurrent_ = false;
01085     schurCurrent_ = false;
01086 
01087     // the solver is initialized
01088     initialized_ = true;
01089 
01090     if (om_->isVerbosity( Debug ) ) {
01091       // Check almost everything here
01092       CheckList chk;
01093       chk.checkV = true;
01094       chk.checkArn = true;
01095       chk.checkAux = true;
01096       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
01097     }
01098 
01099     // Print information on current status
01100     if (om_->isVerbosity(Debug)) {
01101       currentStatus( om_->stream(Debug) );
01102     }
01103     else if (om_->isVerbosity(IterationDetails)) {
01104       currentStatus( om_->stream(IterationDetails) );
01105     }
01106   }
01107 
01108 
01110   // initialize the solver with default state
01111   template <class ScalarType, class MV, class OP>
01112   void BlockKrylovSchur<ScalarType,MV,OP>::initialize()
01113   {
01114     BlockKrylovSchurState<ScalarType,MV> empty;
01115     initialize(empty);
01116   }
01117 
01118 
01120   // Perform BlockKrylovSchur iterations until the StatusTest tells us to stop.
01121   template <class ScalarType, class MV, class OP>
01122   void BlockKrylovSchur<ScalarType,MV,OP>::iterate()
01123   {
01124     //
01125     // Allocate/initialize data structures
01126     //
01127     if (initialized_ == false) {
01128       initialize();
01129     }
01130 
01131     // Compute the current search dimension. 
01132     // If the problem is non-Hermitian and the blocksize is one, let the solver use the extra vector.
01133     int searchDim = blockSize_*numBlocks_;
01134     if (problem_->isHermitian() == false) {
01135       searchDim++;
01136     } 
01137 
01139     // iterate until the status test tells us to stop.
01140     //
01141     // also break if our basis is full
01142     //
01143     while (tester_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
01144 
01145       iter_++;
01146 
01147       // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
01148       int lclDim = curDim_ + blockSize_; 
01149 
01150       // Get the current part of the basis.
01151       std::vector<int> curind(blockSize_);
01152       for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
01153       Teuchos::RCP<MV> Vnext = MVT::CloneView(*V_,curind);
01154 
01155       // Get a view of the previous vectors
01156       // this is used for orthogonalization and for computing V^H K H
01157       for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
01158       Teuchos::RCP<MV> Vprev = MVT::CloneView(*V_,curind);
01159 
01160       // Compute the next vector in the Krylov basis:  Vnext = Op*Vprev
01161       {
01162         Teuchos::TimeMonitor lcltimer( *timerOp_ );
01163         OPT::Apply(*Op_,*Vprev,*Vnext);
01164         count_ApplyOp_ += blockSize_;
01165       }
01166       Vprev = Teuchos::null;
01167       
01168       // Remove all previous Krylov-Schur basis vectors and auxVecs from Vnext
01169       {
01170         Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01171         
01172         // Get a view of all the previous vectors
01173         std::vector<int> prevind(lclDim);
01174         for (int i=0; i<lclDim; i++) { prevind[i] = i; }
01175         Vprev = MVT::CloneView(*V_,prevind);
01176         Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
01177         
01178         // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
01179         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
01180           subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
01181                                ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
01182         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
01183         AsubH.append( subH );
01184         
01185         // Add the auxiliary vectors to the current basis vectors if any exist
01186         if (auxVecs_.size() > 0) {
01187           for (unsigned int i=0; i<auxVecs_.size(); i++) {
01188             AVprev.append( auxVecs_[i] );
01189             AsubH.append( Teuchos::null );
01190           }
01191         }
01192         
01193         // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
01194         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
01195           subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
01196                                ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
01197         int rank = orthman_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
01198         TEST_FOR_EXCEPTION(rank != blockSize_,BlockKrylovSchurOrthoFailure,
01199                            "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
01200       }
01201       //
01202       // V has been extended, and H has been extended. 
01203       //
01204       // Update basis dim and release all pointers.
01205       Vnext = Teuchos::null;
01206       curDim_ += blockSize_;
01207       // The Ritz vectors/values and Schur form are no longer current.
01208       ritzVecsCurrent_ = false;
01209       ritzValsCurrent_ = false;
01210       schurCurrent_ = false;
01211       //
01212       // Update Ritz values and residuals if needed
01213       if (!(iter_%stepSize_)) {
01214         computeRitzValues();
01215       }
01216       
01217       // When required, monitor some orthogonalities
01218       if (om_->isVerbosity( Debug ) ) {
01219         // Check almost everything here
01220         CheckList chk;
01221         chk.checkV = true;
01222         chk.checkArn = true;
01223         om_->print( Debug, accuracyCheck(chk, ": after local update") );
01224       }
01225       else if (om_->isVerbosity( OrthoDetails ) ) {
01226         CheckList chk;
01227         chk.checkV = true;
01228         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
01229       }
01230       
01231       // Print information on current iteration
01232       if (om_->isVerbosity(Debug)) {
01233         currentStatus( om_->stream(Debug) );
01234       }
01235       else if (om_->isVerbosity(IterationDetails)) {
01236         currentStatus( om_->stream(IterationDetails) );
01237       }
01238       
01239     } // end while (statusTest == false)
01240     
01241   } // end of iterate()
01242 
01243 
01245   // Check accuracy, orthogonality, and other debugging stuff
01246   // 
01247   // bools specify which tests we want to run (instead of running more than we actually care about)
01248   //
01249   // checkV : V orthonormal
01250   //          orthogonal to auxvecs
01251   // checkAux: check that auxiliary vectors are actually orthonormal
01252   //
01253   // checkArn: check the Arnoldi factorization
01254   //
01255   // NOTE:  This method needs to check the current dimension of the subspace, since it is possible to
01256   //        call this method when curDim_ = 0 (after initialization).
01257   template <class ScalarType, class MV, class OP>
01258   std::string BlockKrylovSchur<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const 
01259   {
01260     std::stringstream os;
01261     os.precision(2);
01262     os.setf(std::ios::scientific, std::ios::floatfield);
01263     MagnitudeType tmp;
01264 
01265     os << " Debugging checks: iteration " << iter_ << where << std::endl;
01266 
01267     // index vectors for V and F
01268     std::vector<int> lclind(curDim_);
01269     for (int i=0; i<curDim_; i++) lclind[i] = i;
01270     std::vector<int> bsind(blockSize_);
01271     for (int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
01272     
01273     Teuchos::RCP<MV> lclV,lclF,lclAV;
01274     if (curDim_)
01275       lclV = MVT::CloneView(*V_,lclind);
01276     lclF = MVT::CloneView(*V_,bsind);
01277 
01278     if (chk.checkV) {
01279       if (curDim_) {
01280         tmp = orthman_->orthonormError(*lclV);
01281         os << " >> Error in V^H M V == I  : " << tmp << std::endl;
01282       }
01283       tmp = orthman_->orthonormError(*lclF);
01284       os << " >> Error in F^H M F == I  : " << tmp << std::endl;
01285       if (curDim_) {
01286         tmp = orthman_->orthogError(*lclV,*lclF);
01287         os << " >> Error in V^H M F == 0  : " << tmp << std::endl;
01288       }
01289       for (unsigned int i=0; i<auxVecs_.size(); i++) {
01290         if (curDim_) {
01291           tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
01292           os << " >> Error in V^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01293         }
01294         tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
01295         os << " >> Error in F^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01296       }
01297     }
01298     
01299     if (chk.checkArn) {
01300 
01301       if (curDim_) {
01302         // Compute AV      
01303         Teuchos::RCP<MV> lclAV = MVT::Clone(*V_,curDim_);
01304         {
01305           Teuchos::TimeMonitor lcltimer( *timerOp_ );
01306           OPT::Apply(*Op_,*lclV,*lclAV);
01307         }
01308         
01309         // Compute AV - VH
01310         Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
01311         MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
01312         
01313         // Compute FB_k^T - (AV-VH)
01314         Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
01315                                                         blockSize_,curDim_, curDim_ );
01316         MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
01317         
01318         // Compute || FE_k^T - (AV-VH) ||
01319         std::vector<MagnitudeType> arnNorms( curDim_ );
01320         orthman_->norm( *lclAV, &arnNorms );
01321         
01322         for (int i=0; i<curDim_; i++) {        
01323         os << " >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
01324         }
01325       }
01326     }
01327 
01328     if (chk.checkAux) {
01329       for (unsigned int i=0; i<auxVecs_.size(); i++) {
01330         tmp = orthman_->orthonormError(*auxVecs_[i]);
01331         os << " >> Error in Aux[" << i << "]^H M Aux[" << i << "] == I : " << tmp << std::endl;
01332         for (unsigned int j=i+1; j<auxVecs_.size(); j++) {
01333           tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01334           os << " >> Error in Aux[" << i << "]^H M Aux[" << j << "] == 0 : " << tmp << std::endl;
01335         }
01336       }
01337     }
01338 
01339     os << std::endl;
01340 
01341     return os.str();
01342   }
01343 
01345   /* Get the current approximate eigenvalues, i.e. Ritz values.
01346    * 
01347    * POST-CONDITIONS:
01348    *
01349    * ritzValues_ contains Ritz w.r.t. V, H
01350    * Q_ contains the Schur vectors w.r.t. H
01351    * schurH_ contains the Schur matrix w.r.t. H
01352    * ritzOrder_ contains the current ordering from sort manager
01353    */
01354 
01355   template <class ScalarType, class MV, class OP>  
01356   void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzValues()
01357   {
01358     // Can only call this if the solver is initialized
01359     if (initialized_) {
01360 
01361       // This just updates the Ritz values and residuals.
01362       // --> ritzValsCurrent_ will be set to 'true' by this method.
01363       if (!ritzValsCurrent_) {
01364         // Compute the current Ritz values, through computing the Schur form
01365         //   without updating the current projection matrix or sorting the Schur form.
01366         computeSchurForm( false );
01367       }
01368     }
01369   }
01370 
01372   /* Get the current approximate eigenvectors, i.e. Ritz vectors.
01373    * 
01374    * POST-CONDITIONS:
01375    *
01376    * ritzValues_ contains Ritz w.r.t. V, H
01377    * ritzVectors_ is first blockSize_ Ritz vectors w.r.t. V, H
01378    * Q_ contains the Schur vectors w.r.t. H
01379    * schurH_ contains the Schur matrix w.r.t. H
01380    * ritzOrder_ contains the current ordering from sort manager
01381    */
01382 
01383   template <class ScalarType, class MV, class OP>  
01384   void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzVectors()
01385   {
01386     Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
01387 
01388     TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
01389                        "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
01390 
01391     TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
01392                        "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
01393 
01394 
01395     // Check to see if the current subspace dimension is non-trivial and the solver is initialized
01396     if (curDim_ && initialized_) {
01397 
01398       // Check to see if the Ritz vectors are current.
01399       if (!ritzVecsCurrent_) {
01400         
01401         // Check to see if the Schur factorization of H (schurH_, Q) is current and sorted.
01402         if (!schurCurrent_) {
01403           // Compute the Schur factorization of the current H, which will not directly change H,
01404           // the factorization will be sorted and placed in (schurH_, Q)
01405           computeSchurForm( true );
01406         }
01407         
01408         // After the Schur form is computed, then the Ritz values are current.
01409         // Thus, I can check the Ritz index vector to see if I have enough space for the Ritz vectors requested.
01410         TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
01411                            "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
01412 
01413         Teuchos::LAPACK<int,ScalarType> lapack;
01414         Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01415 
01416         // Compute the Ritz vectors.
01417         // --> For a Hermitian problem this is simply the current basis times the first numRitzVecs_ Schur vectors
01418         //     
01419         // --> For a non-Hermitian problem, this involves solving the projected eigenproblem, then
01420         //     placing the product of the current basis times the first numRitzVecs_ Schur vectors times the
01421         //     eigenvectors of interest into the Ritz vectors.
01422 
01423         // Get a view of the current Krylov-Schur basis vectors and Schur vectors
01424         std::vector<int> curind( curDim_ );
01425         for (int i=0; i<curDim_; i++) { curind[i] = i; }
01426         Teuchos::RCP<MV> Vtemp = MVT::CloneView( *V_, curind );
01427         if (problem_->isHermitian()) {
01428           // Get a view into the current Schur vectors
01429           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
01430 
01431           // Compute the current Ritz vectors      
01432           MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
01433           
01434         } else {
01435 
01436           // Get a view into the current Schur vectors.
01437           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01438           
01439           // Get a set of work vectors to hold the current Ritz vectors.
01440           Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
01441 
01442           // Compute the current Krylov-Schur vectors.
01443           MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );          
01444 
01445           //  Now compute the eigenvectors of the Schur form
01446           //  Reset the dense matrix and compute the eigenvalues of the Schur form.
01447           //
01448           // Allocate the work space. This space will be used below for calls to:
01449           // * TREVC (requires 3*N for real, 2*N for complex) 
01450 
01451           int lwork = 3*curDim_;
01452           std::vector<ScalarType> work( lwork );
01453           std::vector<MagnitudeType> rwork( curDim_ );
01454           char side = 'R';
01455           int mm, info = 0; 
01456           const int ldvl = 1;
01457           ScalarType vl[ ldvl ];
01458           Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
01459           lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01460                         copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01461           TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01462                              "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC returned info != 0.");
01463 
01464           // Get a view into the eigenvectors of the Schur form
01465           Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
01466           
01467           // Convert back to Ritz vectors of the operator.
01468           std::vector<int> curind( (numRitzVecs_) );
01469           for (int i=0; i<(int)curind.size(); i++) { curind[i] = i; }
01470 
01471           Teuchos::RCP<MV> view_ritzVectors = MVT::CloneView( *ritzVectors_, curind );
01472           MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
01473 
01474           // Compute the norm of the new Ritz vectors
01475           std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
01476           MVT::MvNorm( *view_ritzVectors, &ritzNrm );
01477 
01478           // Release memory used to compute Ritz vectors before scaling the current vectors.
01479           tmpritzVectors_ = Teuchos::null;
01480           view_ritzVectors = Teuchos::null;
01481           
01482           // Scale the Ritz vectors to have Euclidean norm.
01483           ScalarType ritzScale = ST_ONE;
01484           for (int i=0; i<numRitzVecs_; i++) {
01485             
01486             // If this is a conjugate pair then normalize by the real and imaginary parts.
01487             if (ritzIndex_[i] == 1 ) {
01488               ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
01489               std::vector<int> newind(2);
01490               newind[0] = i; newind[1] = i+1;
01491               tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01492               view_ritzVectors = MVT::CloneView( *ritzVectors_, newind );
01493               MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01494 
01495               // Increment counter for imaginary part
01496               i++;
01497             } else {
01498 
01499               // This is a real Ritz value, normalize the vector
01500               std::vector<int> newind(1);
01501               newind[0] = i;
01502               tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01503               view_ritzVectors = MVT::CloneView( *ritzVectors_, newind );
01504               MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01505             }              
01506           }
01507           
01508         } // if (problem_->isHermitian()) 
01509         
01510         // The current Ritz vectors have been computed.
01511         ritzVecsCurrent_ = true;
01512         
01513       } // if (!ritzVecsCurrent_)      
01514     } // if (curDim_)    
01515   } // computeRitzVectors()
01516   
01517 
01519   /* Get the current approximate eigenvalues, i.e. Ritz values.
01520    * 
01521    * POST-CONDITIONS:
01522    *
01523    * ritzValues_ contains Ritz w.r.t. V, H
01524    * Q_ contains the Schur vectors w.r.t. H
01525    * schurH_ contains the Schur matrix w.r.t. H
01526    * ritzOrder_ contains the current ordering from sort manager
01527    * schurCurrent_ = true if sort = true; i.e. the Schur form is sorted according to the index
01528    *  vector returned by the sort manager.
01529    */
01530   template <class ScalarType, class MV, class OP>
01531   void BlockKrylovSchur<ScalarType,MV,OP>::computeSchurForm( const bool sort )
01532   {
01533     // local timer
01534     Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
01535 
01536     // Check to see if the dimension of the factorization is greater than zero.
01537     if (curDim_) {
01538 
01539       // Check to see if the Schur factorization is current.
01540       if (!schurCurrent_) {
01541         
01542         // Check to see if the Ritz values are current
01543         // --> If they are then the Schur factorization is current but not sorted.
01544         if (!ritzValsCurrent_) {
01545           Teuchos::LAPACK<int,ScalarType> lapack; 
01546           Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01547           Teuchos::BLAS<int,ScalarType> blas;
01548           Teuchos::BLAS<int,MagnitudeType> blas_mag;
01549           
01550           // Get a view into Q, the storage for H's Schur vectors.
01551           Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01552           
01553           // Get a copy of H to compute/sort the Schur form.
01554           schurH_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
01555           //
01556           //---------------------------------------------------
01557           // Compute the Schur factorization of subH
01558           // ---> Use driver GEES to first reduce to upper Hessenberg 
01559           //         form and then compute Schur form, outputting Ritz values
01560           //---------------------------------------------------
01561           //
01562           // Allocate the work space. This space will be used below for calls to:
01563           // * GEES  (requires 3*N for real, 2*N for complex)
01564           // * TREVC (requires 3*N for real, 2*N for complex) 
01565           // * TREXC (requires N for real, none for complex)
01566           // Furthermore, GEES requires a real array of length curDim_ (for complex datatypes)
01567           //
01568           int lwork = 3*curDim_;
01569           std::vector<ScalarType> work( lwork );
01570           std::vector<MagnitudeType> rwork( curDim_ );
01571           std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
01572           std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
01573           std::vector<int> bwork( curDim_ );
01574           int info = 0, sdim = 0; 
01575           char jobvs = 'V';
01576           lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
01577                        &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork, 
01578                        &rwork[0], &bwork[0], &info );
01579           
01580           TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01581                              "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES returned info != 0.");
01582           //
01583           //---------------------------------------------------
01584           // Use the Krylov-Schur factorization to compute the current Ritz residuals 
01585           // for ALL the eigenvalues estimates (Ritz values)
01586           //           || Ax - x\theta || = || U_m+1*B_m+1^H*Q*s || 
01587           //                              = || B_m+1^H*Q*s ||
01588           //
01589           // where U_m+1 is the current Krylov-Schur basis, Q are the Schur vectors, and x = U_m+1*Q*s
01590           // NOTE: This means that s = e_i if the problem is hermitian, else the eigenvectors
01591           //       of the Schur form need to be computed.
01592           //
01593           // First compute H_{m+1,m}*B_m^T, then determine what 's' is.
01594           //---------------------------------------------------
01595           //
01596           // Get current B_m+1
01597           Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View, *H_,
01598                                                           blockSize_, curDim_, curDim_ );
01599           //
01600           // Compute B_m+1^H*Q
01601           Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
01602           blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE, 
01603                      curB.values(), curB.stride(), subQ.values(), subQ.stride(), 
01604                      ST_ZERO, subB.values(), subB.stride() );
01605           //
01606           // Determine what 's' is and compute Ritz residuals.
01607           //
01608           ScalarType* b_ptr = subB.values();
01609           if (problem_->isHermitian()) {
01610             //
01611             // 's' is the i-th canonical basis vector.
01612             //
01613             for (int i=0; i<curDim_ ; i++) {
01614               ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
01615             }   
01616           } else {
01617             //
01618             //  Compute S: the eigenvectors of the block upper triangular, Schur matrix.
01619             //
01620             char side = 'R';
01621             int mm;
01622             const int ldvl = 1;
01623             ScalarType vl[ ldvl ];
01624             Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
01625             lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01626                           S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01627             
01628             TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01629                                "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC returned info != 0.");
01630             //
01631             // Scale the eigenvectors so that their Euclidean norms are all one.
01632             //
01633             scaleRitzVectors( tmp_iRitzValues, &S );
01634             //
01635             // Compute ritzRes = *B_m+1^H*Q*S where the i-th column of S is 's' for the i-th Ritz-value
01636             //
01637             Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
01638             blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE, 
01639                        subB.values(), subB.stride(), S.values(), S.stride(), 
01640                        ST_ZERO, ritzRes.values(), ritzRes.stride() );
01641 
01642             /* TO DO:  There's be an incorrect assumption made in the computation of the Ritz residuals.
01643                        This assumption is that the next vector in the Krylov subspace is Euclidean orthonormal.
01644                        It may not be normalized using Euclidean norm.
01645             Teuchos::RCP<MV> ritzResVecs = MVT::Clone( *V_, curDim_ );
01646             std::vector<int> curind(blockSize_);
01647             for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
01648             Teuchos::RCP<MV> Vtemp = MVT::CloneView(*V_,curind);     
01649             
01650             MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, ritzRes, ST_ZERO, *ritzResVecs );
01651             std::vector<MagnitudeType> ritzResNrms(curDim_);
01652             MVT::MvNorm( *ritzResVecs, &ritzResNrms );
01653             i = 0;
01654             while( i < curDim_ ) {
01655               if ( tmp_ritzValues[curDim_+i] != MT_ZERO ) {
01656                 ritzResiduals_[i] = lapack_mag.LAPY2( ritzResNrms[i], ritzResNrms[i+1] );
01657                 ritzResiduals_[i+1] = ritzResiduals_[i];
01658                 i = i+2;
01659               } else {
01660                 ritzResiduals_[i] = ritzResNrms[i];
01661                 i++;
01662               }
01663             }
01664             */
01665             //
01666             // Compute the Ritz residuals for each Ritz value.
01667             // 
01668             computeRitzResiduals( tmp_iRitzValues, ritzRes, &ritzResiduals_ );
01669           }
01670           //
01671           // Sort the Ritz values.
01672           //
01673           {
01674             Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
01675             int i=0;
01676             if (problem_->isHermitian()) {
01677               //
01678               // Sort using just the real part of the Ritz values.
01679               sm_->sort( this, curDim_, tmp_rRitzValues, &ritzOrder_ ); // don't catch exception
01680               ritzIndex_.clear();
01681               while ( i < curDim_ ) {
01682                 // The Ritz value is not complex.
01683                 ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
01684                 ritzIndex_.push_back(0);
01685                 i++;
01686               }
01687             }
01688             else {
01689               //
01690               // Sort using both the real and imaginary parts of the Ritz values.
01691               sm_->sort( this, curDim_, tmp_rRitzValues, tmp_iRitzValues, &ritzOrder_ );
01692               sortRitzValues( tmp_rRitzValues, tmp_iRitzValues, &ritzValues_, &ritzOrder_, &ritzIndex_ );
01693             }
01694             //
01695             // Sort the ritzResiduals_ based on the ordering from the Sort Manager.
01696             std::vector<MagnitudeType> ritz2( curDim_ );
01697             for (int i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
01698             blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
01699             
01700             // The Ritz values have now been updated.
01701             ritzValsCurrent_ = true;
01702           }
01703 
01704         } // if (!ritzValsCurrent_) ...
01705         // 
01706         //---------------------------------------------------
01707         // * The Ritz values and residuals have been updated at this point.
01708         // 
01709         // * The Schur factorization of the projected matrix has been computed,
01710         //   and is stored in (schurH_, Q_).
01711         //
01712         // Now the Schur factorization needs to be sorted.
01713         //---------------------------------------------------
01714         //
01715         // Sort the Schur form using the ordering from the Sort Manager.
01716         if (sort) {
01717           sortSchurForm( *schurH_, *Q_, ritzOrder_ );    
01718           //
01719           // Indicate the Schur form in (schurH_, Q_) is current and sorted
01720           schurCurrent_ = true;
01721         }
01722       } // if (!schurCurrent_) ...
01723   
01724     } // if (curDim_) ...
01725   
01726   } // computeSchurForm( ... )
01727   
01728 
01730   // Sort the Schur form of H stored in (H,Q) using the ordering vector.
01731   template <class ScalarType, class MV, class OP>
01732   void BlockKrylovSchur<ScalarType,MV,OP>::sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
01733                                                           Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
01734                                                           std::vector<int>& order ) 
01735   {
01736     // local timer
01737     Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
01738     //
01739     //---------------------------------------------------
01740     // Reorder real Schur factorization, remember to add one to the indices for the
01741     // fortran call and determine offset.  The offset is necessary since the TREXC
01742     // method reorders in a nonsymmetric fashion, thus we use the reordering in
01743     // a stack-like fashion.  Also take into account conjugate pairs, which may mess
01744     // up the reordering, since the pair is moved if one of the pair is moved.
01745     //---------------------------------------------------
01746     //
01747     int i = 0, nevtemp = 0;
01748     char compq = 'V';
01749     std::vector<int> offset2( curDim_ );
01750     std::vector<int> order2( curDim_ );
01751 
01752     // LAPACK objects.
01753     Teuchos::LAPACK<int,ScalarType> lapack; 
01754     int lwork = 3*curDim_;
01755     std::vector<ScalarType> work( lwork );
01756 
01757     while (i < curDim_) {
01758       if ( ritzIndex_[i] != 0 ) { // This is the first value of a complex conjugate pair
01759         offset2[nevtemp] = 0;
01760         for (int j=i; j<curDim_; j++) {
01761           if (order[j] > order[i]) { offset2[nevtemp]++; }
01762         }
01763         order2[nevtemp] = order[i];
01764         i = i+2;
01765       } else {
01766         offset2[nevtemp] = 0;
01767         for (int j=i; j<curDim_; j++) {
01768           if (order[j] > order[i]) { offset2[nevtemp]++; }
01769         }
01770         order2[nevtemp] = order[i];
01771         i++;
01772       }
01773       nevtemp++;
01774     }
01775     ScalarType *ptr_h = H.values();
01776     ScalarType *ptr_q = Q.values();
01777     int ldh = H.stride(), ldq = Q.stride();
01778     int info = 0;
01779     for (i=nevtemp-1; i>=0; i--) {
01780       lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, order2[i]+1+offset2[i], 
01781                     1, &work[0], &info );
01782       TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01783                          "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC returned info != 0.");
01784     }
01785   }
01786 
01788   // Print the current status of the solver
01789   template <class ScalarType, class MV, class OP>
01790   void BlockKrylovSchur<ScalarType,MV,OP>::currentStatus(std::ostream &os) 
01791   {
01792     using std::endl;
01793 
01794     os.setf(std::ios::scientific, std::ios::floatfield);
01795     os.precision(6);
01796     os <<"================================================================================" << endl;
01797     os << endl;
01798     os <<"                         BlockKrylovSchur Solver Status" << endl;
01799     os << endl;
01800     os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01801     os <<"The number of iterations performed is " <<iter_<<endl;
01802     os <<"The block size is         " << blockSize_<<endl;
01803     os <<"The number of blocks is   " << numBlocks_<<endl;
01804     os <<"The current basis size is " << curDim_<<endl;
01805     os <<"The number of auxiliary vectors is    " << numAuxVecs_ << endl;
01806     os <<"The number of operations Op*x   is "<<count_ApplyOp_<<endl;
01807 
01808     os.setf(std::ios_base::right, std::ios_base::adjustfield);
01809 
01810     os << endl;
01811     if (initialized_) {
01812       os <<"CURRENT RITZ VALUES             "<<endl;
01813       if (ritzIndex_.size() != 0) {
01814         int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
01815         if (problem_->isHermitian()) {
01816           os << std::setw(20) << "Ritz Value" 
01817              << std::setw(20) << "Ritz Residual"
01818              << endl;
01819           os <<"--------------------------------------------------------------------------------"<<endl;
01820           for (int i=0; i<numPrint; i++) {
01821             os << std::setw(20) << ritzValues_[i].realpart 
01822                << std::setw(20) << ritzResiduals_[i] 
01823                << endl;
01824           }
01825         } else {
01826           os << std::setw(24) << "Ritz Value" 
01827              << std::setw(30) << "Ritz Residual"
01828              << endl;
01829           os <<"--------------------------------------------------------------------------------"<<endl;
01830           for (int i=0; i<numPrint; i++) {
01831             // Print out the real eigenvalue.
01832             os << std::setw(15) << ritzValues_[i].realpart;
01833             if (ritzValues_[i].imagpart < MT_ZERO) {
01834               os << " - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
01835             } else {
01836               os << " + i" << std::setw(15) << ritzValues_[i].imagpart;
01837             }              
01838             os << std::setw(20) << ritzResiduals_[i] << endl;
01839           }
01840         }
01841       } else {
01842         os << std::setw(20) << "[ NONE COMPUTED ]" << endl;
01843       }
01844     }
01845     os << endl;
01846     os <<"================================================================================" << endl;
01847     os << endl;
01848   }
01849   
01850 } // End of namespace Anasazi
01851 
01852 #endif
01853 
01854 // End of file AnasaziBlockKrylovSchur.hpp

Generated on Tue Oct 20 12:45:18 2009 for Anasazi by doxygen 1.4.7