29 #ifndef ANASAZI_GENERALIZED_DAVIDSON_HPP 30 #define ANASAZI_GENERALIZED_DAVIDSON_HPP 38 #include "Teuchos_RCPDecl.hpp" 39 #include "Teuchos_ParameterList.hpp" 40 #include "Teuchos_SerialDenseMatrix.hpp" 41 #include "Teuchos_SerialDenseVector.hpp" 42 #include "Teuchos_Array.hpp" 43 #include "Teuchos_BLAS.hpp" 44 #include "Teuchos_LAPACK.hpp" 62 template <
class ScalarType,
class MV>
77 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
VAV;
80 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
VBV;
83 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
S;
86 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
T;
89 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
Q;
92 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >
Z;
95 std::vector< Value<ScalarType> >
eVals;
98 BV(Teuchos::null), VAV(Teuchos::null),
99 VBV(Teuchos::null), S(Teuchos::null),
100 T(Teuchos::null), Q(Teuchos::null),
101 Z(Teuchos::null), eVals(0) {}
126 template <
class ScalarType,
class MV,
class OP>
133 typedef Teuchos::ScalarTraits<ScalarType> ST;
134 typedef typename ST::magnitudeType MagnitudeType;
135 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
164 Teuchos::ParameterList &pl);
202 if( !d_ritzVectorsValid )
203 computeRitzVectors();
210 std::vector< Value<ScalarType> > getRitzValues();
217 if( !d_ritzIndexValid )
229 return d_expansionIndices;
235 std::vector<MagnitudeType> getResNorms();
240 std::vector<MagnitudeType> getResNorms(
int numWanted);
273 RCP<StatusTest<ScalarType,MV,OP> >
getStatusTest()
const {
return d_tester; }
288 void setBlockSize(
int blockSize);
293 void setSize(
int blockSize,
int maxSubDim);
298 Teuchos::Array< RCP<const MV> >
getAuxVecs()
const {
return d_auxVecs; }
307 void setAuxVecs(
const Teuchos::Array< RCP<const MV> > &auxVecs );
317 void currentStatus( std::ostream &myout );
327 void sortProblem(
int numWanted );
332 void expandSearchSpace();
335 void applyOperators();
338 void updateProjections();
341 void solveProjectedEigenproblem();
344 void computeProjectedEigenvectors( Teuchos::SerialDenseMatrix<int,ScalarType> &X );
347 void scaleEigenvectors(
const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
348 Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
349 Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta );
352 void sortValues( std::vector<MagnitudeType> &realParts,
353 std::vector<MagnitudeType> &imagParts,
354 std::vector<int> &permVec,
358 void computeResidual();
361 void computeRitzIndex();
364 void computeRitzVectors();
367 RCP<Eigenproblem<ScalarType,MV,OP> > d_problem;
368 Teuchos::ParameterList d_pl;
377 int d_maxSubspaceDim;
380 std::string d_initType;
382 bool d_relativeConvergence;
385 RCP<OutputManager<ScalarType> > d_outputMan;
386 RCP<OrthoManager<ScalarType,MV> > d_orthoMan;
387 RCP<SortManager<MagnitudeType> > d_sortMan;
388 RCP<StatusTest<ScalarType,MV,OP> > d_tester;
391 std::vector< Value<ScalarType> > d_eigenvalues;
395 std::vector<MagnitudeType> d_resNorms;
401 RCP<MV> d_ritzVecSpace;
403 Teuchos::Array< RCP<const MV> > d_auxVecs;
406 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VAV;
407 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VBV;
408 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_S;
409 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_T;
410 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Q;
411 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Z;
414 std::vector<MagnitudeType> d_alphar;
415 std::vector<MagnitudeType> d_alphai;
416 std::vector<MagnitudeType> d_betar;
417 std::vector<int> d_ritzIndex;
418 std::vector<int> d_convergedIndices;
419 std::vector<int> d_expansionIndices;
432 std::string d_testSpace;
440 bool d_ritzIndexValid;
441 bool d_ritzVectorsValid;
448 template <
class MagnitudeType,
class MV,
class OP>
453 typedef std::complex<MagnitudeType> ScalarType;
460 Teuchos::ParameterList &pl)
463 MagnitudeType::this_class_is_missing_a_specialization();
474 template <
class ScalarType,
class MV,
class OP>
481 Teuchos::ParameterList &pl )
483 TEUCHOS_TEST_FOR_EXCEPTION( problem == Teuchos::null, std::invalid_argument,
"No Eigenproblem given to solver." );
484 TEUCHOS_TEST_FOR_EXCEPTION( outputman == Teuchos::null, std::invalid_argument,
"No OutputManager given to solver." );
485 TEUCHOS_TEST_FOR_EXCEPTION( orthoman == Teuchos::null, std::invalid_argument,
"No OrthoManager given to solver." );
486 TEUCHOS_TEST_FOR_EXCEPTION( sortman == Teuchos::null, std::invalid_argument,
"No SortManager given to solver." );
487 TEUCHOS_TEST_FOR_EXCEPTION( tester == Teuchos::null, std::invalid_argument,
"No StatusTest given to solver." );
488 TEUCHOS_TEST_FOR_EXCEPTION( !problem->isProblemSet(), std::invalid_argument,
"Problem has not been set." );
492 TEUCHOS_TEST_FOR_EXCEPTION( problem->getA()==Teuchos::null && problem->getOperator()==Teuchos::null,
493 std::invalid_argument,
"Either A or Operator must be non-null on Eigenproblem");
494 d_A = problem->getA()!=Teuchos::null ? problem->getA() : problem->getOperator();
495 d_B = problem->getM();
496 d_P = problem->getPrec();
498 d_outputMan = outputman;
500 d_orthoMan = orthoman;
503 d_NEV = d_problem->getNEV();
504 d_initType = d_pl.get<std::string>(
"Initial Guess",
"Random");
505 d_numToPrint = d_pl.get<
int>(
"Print Number of Ritz Values",-1);
506 d_useLeading = d_pl.get<
bool>(
"Use Leading Vectors",
false);
508 if( d_B != Teuchos::null )
513 if( d_P != Teuchos::null )
518 d_testSpace = d_pl.get<std::string>(
"Test Space",
"V");
519 TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace!=
"V" && d_testSpace!=
"AV" && d_testSpace!=
"BV", std::invalid_argument,
520 "Anasazi::GeneralizedDavidson: Test Space must be V, AV, or BV" );
521 TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace==
"V" ?
false : !d_haveB, std::invalid_argument,
522 "Anasazi::GeneralizedDavidson: Test Space must be V for standard eigenvalue problem" );
525 int blockSize = d_pl.get<
int>(
"Block Size",1);
526 int maxSubDim = d_pl.get<
int>(
"Maximum Subspace Dimension",3*d_NEV*blockSize);
528 d_maxSubspaceDim = -1;
529 setSize( blockSize, maxSubDim );
530 d_relativeConvergence = d_pl.get<
bool>(
"Relative Convergence Tolerance",
false);
533 TEUCHOS_TEST_FOR_EXCEPTION( d_blockSize <= 0, std::invalid_argument,
"Block size must be positive");
534 TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim <= 0, std::invalid_argument,
"Maximum Subspace Dimension must be positive" );
535 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV()+2 > pl.get<
int>(
"Maximum Subspace Dimension"),
536 std::invalid_argument,
"Maximum Subspace Dimension must be strictly greater than NEV");
537 TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim > MVT::GetGlobalLength(*problem->getInitVec()),
538 std::invalid_argument,
"Maximum Subspace Dimension cannot exceed problem size");
544 d_ritzIndexValid =
false;
545 d_ritzVectorsValid =
false;
552 template <
class ScalarType,
class MV,
class OP>
558 d_outputMan->stream(
Warnings) <<
"WARNING: GeneralizedDavidson::iterate called without first calling initialize" << std::endl;
559 d_outputMan->stream(
Warnings) <<
" Default initialization will be performed" << std::endl;
564 if( d_outputMan->isVerbosity(
Debug) )
566 currentStatus( d_outputMan->stream(
Debug) );
573 while( d_tester->getStatus() !=
Passed && d_curDim+d_expansionSize <= d_maxSubspaceDim )
583 solveProjectedEigenproblem();
588 int numToSort = std::max(d_blockSize,d_NEV);
589 numToSort = std::min(numToSort,d_curDim);
590 sortProblem( numToSort );
595 if( d_outputMan->isVerbosity(
Debug) )
597 currentStatus( d_outputMan->stream(
Debug) );
609 template <
class ScalarType,
class MV,
class OP>
623 state.
eVals = getRitzValues();
630 template <
class ScalarType,
class MV,
class OP>
633 setSize(blockSize,d_maxSubspaceDim);
639 template <
class ScalarType,
class MV,
class OP>
642 if( blockSize != d_blockSize || maxSubDim != d_maxSubspaceDim )
644 d_blockSize = blockSize;
645 d_maxSubspaceDim = maxSubDim;
646 d_initialized =
false;
648 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Allocating eigenproblem" 649 <<
" state with block size of " << d_blockSize
650 <<
" and maximum subspace dimension of " << d_maxSubspaceDim << std::endl;
653 d_alphar.resize(d_maxSubspaceDim);
654 d_alphai.resize(d_maxSubspaceDim);
655 d_betar.resize(d_maxSubspaceDim);
658 int msd = d_maxSubspaceDim;
661 RCP<const MV> initVec = d_problem->getInitVec();
664 d_V = MVT::Clone(*initVec, msd);
665 d_AV = MVT::Clone(*initVec, msd);
668 d_VAV = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
669 d_S = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
670 d_Q = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
671 d_Z = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
676 d_BV = MVT::Clone(*initVec, msd);
677 d_VBV = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
678 d_T = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
689 d_R = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
690 d_ritzVecSpace = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
705 template <
class ScalarType,
class MV,
class OP>
710 d_curDim = (state.
curDim > 0 ? state.
curDim : d_blockSize );
712 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Initializing state" 713 <<
" with subspace dimension " << d_curDim << std::endl;
716 std::vector<int> initInds(d_curDim);
717 for(
int i=0; i<d_curDim; ++i )
721 RCP<MV> V1 = MVT::CloneViewNonConst(*d_V,initInds);
724 bool reset_V =
false;;
725 if( state.
curDim > 0 && state.
V != Teuchos::null && MVT::GetNumberVecs(*state.
V) >= d_curDim )
728 MVT::SetBlock(*state.
V,initInds,*V1);
732 else if( MVT::GetNumberVecs(*d_problem->getInitVec()) < d_blockSize || d_initType ==
"Random" )
740 RCP<const MV> initVec = MVT::CloneView(*d_problem->getInitVec(),initInds);
741 MVT::SetBlock(*initVec,initInds,*V1);
748 int rank = d_orthoMan->projectAndNormalize( *V1, d_auxVecs );
749 TEUCHOS_TEST_FOR_EXCEPTION( rank < d_blockSize, std::logic_error,
750 "Anasazi::GeneralizedDavidson::initialize(): Error generating initial orthonormal basis" );
753 if( d_outputMan->isVerbosity(
Debug) )
755 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Error in V^T V == I: " 756 << d_orthoMan->orthonormError( *V1 ) << std::endl;
760 RCP<MV> AV1 = MVT::CloneViewNonConst(*d_AV,initInds);
764 if( !reset_V && state.
AV != Teuchos::null && MVT::GetNumberVecs(*state.
AV) >= d_curDim )
766 if( state.
AV != d_AV )
767 MVT::SetBlock(*state.
AV,initInds,*AV1);
772 OPT::Apply( *d_A, *V1, *AV1 );
773 d_opApplies += MVT::GetNumberVecs( *V1 );
777 Teuchos::SerialDenseMatrix<int,ScalarType> VAV1( Teuchos::View, *d_VAV, d_curDim, d_curDim );
780 if( !reset_V && state.
VAV != Teuchos::null && state.
VAV->numRows() >= d_curDim && state.
VAV->numCols() >= d_curDim )
782 if( state.
VAV != d_VAV )
784 Teuchos::SerialDenseMatrix<int,ScalarType> state_VAV( Teuchos::View, *state.
VAV, d_curDim, d_curDim );
785 VAV1.assign( state_VAV );
791 if( d_testSpace ==
"V" )
793 MVT::MvTransMv( ST::one(), *V1, *AV1, VAV1 );
795 else if( d_testSpace ==
"AV" )
797 MVT::MvTransMv( ST::one(), *AV1, *AV1, VAV1 );
799 else if( d_testSpace ==
"BV" )
801 RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
802 MVT::MvTransMv( ST::one(), *BV1, *AV1, VAV1 );
809 RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
813 if( !reset_V && state.
BV != Teuchos::null && MVT::GetNumberVecs(*state.
BV) >= d_curDim )
815 if( state.
BV != d_BV )
816 MVT::SetBlock(*state.
BV,initInds,*BV1);
821 OPT::Apply( *d_B, *V1, *BV1 );
825 Teuchos::SerialDenseMatrix<int,ScalarType> VBV1( Teuchos::View, *d_VBV, d_curDim, d_curDim );
828 if( !reset_V && state.
VBV != Teuchos::null && state.
VBV->numRows() >= d_curDim && state.
VBV->numCols() >= d_curDim )
830 if( state.
VBV != d_VBV )
832 Teuchos::SerialDenseMatrix<int,ScalarType> state_VBV( Teuchos::View, *state.
VBV, d_curDim, d_curDim );
833 VBV1.assign( state_VBV );
839 if( d_testSpace ==
"V" )
841 MVT::MvTransMv( ST::one(), *V1, *BV1, VBV1 );
843 else if( d_testSpace ==
"AV" )
845 MVT::MvTransMv( ST::one(), *AV1, *BV1, VBV1 );
847 else if( d_testSpace ==
"BV" )
849 MVT::MvTransMv( ST::one(), *BV1, *BV1, VBV1 );
855 solveProjectedEigenproblem();
858 int numToSort = std::max(d_blockSize,d_NEV);
859 numToSort = std::min(numToSort,d_curDim);
860 sortProblem( numToSort );
866 d_initialized =
true;
872 template <
class ScalarType,
class MV,
class OP>
882 template <
class ScalarType,
class MV,
class OP>
883 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
886 return getResNorms(d_residualSize);
892 template <
class ScalarType,
class MV,
class OP>
893 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
896 std::vector<int> resIndices(numWanted);
897 for(
int i=0; i<numWanted; ++i )
900 RCP<const MV> R_checked = MVT::CloneView( *d_R, resIndices );
902 std::vector<MagnitudeType> resNorms;
903 d_orthoMan->norm( *R_checked, resNorms );
911 template <
class ScalarType,
class MV,
class OP>
914 std::vector< Value<ScalarType> > ritzValues;
915 for(
int ival=0; ival<d_curDim; ++ival )
918 thisVal.
realpart = d_alphar[ival] / d_betar[ival];
919 if( d_betar[ival] != MT::zero() )
920 thisVal.
imagpart = d_alphai[ival] / d_betar[ival];
924 ritzValues.push_back( thisVal );
939 template <
class ScalarType,
class MV,
class OP>
941 const Teuchos::Array< RCP<const MV> > &auxVecs )
946 typename Teuchos::Array< RCP<const MV> >::const_iterator arrItr;
948 for( arrItr=auxVecs.begin(); arrItr!=auxVecs.end(); ++arrItr )
950 numAuxVecs += MVT::GetNumberVecs( *(*arrItr) );
953 d_initialized =
false;
959 template <
class ScalarType,
class MV,
class OP>
963 std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
964 std::vector< Value<ScalarType> > ritzVals = getRitzValues();
965 for(
int i=0; i<d_curDim; ++i )
967 realRitz[i] = ritzVals[i].realpart;
968 imagRitz[i] = ritzVals[i].imagpart;
971 std::vector<int> permVec;
972 sortValues( realRitz, imagRitz, permVec, d_curDim );
974 std::vector<int> sel(d_curDim,0);
975 for(
int ii=0; ii<numWanted; ++ii )
976 sel[ permVec[ii] ]=1;
983 int work_size=10*d_maxSubspaceDim+16;
984 std::vector<ScalarType> work(work_size);
990 Teuchos::LAPACK<int,ScalarType> lapack;
991 lapack.TGSEN( ijob, wantq, wantz, &sel[0], d_curDim, d_S->values(), d_S->stride(), d_T->values(), d_T->stride(),
992 &d_alphar[0], &d_alphai[0], &d_betar[0], d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(),
993 &sdim, 0, 0, 0, &work[0], work_size, &iwork, iwork_size, &info );
995 d_ritzIndexValid =
false;
996 d_ritzVectorsValid =
false;
998 std::stringstream ss;
999 ss <<
"Anasazi::GeneralizedDavidson: TGSEN returned error code " << info << std::endl;
1000 TEUCHOS_TEST_FOR_EXCEPTION( info<0, std::runtime_error, ss.str() );
1006 d_outputMan->stream(
Warnings) <<
"WARNING: " << ss.str() << std::endl;
1007 d_outputMan->stream(
Warnings) <<
" Problem may not be correctly sorted" << std::endl;
1011 char getCondNum =
'N';
1014 int work_size = d_curDim;
1015 std::vector<ScalarType> work(work_size);
1020 Teuchos::LAPACK<int,ScalarType> lapack;
1021 lapack.TRSEN (getCondNum, getQ, &sel[0], d_curDim, d_S->values (),
1022 d_S->stride (), d_Z->values (), d_Z->stride (),
1023 &d_alphar[0], &d_alphai[0], &subDim, 0, 0, &work[0],
1024 work_size, &iwork, iwork_size, &info );
1026 std::fill( d_betar.begin(), d_betar.end(), ST::one() );
1028 d_ritzIndexValid =
false;
1029 d_ritzVectorsValid =
false;
1031 TEUCHOS_TEST_FOR_EXCEPTION(
1032 info < 0, std::runtime_error,
"Anasazi::GeneralizedDavidson::" 1033 "sortProblem: LAPACK routine TRSEN returned error code INFO = " 1034 << info <<
" < 0. This indicates that argument " << -info
1035 <<
" (counting starts with one) of TRSEN had an illegal value.");
1041 TEUCHOS_TEST_FOR_EXCEPTION(
1042 info == 1, std::runtime_error,
"Anasazi::GeneralizedDavidson::" 1043 "sortProblem: LAPACK routine TRSEN returned error code INFO = 1. " 1044 "This indicates that the reordering failed because some eigenvalues " 1045 "are too close to separate (the problem is very ill-conditioned).");
1047 TEUCHOS_TEST_FOR_EXCEPTION(
1048 info > 1, std::logic_error,
"Anasazi::GeneralizedDavidson::" 1049 "sortProblem: LAPACK routine TRSEN returned error code INFO = " 1050 << info <<
" > 1. This should not be possible. It may indicate an " 1051 "error either in LAPACK itself, or in Teuchos' LAPACK wrapper.");
1063 template <
class ScalarType,
class MV,
class OP>
1068 std::vector<int> newIndices(d_expansionSize);
1069 for(
int i=0; i<d_expansionSize; ++i )
1071 newIndices[i] = d_curDim+i;
1075 std::vector<int> curIndices(d_curDim);
1076 for(
int i=0; i<d_curDim; ++i )
1080 RCP<MV> V_new = MVT::CloneViewNonConst( *d_V, newIndices);
1081 RCP<const MV> V_cur = MVT::CloneView( *d_V, curIndices);
1082 RCP<const MV> R_active = MVT::CloneView( *d_R, d_expansionIndices);
1087 OPT::Apply( *d_P, *R_active, *V_new );
1092 MVT::SetBlock( *R_active, newIndices, *d_V );
1096 Teuchos::Array< RCP<const MV> > against = d_auxVecs;
1097 against.push_back( V_cur );
1098 int rank = d_orthoMan->projectAndNormalize(*V_new,against);
1100 if( d_outputMan->isVerbosity(
Debug) )
1102 std::vector<int> allIndices(d_curDim+d_expansionSize);
1103 for(
int i=0; i<d_curDim+d_expansionSize; ++i )
1106 RCP<const MV> V_all = MVT::CloneView( *d_V, allIndices );
1108 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Error in V^T V == I: " 1109 << d_orthoMan->orthonormError( *V_all ) << std::endl;
1112 TEUCHOS_TEST_FOR_EXCEPTION( rank != d_expansionSize, std::runtime_error,
1113 "Anasazi::GeneralizedDavidson::ExpandSearchSpace(): Orthonormalization of new vectors failed" );
1120 template <
class ScalarType,
class MV,
class OP>
1124 std::vector<int> newIndices(d_expansionSize);
1125 for(
int i=0; i<d_expansionSize; ++i )
1126 newIndices[i] = d_curDim+i;
1129 RCP<const MV> V_new = MVT::CloneView( *d_V, newIndices );
1130 RCP<MV> AV_new = MVT::CloneViewNonConst( *d_AV, newIndices );
1133 OPT::Apply( *d_A, *V_new, *AV_new );
1134 d_opApplies += MVT::GetNumberVecs( *V_new );
1139 RCP<MV> BV_new = MVT::CloneViewNonConst( *d_BV, newIndices );
1140 OPT::Apply( *d_B, *V_new, *BV_new );
1147 template <
class ScalarType,
class MV,
class OP>
1151 std::vector<int> newIndices(d_expansionSize);
1152 for(
int i=0; i<d_expansionSize; ++i )
1153 newIndices[i] = d_curDim+i;
1155 std::vector<int> curIndices(d_curDim);
1156 for(
int i=0; i<d_curDim; ++i )
1159 std::vector<int> allIndices(d_curDim+d_expansionSize);
1160 for(
int i=0; i<d_curDim+d_expansionSize; ++i )
1164 RCP<const MV> W_new, W_all;
1165 if( d_testSpace ==
"V" )
1167 W_new = MVT::CloneView(*d_V, newIndices );
1168 W_all = MVT::CloneView(*d_V, allIndices );
1170 else if( d_testSpace ==
"AV" )
1172 W_new = MVT::CloneView(*d_AV, newIndices );
1173 W_all = MVT::CloneView(*d_AV, allIndices );
1175 else if( d_testSpace ==
"BV" )
1177 W_new = MVT::CloneView(*d_BV, newIndices );
1178 W_all = MVT::CloneView(*d_BV, allIndices );
1182 RCP<const MV> AV_new = MVT::CloneView(*d_AV, newIndices);
1183 RCP<const MV> AV_current = MVT::CloneView(*d_AV, curIndices);
1186 Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastrow( Teuchos::View, *d_VAV, d_expansionSize, d_curDim, d_curDim, 0 );
1187 MVT::MvTransMv( ST::one(), *W_new, *AV_current, VAV_lastrow );
1190 Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastcol( Teuchos::View, *d_VAV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
1191 MVT::MvTransMv( ST::one(), *W_all, *AV_new, VAV_lastcol );
1196 RCP<const MV> BV_new = MVT::CloneView(*d_BV, newIndices);
1197 RCP<const MV> BV_current = MVT::CloneView(*d_BV, curIndices);
1200 Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastrow( Teuchos::View, *d_VBV, d_expansionSize, d_curDim, d_curDim, 0 );
1201 MVT::MvTransMv( ST::one(), *W_new, *BV_current, VBV_lastrow );
1204 Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastcol( Teuchos::View, *d_VBV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
1205 MVT::MvTransMv( ST::one(), *W_all, *BV_new, VBV_lastcol );
1209 d_curDim += d_expansionSize;
1211 d_ritzIndexValid =
false;
1212 d_ritzVectorsValid =
false;
1218 template <
class ScalarType,
class MV,
class OP>
1226 d_S->assign(*d_VAV);
1227 d_T->assign(*d_VBV);
1230 char leftVecs =
'V';
1231 char rightVecs =
'V';
1232 char sortVals =
'N';
1235 Teuchos::LAPACK<int,ScalarType> lapack;
1239 std::vector<ScalarType> work(1);
1240 lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
1241 d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
1242 d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
1244 work_size = work[0];
1245 work.resize(work_size);
1246 lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
1247 d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
1248 d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
1250 d_ritzIndexValid =
false;
1251 d_ritzVectorsValid =
false;
1253 std::stringstream ss;
1254 ss <<
"Anasazi::GeneralizedDavidson: GGES returned error code " << info << std::endl;
1255 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1261 d_S->assign(*d_VAV);
1266 int work_size = 3*d_curDim;
1267 std::vector<ScalarType> work(work_size);
1270 Teuchos::LAPACK<int,ScalarType> lapack;
1271 lapack.GEES( vecs, d_curDim, d_S->values(), d_S->stride(), &sdim, &d_alphar[0], &d_alphai[0],
1272 d_Z->values(), d_Z->stride(), &work[0], work_size, 0, 0, &info);
1274 std::fill( d_betar.begin(), d_betar.end(), ST::one() );
1276 d_ritzIndexValid =
false;
1277 d_ritzVectorsValid =
false;
1279 std::stringstream ss;
1280 ss <<
"Anasazi::GeneralizedDavidson: GEES returned error code " << info << std::endl;
1281 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1293 template <
class ScalarType,
class MV,
class OP>
1296 if( d_ritzIndexValid )
1299 d_ritzIndex.resize( d_curDim );
1301 while( i < d_curDim )
1303 if( d_alphai[i] == ST::zero() )
1311 d_ritzIndex[i+1] = -1;
1315 d_ritzIndexValid =
true;
1326 template <
class ScalarType,
class MV,
class OP>
1329 if( d_ritzVectorsValid )
1336 std::vector<int> checkedIndices(d_residualSize);
1337 for(
int ii=0; ii<d_residualSize; ++ii )
1338 checkedIndices[ii] = ii;
1341 Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
1342 computeProjectedEigenvectors( X );
1345 Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View,X,d_curDim,d_residualSize);
1348 d_ritzVecs = MVT::CloneViewNonConst( *d_ritzVecSpace, checkedIndices );
1350 std::vector<int> curIndices(d_curDim);
1351 for(
int i=0; i<d_curDim; ++i )
1354 RCP<const MV> V_current = MVT::CloneView( *d_V, curIndices );
1357 MVT::MvTimesMatAddMv(ST::one(),*V_current,X_wanted,ST::zero(),*d_ritzVecs);
1360 std::vector<MagnitudeType> scale(d_residualSize);
1361 MVT::MvNorm( *d_ritzVecs, scale );
1362 Teuchos::LAPACK<int,ScalarType> lapack;
1363 for(
int i=0; i<d_residualSize; ++i )
1365 if( d_ritzIndex[i] == 0 )
1367 scale[i] = 1.0/scale[i];
1369 else if( d_ritzIndex[i] == 1 )
1371 MagnitudeType nrm = lapack.LAPY2(scale[i],scale[i+1]);
1373 scale[i+1] = 1.0/nrm;
1376 MVT::MvScale( *d_ritzVecs, scale );
1378 d_ritzVectorsValid =
true;
1385 template <
class ScalarType,
class MV,
class OP>
1387 std::vector<MagnitudeType> &imagParts,
1388 std::vector<int> &permVec,
1393 TEUCHOS_TEST_FOR_EXCEPTION( (
int) realParts.size()<N, std::runtime_error,
1394 "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1395 TEUCHOS_TEST_FOR_EXCEPTION( (
int) imagParts.size()<N, std::runtime_error,
1396 "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1398 RCP< std::vector<int> > rcpPermVec = Teuchos::rcpFromRef(permVec);
1400 d_sortMan->sort( realParts, imagParts, rcpPermVec, N );
1402 d_ritzIndexValid =
false;
1403 d_ritzVectorsValid =
false;
1417 template <
class ScalarType,
class MV,
class OP>
1419 Teuchos::SerialDenseMatrix<int,ScalarType> &X )
1421 int N = X.numRows();
1424 Teuchos::SerialDenseMatrix<int,ScalarType>
S(Teuchos::Copy, *d_S, N, N);
1425 Teuchos::SerialDenseMatrix<int,ScalarType>
T(Teuchos::Copy, *d_T, N, N);
1426 Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Q, N, N);
1428 char whichVecs =
'R';
1430 int work_size = 6*d_maxSubspaceDim;
1431 std::vector<ScalarType> work(work_size,ST::zero());
1435 Teuchos::LAPACK<int,ScalarType> lapack;
1436 lapack.TGEVC( whichVecs, howMany, 0, N, S.values(), S.stride(), T.values(), T.stride(),
1437 VL.values(), VL.stride(), X.values(), X.stride(), N, &M, &work[0], &info );
1439 std::stringstream ss;
1440 ss <<
"Anasazi::GeneralizedDavidson: TGEVC returned error code " << info << std::endl;
1441 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1445 Teuchos::SerialDenseMatrix<int,ScalarType>
S(Teuchos::Copy, *d_S, N, N);
1446 Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Z, N, N);
1448 char whichVecs =
'R';
1451 std::vector<ScalarType> work(3*N);
1455 Teuchos::LAPACK<int,ScalarType> lapack;
1457 lapack.TREVC( whichVecs, howMany, &sel, N, S.values(), S.stride(), VL.values(), VL.stride(),
1458 X.values(), X.stride(), N, &m, &work[0], &info );
1460 std::stringstream ss;
1461 ss <<
"Anasazi::GeneralizedDavidson: TREVC returned error code " << info << std::endl;
1462 TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1469 template <
class ScalarType,
class MV,
class OP>
1471 const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
1472 Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
1473 Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta )
1475 int Nr = X.numRows();
1476 int Nc = X.numCols();
1478 TEUCHOS_TEST_FOR_EXCEPTION( Nr>d_curDim, std::logic_error,
1479 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1480 TEUCHOS_TEST_FOR_EXCEPTION( Nc>d_curDim, std::logic_error,
1481 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1482 TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numRows()!=Nr, std::logic_error,
1483 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xalpha does not match X");
1484 TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numCols()!=Nc, std::logic_error,
1485 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xalpha does not match X");
1486 TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numRows()!=Nr, std::logic_error,
1487 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xbeta does not match X");
1488 TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numCols()!=Nc, std::logic_error,
1489 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xbeta does not match X");
1493 Teuchos::SerialDenseMatrix<int,ScalarType> Alpha(Nc,Nc,
true);
1494 Teuchos::SerialDenseMatrix<int,ScalarType> Beta(Nc,Nc,
true);
1498 for(
int i=0; i<Nc; ++i )
1500 if( d_ritzIndex[i] == 0 )
1502 Alpha(i,i) = d_alphar[i];
1503 Beta(i,i) = d_betar[i];
1505 else if( d_ritzIndex[i] == 1 )
1507 Alpha(i,i) = d_alphar[i];
1508 Alpha(i,i+1) = d_alphai[i];
1509 Beta(i,i) = d_betar[i];
1513 Alpha(i,i-1) = d_alphai[i];
1514 Alpha(i,i) = d_alphar[i];
1515 Beta(i,i) = d_betar[i];
1522 err = X_alpha.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Alpha, ST::zero() );
1523 std::stringstream astream;
1524 astream <<
"GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1525 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, astream.str() );
1528 err = X_beta.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Beta, ST::zero() );
1529 std::stringstream bstream;
1530 bstream <<
"GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1531 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, bstream.str() );
1537 template <
class ScalarType,
class MV,
class OP>
1543 d_residualSize = std::max( d_blockSize, d_NEV );
1544 if( d_curDim < d_residualSize )
1546 d_residualSize = d_curDim;
1548 else if( d_ritzIndex[d_residualSize-1] == 1 )
1554 std::vector<int> residualIndices(d_residualSize);
1555 for(
int i=0; i<d_residualSize; ++i )
1556 residualIndices[i] = i;
1559 Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
1562 computeProjectedEigenvectors( X );
1565 Teuchos::SerialDenseMatrix<int,ScalarType> X_alpha(d_curDim,d_residualSize);
1566 Teuchos::SerialDenseMatrix<int,ScalarType> X_beta(d_curDim,d_residualSize);
1569 Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View, X, d_curDim, d_residualSize);
1572 scaleEigenvectors( X_wanted, X_alpha, X_beta );
1575 RCP<MV> R_active = MVT::CloneViewNonConst( *d_R, residualIndices );
1578 std::vector<int> activeIndices(d_curDim);
1579 for(
int i=0; i<d_curDim; ++i )
1583 RCP<const MV> AV_active = MVT::CloneView( *d_AV, activeIndices );
1584 MVT::MvTimesMatAddMv(ST::one(),*AV_active, X_beta, ST::zero(),*R_active);
1588 RCP<const MV> BV_active = MVT::CloneView( *d_BV, activeIndices );
1589 MVT::MvTimesMatAddMv(ST::one(),*BV_active, X_alpha,-ST::one(), *R_active);
1593 RCP<const MV> V_active = MVT::CloneView( *d_V, activeIndices );
1594 MVT::MvTimesMatAddMv(ST::one(),*V_active, X_alpha,-ST::one(), *R_active);
1609 Teuchos::LAPACK<int,ScalarType> lapack;
1610 Teuchos::BLAS<int,ScalarType> blas;
1611 std::vector<MagnitudeType> resScaling(d_residualSize);
1612 for(
int icol=0; icol<d_residualSize; ++icol )
1614 if( d_ritzIndex[icol] == 0 )
1616 MagnitudeType Xnrm = blas.NRM2( d_curDim, X_wanted[icol], 1);
1617 MagnitudeType ABscaling = d_relativeConvergence ? d_alphar[icol] : d_betar[icol];
1618 resScaling[icol] = MT::one() / (Xnrm * ABscaling);
1620 else if( d_ritzIndex[icol] == 1 )
1622 MagnitudeType Xnrm1 = blas.NRM2( d_curDim, X_wanted[icol], 1 );
1623 MagnitudeType Xnrm2 = blas.NRM2( d_curDim, X_wanted[icol+1], 1 );
1624 MagnitudeType Xnrm = lapack.LAPY2(Xnrm1,Xnrm2);
1625 MagnitudeType ABscaling = d_relativeConvergence ? lapack.LAPY2(d_alphar[icol],d_alphai[icol])
1627 resScaling[icol] = MT::one() / (Xnrm * ABscaling);
1628 resScaling[icol+1] = MT::one() / (Xnrm * ABscaling);
1631 MVT::MvScale( *R_active, resScaling );
1634 d_resNorms.resize(d_residualSize);
1635 MVT::MvNorm(*R_active,d_resNorms);
1642 for(
int i=0; i<d_residualSize; ++i )
1644 if( d_ritzIndex[i] == 1 )
1646 MagnitudeType nrm = lapack.LAPY2(d_resNorms[i],d_resNorms[i+1]);
1647 d_resNorms[i] = nrm;
1648 d_resNorms[i+1] = nrm;
1653 d_tester->checkStatus(
this);
1656 if( d_useLeading || d_blockSize >= d_NEV )
1658 d_expansionSize=d_blockSize;
1659 if( d_ritzIndex[d_blockSize-1]==1 )
1662 d_expansionIndices.resize(d_expansionSize);
1663 for(
int i=0; i<d_expansionSize; ++i )
1664 d_expansionIndices[i] = i;
1668 std::vector<int> convergedVectors = d_tester->whichVecs();
1672 for( startVec=0; startVec<d_residualSize; ++startVec )
1674 if( std::find(convergedVectors.begin(),convergedVectors.end(),startVec)==convergedVectors.end() )
1680 int endVec = startVec + d_blockSize - 1;
1681 if( endVec > (d_residualSize-1) )
1683 endVec = d_residualSize-1;
1684 startVec = d_residualSize-d_blockSize;
1688 if( d_ritzIndex[startVec]==-1 )
1694 if( d_ritzIndex[endVec]==1 )
1697 d_expansionSize = 1+endVec-startVec;
1698 d_expansionIndices.resize(d_expansionSize);
1699 for(
int i=0; i<d_expansionSize; ++i )
1700 d_expansionIndices[i] = startVec+i;
1707 template <
class ScalarType,
class MV,
class OP>
1712 myout.setf(std::ios::scientific, std::ios::floatfield);
1715 myout <<
"================================================================================" << endl;
1717 myout <<
" GeneralizedDavidson Solver Solver Status" << endl;
1719 myout <<
"The solver is "<<(d_initialized ?
"initialized." :
"not initialized.") << endl;
1720 myout <<
"The number of iterations performed is " << d_iteration << endl;
1721 myout <<
"The number of operator applies performed is " << d_opApplies << endl;
1722 myout <<
"The block size is " << d_expansionSize << endl;
1723 myout <<
"The current basis size is " << d_curDim << endl;
1724 myout <<
"The number of requested eigenvalues is " << d_NEV << endl;
1725 myout <<
"The number of converged values is " << d_tester->howMany() << endl;
1728 myout.setf(std::ios_base::right, std::ios_base::adjustfield);
1732 myout <<
"CURRENT RITZ VALUES" << endl;
1734 myout << std::setw(24) <<
"Ritz Value" 1735 << std::setw(30) <<
"Residual Norm" << endl;
1736 myout <<
"--------------------------------------------------------------------------------" << endl;
1737 if( d_residualSize > 0 )
1739 std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
1740 std::vector< Value<ScalarType> > ritzVals = getRitzValues();
1741 for(
int i=0; i<d_curDim; ++i )
1743 realRitz[i] = ritzVals[i].realpart;
1744 imagRitz[i] = ritzVals[i].imagpart;
1746 std::vector<int> permvec;
1747 sortValues( realRitz, imagRitz, permvec, d_curDim );
1749 int numToPrint = std::max( d_numToPrint, d_NEV );
1750 numToPrint = std::min( d_curDim, numToPrint );
1757 if( d_ritzIndex[permvec[numToPrint-1]] != 0 )
1761 while( i<numToPrint )
1763 if( imagRitz[i] == ST::zero() )
1765 myout << std::setw(15) << realRitz[i];
1766 myout <<
" + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1767 if( i < d_residualSize )
1768 myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1770 myout <<
" Not Computed" << endl;
1777 myout << std::setw(15) << realRitz[i];
1778 myout <<
" + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1779 if( i < d_residualSize )
1780 myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1782 myout <<
" Not Computed" << endl;
1785 myout << std::setw(15) << realRitz[i];
1786 myout <<
" - i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1787 if( i < d_residualSize )
1788 myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1790 myout <<
" Not Computed" << endl;
1798 myout << std::setw(20) <<
"[ NONE COMPUTED ]" << endl;
1802 myout <<
"================================================================================" << endl;
1808 #endif // ANASAZI_GENERALIZED_DAVIDSON_HPP void sortProblem(int numWanted)
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxilliary vectors.
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > tester)
Set status test.
RCP< MV > V
Orthonormal basis for search subspace.
RCP< MV > AV
Image of V under A.
std::vector< MagnitudeType > getRitzRes2Norms()
Get the current Ritz residual norms (2-norm)
Teuchos::ScalarTraits< ScalarType >::magnitudeType imagpart
The imaginary component of the eigenvalue.
std::vector< int > getBlockIndex() const
Get indices of current block.
This class defines the interface required by an eigensolver and status test class to compute solution...
Solves eigenvalue problem using generalized Davidson method.
std::vector< Value< ScalarType > > getRitzValues()
Get the current Ritz values.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > T
Right quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
Virtual base class which defines basic traits for the operator type.
GeneralizedDavidsonState< ScalarType, MV > getState()
Get the current state of the eigensolver.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VAV
Projection of A onto V.
int getMaxSubspaceDim() const
Get maximum subspace dimension.
void setBlockSize(int blockSize)
Set block size.
int curDim
The current subspace dimension.
std::vector< int > getRitzIndex()
Get the current Ritz index vector.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get eigenproblem.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
void iterate()
Solves the eigenvalue problem.
std::vector< MagnitudeType > getRes2Norms()
Get the current residual norms (2-norm)
RCP< MV > BV
Image of V under B.
Structure to contain pointers to GeneralizedDavidson state variables.
Abstract class definition for Anasazi Output Managers.
std::vector< Value< ScalarType > > eVals
Vector of generalized eigenvalues.
void setSize(int blockSize, int maxSubDim)
Set problem size.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Output managers remove the need for the eigensolver to know any information about the required output...
std::vector< MagnitudeType > getResNorms()
Get the current residual norms (w.r.t. norm defined by OrthoManager)
Templated virtual class for providing orthogonalization/orthonormalization methods.
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get status test.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Q
Left generalized Schur vectors from QZ decomposition of (VAV,VBV)
Traits class which defines basic operations on multivectors.
This struct is used for storing eigenvalues and Ritz values, as a pair of real values.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Virtual base class which defines the interface between an eigensolver and a class whose job is the so...
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.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > S
Left quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VBV
Projection of B onto V.
void initialize()
Initialize the eigenvalue problem.
bool isInitialized() const
Query if the solver is in an initialized state.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Z
Right generalized Schur vectors from QZ decomposition of (VAV,VBV)
void setAuxVecs(const Teuchos::Array< RCP< const MV > > &auxVecs)
Set auxilliary vectors.
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
int getBlockSize() const
Get block size.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
RCP< const MV > getRitzVectors()
Get the current Ritz vectors.
void currentStatus(std::ostream &myout)
Print current status of solver.
Teuchos::ScalarTraits< ScalarType >::magnitudeType realpart
The real component of the eigenvalue.
Common interface of stopping criteria for Anasazi's solvers.
int getNumIters() const
Get number of iterations.
int getCurSubspaceDim() const
Get current subspace dimension.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Declaration and definition of Anasazi::StatusTest.
void resetNumIters()
Reset the number of iterations.