Anasazi Version of the Day
AnasaziGeneralizedDavidson.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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef ANASAZI_GENERALIZED_DAVIDSON_HPP
00030 #define ANASAZI_GENERALIZED_DAVIDSON_HPP
00031 
00038 #include "Teuchos_RCPDecl.hpp"
00039 #include "Teuchos_ParameterList.hpp"
00040 #include "Teuchos_SerialDenseMatrix.hpp"
00041 #include "Teuchos_SerialDenseVector.hpp"
00042 #include "Teuchos_Array.hpp"
00043 #include "Teuchos_BLAS.hpp"
00044 #include "Teuchos_LAPACK.hpp"
00045 
00046 #include "AnasaziConfigDefs.hpp"
00047 #include "AnasaziTypes.hpp"
00048 #include "AnasaziEigenproblem.hpp"
00049 #include "AnasaziEigensolver.hpp"
00050 #include "AnasaziOrthoManager.hpp"
00051 #include "AnasaziOutputManager.hpp"
00052 #include "AnasaziSortManager.hpp"
00053 #include "AnasaziStatusTest.hpp"
00054 
00055 using Teuchos::RCP;
00056 
00057 namespace Anasazi {
00058 
00062 template <class ScalarType, class MV>
00063 struct GeneralizedDavidsonState {
00065     int curDim;
00066 
00068     RCP<MV> V;
00069 
00071     RCP<MV> AV;
00072 
00074     RCP<MV> BV;
00075 
00077     RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > VAV;
00078 
00080     RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > VBV;
00081 
00083     RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > S;
00084 
00086     RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > T;
00087 
00089     RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > Q;
00090 
00092     RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > Z;
00093 
00095     std::vector< Value<ScalarType> > eVals;
00096 
00097     GeneralizedDavidsonState() : curDim(0), V(Teuchos::null), AV(Teuchos::null),
00098                                  BV(Teuchos::null), VAV(Teuchos::null),
00099                                  VBV(Teuchos::null), S(Teuchos::null),
00100                                  T(Teuchos::null), Q(Teuchos::null),
00101                                  Z(Teuchos::null), eVals(0) {}
00102 
00103 };
00104 
00105 
00126 template <class ScalarType, class MV, class OP>
00127 class GeneralizedDavidson : public Eigensolver<ScalarType,MV,OP>
00128 {
00129   private:
00130     // Convenience Typedefs
00131     typedef MultiVecTraits<ScalarType,MV>            MVT;
00132     typedef OperatorTraits<ScalarType,MV,OP>         OPT;
00133     typedef Teuchos::ScalarTraits<ScalarType>        ST;
00134     typedef typename ST::magnitudeType               MagnitudeType;
00135     typedef Teuchos::ScalarTraits<MagnitudeType>     MT;
00136 
00137   public:
00138 
00159     GeneralizedDavidson(const RCP<Eigenproblem<ScalarType,MV,OP> >  &problem,
00160                         const RCP<SortManager<MagnitudeType> >      &sortman,
00161                         const RCP<OutputManager<ScalarType> >       &outputman,
00162                         const RCP<StatusTest<ScalarType,MV,OP> >    &tester,
00163                         const RCP<OrthoManager<ScalarType,MV> >     &orthoman,
00164                         Teuchos::ParameterList                      &pl);
00165 
00169     void iterate();
00170 
00180     void initialize();
00181 
00185     void initialize( GeneralizedDavidsonState<ScalarType,MV> state );
00186 
00190     int getNumIters() const { return d_iteration; }
00191 
00195     void resetNumIters() { d_iteration=0; d_opApplies=0; }
00196 
00200     RCP<const MV> getRitzVectors()
00201     {
00202         if( !d_ritzVectorsValid )
00203             computeRitzVectors();
00204         return d_ritzVecs;
00205     }
00206 
00210     std::vector< Value<ScalarType> > getRitzValues();
00211 
00215     std::vector<int> getRitzIndex()
00216     {
00217         if( !d_ritzIndexValid )
00218             computeRitzIndex();
00219         return d_ritzIndex;
00220     }
00221 
00227     std::vector<int> getBlockIndex() const
00228     {
00229         return d_expansionIndices;
00230     }
00231 
00235     std::vector<MagnitudeType> getResNorms();
00236 
00240     std::vector<MagnitudeType> getResNorms(int numWanted);
00241 
00245     std::vector<MagnitudeType> getRes2Norms() { return d_resNorms; }
00246 
00253     std::vector<MagnitudeType> getRitzRes2Norms() { return d_resNorms; }
00254 
00258     int getCurSubspaceDim() const { return d_curDim; }
00259 
00263     int getMaxSubspaceDim() const { return d_maxSubspaceDim; }
00264 
00268     void setStatusTest( RCP<StatusTest<ScalarType,MV,OP> > tester) { d_tester = tester; }
00269 
00273     RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const { return d_tester; }
00274 
00278     const Eigenproblem<ScalarType,MV,OP> & getProblem() const { return *d_problem; }
00279 
00283     int getBlockSize() const { return d_expansionSize; }
00284 
00288     void setBlockSize(int blockSize);
00289 
00293     void setSize(int blockSize, int maxSubDim);
00294 
00298     Teuchos::Array< RCP<const MV> > getAuxVecs() const { return d_auxVecs; }
00299 
00307     void setAuxVecs( const Teuchos::Array< RCP<const MV> > &auxVecs );
00308 
00312     bool isInitialized() const { return d_initialized; }
00313 
00317     void currentStatus( std::ostream &myout );
00318 
00322     GeneralizedDavidsonState<ScalarType,MV> getState();
00323 
00327     void sortProblem( int numWanted );
00328 
00329   private:
00330 
00331     // Expand subspace
00332     void expandSearchSpace();
00333 
00334     // Apply Operators
00335     void applyOperators();
00336 
00337     // Update projections
00338     void updateProjections();
00339 
00340     // Solve projected eigenproblem
00341     void solveProjectedEigenproblem();
00342 
00343     // Compute eigenvectors of matrix pair
00344     void computeProjectedEigenvectors( Teuchos::SerialDenseMatrix<int,ScalarType> &X );
00345 
00346     // Scale projected eigenvectors by alpha/beta
00347     void scaleEigenvectors( const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
00348                                   Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
00349                                   Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta );
00350 
00351     // Sort vectors of pairs
00352     void sortValues( std::vector<MagnitudeType> &realParts,
00353                      std::vector<MagnitudeType> &imagParts,
00354                      std::vector<int>    &permVec,
00355                      int N);
00356 
00357     // Compute Residual
00358     void computeResidual();
00359 
00360     // Update the current Ritz index vector
00361     void computeRitzIndex();
00362 
00363     // Compute the current Ritz vectors
00364     void computeRitzVectors();
00365 
00366     // Operators
00367     RCP<Eigenproblem<ScalarType,MV,OP> > d_problem;
00368     Teuchos::ParameterList d_pl;
00369     RCP<const OP> d_A;
00370     RCP<const OP> d_B;
00371     RCP<const OP> d_P;
00372     bool d_haveB;
00373     bool d_haveP;
00374 
00375     // Parameters
00376     int d_blockSize;
00377     int d_maxSubspaceDim;
00378     int d_NEV;
00379     int d_numToPrint;
00380     std::string d_initType;
00381     int d_verbosity;
00382     bool d_relativeConvergence;
00383 
00384     // Managers
00385     RCP<OutputManager<ScalarType> >     d_outputMan;
00386     RCP<OrthoManager<ScalarType,MV> >   d_orthoMan;
00387     RCP<SortManager<MagnitudeType> >    d_sortMan;
00388     RCP<StatusTest<ScalarType,MV,OP> >  d_tester;
00389 
00390     // Eigenvalues
00391     std::vector< Value<ScalarType> > d_eigenvalues;
00392 
00393     // Residual Vector
00394     RCP<MV> d_R;
00395     std::vector<MagnitudeType> d_resNorms;
00396 
00397     // Subspace Vectors
00398     RCP<MV> d_V;
00399     RCP<MV> d_AV;
00400     RCP<MV> d_BV;
00401     RCP<MV> d_ritzVecSpace;
00402     RCP<MV> d_ritzVecs;
00403     Teuchos::Array< RCP<const MV> > d_auxVecs;
00404 
00405     // Serial Matrices
00406     RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VAV;
00407     RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VBV;
00408     RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_S;
00409     RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_T;
00410     RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Q;
00411     RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Z;
00412 
00413     // Arrays for holding Ritz values
00414     std::vector<MagnitudeType> d_alphar;
00415     std::vector<MagnitudeType> d_alphai;
00416     std::vector<MagnitudeType> d_betar;
00417     std::vector<int>    d_ritzIndex;
00418     std::vector<int>    d_convergedIndices;
00419     std::vector<int>    d_expansionIndices;
00420 
00421     // Current subspace dimension
00422     int d_curDim;
00423 
00424     // How many vectors are to be added to the subspace
00425     int d_expansionSize;
00426 
00427     // Should subspace expansion use leading vectors
00428     //  (if false, will use leading unconverged vectors)
00429     bool d_useLeading;
00430 
00431     // What should be used for test subspace (V, AV, or BV)
00432     std::string d_testSpace;
00433 
00434     // How many residual vectors are valid
00435     int d_residualSize;
00436 
00437     int  d_iteration;
00438     int  d_opApplies;
00439     bool d_initialized;
00440     bool d_ritzIndexValid;
00441     bool d_ritzVectorsValid;
00442 
00443 };
00444 
00445 //---------------------------------------------------------------------------//
00446 // Prevent instantiation on complex scalar type
00447 //---------------------------------------------------------------------------//
00448 template <class MagnitudeType, class MV, class OP>
00449 class GeneralizedDavidson<std::complex<MagnitudeType>,MV,OP>
00450 {
00451   public:
00452 
00453     typedef std::complex<MagnitudeType> ScalarType;
00454     GeneralizedDavidson(
00455         const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00456         const RCP<SortManager<MagnitudeType> >     &sortman,
00457         const RCP<OutputManager<ScalarType> >      &outputman,
00458         const RCP<StatusTest<ScalarType,MV,OP> >   &tester,
00459         const RCP<OrthoManager<ScalarType,MV> >    &orthoman,
00460         Teuchos::ParameterList                     &pl)
00461     {
00462         // Provide a compile error when attempting to instantiate on complex type
00463         MagnitudeType::this_class_is_missing_a_specialization();
00464     }
00465 };
00466 
00467 //---------------------------------------------------------------------------//
00468 // PUBLIC METHODS
00469 //---------------------------------------------------------------------------//
00470 
00471 //---------------------------------------------------------------------------//
00472 // Constructor
00473 //---------------------------------------------------------------------------//
00474 template <class ScalarType, class MV, class OP>
00475 GeneralizedDavidson<ScalarType,MV,OP>::GeneralizedDavidson(
00476         const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00477         const RCP<SortManager<MagnitudeType> >     &sortman,
00478         const RCP<OutputManager<ScalarType> >      &outputman,
00479         const RCP<StatusTest<ScalarType,MV,OP> >   &tester,
00480         const RCP<OrthoManager<ScalarType,MV> >    &orthoman,
00481         Teuchos::ParameterList                     &pl )
00482 {
00483     TEUCHOS_TEST_FOR_EXCEPTION(   problem == Teuchos::null, std::invalid_argument, "No Eigenproblem given to solver." );
00484     TEUCHOS_TEST_FOR_EXCEPTION( outputman == Teuchos::null, std::invalid_argument, "No OutputManager given to solver." );
00485     TEUCHOS_TEST_FOR_EXCEPTION(  orthoman == Teuchos::null, std::invalid_argument, "No OrthoManager given to solver." );
00486     TEUCHOS_TEST_FOR_EXCEPTION(   sortman == Teuchos::null, std::invalid_argument, "No SortManager given to solver." );
00487     TEUCHOS_TEST_FOR_EXCEPTION(    tester == Teuchos::null, std::invalid_argument, "No StatusTest given to solver." );
00488     TEUCHOS_TEST_FOR_EXCEPTION(   !problem->isProblemSet(), std::invalid_argument, "Problem has not been set." );
00489 
00490     d_problem = problem;
00491     d_pl = pl;
00492     TEUCHOS_TEST_FOR_EXCEPTION( problem->getA()==Teuchos::null && problem->getOperator()==Teuchos::null,
00493                                 std::invalid_argument, "Either A or Operator must be non-null on Eigenproblem");
00494     d_A = problem->getA()!=Teuchos::null ? problem->getA() : problem->getOperator();
00495     d_B = problem->getM();
00496     d_P = problem->getPrec();
00497     d_sortMan = sortman;
00498     d_outputMan = outputman;
00499     d_tester = tester;
00500     d_orthoMan = orthoman;
00501 
00502     // Pull entries from the ParameterList and Eigenproblem
00503     d_NEV        = d_problem->getNEV();
00504     d_initType   = d_pl.get<std::string>("Initial Guess","Random");
00505     d_numToPrint = d_pl.get<int>("Print Number of Ritz Values",-1);
00506     d_useLeading = d_pl.get<bool>("Use Leading Vectors",false);
00507 
00508     if( d_B != Teuchos::null )
00509         d_haveB = true;
00510     else
00511         d_haveB = false;
00512 
00513     if( d_P != Teuchos::null )
00514         d_haveP = true;
00515     else
00516         d_haveP = false;
00517 
00518     d_testSpace = d_pl.get<std::string>("Test Space","V");
00519     TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace!="V" && d_testSpace!="AV" && d_testSpace!="BV", std::invalid_argument,
00520         "Anasazi::GeneralizedDavidson: Test Space must be V, AV, or BV" );
00521     TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace=="V" ? false : !d_haveB, std::invalid_argument,
00522         "Anasazi::GeneralizedDavidson: Test Space must be V for standard eigenvalue problem" );
00523 
00524     // Allocate space for subspace vectors, projected matrices
00525     int blockSize  = d_pl.get<int>("Block Size",1);
00526     int maxSubDim  = d_pl.get<int>("Maximum Subspace Dimension",3*d_NEV*blockSize);
00527     d_blockSize      = -1;
00528     d_maxSubspaceDim = -1;
00529     setSize( blockSize, maxSubDim );
00530     d_relativeConvergence = d_pl.get<bool>("Relative Convergence Tolerance",false);
00531 
00532     // Make sure subspace size is consistent with requested eigenvalues
00533     TEUCHOS_TEST_FOR_EXCEPTION( d_blockSize <= 0, std::invalid_argument, "Block size must be positive");
00534     TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim <= 0, std::invalid_argument, "Maximum Subspace Dimension must be positive" );
00535     TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV()+2 > pl.get<int>("Maximum Subspace Dimension"),
00536                                 std::invalid_argument, "Maximum Subspace Dimension must be strictly greater than NEV");
00537     TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim > MVT::GetVecLength(*problem->getInitVec()),
00538                                 std::invalid_argument, "Maximum Subspace Dimension cannot exceed problem size");
00539 
00540 
00541     d_curDim = 0;
00542     d_iteration = 0;
00543     d_opApplies = 0;
00544     d_ritzIndexValid = false;
00545     d_ritzVectorsValid = false;
00546 }
00547 
00548 
00549 //---------------------------------------------------------------------------//
00550 // Iterate
00551 //---------------------------------------------------------------------------//
00552 template <class ScalarType, class MV, class OP>
00553 void GeneralizedDavidson<ScalarType,MV,OP>::iterate()
00554 {
00555     // Initialize Problem
00556     if( !d_initialized )
00557     {
00558         d_outputMan->stream(Warnings) << "WARNING: GeneralizedDavidson::iterate called without first calling initialize" << std::endl;
00559         d_outputMan->stream(Warnings) << "         Default initialization will be performed" << std::endl;
00560         initialize();
00561     }
00562 
00563     // Print current status
00564     if( d_outputMan->isVerbosity(Debug) )
00565     {
00566         currentStatus( d_outputMan->stream(Debug) );
00567     }
00568     else if( d_outputMan->isVerbosity(IterationDetails) )
00569     {
00570         currentStatus( d_outputMan->stream(IterationDetails) );
00571     }
00572 
00573     while( d_tester->getStatus() != Passed && d_curDim+d_expansionSize <= d_maxSubspaceDim )
00574     {
00575         d_iteration++;
00576 
00577         expandSearchSpace();
00578 
00579         applyOperators();
00580 
00581         updateProjections();
00582 
00583         solveProjectedEigenproblem();
00584 
00585         // Make sure the most significant Ritz values are in front
00586         // We want the greater of the block size and the number of
00587         //  requested values, but can't exceed the current dimension
00588         int numToSort = std::max(d_blockSize,d_NEV);
00589         numToSort = std::min(numToSort,d_curDim);
00590         sortProblem( numToSort );
00591 
00592         computeResidual();
00593 
00594         // Print current status
00595         if( d_outputMan->isVerbosity(Debug) )
00596         {
00597             currentStatus( d_outputMan->stream(Debug) );
00598         }
00599         else if( d_outputMan->isVerbosity(IterationDetails) )
00600         {
00601             currentStatus( d_outputMan->stream(IterationDetails) );
00602         }
00603     }
00604 }
00605 
00606 //---------------------------------------------------------------------------//
00607 // Return the current state struct
00608 //---------------------------------------------------------------------------//
00609 template <class ScalarType, class MV, class OP>
00610 GeneralizedDavidsonState<ScalarType,MV> GeneralizedDavidson<ScalarType,MV,OP>::getState()
00611 {
00612     GeneralizedDavidsonState<ScalarType,MV> state;
00613     state.curDim = d_curDim;
00614     state.V      = d_V;
00615     state.AV     = d_AV;
00616     state.BV     = d_BV;
00617     state.VAV    = d_VAV;
00618     state.VBV    = d_VBV;
00619     state.S      = d_S;
00620     state.T      = d_T;
00621     state.Q      = d_Q;
00622     state.Z      = d_Z;
00623     state.eVals  = getRitzValues();
00624     return state;
00625 }
00626 
00627 //---------------------------------------------------------------------------//
00628 // Set block size
00629 //---------------------------------------------------------------------------//
00630 template <class ScalarType, class MV, class OP>
00631 void GeneralizedDavidson<ScalarType,MV,OP>::setBlockSize(int blockSize)
00632 {
00633     setSize(blockSize,d_maxSubspaceDim);
00634 }
00635 
00636 //---------------------------------------------------------------------------//
00637 // Set block size and maximum subspace dimension.
00638 //---------------------------------------------------------------------------//
00639 template <class ScalarType, class MV, class OP>
00640 void GeneralizedDavidson<ScalarType,MV,OP>::setSize(int blockSize, int maxSubDim )
00641 {
00642     if( blockSize != d_blockSize || maxSubDim != d_maxSubspaceDim )
00643     {
00644         d_blockSize = blockSize;
00645         d_maxSubspaceDim = maxSubDim;
00646         d_initialized = false;
00647 
00648         d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Allocating eigenproblem"
00649             << " state with block size of " << d_blockSize
00650             << " and maximum subspace dimension of " << d_maxSubspaceDim << std::endl;
00651 
00652         // Resize arrays for Ritz values
00653         d_alphar.resize(d_maxSubspaceDim);
00654         d_alphai.resize(d_maxSubspaceDim);
00655         d_betar.resize(d_maxSubspaceDim);
00656 
00657         // Shorten for convenience here
00658         int msd = d_maxSubspaceDim;
00659 
00660         // Temporarily save initialization vector to clone needed vectors
00661         RCP<const MV> initVec = d_problem->getInitVec();
00662 
00663         // Allocate subspace vectors
00664         d_V            = MVT::Clone(*initVec, msd);
00665         d_AV           = MVT::Clone(*initVec, msd);
00666 
00667         // Allocate serial matrices
00668         d_VAV = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
00669         d_S   = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
00670         d_Q   = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
00671         d_Z   = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
00672 
00673         // If this is generalized eigenproblem, allocate B components
00674         if( d_haveB )
00675         {
00676             d_BV  = MVT::Clone(*initVec, msd);
00677             d_VBV = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
00678             d_T   = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
00679         }
00680 
00681         /* Allocate space for residual and Ritz vectors
00682          * The residual serves two purposes in the Davidson algorithm --
00683          *  subspace expansion (via the preconditioner) and convergence checking.
00684          * We need "Block Size" vectors for subspace expantion and NEV vectors
00685          *  for convergence checking.  Allocate space for max of these, one
00686          *  extra to avoid splitting conjugate pairs
00687          * Allocate one more than "Block Size" to avoid splitting a conjugate pair
00688          */
00689         d_R = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
00690         d_ritzVecSpace = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
00691     }
00692 }
00693 
00694 //---------------------------------------------------------------------------//
00695 /*
00696  * Initialize the eigenvalue problem
00697  *
00698  * Anything on the state that is not null is assumed to be valid.
00699  * Anything not present on the state will be generated
00700  * Very limited error checking can be performed to ensure the validity of
00701  * state components (e.g. we cannot verify that state.AV actually corresponds
00702  * to A*state.V), so this function should be used carefully.
00703  */
00704 //---------------------------------------------------------------------------//
00705 template <class ScalarType, class MV, class OP>
00706 void GeneralizedDavidson<ScalarType,MV,OP>::initialize( GeneralizedDavidsonState<ScalarType,MV> state )
00707 {
00708     // If state has nonzero dimension, we initialize from that, otherwise
00709     //  we'll pick d_blockSize vectors to start with
00710     d_curDim = (state.curDim > 0 ? state.curDim : d_blockSize );
00711 
00712     d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Initializing state"
00713         << " with subspace dimension " << d_curDim << std::endl;
00714 
00715     // Index for 1st block_size vectors
00716     std::vector<int> initInds(d_curDim);
00717     for( int i=0; i<d_curDim; ++i )
00718         initInds[i] = i;
00719 
00720     // View of vectors that need to be initialized
00721     RCP<MV>  V1 = MVT::CloneViewNonConst(*d_V,initInds);
00722 
00723     // If state's dimension is large enough, use state.V to initialize
00724     bool reset_V = false;;
00725     if( state.curDim > 0 && state.V != Teuchos::null && MVT::GetNumberVecs(*state.V) >= d_curDim )
00726     {
00727         if( state.V != d_V )
00728             MVT::SetBlock(*state.V,initInds,*V1);
00729     }
00730     // If there aren't enough vectors in problem->getInitVec() or the user specifically
00731     //  wants to use random data, set V to random
00732     else if( MVT::GetNumberVecs(*d_problem->getInitVec()) < d_blockSize || d_initType == "Random" )
00733     {
00734         MVT::MvRandom(*V1);
00735         reset_V = true;
00736     }
00737     // Use vectors in problem->getInitVec()
00738     else
00739     {
00740         RCP<const MV> initVec = MVT::CloneView(*d_problem->getInitVec(),initInds);
00741         MVT::SetBlock(*initVec,initInds,*V1);
00742         reset_V = true;
00743     }
00744 
00745     // If we reset V, it needs to be orthonormalized
00746     if( reset_V )
00747     {
00748         int rank = d_orthoMan->projectAndNormalize( *V1, d_auxVecs );
00749         TEUCHOS_TEST_FOR_EXCEPTION( rank < d_blockSize, std::logic_error,
00750             "Anasazi::GeneralizedDavidson::initialize(): Error generating initial orthonormal basis" );
00751     }
00752 
00753     if( d_outputMan->isVerbosity(Debug) )
00754     {
00755         d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Error in V^T V == I: "
00756             << d_orthoMan->orthonormError( *V1 ) << std::endl;
00757     }
00758 
00759     // Now process AV
00760     RCP<MV> AV1 = MVT::CloneViewNonConst(*d_AV,initInds);
00761 
00762     // If AV in the state is valid and of appropriate size, use it
00763     // We have no way to check that AV is actually A*V
00764     if( !reset_V && state.AV != Teuchos::null && MVT::GetNumberVecs(*state.AV) >= d_curDim )
00765     {
00766         if( state.AV != d_AV )
00767             MVT::SetBlock(*state.AV,initInds,*AV1);
00768     }
00769     // Otherwise apply A to V
00770     else
00771     {
00772         OPT::Apply( *d_A, *V1, *AV1 );
00773         d_opApplies += MVT::GetNumberVecs( *V1 );
00774     }
00775 
00776     // Views of matrix to be updated
00777     Teuchos::SerialDenseMatrix<int,ScalarType> VAV1( Teuchos::View, *d_VAV, d_curDim, d_curDim );
00778 
00779     // If the state has a valid VAV, use it
00780     if( !reset_V && state.VAV != Teuchos::null && state.VAV->numRows() >= d_curDim && state.VAV->numCols() >= d_curDim )
00781     {
00782         if( state.VAV != d_VAV )
00783         {
00784             Teuchos::SerialDenseMatrix<int,ScalarType> state_VAV( Teuchos::View, *state.VAV, d_curDim, d_curDim );
00785             VAV1.assign( state_VAV );
00786         }
00787     }
00788     // Otherwise compute VAV from V,AV
00789     else
00790     {
00791         if( d_testSpace == "V" )
00792         {
00793             MVT::MvTransMv( ST::one(),  *V1, *AV1, VAV1 );
00794         }
00795         else if( d_testSpace == "AV" )
00796         {
00797             MVT::MvTransMv( ST::one(), *AV1, *AV1, VAV1 );
00798         }
00799         else if( d_testSpace == "BV" )
00800         {
00801             RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
00802             MVT::MvTransMv( ST::one(), *BV1, *AV1, VAV1 );
00803         }
00804     }
00805 
00806     // Process BV if we have it
00807     if( d_haveB )
00808     {
00809         RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
00810 
00811         // If BV in the state is valid and of appropriate size, use it
00812         // We have no way to check that BV is actually B*V
00813         if( !reset_V && state.BV != Teuchos::null && MVT::GetNumberVecs(*state.BV) >= d_curDim )
00814         {
00815             if( state.BV != d_BV )
00816                 MVT::SetBlock(*state.BV,initInds,*BV1);
00817         }
00818         // Otherwise apply B to V
00819         else
00820         {
00821             OPT::Apply( *d_B, *V1, *BV1 );
00822         }
00823 
00824         // Views of matrix to be updated
00825         Teuchos::SerialDenseMatrix<int,ScalarType> VBV1( Teuchos::View, *d_VBV, d_curDim, d_curDim );
00826 
00827         // If the state has a valid VBV, use it
00828         if( !reset_V && state.VBV != Teuchos::null && state.VBV->numRows() >= d_curDim && state.VBV->numCols() >= d_curDim )
00829         {
00830             if( state.VBV != d_VBV )
00831             {
00832                 Teuchos::SerialDenseMatrix<int,ScalarType> state_VBV( Teuchos::View, *state.VBV, d_curDim, d_curDim );
00833                 VBV1.assign( state_VBV );
00834             }
00835         }
00836         // Otherwise compute VBV from V,BV
00837         else
00838         {
00839             if( d_testSpace == "V" )
00840             {
00841                 MVT::MvTransMv( ST::one(),  *V1, *BV1, VBV1 );
00842             }
00843             else if( d_testSpace == "AV" )
00844             {
00845                 MVT::MvTransMv( ST::one(), *AV1, *BV1, VBV1 );
00846             }
00847             else if( d_testSpace == "BV" )
00848             {
00849                 MVT::MvTransMv( ST::one(), *BV1, *BV1, VBV1 );
00850             }
00851         }
00852     }
00853 
00854     // Update Ritz values
00855     solveProjectedEigenproblem();
00856 
00857     // Sort
00858     int numToSort = std::max(d_blockSize,d_NEV);
00859     numToSort = std::min(numToSort,d_curDim);
00860     sortProblem( numToSort );
00861 
00862     // Get valid residual
00863     computeResidual();
00864 
00865     // Set solver to initialized
00866     d_initialized = true;
00867 }
00868 
00869 //---------------------------------------------------------------------------//
00870 // Initialize the eigenvalue problem with empty state
00871 //---------------------------------------------------------------------------//
00872 template <class ScalarType, class MV, class OP>
00873 void GeneralizedDavidson<ScalarType,MV,OP>::initialize()
00874 {
00875     GeneralizedDavidsonState<ScalarType,MV> empty;
00876     initialize( empty );
00877 }
00878 
00879 //---------------------------------------------------------------------------//
00880 // Get current residual norms
00881 //---------------------------------------------------------------------------//
00882 template <class ScalarType, class MV, class OP>
00883 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
00884     GeneralizedDavidson<ScalarType,MV,OP>::getResNorms()
00885 {
00886     return getResNorms(d_residualSize);
00887 }
00888 
00889 //---------------------------------------------------------------------------//
00890 // Get current residual norms
00891 //---------------------------------------------------------------------------//
00892 template <class ScalarType, class MV, class OP>
00893 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
00894     GeneralizedDavidson<ScalarType,MV,OP>::getResNorms(int numWanted)
00895 {
00896     std::vector<int> resIndices(numWanted);
00897     for( int i=0; i<numWanted; ++i )
00898         resIndices[i]=i;
00899 
00900     RCP<const MV> R_checked = MVT::CloneView( *d_R, resIndices );
00901 
00902     std::vector<MagnitudeType> resNorms;
00903     d_orthoMan->norm( *R_checked, resNorms );
00904 
00905     return resNorms;
00906 }
00907 
00908 //---------------------------------------------------------------------------//
00909 // Get current Ritz values
00910 //---------------------------------------------------------------------------//
00911 template <class ScalarType, class MV, class OP>
00912 std::vector< Value<ScalarType> > GeneralizedDavidson<ScalarType,MV,OP>::getRitzValues()
00913 {
00914     std::vector< Value<ScalarType> > ritzValues;
00915     for( int ival=0; ival<d_curDim; ++ival )
00916     {
00917         Value<ScalarType> thisVal;
00918         thisVal.realpart = d_alphar[ival] / d_betar[ival];
00919         if( d_betar[ival] != MT::zero() )
00920             thisVal.imagpart = d_alphai[ival] / d_betar[ival];
00921         else
00922             thisVal.imagpart = MT::zero();
00923 
00924         ritzValues.push_back( thisVal );
00925     }
00926 
00927     return ritzValues;
00928 }
00929 
00930 //---------------------------------------------------------------------------//
00931 /*
00932  * Set auxilliary vectors
00933  *
00934  * Manually setting the auxilliary vectors invalidates the current state
00935  * of the solver.  Reuse of any components of the solver requires extracting
00936  * the state, orthogonalizing V against the aux vecs and reinitializing.
00937  */
00938 //---------------------------------------------------------------------------//
00939 template <class ScalarType, class MV, class OP>
00940 void GeneralizedDavidson<ScalarType,MV,OP>::setAuxVecs(
00941         const Teuchos::Array< RCP<const MV> > &auxVecs )
00942 {
00943     d_auxVecs = auxVecs;
00944 
00945     // Set state to uninitialized if any vectors were set here
00946     typename Teuchos::Array< RCP<const MV> >::const_iterator arrItr;
00947     int numAuxVecs=0;
00948     for( arrItr=auxVecs.begin(); arrItr!=auxVecs.end(); ++arrItr )
00949     {
00950         numAuxVecs += MVT::GetNumberVecs( *(*arrItr) );
00951     }
00952     if( numAuxVecs > 0 )
00953         d_initialized = false;
00954 }
00955 
00956 //---------------------------------------------------------------------------//
00957 // Reorder Schur form, bringing wanted values to front
00958 //---------------------------------------------------------------------------//
00959 template <class ScalarType, class MV, class OP>
00960 void GeneralizedDavidson<ScalarType,MV,OP>::sortProblem( int numWanted )
00961 {
00962     // Get permutation vector
00963     std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
00964     std::vector< Value<ScalarType> > ritzVals = getRitzValues();
00965     for( int i=0; i<d_curDim; ++i )
00966     {
00967         realRitz[i] = ritzVals[i].realpart;
00968         imagRitz[i] = ritzVals[i].imagpart;
00969     }
00970 
00971     std::vector<int> permVec;
00972     sortValues( realRitz, imagRitz, permVec, d_curDim );
00973 
00974     std::vector<int> sel(d_curDim,0);
00975     for( int ii=0; ii<numWanted; ++ii )
00976         sel[ permVec[ii] ]=1;
00977 
00978     if( d_haveB )
00979     {
00980         int ijob  = 0; // reorder only, no condition number estimates
00981         int wantq = 1; // keep left Schur vectors
00982         int wantz = 1; // keep right Schur vectors
00983         int work_size=10*d_maxSubspaceDim+16;
00984         std::vector<ScalarType> work(work_size);
00985         int sdim   = 0;
00986         int iwork_size = 1;
00987         int iwork;
00988         int info   = 0;
00989 
00990         Teuchos::LAPACK<int,ScalarType> lapack;
00991         lapack.TGSEN( ijob, wantq, wantz, &sel[0], d_curDim, d_S->values(), d_S->stride(), d_T->values(), d_T->stride(),
00992                       &d_alphar[0], &d_alphai[0], &d_betar[0], d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(),
00993                       &sdim, 0, 0, 0, &work[0], work_size, &iwork, iwork_size, &info );
00994 
00995         d_ritzIndexValid   = false;
00996         d_ritzVectorsValid = false;
00997 
00998         std::stringstream ss;
00999         ss << "Anasazi::GeneralizedDavidson: TGSEN returned error code " << info << std::endl;
01000         TEUCHOS_TEST_FOR_EXCEPTION( info<0, std::runtime_error, ss.str() );
01001         if( info > 0 )
01002         {
01003             // Only issue a warning for positive error code, this usually indicates
01004             //  that the system has not been fully reordered, presumably due to ill-conditioning.
01005             // This is usually not detrimental to the calculation.
01006             d_outputMan->stream(Warnings) << "WARNING: " << ss.str() << std::endl;
01007             d_outputMan->stream(Warnings) << "  Problem may not be correctly sorted" << std::endl;
01008         }
01009     }
01010     else {
01011       char getCondNum = 'N'; // no condition number estimates
01012       char getQ = 'V';       // keep Schur vectors
01013       int subDim = 0;
01014       int work_size = d_curDim;
01015       std::vector<ScalarType> work(work_size);
01016       int iwork_size = 1;
01017       int iwork = 0;
01018       int info = 0;
01019 
01020       Teuchos::LAPACK<int,ScalarType> lapack;
01021       lapack.TRSEN (getCondNum, getQ, &sel[0], d_curDim, d_S->values (),
01022                     d_S->stride (), d_Z->values (), d_Z->stride (),
01023                     &d_alphar[0], &d_alphai[0], &subDim, 0, 0, &work[0],
01024                     work_size, &iwork, iwork_size, &info );
01025 
01026       std::fill( d_betar.begin(), d_betar.end(), ST::one() );
01027 
01028       d_ritzIndexValid = false;
01029       d_ritzVectorsValid = false;
01030 
01031       TEUCHOS_TEST_FOR_EXCEPTION(
01032         info < 0, std::runtime_error, "Anasazi::GeneralizedDavidson::"
01033         "sortProblem: LAPACK routine TRSEN returned error code INFO = "
01034         << info << " < 0.  This indicates that argument " << -info
01035         << " (counting starts with one) of TRSEN had an illegal value.");
01036 
01037       // LAPACK's documentation suggests that this should only happen
01038       // in the real-arithmetic case, because I only see INFO == 1
01039       // possible for DTRSEN, not for ZTRSEN.  Nevertheless, it's
01040       // harmless to check regardless.
01041       TEUCHOS_TEST_FOR_EXCEPTION(
01042         info == 1, std::runtime_error, "Anasazi::GeneralizedDavidson::"
01043         "sortProblem: LAPACK routine TRSEN returned error code INFO = 1.  "
01044         "This indicates that the reordering failed because some eigenvalues "
01045         "are too close to separate (the problem is very ill-conditioned).");
01046 
01047       TEUCHOS_TEST_FOR_EXCEPTION(
01048         info > 1, std::logic_error, "Anasazi::GeneralizedDavidson::"
01049         "sortProblem: LAPACK routine TRSEN returned error code INFO = "
01050         << info << " > 1.  This should not be possible.  It may indicate an "
01051         "error either in LAPACK itself, or in Teuchos' LAPACK wrapper.");
01052     }
01053 }
01054 
01055 
01056 //---------------------------------------------------------------------------//
01057 // PRIVATE IMPLEMENTATION
01058 //---------------------------------------------------------------------------//
01059 
01060 //---------------------------------------------------------------------------//
01061 // Expand subspace using preconditioner and orthogonalize
01062 //---------------------------------------------------------------------------//
01063 template <class ScalarType, class MV, class OP>
01064 void GeneralizedDavidson<ScalarType,MV,OP>::expandSearchSpace()
01065 {
01066     // Get indices into relevant portion of residual and
01067     //  location to be added to search space
01068     std::vector<int> newIndices(d_expansionSize);
01069     for( int i=0; i<d_expansionSize; ++i )
01070     {
01071         newIndices[i] = d_curDim+i;
01072     }
01073 
01074     // Get indices into pre-existing search space
01075     std::vector<int> curIndices(d_curDim);
01076     for( int i=0; i<d_curDim; ++i )
01077         curIndices[i] = i;
01078 
01079     // Get View of vectors
01080     RCP<MV>       V_new    = MVT::CloneViewNonConst( *d_V, newIndices);
01081     RCP<const MV> V_cur    = MVT::CloneView(         *d_V, curIndices);
01082     RCP<const MV> R_active = MVT::CloneView(         *d_R, d_expansionIndices);
01083 
01084     if( d_haveP )
01085     {
01086         // Apply Preconditioner to Residual
01087         OPT::Apply( *d_P, *R_active, *V_new );
01088     }
01089     else
01090     {
01091         // Just copy the residual
01092         MVT::SetBlock( *R_active, newIndices, *d_V );
01093     }
01094 
01095     // Normalize new vector against existing vectors in V plus auxVecs
01096     Teuchos::Array< RCP<const MV> > against = d_auxVecs;
01097     against.push_back( V_cur );
01098     int rank = d_orthoMan->projectAndNormalize(*V_new,against);
01099 
01100     if( d_outputMan->isVerbosity(Debug) )
01101     {
01102         std::vector<int> allIndices(d_curDim+d_expansionSize);
01103         for( int i=0; i<d_curDim+d_expansionSize; ++i )
01104             allIndices[i]=i;
01105 
01106         RCP<const MV> V_all = MVT::CloneView( *d_V, allIndices );
01107 
01108         d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Error in V^T V == I: "
01109             << d_orthoMan->orthonormError( *V_all ) << std::endl;
01110     }
01111 
01112     TEUCHOS_TEST_FOR_EXCEPTION( rank != d_expansionSize, std::runtime_error,
01113         "Anasazi::GeneralizedDavidson::ExpandSearchSpace(): Orthonormalization of new vectors failed" );
01114 
01115 }
01116 
01117 //---------------------------------------------------------------------------//
01118 // Apply operators
01119 //---------------------------------------------------------------------------//
01120 template <class ScalarType, class MV, class OP>
01121 void GeneralizedDavidson<ScalarType,MV,OP>::applyOperators()
01122 {
01123     // Get indices for different components
01124     std::vector<int> newIndices(d_expansionSize);
01125     for( int i=0; i<d_expansionSize; ++i )
01126         newIndices[i] = d_curDim+i;
01127 
01128     // Get Views
01129     RCP<const MV>  V_new = MVT::CloneView(         *d_V,  newIndices );
01130     RCP<MV>       AV_new = MVT::CloneViewNonConst( *d_AV, newIndices );
01131 
01132     // Multiply by A
01133     OPT::Apply( *d_A, *V_new, *AV_new );
01134     d_opApplies += MVT::GetNumberVecs( *V_new );
01135 
01136     // Multiply by B
01137     if( d_haveB )
01138     {
01139         RCP<MV>       BV_new = MVT::CloneViewNonConst( *d_BV, newIndices );
01140         OPT::Apply( *d_B, *V_new, *BV_new );
01141     }
01142 }
01143 
01144 //---------------------------------------------------------------------------//
01145 // Update projected matrices.
01146 //---------------------------------------------------------------------------//
01147 template <class ScalarType, class MV, class OP>
01148 void GeneralizedDavidson<ScalarType,MV,OP>::updateProjections()
01149 {
01150     // Get indices for different components
01151     std::vector<int> newIndices(d_expansionSize);
01152     for( int i=0; i<d_expansionSize; ++i )
01153         newIndices[i] = d_curDim+i;
01154 
01155     std::vector<int> curIndices(d_curDim);
01156     for( int i=0; i<d_curDim; ++i )
01157         curIndices[i] = i;
01158 
01159     std::vector<int> allIndices(d_curDim+d_expansionSize);
01160     for( int i=0; i<d_curDim+d_expansionSize; ++i )
01161         allIndices[i] = i;
01162 
01163     // Test subspace can be V, AV, or BV
01164     RCP<const MV> W_new, W_all;
01165     if( d_testSpace == "V" )
01166     {
01167         W_new = MVT::CloneView(*d_V, newIndices );
01168         W_all = MVT::CloneView(*d_V, allIndices );
01169     }
01170     else if( d_testSpace == "AV" )
01171     {
01172         W_new = MVT::CloneView(*d_AV, newIndices );
01173         W_all = MVT::CloneView(*d_AV, allIndices );
01174     }
01175     else if( d_testSpace == "BV" )
01176     {
01177         W_new = MVT::CloneView(*d_BV, newIndices );
01178         W_all = MVT::CloneView(*d_BV, allIndices );
01179     }
01180 
01181     // Get views of AV
01182     RCP<const MV>     AV_new = MVT::CloneView(*d_AV, newIndices);
01183     RCP<const MV> AV_current = MVT::CloneView(*d_AV, curIndices);
01184 
01185     // Last block_size rows of VAV (minus final entry)
01186     Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastrow( Teuchos::View, *d_VAV, d_expansionSize, d_curDim, d_curDim, 0 );
01187     MVT::MvTransMv( ST::one(), *W_new, *AV_current, VAV_lastrow );
01188 
01189     // Last block_size columns of VAV
01190     Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastcol( Teuchos::View, *d_VAV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
01191     MVT::MvTransMv( ST::one(), *W_all, *AV_new, VAV_lastcol );
01192 
01193     if( d_haveB )
01194     {
01195         // Get views of BV
01196         RCP<const MV>     BV_new = MVT::CloneView(*d_BV, newIndices);
01197         RCP<const MV> BV_current = MVT::CloneView(*d_BV, curIndices);
01198 
01199         // Last block_size rows of VBV (minus final entry)
01200         Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastrow( Teuchos::View, *d_VBV, d_expansionSize, d_curDim, d_curDim, 0 );
01201         MVT::MvTransMv( ST::one(), *W_new, *BV_current, VBV_lastrow );
01202 
01203         // Last block_size columns of VBV
01204         Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastcol( Teuchos::View, *d_VBV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
01205         MVT::MvTransMv( ST::one(), *W_all, *BV_new, VBV_lastcol );
01206     }
01207 
01208     // All bases are expanded, increase current subspace dimension
01209     d_curDim += d_expansionSize;
01210 
01211     d_ritzIndexValid   = false;
01212     d_ritzVectorsValid = false;
01213 }
01214 
01215 //---------------------------------------------------------------------------//
01216 // Solve low dimensional eigenproblem using LAPACK
01217 //---------------------------------------------------------------------------//
01218 template <class ScalarType, class MV, class OP>
01219 void GeneralizedDavidson<ScalarType,MV,OP>::solveProjectedEigenproblem()
01220 {
01221     if( d_haveB )
01222     {
01223         // VAV and VBV need to stay unchanged, GGES will overwrite
01224         //  S and T with the triangular matrices from the generalized
01225         //  Schur form
01226         d_S->assign(*d_VAV);
01227         d_T->assign(*d_VBV);
01228 
01229         // Get QZ Decomposition of Projected Problem
01230         char leftVecs  = 'V'; // compute left vectors
01231         char rightVecs = 'V'; // compute right vectors
01232         char sortVals  = 'N'; // don't sort
01233         int sdim;
01234         // int work_size = 10*d_curDim+16;
01235         Teuchos::LAPACK<int,ScalarType> lapack;
01236         int info;
01237         // workspace query
01238         int work_size = -1;
01239         std::vector<ScalarType> work(1);
01240         lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
01241                      d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
01242                      d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
01243         // actual call
01244         work_size = work[0];
01245         work.resize(work_size);
01246         lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
01247                      d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
01248                      d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
01249 
01250         d_ritzIndexValid   = false;
01251         d_ritzVectorsValid = false;
01252 
01253         std::stringstream ss;
01254         ss << "Anasazi::GeneralizedDavidson: GGES returned error code " << info << std::endl;
01255         TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
01256     }
01257     else
01258     {
01259         // VAV needs to stay unchanged, GGES will overwrite
01260         //  S with the triangular matrix from the Schur form
01261         d_S->assign(*d_VAV);
01262 
01263         // Get QR Decomposition of Projected Problem
01264         char vecs = 'V'; // compute Schur vectors
01265         int sdim;
01266         int work_size = 3*d_curDim;
01267         std::vector<ScalarType>  work(work_size);
01268         int info;
01269 
01270         Teuchos::LAPACK<int,ScalarType> lapack;
01271         lapack.GEES( vecs, d_curDim, d_S->values(), d_S->stride(), &sdim, &d_alphar[0], &d_alphai[0],
01272                      d_Z->values(), d_Z->stride(), &work[0], work_size, 0, 0, &info);
01273 
01274         std::fill( d_betar.begin(), d_betar.end(), ST::one() );
01275 
01276         d_ritzIndexValid   = false;
01277         d_ritzVectorsValid = false;
01278 
01279         std::stringstream ss;
01280         ss << "Anasazi::GeneralizedDavidson: GEES returned error code " << info << std::endl;
01281         TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
01282     }
01283 }
01284 
01285 //---------------------------------------------------------------------------//
01286 /*
01287  * Get index vector into current Ritz values/vectors
01288  *
01289  * The current ordering of d_alphar, d_alphai, d_betar will be used.
01290  * Reordering those vectors will invalidate the vector returned here.
01291  */
01292 //---------------------------------------------------------------------------//
01293 template <class ScalarType, class MV, class OP>
01294 void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzIndex()
01295 {
01296     if( d_ritzIndexValid )
01297         return;
01298 
01299     d_ritzIndex.resize( d_curDim );
01300     int i=0;
01301     while( i < d_curDim )
01302     {
01303         if( d_alphai[i] == ST::zero() )
01304         {
01305             d_ritzIndex[i] = 0;
01306             i++;
01307         }
01308         else
01309         {
01310             d_ritzIndex[i]   =  1;
01311             d_ritzIndex[i+1] = -1;
01312             i+=2;
01313         }
01314     }
01315     d_ritzIndexValid = true;
01316 }
01317 
01318 //---------------------------------------------------------------------------//
01319 /*
01320  * Compute current Ritz vectors
01321  *
01322  * The current ordering of d_alphar, d_alphai, d_betar will be used.
01323  * Reordering those vectors will invalidate the vector returned here.
01324  */
01325 //---------------------------------------------------------------------------//
01326 template <class ScalarType, class MV, class OP>
01327 void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzVectors()
01328 {
01329     if( d_ritzVectorsValid )
01330         return;
01331 
01332     // Make Ritz indices current
01333     computeRitzIndex();
01334 
01335     // Get indices of converged vector
01336     std::vector<int> checkedIndices(d_residualSize);
01337     for( int ii=0; ii<d_residualSize; ++ii )
01338         checkedIndices[ii] = ii;
01339 
01340     // Get eigenvectors of projected system
01341     Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
01342     computeProjectedEigenvectors( X );
01343 
01344     // Get view of wanted vectors
01345     Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View,X,d_curDim,d_residualSize);
01346 
01347     // Get views of relevant portion of V, evecs
01348     d_ritzVecs = MVT::CloneViewNonConst( *d_ritzVecSpace, checkedIndices );
01349 
01350     std::vector<int> curIndices(d_curDim);
01351     for( int i=0; i<d_curDim; ++i )
01352         curIndices[i] = i;
01353 
01354     RCP<const MV> V_current = MVT::CloneView( *d_V, curIndices );
01355 
01356     // Now form Ritz vector
01357     MVT::MvTimesMatAddMv(ST::one(),*V_current,X_wanted,ST::zero(),*d_ritzVecs);
01358 
01359     // Normalize vectors, conjugate pairs get normalized together
01360     std::vector<MagnitudeType> scale(d_residualSize);
01361     MVT::MvNorm( *d_ritzVecs, scale );
01362     Teuchos::LAPACK<int,ScalarType> lapack;
01363     for( int i=0; i<d_residualSize; ++i )
01364     {
01365         if( d_ritzIndex[i] == 0 )
01366         {
01367             scale[i] = 1.0/scale[i];
01368         }
01369         else if( d_ritzIndex[i] == 1 )
01370         {
01371             MagnitudeType nrm = lapack.LAPY2(scale[i],scale[i+1]);
01372             scale[i]   = 1.0/nrm;
01373             scale[i+1] = 1.0/nrm;
01374         }
01375     }
01376     MVT::MvScale( *d_ritzVecs, scale );
01377 
01378     d_ritzVectorsValid = true;
01379 
01380 }
01381 
01382 //---------------------------------------------------------------------------//
01383 // Use sort manager to sort generalized eigenvalues
01384 //---------------------------------------------------------------------------//
01385 template <class ScalarType, class MV, class OP>
01386 void GeneralizedDavidson<ScalarType,MV,OP>::sortValues( std::vector<MagnitudeType> &realParts,
01387                                              std::vector<MagnitudeType> &imagParts,
01388                                              std::vector<int>    &permVec,
01389                                              int N)
01390 {
01391     permVec.resize(N);
01392 
01393     TEUCHOS_TEST_FOR_EXCEPTION( (int) realParts.size()<N, std::runtime_error,
01394         "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
01395     TEUCHOS_TEST_FOR_EXCEPTION( (int) imagParts.size()<N, std::runtime_error,
01396         "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
01397 
01398     RCP< std::vector<int> > rcpPermVec = Teuchos::rcpFromRef(permVec);
01399 
01400     d_sortMan->sort( realParts, imagParts, rcpPermVec, N );
01401 
01402     d_ritzIndexValid = false;
01403     d_ritzVectorsValid = false;
01404 }
01405 
01406 //---------------------------------------------------------------------------//
01407 /*
01408  * Compute (right) scaled eigenvectors of a pair of dense matrices
01409  *
01410  * This routine computes the eigenvectors for the generalized eigenvalue
01411  * problem \f$ \beta A x = \alpha B x \f$.  The input matrices are the upper
01412  * quasi-triangular matrices S and T from a real QZ decomposition,
01413  * the routine dtgevc will back-calculate the eigenvectors of the original
01414  * pencil (A,B) using the orthogonal matrices Q and Z.
01415  */
01416 //---------------------------------------------------------------------------//
01417 template <class ScalarType, class MV, class OP>
01418 void GeneralizedDavidson<ScalarType,MV,OP>::computeProjectedEigenvectors(
01419         Teuchos::SerialDenseMatrix<int,ScalarType> &X )
01420 {
01421     int N = X.numRows();
01422     if( d_haveB )
01423     {
01424         Teuchos::SerialDenseMatrix<int,ScalarType>  S(Teuchos::Copy, *d_S, N, N);
01425         Teuchos::SerialDenseMatrix<int,ScalarType>  T(Teuchos::Copy, *d_T, N, N);
01426         Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Q, N, N);
01427 
01428         char whichVecs = 'R'; // only need right eigenvectors
01429         char howMany   = 'B'; // back-compute eigenvectors of original A,B (we have Schur decomposition here)
01430         int work_size = 6*d_maxSubspaceDim;
01431         std::vector<ScalarType> work(work_size,ST::zero());
01432         int info;
01433         int M;
01434 
01435         Teuchos::LAPACK<int,ScalarType> lapack;
01436         lapack.TGEVC( whichVecs, howMany, 0, N, S.values(), S.stride(), T.values(), T.stride(),
01437                       VL.values(), VL.stride(), X.values(), X.stride(), N, &M, &work[0], &info );
01438 
01439         std::stringstream ss;
01440         ss << "Anasazi::GeneralizedDavidson: TGEVC returned error code " << info << std::endl;
01441         TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
01442     }
01443     else
01444     {
01445         Teuchos::SerialDenseMatrix<int,ScalarType>  S(Teuchos::Copy, *d_S, N, N);
01446         Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Z, N, N);
01447 
01448         char whichVecs = 'R'; // only need right eigenvectors
01449         char howMany   = 'B'; // back-compute eigenvectors of original A (we have Schur decomposition here)
01450         int sel = 0;
01451         std::vector<ScalarType> work(3*N);
01452         int m;
01453         int info;
01454 
01455         Teuchos::LAPACK<int,ScalarType> lapack;
01456 
01457         lapack.TREVC( whichVecs, howMany, &sel, N, S.values(), S.stride(), VL.values(), VL.stride(),
01458                       X.values(), X.stride(), N, &m, &work[0], &info );
01459 
01460         std::stringstream ss;
01461         ss << "Anasazi::GeneralizedDavidson: TREVC returned error code " << info << std::endl;
01462         TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
01463     }
01464 }
01465 
01466 //---------------------------------------------------------------------------//
01467 // Scale eigenvectors by quasi-diagonal matrices alpha and beta
01468 //---------------------------------------------------------------------------//
01469 template <class ScalarType, class MV, class OP>
01470 void GeneralizedDavidson<ScalarType,MV,OP>::scaleEigenvectors(
01471         const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
01472               Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
01473               Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta )
01474 {
01475     int Nr = X.numRows();
01476     int Nc = X.numCols();
01477 
01478     TEUCHOS_TEST_FOR_EXCEPTION( Nr>d_curDim, std::logic_error,
01479         "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
01480     TEUCHOS_TEST_FOR_EXCEPTION( Nc>d_curDim, std::logic_error,
01481         "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
01482     TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numRows()!=Nr, std::logic_error,
01483         "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xalpha does not match X");
01484     TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numCols()!=Nc, std::logic_error,
01485         "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xalpha does not match X");
01486     TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numRows()!=Nr, std::logic_error,
01487         "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xbeta does not match X");
01488     TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numCols()!=Nc, std::logic_error,
01489         "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xbeta does not match X");
01490 
01491     // Now form quasi-diagonal matrices
01492     //  containing alpha and beta
01493     Teuchos::SerialDenseMatrix<int,ScalarType> Alpha(Nc,Nc,true);
01494     Teuchos::SerialDenseMatrix<int,ScalarType> Beta(Nc,Nc,true);
01495 
01496     computeRitzIndex();
01497 
01498     for( int i=0; i<Nc; ++i )
01499     {
01500         if( d_ritzIndex[i] == 0 )
01501         {
01502             Alpha(i,i) = d_alphar[i];
01503             Beta(i,i)  = d_betar[i];
01504         }
01505         else if( d_ritzIndex[i] == 1 )
01506         {
01507             Alpha(i,i)   = d_alphar[i];
01508             Alpha(i,i+1) = d_alphai[i];
01509             Beta(i,i)    = d_betar[i];
01510         }
01511         else
01512         {
01513             Alpha(i,i-1) = d_alphai[i];
01514             Alpha(i,i)   = d_alphar[i];
01515             Beta(i,i)    = d_betar[i];
01516         }
01517     }
01518 
01519     int err;
01520 
01521     // Multiply the eigenvectors by alpha
01522     err = X_alpha.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Alpha, ST::zero() );
01523     std::stringstream astream;
01524     astream << "GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
01525     TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, astream.str() );
01526 
01527     // Multiply the eigenvectors by beta
01528     err = X_beta.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Beta, ST::zero() );
01529     std::stringstream bstream;
01530     bstream << "GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
01531     TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, bstream.str() );
01532 }
01533 
01534 //---------------------------------------------------------------------------//
01535 // Compute residual
01536 //---------------------------------------------------------------------------//
01537 template <class ScalarType, class MV, class OP>
01538 void GeneralizedDavidson<ScalarType,MV,OP>::computeResidual()
01539 {
01540     computeRitzIndex();
01541 
01542     // Determine how many residual vectors need to be computed
01543     d_residualSize = std::max( d_blockSize, d_NEV );
01544     if( d_curDim < d_residualSize )
01545     {
01546         d_residualSize = d_curDim;
01547     }
01548     else if( d_ritzIndex[d_residualSize-1] == 1 )
01549     {
01550         d_residualSize++;
01551     }
01552 
01553     // Get indices of all valid residual vectors
01554     std::vector<int> residualIndices(d_residualSize);
01555     for( int i=0; i<d_residualSize; ++i )
01556         residualIndices[i] = i;
01557 
01558     // X will store (right) eigenvectors of projected system
01559     Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
01560 
01561     // Get eigenvectors of projected problem -- computed from previous Schur decomposition
01562     computeProjectedEigenvectors( X );
01563 
01564     // X_alpha and X_beta will be eigenvectors right-multiplied by alpha, beta (which are quasi-diagonal portions of S,T)
01565     Teuchos::SerialDenseMatrix<int,ScalarType> X_alpha(d_curDim,d_residualSize);
01566     Teuchos::SerialDenseMatrix<int,ScalarType>  X_beta(d_curDim,d_residualSize);
01567 
01568     // X_wanted is the wanted portion of X
01569     Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View, X, d_curDim, d_residualSize);
01570 
01571     // Scale Eigenvectors by alpha or beta
01572     scaleEigenvectors( X_wanted, X_alpha, X_beta );
01573 
01574     // Get view of residual vector(s)
01575     RCP<MV> R_active = MVT::CloneViewNonConst( *d_R, residualIndices );
01576 
01577     // View of active portion of AV,BV
01578     std::vector<int> activeIndices(d_curDim);
01579     for( int i=0; i<d_curDim; ++i )
01580         activeIndices[i]=i;
01581 
01582     // Compute residual
01583     RCP<const MV> AV_active = MVT::CloneView( *d_AV, activeIndices );
01584     MVT::MvTimesMatAddMv(ST::one(),*AV_active, X_beta,  ST::zero(),*R_active);
01585 
01586     if( d_haveB )
01587     {
01588         RCP<const MV> BV_active = MVT::CloneView( *d_BV, activeIndices );
01589         MVT::MvTimesMatAddMv(ST::one(),*BV_active, X_alpha,-ST::one(), *R_active);
01590     }
01591     else
01592     {
01593         RCP<const MV> V_active = MVT::CloneView( *d_V, activeIndices );
01594         MVT::MvTimesMatAddMv(ST::one(),*V_active, X_alpha,-ST::one(), *R_active);
01595     }
01596 
01597     /* Apply a scaling to the residual
01598      * For generalized eigenvalue problems, LAPACK scales eigenvectors
01599      *  to have unit length in the infinity norm, we want them to have unit
01600      *  length in the 2-norm.  For conjugate pairs, the scaling is such that
01601      *  |xr|^2 + |xi|^2 = 1
01602      * Additionally, the residual is currently computed as r=beta*A*x-alpha*B*x
01603      *  but the "standard" residual is r=A*x-(alpha/beta)*B*x, or if we want
01604      *  to scale the residual by the Ritz value then it is r=(beta/alpha)*A*x-B*x
01605      *  Performing the scaling this way allows us to avoid the possibility of
01606      *  diving by infinity or zero if the StatusTest were allowed to handle the
01607      *  scaling.
01608      */
01609     Teuchos::LAPACK<int,ScalarType> lapack;
01610     Teuchos::BLAS<int,ScalarType> blas;
01611     std::vector<MagnitudeType> resScaling(d_residualSize);
01612     for( int icol=0; icol<d_residualSize; ++icol )
01613     {
01614         if( d_ritzIndex[icol] == 0 )
01615         {
01616             MagnitudeType Xnrm = blas.NRM2( d_curDim, X_wanted[icol], 1);
01617             MagnitudeType ABscaling = d_relativeConvergence ? d_alphar[icol] : d_betar[icol];
01618             resScaling[icol] = MT::one() / (Xnrm * ABscaling);
01619         }
01620         else if( d_ritzIndex[icol] == 1 )
01621         {
01622             MagnitudeType Xnrm1 = blas.NRM2( d_curDim, X_wanted[icol],   1 );
01623             MagnitudeType Xnrm2 = blas.NRM2( d_curDim, X_wanted[icol+1], 1 );
01624             MagnitudeType Xnrm  = lapack.LAPY2(Xnrm1,Xnrm2);
01625             MagnitudeType ABscaling = d_relativeConvergence ? lapack.LAPY2(d_alphar[icol],d_alphai[icol])
01626                                                             : d_betar[icol];
01627             resScaling[icol]   = MT::one() / (Xnrm * ABscaling);
01628             resScaling[icol+1] = MT::one() / (Xnrm * ABscaling);
01629         }
01630     }
01631     MVT::MvScale( *R_active, resScaling );
01632 
01633     // Compute residual norms
01634     d_resNorms.resize(d_residualSize);
01635     MVT::MvNorm(*R_active,d_resNorms);
01636 
01637     // If Ritz value i is real, then the corresponding residual vector
01638     //  is the true residual
01639     // If Ritz values i and i+1 form a conjugate pair, then the
01640     //  corresponding residual vectors are the real and imaginary components
01641     //  of the residual.  Adjust the residual norms appropriately...
01642     for( int i=0; i<d_residualSize; ++i )
01643     {
01644         if( d_ritzIndex[i] == 1 )
01645         {
01646             MagnitudeType nrm = lapack.LAPY2(d_resNorms[i],d_resNorms[i+1]);
01647             d_resNorms[i]   = nrm;
01648             d_resNorms[i+1] = nrm;
01649         }
01650     }
01651 
01652     // Evaluate with status test
01653     d_tester->checkStatus(this);
01654 
01655     // Determine which residual vectors should be used for subspace expansion
01656     if( d_useLeading || d_blockSize >= d_NEV )
01657     {
01658         d_expansionSize=d_blockSize;
01659         if( d_ritzIndex[d_blockSize-1]==1 )
01660             d_expansionSize++;
01661 
01662         d_expansionIndices.resize(d_expansionSize);
01663         for( int i=0; i<d_expansionSize; ++i )
01664             d_expansionIndices[i] = i;
01665     }
01666     else
01667     {
01668         std::vector<int> convergedVectors = d_tester->whichVecs();
01669 
01670         // Get index of first unconverged vector
01671         int startVec;
01672         for( startVec=0; startVec<d_residualSize; ++startVec )
01673         {
01674             if( std::find(convergedVectors.begin(),convergedVectors.end(),startVec)==convergedVectors.end() )
01675                 break;
01676         }
01677 
01678         // Now get a contiguous block of indices starting at startVec
01679         // If this crosses the end of our residual vectors, take the final d_blockSize vectors
01680         int endVec = startVec + d_blockSize - 1;
01681         if( endVec > (d_residualSize-1) )
01682         {
01683             endVec   = d_residualSize-1;
01684             startVec = d_residualSize-d_blockSize;
01685         }
01686 
01687         // Don't split conjugate pairs on either end of the range
01688         if( d_ritzIndex[startVec]==-1 )
01689         {
01690             startVec--;
01691             endVec--;
01692         }
01693 
01694         if( d_ritzIndex[endVec]==1 )
01695             endVec++;
01696 
01697         d_expansionSize = 1+endVec-startVec;
01698         d_expansionIndices.resize(d_expansionSize);
01699         for( int i=0; i<d_expansionSize; ++i )
01700             d_expansionIndices[i] = startVec+i;
01701     }
01702 }
01703 
01704 //---------------------------------------------------------------------------//
01705 // Print current status.
01706 //---------------------------------------------------------------------------//
01707 template <class ScalarType, class MV, class OP>
01708 void GeneralizedDavidson<ScalarType,MV,OP>::currentStatus( std::ostream &myout )
01709 {
01710     using std::endl;
01711 
01712     myout.setf(std::ios::scientific, std::ios::floatfield);
01713     myout.precision(6);
01714     myout <<endl;
01715     myout <<"================================================================================" << endl;
01716     myout << endl;
01717     myout <<"                    GeneralizedDavidson Solver Solver Status" << endl;
01718     myout << endl;
01719     myout <<"The solver is "<<(d_initialized ? "initialized." : "not initialized.") << endl;
01720     myout <<"The number of iterations performed is " << d_iteration << endl;
01721     myout <<"The number of operator applies performed is " << d_opApplies << endl;
01722     myout <<"The block size is         " << d_expansionSize << endl;
01723     myout <<"The current basis size is " << d_curDim << endl;
01724     myout <<"The number of requested eigenvalues is " << d_NEV << endl;
01725     myout <<"The number of converged values is " << d_tester->howMany() << endl;
01726     myout << endl;
01727 
01728     myout.setf(std::ios_base::right, std::ios_base::adjustfield);
01729 
01730     if( d_initialized )
01731     {
01732         myout << "CURRENT RITZ VALUES" << endl;
01733 
01734         myout << std::setw(24) << "Ritz Value"
01735               << std::setw(30) << "Residual Norm" << endl;
01736         myout << "--------------------------------------------------------------------------------" << endl;
01737         if( d_residualSize > 0 )
01738         {
01739             std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
01740             std::vector< Value<ScalarType> > ritzVals = getRitzValues();
01741             for( int i=0; i<d_curDim; ++i )
01742             {
01743                 realRitz[i] = ritzVals[i].realpart;
01744                 imagRitz[i] = ritzVals[i].imagpart;
01745             }
01746             std::vector<int> permvec;
01747             sortValues( realRitz, imagRitz, permvec, d_curDim );
01748 
01749             int numToPrint = std::max( d_numToPrint, d_NEV );
01750             numToPrint = std::min( d_curDim, numToPrint );
01751 
01752             // Because the sort manager does not use a stable sort, occasionally
01753             //  the portion of a conjugate pair with negative imaginary part will be placed
01754             //  first...in that case the following will not give the usual expected behavior
01755             //  and an extra value will be printed.  This is only an issue with the output
01756             //  format because the actually sorting of Schur forms is guaranteed to be stable.
01757             if( d_ritzIndex[permvec[numToPrint-1]] != 0 )
01758                 numToPrint++;
01759 
01760             int i=0;
01761             while( i<numToPrint )
01762             {
01763                 if( imagRitz[i] == ST::zero() )
01764                 {
01765                     myout << std::setw(15) << realRitz[i];
01766                     myout << " + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
01767                     if( i < d_residualSize )
01768                         myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
01769                     else
01770                         myout << "        Not Computed" << endl;
01771 
01772                     i++;
01773                 }
01774                 else
01775                 {
01776                     // Positive imaginary part
01777                     myout << std::setw(15) << realRitz[i];
01778                     myout << " + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
01779                     if( i < d_residualSize )
01780                         myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
01781                     else
01782                         myout << "        Not Computed" << endl;
01783 
01784                     // Negative imaginary part
01785                     myout << std::setw(15) << realRitz[i];
01786                     myout << " - i" << std::setw(15) << ST::magnitude( imagRitz[i] );
01787                     if( i < d_residualSize )
01788                         myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
01789                     else
01790                         myout << "        Not Computed" << endl;
01791 
01792                     i+=2;
01793                 }
01794             }
01795         }
01796         else
01797         {
01798             myout << std::setw(20) << "[ NONE COMPUTED ]" << endl;
01799         }
01800     }
01801     myout << endl;
01802     myout << "================================================================================" << endl;
01803     myout << endl;
01804 }
01805 
01806 } // namespace Anasazi
01807 
01808 #endif // ANASAZI_GENERALIZED_DAVIDSON_HPP
01809 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends