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

Generated on Thu Sep 18 12:31:35 2008 for Anasazi by doxygen 1.3.9.1