29 #ifndef ANASAZI_GENERALIZED_DAVIDSON_HPP 30 #define ANASAZI_GENERALIZED_DAVIDSON_HPP 62 template <
class ScalarType,
class MV>
95 std::vector< Value<ScalarType> >
eVals;
126 template <
class ScalarType,
class MV,
class OP>
202 if( !d_ritzVectorsValid )
203 computeRitzVectors();
217 if( !d_ritzIndexValid )
229 return d_expansionIndices;
240 std::vector<MagnitudeType>
getResNorms(
int numWanted);
293 void setSize(
int blockSize,
int maxSubDim);
332 void expandSearchSpace();
335 void applyOperators();
338 void updateProjections();
341 void solveProjectedEigenproblem();
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();
377 int d_maxSubspaceDim;
380 std::string d_initType;
382 bool d_relativeConvergence;
391 std::vector< Value<ScalarType> > d_eigenvalues;
395 std::vector<MagnitudeType> d_resNorms;
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>
449 class GeneralizedDavidson<std::complex<MagnitudeType>,MV,OP>
453 typedef std::complex<MagnitudeType> ScalarType;
455 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
456 const RCP<SortManager<MagnitudeType> > &sortman,
457 const RCP<OutputManager<ScalarType> > &outputman,
458 const RCP<StatusTest<ScalarType,MV,OP> > &tester,
459 const RCP<OrthoManager<ScalarType,MV> > &orthoman,
463 MagnitudeType::this_class_is_missing_a_specialization();
474 template <
class ScalarType,
class MV,
class OP>
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");
520 "Anasazi::GeneralizedDavidson: Test Space must be V, AV, or BV" );
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);
536 std::invalid_argument,
"Maximum Subspace Dimension must be strictly greater than NEV");
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;
664 d_V = MVT::Clone(*initVec, msd);
665 d_AV = MVT::Clone(*initVec, msd);
676 d_BV = MVT::Clone(*initVec, 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 );
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 );
780 if( !reset_V && state.
VAV != Teuchos::null && state.
VAV->numRows() >= d_curDim && state.
VAV->numCols() >= d_curDim )
782 if( state.
VAV != d_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 );
828 if( !reset_V && state.
VBV != Teuchos::null && state.
VBV->numRows() >= d_curDim && state.
VBV->numCols() >= d_curDim )
830 if( state.
VBV != d_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>
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);
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;
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);
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;
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.");
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).");
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);
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 );
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 )
1108 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidson: Error in V^T V == I: " 1109 << d_orthoMan->orthonormError( *V_all ) << std::endl;
1113 "Anasazi::GeneralizedDavidson::ExpandSearchSpace(): Orthonormalization of new vectors failed" );
1120 template <
class ScalarType,
class MV,
class OP>
1121 void GeneralizedDavidson<ScalarType,MV,OP>::applyOperators()
1124 std::vector<int> newIndices(d_expansionSize);
1125 for(
int i=0; i<d_expansionSize; ++i )
1126 newIndices[i] = d_curDim+i;
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>
1148 void GeneralizedDavidson<ScalarType,MV,OP>::updateProjections()
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 )
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 );
1183 RCP<const MV> AV_current = MVT::CloneView(*d_AV, curIndices);
1187 MVT::MvTransMv( ST::one(), *W_new, *AV_current, VAV_lastrow );
1191 MVT::MvTransMv( ST::one(), *W_all, *AV_new, VAV_lastcol );
1197 RCP<const MV> BV_current = MVT::CloneView(*d_BV, curIndices);
1201 MVT::MvTransMv( ST::one(), *W_new, *BV_current, VBV_lastrow );
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>
1219 void GeneralizedDavidson<ScalarType,MV,OP>::solveProjectedEigenproblem()
1226 d_S->assign(*d_VAV);
1227 d_T->assign(*d_VBV);
1230 char leftVecs =
'V';
1231 char rightVecs =
'V';
1232 char sortVals =
'N';
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;
1261 d_S->assign(*d_VAV);
1266 int work_size = 3*d_curDim;
1267 std::vector<ScalarType> work(work_size);
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;
1293 template <
class ScalarType,
class MV,
class OP>
1294 void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzIndex()
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>
1327 void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzVectors()
1329 if( d_ritzVectorsValid )
1336 std::vector<int> checkedIndices(d_residualSize);
1337 for(
int ii=0; ii<d_residualSize; ++ii )
1338 checkedIndices[ii] = ii;
1342 computeProjectedEigenvectors( X );
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 );
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>
1386 void GeneralizedDavidson<ScalarType,MV,OP>::sortValues( std::vector<MagnitudeType> &realParts,
1387 std::vector<MagnitudeType> &imagParts,
1388 std::vector<int> &permVec,
1394 "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1396 "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1400 d_sortMan->sort( realParts, imagParts, rcpPermVec, N );
1402 d_ritzIndexValid =
false;
1403 d_ritzVectorsValid =
false;
1417 template <
class ScalarType,
class MV,
class OP>
1418 void GeneralizedDavidson<ScalarType,MV,OP>::computeProjectedEigenvectors(
1428 char whichVecs =
'R';
1430 int work_size = 6*d_maxSubspaceDim;
1431 std::vector<ScalarType> work(work_size,ST::zero());
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;
1448 char whichVecs =
'R';
1451 std::vector<ScalarType> work(3*N);
1457 lapack.
TREVC( whichVecs, howMany, &sel, N, S.values(), S.stride(), VL.values(), VL.stride(),
1460 std::stringstream ss;
1461 ss <<
"Anasazi::GeneralizedDavidson: TREVC returned error code " << info << std::endl;
1469 template <
class ScalarType,
class MV,
class OP>
1470 void GeneralizedDavidson<ScalarType,MV,OP>::scaleEigenvectors(
1479 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1481 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1483 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xalpha does not match X");
1485 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xalpha does not match X");
1487 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xbeta does not match X");
1489 "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xbeta does not match X");
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];
1523 std::stringstream astream;
1524 astream <<
"GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1529 std::stringstream bstream;
1530 bstream <<
"GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1537 template <
class ScalarType,
class MV,
class OP>
1538 void GeneralizedDavidson<ScalarType,MV,OP>::computeResidual()
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;
1562 computeProjectedEigenvectors( X );
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);
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)
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.
ScalarType LAPY2(const ScalarType x, const ScalarType y) const
std::vector< MagnitudeType > getRitzRes2Norms()
Get the current Ritz residual norms (2-norm)
ScalarType * values() const
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get status test.
Teuchos::ScalarTraits< ScalarType >::magnitudeType imagpart
The imaginary component of the eigenvalue.
This class defines the interface required by an eigensolver and status test class to compute solution...
int getBlockSize() const
Get block size.
void TGEVC(const char SIDE, const char HOWMNY, const OrdinalType *SELECT, const OrdinalType n, ScalarType *S, const OrdinalType lds, ScalarType *P, const OrdinalType ldp, ScalarType *VL, const OrdinalType ldvl, ScalarType *VR, const OrdinalType ldvr, const OrdinalType mm, OrdinalType *M, ScalarType *WORK, OrdinalType *info) const
Solves eigenvalue problem using generalized Davidson method.
std::vector< Value< ScalarType > > getRitzValues()
Get the current Ritz values.
T & get(const std::string &name, T def_value)
void TRSEN(const char JOB, const char COMPQ, const OrdinalType *SELECT, const OrdinalType n, ScalarType *T, const OrdinalType ldt, ScalarType *Q, const OrdinalType ldq, MagnitudeType *WR, MagnitudeType *WI, OrdinalType *M, ScalarType *S, MagnitudeType *SEP, ScalarType *WORK, const OrdinalType lwork, OrdinalType *IWORK, const OrdinalType liwork, OrdinalType *info) const
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > T
Right quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
void setBlockSize(int blockSize)
Set block size.
int curDim
The current subspace dimension.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType n, const ScalarType *x, const OrdinalType incx) const
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.
void GEES(const char JOBVS, const char SORT, OrdinalType(*ptr2func)(ScalarType *, ScalarType *), const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *sdim, ScalarType *WR, ScalarType *WI, ScalarType *VS, const OrdinalType ldvs, ScalarType *WORK, const OrdinalType lwork, OrdinalType *BWORK, OrdinalType *info) const
Abstract class definition for Anasazi Output Managers.
std::vector< Value< ScalarType > > eVals
Vector of generalized eigenvalues.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setSize(int blockSize, int maxSubDim)
Set problem size.
OrdinalType numRows() const
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...
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxilliary vectors.
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< Teuchos::SerialDenseMatrix< int, ScalarType > > Q
Left generalized Schur vectors from QZ decomposition of (VAV,VBV)
Traits class which defines basic operations on multivectors.
void GGES(const char JOBVL, const char JOBVR, const char SORT, OrdinalType(*ptr2func)(ScalarType *, ScalarType *, ScalarType *), const OrdinalType n, ScalarType *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb, OrdinalType *sdim, MagnitudeType *ALPHAR, MagnitudeType *ALPHAI, MagnitudeType *BETA, ScalarType *VL, const OrdinalType ldvl, ScalarType *VR, const OrdinalType ldvr, ScalarType *WORK, const OrdinalType lwork, OrdinalType *BWORK, OrdinalType *info) const
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...
bool isInitialized() const
Query if the solver is in an initialized state.
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.
int getCurSubspaceDim() const
Get current subspace dimension.
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.
void push_back(const value_type &x)
int getNumIters() const
Get number of iterations.
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.
void TGSEN(const OrdinalType ijob, const OrdinalType wantq, const OrdinalType wantz, const OrdinalType *SELECT, const OrdinalType n, ScalarType *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb, MagnitudeType *ALPHAR, MagnitudeType *ALPHAI, MagnitudeType *BETA, ScalarType *Q, const OrdinalType ldq, ScalarType *Z, const OrdinalType ldz, OrdinalType *M, MagnitudeType *PL, MagnitudeType *PR, MagnitudeType *DIF, ScalarType *WORK, const OrdinalType lwork, OrdinalType *IWORK, const OrdinalType liwork, OrdinalType *info) const
void TREVC(const char SIDE, const char HOWMNY, OrdinalType *select, const OrdinalType n, const ScalarType *T, const OrdinalType ldt, ScalarType *VL, const OrdinalType ldvl, ScalarType *VR, const OrdinalType ldvr, const OrdinalType mm, OrdinalType *m, ScalarType *WORK, OrdinalType *info) const
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
OrdinalType stride() const
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.
OrdinalType numCols() const
Common interface of stopping criteria for Anasazi's solvers.
std::vector< int > getBlockIndex() const
Get indices of current block.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
int getMaxSubspaceDim() const
Get maximum subspace dimension.
Declaration and definition of Anasazi::StatusTest.
void resetNumIters()
Reset the number of iterations.