42#ifndef ANASAZI_GENERALIZED_DAVIDSON_HPP
43#define ANASAZI_GENERALIZED_DAVIDSON_HPP
51#include "Teuchos_RCPDecl.hpp"
52#include "Teuchos_ParameterList.hpp"
53#include "Teuchos_SerialDenseMatrix.hpp"
54#include "Teuchos_SerialDenseVector.hpp"
55#include "Teuchos_Array.hpp"
56#include "Teuchos_BLAS.hpp"
57#include "Teuchos_LAPACK.hpp"
58#include "Teuchos_FancyOStream.hpp"
76template <
class ScalarType,
class MV>
91 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
VAV;
94 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
VBV;
97 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
S;
100 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
T;
103 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
Q;
106 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
Z;
109 std::vector< Value<ScalarType> >
eVals;
112 BV(Teuchos::null),
VAV(Teuchos::null),
113 VBV(Teuchos::null),
S(Teuchos::null),
114 T(Teuchos::null),
Q(Teuchos::null),
115 Z(Teuchos::null),
eVals(0) {}
140template <
class ScalarType,
class MV,
class OP>
145 typedef MultiVecTraits<ScalarType,MV> MVT;
146 typedef OperatorTraits<ScalarType,MV,OP> OPT;
147 typedef Teuchos::ScalarTraits<ScalarType> ST;
148 typedef typename ST::magnitudeType MagnitudeType;
149 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
174 const RCP<SortManager<MagnitudeType> > &sortman,
175 const RCP<OutputManager<ScalarType> > &outputman,
176 const RCP<StatusTest<ScalarType,MV,OP> > &tester,
177 const RCP<OrthoManager<ScalarType,MV> > &orthoman,
178 Teuchos::ParameterList &pl);
199 void initialize( GeneralizedDavidsonState<ScalarType,MV>& state );
216 if( !d_ritzVectorsValid )
217 computeRitzVectors();
231 if( !d_ritzIndexValid )
243 return d_expansionIndices;
254 std::vector<MagnitudeType>
getResNorms(
int numWanted);
282 void setStatusTest( RCP<StatusTest<ScalarType,MV,OP> > tester) { d_tester = tester; }
287 RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const {
return d_tester; }
292 const Eigenproblem<ScalarType,MV,OP> &
getProblem()
const {
return *d_problem; }
307 void setSize(
int blockSize,
int maxSubDim);
312 Teuchos::Array< RCP<const MV> >
getAuxVecs()
const {
return d_auxVecs; }
321 void setAuxVecs(
const Teuchos::Array< RCP<const MV> > &auxVecs );
336 GeneralizedDavidsonState<ScalarType,MV>
getState();
346 void expandSearchSpace();
349 void applyOperators();
352 void updateProjections();
355 void solveProjectedEigenproblem();
358 void computeProjectedEigenvectors( Teuchos::SerialDenseMatrix<int,ScalarType> &X );
361 void scaleEigenvectors(
const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
362 Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
363 Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta );
366 void sortValues( std::vector<MagnitudeType> &realParts,
367 std::vector<MagnitudeType> &imagParts,
368 std::vector<int> &permVec,
372 void computeResidual();
375 void computeRitzIndex();
378 void computeRitzVectors();
381 RCP<Eigenproblem<ScalarType,MV,OP> > d_problem;
382 Teuchos::ParameterList d_pl;
391 int d_maxSubspaceDim;
394 std::string d_initType;
396 bool d_relativeConvergence;
399 RCP<OutputManager<ScalarType> > d_outputMan;
400 RCP<OrthoManager<ScalarType,MV> > d_orthoMan;
401 RCP<SortManager<MagnitudeType> > d_sortMan;
402 RCP<StatusTest<ScalarType,MV,OP> > d_tester;
405 std::vector< Value<ScalarType> > d_eigenvalues;
409 std::vector<MagnitudeType> d_resNorms;
415 RCP<MV> d_ritzVecSpace;
417 Teuchos::Array< RCP<const MV> > d_auxVecs;
420 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VAV;
421 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VBV;
422 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_S;
423 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_T;
424 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Q;
425 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Z;
428 std::vector<MagnitudeType> d_alphar;
429 std::vector<MagnitudeType> d_alphai;
430 std::vector<MagnitudeType> d_betar;
431 std::vector<int> d_ritzIndex;
432 std::vector<int> d_convergedIndices;
433 std::vector<int> d_expansionIndices;
446 std::string d_testSpace;
454 bool d_ritzIndexValid;
455 bool d_ritzVectorsValid;
462template <
class MagnitudeType,
class MV,
class OP>
463class GeneralizedDavidson<std::complex<MagnitudeType>,MV,OP>
467 typedef std::complex<MagnitudeType> ScalarType;
469 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
470 const RCP<SortManager<MagnitudeType> > &sortman,
471 const RCP<OutputManager<ScalarType> > &outputman,
472 const RCP<StatusTest<ScalarType,MV,OP> > &tester,
473 const RCP<OrthoManager<ScalarType,MV> > &orthoman,
474 Teuchos::ParameterList &pl)
477 MagnitudeType::this_class_is_missing_a_specialization();
488template <
class ScalarType,
class MV,
class OP>
490 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
491 const RCP<SortManager<MagnitudeType> > &sortman,
492 const RCP<OutputManager<ScalarType> > &outputman,
493 const RCP<StatusTest<ScalarType,MV,OP> > &tester,
494 const RCP<OrthoManager<ScalarType,MV> > &orthoman,
495 Teuchos::ParameterList &pl )
497 TEUCHOS_TEST_FOR_EXCEPTION( problem == Teuchos::null, std::invalid_argument,
"No Eigenproblem given to solver." );
498 TEUCHOS_TEST_FOR_EXCEPTION( outputman == Teuchos::null, std::invalid_argument,
"No OutputManager given to solver." );
499 TEUCHOS_TEST_FOR_EXCEPTION( orthoman == Teuchos::null, std::invalid_argument,
"No OrthoManager given to solver." );
500 TEUCHOS_TEST_FOR_EXCEPTION( sortman == Teuchos::null, std::invalid_argument,
"No SortManager given to solver." );
501 TEUCHOS_TEST_FOR_EXCEPTION( tester == Teuchos::null, std::invalid_argument,
"No StatusTest given to solver." );
502 TEUCHOS_TEST_FOR_EXCEPTION( !problem->isProblemSet(), std::invalid_argument,
"Problem has not been set." );
506 TEUCHOS_TEST_FOR_EXCEPTION( problem->getA()==Teuchos::null && problem->getOperator()==Teuchos::null,
507 std::invalid_argument,
"Either A or Operator must be non-null on Eigenproblem");
508 d_A = problem->getA()!=Teuchos::null ? problem->getA() : problem->getOperator();
509 d_B = problem->getM();
510 d_P = problem->getPrec();
512 d_outputMan = outputman;
514 d_orthoMan = orthoman;
517 d_NEV = d_problem->getNEV();
518 d_initType = d_pl.get<std::string>(
"Initial Guess",
"Random");
519 d_numToPrint = d_pl.get<
int>(
"Print Number of Ritz Values",-1);
520 d_useLeading = d_pl.get<
bool>(
"Use Leading Vectors",
false);
522 if( d_B != Teuchos::null )
527 if( d_P != Teuchos::null )
532 d_testSpace = d_pl.get<std::string>(
"Test Space",
"V");
533 TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace!=
"V" && d_testSpace!=
"AV" && d_testSpace!=
"BV", std::invalid_argument,
534 "Anasazi::GeneralizedDavidson: Test Space must be V, AV, or BV" );
535 TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace==
"V" ?
false : !d_haveB, std::invalid_argument,
536 "Anasazi::GeneralizedDavidson: Test Space must be V for standard eigenvalue problem" );
539 int blockSize = d_pl.get<
int>(
"Block Size",1);
540 int maxSubDim = d_pl.get<
int>(
"Maximum Subspace Dimension",3*d_NEV*blockSize);
542 d_maxSubspaceDim = -1;
543 setSize( blockSize, maxSubDim );
544 d_relativeConvergence = d_pl.get<
bool>(
"Relative Convergence Tolerance",
false);
547 TEUCHOS_TEST_FOR_EXCEPTION( d_blockSize <= 0, std::invalid_argument,
"Block size must be positive");
548 TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim <= 0, std::invalid_argument,
"Maximum Subspace Dimension must be positive" );
549 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV()+2 > pl.get<
int>(
"Maximum Subspace Dimension"),
550 std::invalid_argument,
"Maximum Subspace Dimension must be strictly greater than NEV");
551 TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim > MVT::GetGlobalLength(*problem->getInitVec()),
552 std::invalid_argument,
"Maximum Subspace Dimension cannot exceed problem size");
558 d_ritzIndexValid =
false;
559 d_ritzVectorsValid =
false;
566template <
class ScalarType,
class MV,
class OP>
572 d_outputMan->stream(
Warnings) <<
"WARNING: GeneralizedDavidson::iterate called without first calling initialize" << std::endl;
573 d_outputMan->stream(
Warnings) <<
" Default initialization will be performed" << std::endl;
578 if( d_outputMan->isVerbosity(
Debug) )
580 currentStatus( d_outputMan->stream(
Debug) );
587 while( d_tester->getStatus() !=
Passed && d_curDim+d_expansionSize <= d_maxSubspaceDim )
597 solveProjectedEigenproblem();
602 int numToSort = std::max(d_blockSize,d_NEV);
603 numToSort = std::min(numToSort,d_curDim);
604 sortProblem( numToSort );
609 if( d_outputMan->isVerbosity(
Debug) )
611 currentStatus( d_outputMan->stream(
Debug) );
623template <
class ScalarType,
class MV,
class OP>
626 GeneralizedDavidsonState<ScalarType,MV> state;
627 state.curDim = d_curDim;
637 state.eVals = getRitzValues();
644template <
class ScalarType,
class MV,
class OP>
647 setSize(blockSize,d_maxSubspaceDim);
653template <
class ScalarType,
class MV,
class OP>
656 if( blockSize != d_blockSize || maxSubDim != d_maxSubspaceDim )
658 d_blockSize = blockSize;
659 d_maxSubspaceDim = maxSubDim;
660 d_initialized =
false;
662 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Allocating eigenproblem"
663 <<
" state with block size of " << d_blockSize
664 <<
" and maximum subspace dimension of " << d_maxSubspaceDim << std::endl;
667 d_alphar.resize(d_maxSubspaceDim);
668 d_alphai.resize(d_maxSubspaceDim);
669 d_betar.resize(d_maxSubspaceDim);
672 int msd = d_maxSubspaceDim;
675 RCP<const MV> initVec = d_problem->getInitVec();
678 d_V = MVT::Clone(*initVec, msd);
679 d_AV = MVT::Clone(*initVec, msd);
682 d_VAV = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
683 d_S = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
684 d_Q = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
685 d_Z = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
690 d_BV = MVT::Clone(*initVec, msd);
691 d_VBV = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
692 d_T = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
703 d_R = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
704 d_ritzVecSpace = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
719template <
class ScalarType,
class MV,
class OP>
724 d_curDim = (state.curDim > 0 ? state.curDim : d_blockSize );
726 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Initializing state"
727 <<
" with subspace dimension " << d_curDim << std::endl;
730 std::vector<int> initInds(d_curDim);
731 for(
int i=0; i<d_curDim; ++i )
735 RCP<MV> V1 = MVT::CloneViewNonConst(*d_V,initInds);
738 bool reset_V =
false;;
739 if( state.curDim > 0 && state.V != Teuchos::null && MVT::GetNumberVecs(*state.V) >= d_curDim )
742 MVT::SetBlock(*state.V,initInds,*V1);
746 else if( MVT::GetNumberVecs(*d_problem->getInitVec()) < d_blockSize || d_initType ==
"Random" )
754 RCP<const MV> initVec = MVT::CloneView(*d_problem->getInitVec(),initInds);
755 MVT::SetBlock(*initVec,initInds,*V1);
762 int rank = d_orthoMan->projectAndNormalize( *V1, d_auxVecs );
763 TEUCHOS_TEST_FOR_EXCEPTION( rank < d_blockSize, std::logic_error,
764 "Anasazi::GeneralizedDavidson::initialize(): Error generating initial orthonormal basis" );
767 if( d_outputMan->isVerbosity(
Debug) )
769 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Error in V^T V == I: "
770 << d_orthoMan->orthonormError( *V1 ) << std::endl;
774 RCP<MV> AV1 = MVT::CloneViewNonConst(*d_AV,initInds);
778 if( !reset_V && state.AV != Teuchos::null && MVT::GetNumberVecs(*state.AV) >= d_curDim )
780 if( state.AV != d_AV )
781 MVT::SetBlock(*state.AV,initInds,*AV1);
786 OPT::Apply( *d_A, *V1, *AV1 );
787 d_opApplies += MVT::GetNumberVecs( *V1 );
791 Teuchos::SerialDenseMatrix<int,ScalarType> VAV1( Teuchos::View, *d_VAV, d_curDim, d_curDim );
794 if( !reset_V && state.VAV != Teuchos::null && state.VAV->numRows() >= d_curDim && state.VAV->numCols() >= d_curDim )
796 if( state.VAV != d_VAV )
798 Teuchos::SerialDenseMatrix<int,ScalarType> state_VAV( Teuchos::View, *state.VAV, d_curDim, d_curDim );
799 VAV1.assign( state_VAV );
805 if( d_testSpace ==
"V" )
807 MVT::MvTransMv( ST::one(), *V1, *AV1, VAV1 );
809 else if( d_testSpace ==
"AV" )
811 MVT::MvTransMv( ST::one(), *AV1, *AV1, VAV1 );
813 else if( d_testSpace ==
"BV" )
815 RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
816 MVT::MvTransMv( ST::one(), *BV1, *AV1, VAV1 );
823 RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
827 if( !reset_V && state.BV != Teuchos::null && MVT::GetNumberVecs(*state.BV) >= d_curDim )
829 if( state.BV != d_BV )
830 MVT::SetBlock(*state.BV,initInds,*BV1);
835 OPT::Apply( *d_B, *V1, *BV1 );
839 Teuchos::SerialDenseMatrix<int,ScalarType> VBV1( Teuchos::View, *d_VBV, d_curDim, d_curDim );
842 if( !reset_V && state.VBV != Teuchos::null && state.VBV->numRows() >= d_curDim && state.VBV->numCols() >= d_curDim )
844 if( state.VBV != d_VBV )
846 Teuchos::SerialDenseMatrix<int,ScalarType> state_VBV( Teuchos::View, *state.VBV, d_curDim, d_curDim );
847 VBV1.assign( state_VBV );
853 if( d_testSpace ==
"V" )
855 MVT::MvTransMv( ST::one(), *V1, *BV1, VBV1 );
857 else if( d_testSpace ==
"AV" )
859 MVT::MvTransMv( ST::one(), *AV1, *BV1, VBV1 );
861 else if( d_testSpace ==
"BV" )
863 MVT::MvTransMv( ST::one(), *BV1, *BV1, VBV1 );
869 solveProjectedEigenproblem();
872 int numToSort = std::max(d_blockSize,d_NEV);
873 numToSort = std::min(numToSort,d_curDim);
874 sortProblem( numToSort );
880 d_initialized =
true;
886template <
class ScalarType,
class MV,
class OP>
889 GeneralizedDavidsonState<ScalarType,MV> empty;
896template <
class ScalarType,
class MV,
class OP>
897std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
900 return getResNorms(d_residualSize);
906template <
class ScalarType,
class MV,
class OP>
907std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
910 std::vector<int> resIndices(numWanted);
911 for(
int i=0; i<numWanted; ++i )
914 RCP<const MV> R_checked = MVT::CloneView( *d_R, resIndices );
916 std::vector<MagnitudeType> resNorms;
917 d_orthoMan->norm( *R_checked, resNorms );
925template <
class ScalarType,
class MV,
class OP>
928 std::vector< Value<ScalarType> > ritzValues;
929 for(
int ival=0; ival<d_curDim; ++ival )
931 Value<ScalarType> thisVal;
932 thisVal.realpart = d_alphar[ival] / d_betar[ival];
933 if( d_betar[ival] != MT::zero() )
934 thisVal.imagpart = d_alphai[ival] / d_betar[ival];
936 thisVal.imagpart = MT::zero();
938 ritzValues.push_back( thisVal );
953template <
class ScalarType,
class MV,
class OP>
955 const Teuchos::Array< RCP<const MV> > &auxVecs )
960 typename Teuchos::Array< RCP<const MV> >::const_iterator arrItr;
962 for( arrItr=auxVecs.begin(); arrItr!=auxVecs.end(); ++arrItr )
964 numAuxVecs += MVT::GetNumberVecs( *(*arrItr) );
967 d_initialized =
false;
973template <
class ScalarType,
class MV,
class OP>
977 std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
978 std::vector< Value<ScalarType> > ritzVals = getRitzValues();
979 for(
int i=0; i<d_curDim; ++i )
981 realRitz[i] = ritzVals[i].realpart;
982 imagRitz[i] = ritzVals[i].imagpart;
985 std::vector<int> permVec;
986 sortValues( realRitz, imagRitz, permVec, d_curDim );
988 std::vector<int> sel(d_curDim,0);
989 for(
int ii=0; ii<numWanted; ++ii )
990 sel[ permVec[ii] ]=1;
997 int work_size=10*d_maxSubspaceDim+16;
998 std::vector<ScalarType> work(work_size);
1004 Teuchos::LAPACK<int,ScalarType> lapack;
1005 lapack.TGSEN( ijob, wantq, wantz, &sel[0], d_curDim, d_S->values(), d_S->stride(), d_T->values(), d_T->stride(),
1006 &d_alphar[0], &d_alphai[0], &d_betar[0], d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(),
1007 &sdim, 0, 0, 0, &work[0], work_size, &iwork, iwork_size, &info );
1009 d_ritzIndexValid =
false;
1010 d_ritzVectorsValid =
false;
1012 std::stringstream ss;
1013 ss <<
"Anasazi::GeneralizedDavidson: TGSEN returned error code " << info << std::endl;
1014 TEUCHOS_TEST_FOR_EXCEPTION( info<0, std::runtime_error, ss.str() );
1020 d_outputMan->stream(
Warnings) <<
"WARNING: " << ss.str() << std::endl;
1021 d_outputMan->stream(
Warnings) <<
" Problem may not be correctly sorted" << std::endl;
1025 char getCondNum =
'N';
1028 int work_size = d_curDim;
1029 std::vector<ScalarType> work(work_size);
1034 Teuchos::LAPACK<int,ScalarType> lapack;
1035 lapack.TRSEN (getCondNum, getQ, &sel[0], d_curDim, d_S->values (),
1036 d_S->stride (), d_Z->values (), d_Z->stride (),
1037 &d_alphar[0], &d_alphai[0], &subDim, 0, 0, &work[0],
1038 work_size, &iwork, iwork_size, &info );
1040 std::fill( d_betar.begin(), d_betar.end(), ST::one() );
1042 d_ritzIndexValid =
false;
1043 d_ritzVectorsValid =
false;
1045 TEUCHOS_TEST_FOR_EXCEPTION(
1046 info < 0, std::runtime_error,
"Anasazi::GeneralizedDavidson::"
1047 "sortProblem: LAPACK routine TRSEN returned error code INFO = "
1048 << info <<
" < 0. This indicates that argument " << -info
1049 <<
" (counting starts with one) of TRSEN had an illegal value.");
1055 TEUCHOS_TEST_FOR_EXCEPTION(
1056 info == 1, std::runtime_error,
"Anasazi::GeneralizedDavidson::"
1057 "sortProblem: LAPACK routine TRSEN returned error code INFO = 1. "
1058 "This indicates that the reordering failed because some eigenvalues "
1059 "are too close to separate (the problem is very ill-conditioned).");
1061 TEUCHOS_TEST_FOR_EXCEPTION(
1062 info > 1, std::logic_error,
"Anasazi::GeneralizedDavidson::"
1063 "sortProblem: LAPACK routine TRSEN returned error code INFO = "
1064 << info <<
" > 1. This should not be possible. It may indicate an "
1065 "error either in LAPACK itself, or in Teuchos' LAPACK wrapper.");
1077template <
class ScalarType,
class MV,
class OP>
1082 std::vector<int> newIndices(d_expansionSize);
1083 for(
int i=0; i<d_expansionSize; ++i )
1085 newIndices[i] = d_curDim+i;
1089 std::vector<int> curIndices(d_curDim);
1090 for(
int i=0; i<d_curDim; ++i )
1094 RCP<MV> V_new = MVT::CloneViewNonConst( *d_V, newIndices);
1095 RCP<const MV> V_cur = MVT::CloneView( *d_V, curIndices);
1096 RCP<const MV> R_active = MVT::CloneView( *d_R, d_expansionIndices);
1101 OPT::Apply( *d_P, *R_active, *V_new );
1106 MVT::SetBlock( *R_active, newIndices, *d_V );
1110 Teuchos::Array< RCP<const MV> > against = d_auxVecs;
1111 against.push_back( V_cur );
1112 int rank = d_orthoMan->projectAndNormalize(*V_new,against);
1114 if( d_outputMan->isVerbosity(
Debug) )
1116 std::vector<int> allIndices(d_curDim+d_expansionSize);
1117 for(
int i=0; i<d_curDim+d_expansionSize; ++i )
1120 RCP<const MV> V_all = MVT::CloneView( *d_V, allIndices );
1122 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Error in V^T V == I: "
1123 << d_orthoMan->orthonormError( *V_all ) << std::endl;
1126 TEUCHOS_TEST_FOR_EXCEPTION( rank != d_expansionSize, std::runtime_error,
1127 "Anasazi::GeneralizedDavidson::ExpandSearchSpace(): Orthonormalization of new vectors failed" );
1134template <
class ScalarType,
class MV,
class OP>
1135void GeneralizedDavidson<ScalarType,MV,OP>::applyOperators()
1138 std::vector<int> newIndices(d_expansionSize);
1139 for(
int i=0; i<d_expansionSize; ++i )
1140 newIndices[i] = d_curDim+i;
1143 RCP<const MV> V_new = MVT::CloneView( *d_V, newIndices );
1144 RCP<MV> AV_new = MVT::CloneViewNonConst( *d_AV, newIndices );
1147 OPT::Apply( *d_A, *V_new, *AV_new );
1148 d_opApplies += MVT::GetNumberVecs( *V_new );
1153 RCP<MV> BV_new = MVT::CloneViewNonConst( *d_BV, newIndices );
1154 OPT::Apply( *d_B, *V_new, *BV_new );
1161template <
class ScalarType,
class MV,
class OP>
1162void GeneralizedDavidson<ScalarType,MV,OP>::updateProjections()
1165 std::vector<int> newIndices(d_expansionSize);
1166 for(
int i=0; i<d_expansionSize; ++i )
1167 newIndices[i] = d_curDim+i;
1169 std::vector<int> curIndices(d_curDim);
1170 for(
int i=0; i<d_curDim; ++i )
1173 std::vector<int> allIndices(d_curDim+d_expansionSize);
1174 for(
int i=0; i<d_curDim+d_expansionSize; ++i )
1178 RCP<const MV> W_new, W_all;
1179 if( d_testSpace ==
"V" )
1181 W_new = MVT::CloneView(*d_V, newIndices );
1182 W_all = MVT::CloneView(*d_V, allIndices );
1184 else if( d_testSpace ==
"AV" )
1186 W_new = MVT::CloneView(*d_AV, newIndices );
1187 W_all = MVT::CloneView(*d_AV, allIndices );
1189 else if( d_testSpace ==
"BV" )
1191 W_new = MVT::CloneView(*d_BV, newIndices );
1192 W_all = MVT::CloneView(*d_BV, allIndices );
1196 RCP<const MV> AV_new = MVT::CloneView(*d_AV, newIndices);
1197 RCP<const MV> AV_current = MVT::CloneView(*d_AV, curIndices);
1200 Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastrow( Teuchos::View, *d_VAV, d_expansionSize, d_curDim, d_curDim, 0 );
1201 MVT::MvTransMv( ST::one(), *W_new, *AV_current, VAV_lastrow );
1204 Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastcol( Teuchos::View, *d_VAV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
1205 MVT::MvTransMv( ST::one(), *W_all, *AV_new, VAV_lastcol );
1210 RCP<const MV> BV_new = MVT::CloneView(*d_BV, newIndices);
1211 RCP<const MV> BV_current = MVT::CloneView(*d_BV, curIndices);
1214 Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastrow( Teuchos::View, *d_VBV, d_expansionSize, d_curDim, d_curDim, 0 );
1215 MVT::MvTransMv( ST::one(), *W_new, *BV_current, VBV_lastrow );
1218 Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastcol( Teuchos::View, *d_VBV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
1219 MVT::MvTransMv( ST::one(), *W_all, *BV_new, VBV_lastcol );
1223 d_curDim += d_expansionSize;
1225 d_ritzIndexValid =
false;
1226 d_ritzVectorsValid =
false;
1232template <
class ScalarType,
class MV,
class OP>
1233void GeneralizedDavidson<ScalarType,MV,OP>::solveProjectedEigenproblem()
1240 d_S->assign(*d_VAV);
1241 d_T->assign(*d_VBV);
1244 char leftVecs =
'V';
1245 char rightVecs =
'V';
1246 char sortVals =
'N';
1249 Teuchos::LAPACK<int,ScalarType> lapack;
1253 std::vector<ScalarType> work(1);
1254 lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
1255 d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
1256 d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
1258 work_size = work[0];
1259 work.resize(work_size);
1260 lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
1261 d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
1262 d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
1264 d_ritzIndexValid =
false;
1265 d_ritzVectorsValid =
false;
1267 std::stringstream ss;
1268 ss <<
"Anasazi::GeneralizedDavidson: GGES returned error code " << info << std::endl;
1269 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1275 d_S->assign(*d_VAV);
1280 int work_size = 3*d_curDim;
1281 std::vector<ScalarType> work(work_size);
1284 Teuchos::LAPACK<int,ScalarType> lapack;
1285 lapack.GEES( vecs, d_curDim, d_S->values(), d_S->stride(), &sdim, &d_alphar[0], &d_alphai[0],
1286 d_Z->values(), d_Z->stride(), &work[0], work_size, 0, 0, &info);
1288 std::fill( d_betar.begin(), d_betar.end(), ST::one() );
1290 d_ritzIndexValid =
false;
1291 d_ritzVectorsValid =
false;
1293 std::stringstream ss;
1294 ss <<
"Anasazi::GeneralizedDavidson: GEES returned error code " << info << std::endl;
1295 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1307template <
class ScalarType,
class MV,
class OP>
1308void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzIndex()
1310 if( d_ritzIndexValid )
1313 d_ritzIndex.resize( d_curDim );
1315 while( i < d_curDim )
1317 if( d_alphai[i] == ST::zero() )
1325 d_ritzIndex[i+1] = -1;
1329 d_ritzIndexValid =
true;
1340template <
class ScalarType,
class MV,
class OP>
1341void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzVectors()
1343 if( d_ritzVectorsValid )
1350 std::vector<int> checkedIndices(d_residualSize);
1351 for(
int ii=0; ii<d_residualSize; ++ii )
1352 checkedIndices[ii] = ii;
1355 Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
1356 computeProjectedEigenvectors( X );
1359 Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View,X,d_curDim,d_residualSize);
1362 d_ritzVecs = MVT::CloneViewNonConst( *d_ritzVecSpace, checkedIndices );
1364 std::vector<int> curIndices(d_curDim);
1365 for(
int i=0; i<d_curDim; ++i )
1368 RCP<const MV> V_current = MVT::CloneView( *d_V, curIndices );
1371 MVT::MvTimesMatAddMv(ST::one(),*V_current,X_wanted,ST::zero(),*d_ritzVecs);
1374 std::vector<MagnitudeType> scale(d_residualSize);
1375 MVT::MvNorm( *d_ritzVecs, scale );
1376 Teuchos::LAPACK<int,ScalarType> lapack;
1377 for(
int i=0; i<d_residualSize; ++i )
1379 if( d_ritzIndex[i] == 0 )
1381 scale[i] = 1.0/scale[i];
1383 else if( d_ritzIndex[i] == 1 )
1385 MagnitudeType nrm = lapack.LAPY2(scale[i],scale[i+1]);
1387 scale[i+1] = 1.0/nrm;
1390 MVT::MvScale( *d_ritzVecs, scale );
1392 d_ritzVectorsValid =
true;
1399template <
class ScalarType,
class MV,
class OP>
1400void GeneralizedDavidson<ScalarType,MV,OP>::sortValues( std::vector<MagnitudeType> &realParts,
1401 std::vector<MagnitudeType> &imagParts,
1402 std::vector<int> &permVec,
1407 TEUCHOS_TEST_FOR_EXCEPTION( (
int) realParts.size()<N, std::runtime_error,
1408 "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1409 TEUCHOS_TEST_FOR_EXCEPTION( (
int) imagParts.size()<N, std::runtime_error,
1410 "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1412 RCP< std::vector<int> > rcpPermVec = Teuchos::rcpFromRef(permVec);
1414 d_sortMan->sort( realParts, imagParts, rcpPermVec, N );
1416 d_ritzIndexValid =
false;
1417 d_ritzVectorsValid =
false;
1431template <
class ScalarType,
class MV,
class OP>
1432void GeneralizedDavidson<ScalarType,MV,OP>::computeProjectedEigenvectors(
1433 Teuchos::SerialDenseMatrix<int,ScalarType> &X )
1435 int N = X.numRows();
1438 Teuchos::SerialDenseMatrix<int,ScalarType> S(Teuchos::Copy, *d_S, N, N);
1439 Teuchos::SerialDenseMatrix<int,ScalarType> T(Teuchos::Copy, *d_T, N, N);
1440 Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Q, N, N);
1442 char whichVecs =
'R';
1444 int work_size = 6*d_maxSubspaceDim;
1445 std::vector<ScalarType> work(work_size,ST::zero());
1449 Teuchos::LAPACK<int,ScalarType> lapack;
1450 lapack.TGEVC( whichVecs, howMany, 0, N, S.values(), S.stride(), T.values(), T.stride(),
1451 VL.values(), VL.stride(), X.values(), X.stride(), N, &M, &work[0], &info );
1453 std::stringstream ss;
1454 ss <<
"Anasazi::GeneralizedDavidson: TGEVC returned error code " << info << std::endl;
1455 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1459 Teuchos::SerialDenseMatrix<int,ScalarType> S(Teuchos::Copy, *d_S, N, N);
1460 Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Z, N, N);
1462 char whichVecs =
'R';
1465 std::vector<ScalarType> work(3*N);
1469 Teuchos::LAPACK<int,ScalarType> lapack;
1471 lapack.TREVC( whichVecs, howMany, &sel, N, S.values(), S.stride(), VL.values(), VL.stride(),
1472 X.values(), X.stride(), N, &m, &work[0], &info );
1474 std::stringstream ss;
1475 ss <<
"Anasazi::GeneralizedDavidson: TREVC returned error code " << info << std::endl;
1476 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1483template <
class ScalarType,
class MV,
class OP>
1484void GeneralizedDavidson<ScalarType,MV,OP>::scaleEigenvectors(
1485 const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
1486 Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
1487 Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta )
1489 int Nr = X.numRows();
1490 int Nc = X.numCols();
1492 TEUCHOS_TEST_FOR_EXCEPTION( Nr>d_curDim, std::logic_error,
1493 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1494 TEUCHOS_TEST_FOR_EXCEPTION( Nc>d_curDim, std::logic_error,
1495 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1496 TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numRows()!=Nr, std::logic_error,
1497 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xalpha does not match X");
1498 TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numCols()!=Nc, std::logic_error,
1499 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xalpha does not match X");
1500 TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numRows()!=Nr, std::logic_error,
1501 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xbeta does not match X");
1502 TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numCols()!=Nc, std::logic_error,
1503 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xbeta does not match X");
1507 Teuchos::SerialDenseMatrix<int,ScalarType> Alpha(Nc,Nc,
true);
1508 Teuchos::SerialDenseMatrix<int,ScalarType> Beta(Nc,Nc,
true);
1512 for(
int i=0; i<Nc; ++i )
1514 if( d_ritzIndex[i] == 0 )
1516 Alpha(i,i) = d_alphar[i];
1517 Beta(i,i) = d_betar[i];
1519 else if( d_ritzIndex[i] == 1 )
1521 Alpha(i,i) = d_alphar[i];
1522 Alpha(i,i+1) = d_alphai[i];
1523 Beta(i,i) = d_betar[i];
1527 Alpha(i,i-1) = d_alphai[i];
1528 Alpha(i,i) = d_alphar[i];
1529 Beta(i,i) = d_betar[i];
1536 err = X_alpha.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Alpha, ST::zero() );
1537 std::stringstream astream;
1538 astream <<
"GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1539 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, astream.str() );
1542 err = X_beta.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Beta, ST::zero() );
1543 std::stringstream bstream;
1544 bstream <<
"GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1545 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, bstream.str() );
1551template <
class ScalarType,
class MV,
class OP>
1552void GeneralizedDavidson<ScalarType,MV,OP>::computeResidual()
1557 d_residualSize = std::max( d_blockSize, d_NEV );
1558 if( d_curDim < d_residualSize )
1560 d_residualSize = d_curDim;
1562 else if( d_ritzIndex[d_residualSize-1] == 1 )
1568 std::vector<int> residualIndices(d_residualSize);
1569 for(
int i=0; i<d_residualSize; ++i )
1570 residualIndices[i] = i;
1573 Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
1576 computeProjectedEigenvectors( X );
1579 Teuchos::SerialDenseMatrix<int,ScalarType> X_alpha(d_curDim,d_residualSize);
1580 Teuchos::SerialDenseMatrix<int,ScalarType> X_beta(d_curDim,d_residualSize);
1583 Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View, X, d_curDim, d_residualSize);
1586 scaleEigenvectors( X_wanted, X_alpha, X_beta );
1589 RCP<MV> R_active = MVT::CloneViewNonConst( *d_R, residualIndices );
1592 std::vector<int> activeIndices(d_curDim);
1593 for(
int i=0; i<d_curDim; ++i )
1597 RCP<const MV> AV_active = MVT::CloneView( *d_AV, activeIndices );
1598 MVT::MvTimesMatAddMv(ST::one(),*AV_active, X_beta, ST::zero(),*R_active);
1602 RCP<const MV> BV_active = MVT::CloneView( *d_BV, activeIndices );
1603 MVT::MvTimesMatAddMv(ST::one(),*BV_active, X_alpha,-ST::one(), *R_active);
1607 RCP<const MV> V_active = MVT::CloneView( *d_V, activeIndices );
1608 MVT::MvTimesMatAddMv(ST::one(),*V_active, X_alpha,-ST::one(), *R_active);
1623 Teuchos::LAPACK<int,ScalarType> lapack;
1624 Teuchos::BLAS<int,ScalarType> blas;
1625 std::vector<MagnitudeType> resScaling(d_residualSize);
1626 for(
int icol=0; icol<d_residualSize; ++icol )
1628 if( d_ritzIndex[icol] == 0 )
1630 MagnitudeType Xnrm = blas.NRM2( d_curDim, X_wanted[icol], 1);
1631 MagnitudeType ABscaling = d_relativeConvergence ? d_alphar[icol] : d_betar[icol];
1632 resScaling[icol] = MT::one() / (Xnrm * ABscaling);
1634 else if( d_ritzIndex[icol] == 1 )
1636 MagnitudeType Xnrm1 = blas.NRM2( d_curDim, X_wanted[icol], 1 );
1637 MagnitudeType Xnrm2 = blas.NRM2( d_curDim, X_wanted[icol+1], 1 );
1638 MagnitudeType Xnrm = lapack.LAPY2(Xnrm1,Xnrm2);
1639 MagnitudeType ABscaling = d_relativeConvergence ? lapack.LAPY2(d_alphar[icol],d_alphai[icol])
1641 resScaling[icol] = MT::one() / (Xnrm * ABscaling);
1642 resScaling[icol+1] = MT::one() / (Xnrm * ABscaling);
1645 MVT::MvScale( *R_active, resScaling );
1648 d_resNorms.resize(d_residualSize);
1649 MVT::MvNorm(*R_active,d_resNorms);
1656 for(
int i=0; i<d_residualSize; ++i )
1658 if( d_ritzIndex[i] == 1 )
1660 MagnitudeType nrm = lapack.LAPY2(d_resNorms[i],d_resNorms[i+1]);
1661 d_resNorms[i] = nrm;
1662 d_resNorms[i+1] = nrm;
1667 d_tester->checkStatus(
this);
1670 if( d_useLeading || d_blockSize >= d_NEV )
1672 d_expansionSize=d_blockSize;
1673 if( d_ritzIndex[d_blockSize-1]==1 )
1676 d_expansionIndices.resize(d_expansionSize);
1677 for(
int i=0; i<d_expansionSize; ++i )
1678 d_expansionIndices[i] = i;
1682 std::vector<int> convergedVectors = d_tester->whichVecs();
1686 for( startVec=0; startVec<d_residualSize; ++startVec )
1688 if( std::find(convergedVectors.begin(),convergedVectors.end(),startVec)==convergedVectors.end() )
1694 int endVec = startVec + d_blockSize - 1;
1695 if( endVec > (d_residualSize-1) )
1697 endVec = d_residualSize-1;
1698 startVec = d_residualSize-d_blockSize;
1702 if( d_ritzIndex[startVec]==-1 )
1708 if( d_ritzIndex[endVec]==1 )
1711 d_expansionSize = 1+endVec-startVec;
1712 d_expansionIndices.resize(d_expansionSize);
1713 for(
int i=0; i<d_expansionSize; ++i )
1714 d_expansionIndices[i] = startVec+i;
1721template <
class ScalarType,
class MV,
class OP>
1726 myout.setf(std::ios::scientific, std::ios::floatfield);
1729 myout <<
"================================================================================" << endl;
1731 myout <<
" GeneralizedDavidson Solver Solver Status" << endl;
1733 myout <<
"The solver is "<<(d_initialized ?
"initialized." :
"not initialized.") << endl;
1734 myout <<
"The number of iterations performed is " << d_iteration << endl;
1735 myout <<
"The number of operator applies performed is " << d_opApplies << endl;
1736 myout <<
"The block size is " << d_expansionSize << endl;
1737 myout <<
"The current basis size is " << d_curDim << endl;
1738 myout <<
"The number of requested eigenvalues is " << d_NEV << endl;
1739 myout <<
"The number of converged values is " << d_tester->howMany() << endl;
1742 myout.setf(std::ios_base::right, std::ios_base::adjustfield);
1746 myout <<
"CURRENT RITZ VALUES" << endl;
1748 myout << std::setw(24) <<
"Ritz Value"
1749 << std::setw(30) <<
"Residual Norm" << endl;
1750 myout <<
"--------------------------------------------------------------------------------" << endl;
1751 if( d_residualSize > 0 )
1753 std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
1754 std::vector< Value<ScalarType> > ritzVals = getRitzValues();
1755 for(
int i=0; i<d_curDim; ++i )
1757 realRitz[i] = ritzVals[i].realpart;
1758 imagRitz[i] = ritzVals[i].imagpart;
1760 std::vector<int> permvec;
1761 sortValues( realRitz, imagRitz, permvec, d_curDim );
1763 int numToPrint = std::max( d_numToPrint, d_NEV );
1764 numToPrint = std::min( d_curDim, numToPrint );
1771 if( d_ritzIndex[permvec[numToPrint-1]] != 0 )
1775 while( i<numToPrint )
1777 if( imagRitz[i] == ST::zero() )
1779 myout << std::setw(15) << realRitz[i];
1780 myout <<
" + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1781 if( i < d_residualSize )
1782 myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1784 myout <<
" Not Computed" << endl;
1791 myout << std::setw(15) << realRitz[i];
1792 myout <<
" + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1793 if( i < d_residualSize )
1794 myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1796 myout <<
" Not Computed" << endl;
1799 myout << std::setw(15) << realRitz[i];
1800 myout <<
" - i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1801 if( i < d_residualSize )
1802 myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1804 myout <<
" Not Computed" << endl;
1812 myout << std::setw(20) <<
"[ NONE COMPUTED ]" << endl;
1816 myout <<
"================================================================================" << endl;
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Abstract class definition for Anasazi Output Managers.
Virtual base class which defines the interface between an eigensolver and a class whose job is the so...
Declaration and definition of Anasazi::StatusTest.
Types and exceptions used within Anasazi solvers and interfaces.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Solves eigenvalue problem using generalized Davidson method.
void sortProblem(int numWanted)
int getCurSubspaceDim() const
Get current subspace dimension.
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get status test.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get eigenproblem.
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > tester)
Set status test.
void currentStatus(std::ostream &myout)
Print current status of solver.
int getMaxSubspaceDim() const
Get maximum subspace dimension.
std::vector< MagnitudeType > getResNorms()
Get the current residual norms (w.r.t. norm defined by OrthoManager)
std::vector< int > getBlockIndex() const
Get indices of current block.
std::vector< int > getRitzIndex()
Get the current Ritz index vector.
GeneralizedDavidson(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const RCP< SortManager< MagnitudeType > > &sortman, const RCP< OutputManager< ScalarType > > &outputman, const RCP< StatusTest< ScalarType, MV, OP > > &tester, const RCP< OrthoManager< ScalarType, MV > > &orthoman, Teuchos::ParameterList &pl)
Constructor.
void initialize()
Initialize the eigenvalue problem.
GeneralizedDavidsonState< ScalarType, MV > getState()
Get the current state of the eigensolver.
bool isInitialized() const
Query if the solver is in an initialized state.
void setAuxVecs(const Teuchos::Array< RCP< const MV > > &auxVecs)
Set auxilliary vectors.
std::vector< Value< ScalarType > > getRitzValues()
Get the current Ritz values.
int getNumIters() const
Get number of iterations.
std::vector< MagnitudeType > getRes2Norms()
Get the current residual norms (2-norm)
void resetNumIters()
Reset the number of iterations.
RCP< const MV > getRitzVectors()
Get the current Ritz vectors.
int getBlockSize() const
Get block size.
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxilliary vectors.
void iterate()
Solves the eigenvalue problem.
std::vector< MagnitudeType > getRitzRes2Norms()
Get the current Ritz residual norms (2-norm)
void setSize(int blockSize, int maxSubDim)
Set problem size.
void setBlockSize(int blockSize)
Set block size.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Structure to contain pointers to GeneralizedDavidson state variables.
RCP< MV > AV
Image of V under A.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > T
Right quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VAV
Projection of A onto V.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Q
Left generalized Schur vectors from QZ decomposition of (VAV,VBV)
int curDim
The current subspace dimension.
std::vector< Value< ScalarType > > eVals
Vector of generalized eigenvalues.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VBV
Projection of B onto V.
RCP< MV > BV
Image of V under B.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Z
Right generalized Schur vectors from QZ decomposition of (VAV,VBV)
RCP< MV > V
Orthonormal basis for search subspace.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > S
Left quasi upper triangular matrix from QZ decomposition of (VAV,VBV)