429 typedef SolverUtils<ScalarType,MV,OP> msutils;
431 const int nev = problem_->getNEV();
434 Teuchos::RCP<Teuchos::FancyOStream>
435 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(
Debug)));
436 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
437 *out <<
"Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
442 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp(
new BasicSort<MagnitudeType>(whch_) );
448 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
449 if (globalTest_ == Teuchos::null) {
450 convtest = Teuchos::rcp(
new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
453 convtest = globalTest_;
455 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
456 = Teuchos::rcp(
new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
458 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
460 if (lockingTest_ == Teuchos::null) {
461 locktest = Teuchos::rcp(
new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
464 locktest = lockingTest_;
468 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
469 alltests.push_back(ordertest);
470 if (locktest != Teuchos::null) alltests.push_back(locktest);
471 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
473 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
476 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
477 if ( printer_->isVerbosity(
Debug) ) {
478 outputtest = Teuchos::rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,
Passed+
Failed+
Undefined ) );
481 outputtest = Teuchos::rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,
Passed ) );
486 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
487 if (ortho_==
"SVQB") {
488 ortho = Teuchos::rcp(
new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
489 }
else if (ortho_==
"DGKS") {
490 ortho = Teuchos::rcp(
new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
492 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS",std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
497 Teuchos::ParameterList plist;
498 plist.set(
"Block Size",blockSize_);
499 plist.set(
"Num Blocks",numBlocks_);
503 Teuchos::RCP<BlockDavidson<ScalarType,MV,OP> > bd_solver
504 = Teuchos::rcp(
new BlockDavidson<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,plist) );
506 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
507 if (probauxvecs != Teuchos::null) {
508 bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
515 int curNumLocked = 0;
516 Teuchos::RCP<MV> lockvecs;
522 if (maxLocked_ > 0) {
523 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
525 std::vector<MagnitudeType> lockvals;
573 Teuchos::RCP<MV> workMV;
574 if (inSituRestart_ ==
false) {
576 if (useLocking_==
true || maxRestarts_ > 0) {
577 workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
581 workMV = Teuchos::null;
590 if (maxRestarts_ > 0) {
591 workMV = MVT::Clone(*problem_->getInitVec(),1);
594 workMV = Teuchos::null;
599 const ScalarType ONE = SCT::one();
600 const ScalarType ZERO = SCT::zero();
601 Teuchos::LAPACK<int,ScalarType> lapack;
602 Teuchos::BLAS<int,ScalarType> blas;
605 Eigensolution<ScalarType,MV> sol;
607 problem_->setSolution(sol);
613#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
614 Teuchos::TimeMonitor slvtimer(*_timerSolve);
620 bd_solver->iterate();
627 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
628 throw AnasaziError(
"Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
635 else if (ordertest->getStatus() ==
Passed ) {
646 else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
648#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
649 Teuchos::TimeMonitor restimer(*_timerRestarting);
652 if ( numRestarts >= maxRestarts_ ) {
657 printer_->stream(
IterationDetails) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
659 BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
660 int curdim = state.curDim;
661 int newdim = numRestartBlocks_*blockSize_;
665 std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
666 *out <<
"Ritz values from solver:\n";
667 for (
unsigned int i=0; i<ritzvalues.size(); i++) {
668 *out << ritzvalues[i].realpart <<
" ";
676 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
677 std::vector<MagnitudeType> theta(curdim);
681 std::stringstream os;
682 os <<
"KK before HEEV...\n"
683 << printMat(*state.KK) <<
"\n";
687 int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
688 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
689 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
690 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
691 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
695 std::stringstream os;
696 *out <<
"Ritz values from heev(KK):\n";
697 for (
unsigned int i=0; i<theta.size(); i++) *out << theta[i] <<
" ";
698 os <<
"\nRitz vectors from heev(KK):\n"
699 << printMat(S) <<
"\n";
707 std::vector<int> order(curdim);
708 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
711 msutils::permuteVectors(order,S);
715 std::stringstream os;
716 *out <<
"Ritz values from heev(KK) after sorting:\n";
717 std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out,
" "));
718 os <<
"\nRitz vectors from heev(KK) after sorting:\n"
719 << printMat(S) <<
"\n";
725 Teuchos::SerialDenseMatrix<int,ScalarType> Sr(Teuchos::View,S,curdim,newdim);
728 std::stringstream os;
729 os <<
"Significant primitive Ritz vectors:\n"
730 << printMat(Sr) <<
"\n";
736 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
738 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim),
739 KKold(Teuchos::View,*state.KK,curdim,curdim);
742 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
743 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
744 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
746 teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
747 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
748 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
750 for (
int j=0; j<newdim-1; ++j) {
751 for (
int i=j+1; i<newdim; ++i) {
752 newKK(i,j) = SCT::conjugate(newKK(j,i));
758 std::stringstream os;
760 << printMat(newKK) <<
"\n";
766 BlockDavidsonState<ScalarType,MV> rstate;
767 rstate.curDim = newdim;
768 rstate.KK = Teuchos::rcpFromRef(newKK);
775 rstate.KX = state.KX;
776 rstate.MX = state.MX;
778 rstate.T = Teuchos::rcp(
new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
780 if (inSituRestart_ ==
true) {
782 *out <<
"Beginning in-situ restart...\n";
786 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
790 std::vector<ScalarType> tau(newdim), work(newdim);
792 lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
793 TEUCHOS_TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
794 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
795 if (printer_->isVerbosity(
Debug)) {
796 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,Sr,newdim,newdim);
797 for (
int j=0; j<newdim; j++) {
798 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
799 for (
int i=j+1; i<newdim; i++) {
803 printer_->stream(
Debug) <<
"||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
810 std::vector<int> curind(curdim);
811 for (
int i=0; i<curdim; i++) curind[i] = i;
812 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
813 msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
819 rstate.V = solverbasis;
823 *out <<
"Beginning ex-situ restart...\n";
826 std::vector<int> curind(curdim), newind(newdim);
827 for (
int i=0; i<curdim; i++) curind[i] = i;
828 for (
int i=0; i<newdim; i++) newind[i] = i;
829 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
830 Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*workMV ,newind);
832 MVT::MvTimesMatAddMv(ONE,*oldV,Sr,ZERO,*newV);
840 bd_solver->initialize(rstate);
847 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
849#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
850 Teuchos::TimeMonitor lcktimer(*_timerLocking);
855 BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
856 const int curdim = state.curDim;
860 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
861 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
862 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (
int)locktest->whichVecs().size(),std::logic_error,
863 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
865 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
866 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
869 std::vector<int> tmp_vector_int;
870 if (curNumLocked + locktest->howMany() > maxLocked_) {
872 tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
875 tmp_vector_int = locktest->whichVecs();
877 const std::vector<int> lockind(tmp_vector_int);
878 const int numNewLocked = lockind.size();
882 const int numUnlocked = curdim-numNewLocked;
883 tmp_vector_int.resize(curdim);
884 for (
int i=0; i<curdim; i++) tmp_vector_int[i] = i;
885 const std::vector<int> curind(tmp_vector_int);
886 tmp_vector_int.resize(numUnlocked);
887 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
888 const std::vector<int> unlockind(tmp_vector_int);
889 tmp_vector_int.clear();
893 if (printer_->isVerbosity(
Debug)) {
894 printer_->print(
Debug,
"Locking vectors: ");
895 for (
unsigned int i=0; i<lockind.size(); i++) {printer_->stream(
Debug) <<
" " << lockind[i];}
896 printer_->print(
Debug,
"\n");
897 printer_->print(
Debug,
"Not locking vectors: ");
898 for (
unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(
Debug) <<
" " << unlockind[i];}
899 printer_->print(
Debug,
"\n");
910 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
911 std::vector<MagnitudeType> theta(curdim);
914 int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
915 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
916 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
917 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
918 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
921 std::vector<int> order(curdim);
922 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
925 msutils::permuteVectors(order,S);
931 Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
932 for (
int i=0; i<numUnlocked; i++) {
933 blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
946 Teuchos::RCP<MV> defV, augV;
947 if (inSituRestart_ ==
true) {
950 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
954 Teuchos::SerialDenseMatrix<int,ScalarType> copySu(Su);
955 std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
957 lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
958 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
959 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
960 if (printer_->isVerbosity(
Debug)) {
961 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
962 for (
int j=0; j<numUnlocked; j++) {
963 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
964 for (
int i=j+1; i<numUnlocked; i++) {
968 printer_->stream(
Debug) <<
"||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
975 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
976 msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
978 std::vector<int> defind(numUnlocked), augind(numNewLocked);
979 for (
int i=0; i<numUnlocked ; i++) defind[i] = i;
980 for (
int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
981 defV = MVT::CloneViewNonConst(*solverbasis,defind);
982 augV = MVT::CloneViewNonConst(*solverbasis,augind);
986 std::vector<int> defind(numUnlocked), augind(numNewLocked);
987 for (
int i=0; i<numUnlocked ; i++) defind[i] = i;
988 for (
int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
989 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
990 defV = MVT::CloneViewNonConst(*workMV,defind);
991 augV = MVT::CloneViewNonConst(*workMV,augind);
993 MVT::MvTimesMatAddMv(ONE,*oldV,Su,ZERO,*defV);
1007 Teuchos::RCP<const MV> curlocked, newLocked;
1008 Teuchos::RCP<MV> augTmp;
1011 if (curNumLocked > 0) {
1012 std::vector<int> curlockind(curNumLocked);
1013 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
1014 curlocked = MVT::CloneView(*lockvecs,curlockind);
1017 curlocked = Teuchos::null;
1020 std::vector<int> augtmpind(numNewLocked);
1021 for (
int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
1022 augTmp = MVT::CloneViewNonConst(*lockvecs,augtmpind);
1024 newLocked = MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
1030 MVT::MvRandom(*augV);
1035 Teuchos::Array<Teuchos::RCP<const MV> > against;
1036 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1037 if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
1038 if (curlocked != Teuchos::null) against.push_back(curlocked);
1039 against.push_back(newLocked);
1040 against.push_back(defV);
1041 if (problem_->getM() != Teuchos::null) {
1042 OPT::Apply(*problem_->getM(),*augV,*augTmp);
1044 ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
1055 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
1057 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked),
1058 KKold(Teuchos::View,*state.KK,curdim,curdim),
1059 KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
1062 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
1063 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1064 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1066 teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
1067 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1068 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1073 OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
1074 Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
1075 KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
1076 MVT::MvTransMv(ONE,*defV,*augTmp,KK12);
1077 MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
1081 defV = Teuchos::null;
1082 augV = Teuchos::null;
1085 for (
int j=0; j<curdim; ++j) {
1086 for (
int i=j+1; i<curdim; ++i) {
1087 newKK(i,j) = SCT::conjugate(newKK(j,i));
1093 augTmp = Teuchos::null;
1095 std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
1096 for (
int i=0; i<numNewLocked; i++) {
1097 lockvals.push_back(allvals[lockind[i]].realpart);
1100 std::vector<int> indlock(numNewLocked);
1101 for (
int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
1102 MVT::SetBlock(*newLocked,indlock,*lockvecs);
1103 newLocked = Teuchos::null;
1105 curNumLocked += numNewLocked;
1106 std::vector<int> curlockind(curNumLocked);
1107 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
1108 curlocked = MVT::CloneView(*lockvecs,curlockind);
1114 ordertest->setAuxVals(lockvals);
1116 Teuchos::Array< Teuchos::RCP<const MV> > aux;
1117 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
1118 aux.push_back(curlocked);
1119 bd_solver->setAuxVecs(aux);
1121 if (curNumLocked == maxLocked_) {
1123 combotest->removeTest(locktest);
1129 BlockDavidsonState<ScalarType,MV> rstate;
1130 rstate.curDim = curdim;
1131 if (inSituRestart_) {
1139 rstate.KK = Teuchos::rcpFromRef(newKK);
1142 bd_solver->initialize(rstate);
1151 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
1156 <<
"Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
1157 << err.what() << std::endl
1158 <<
"Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1164 workMV = Teuchos::null;
1166 sol.numVecs = ordertest->howMany();
1167 if (sol.numVecs > 0) {
1168 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
1169 sol.Espace = sol.Evecs;
1170 sol.Evals.resize(sol.numVecs);
1171 std::vector<MagnitudeType> vals(sol.numVecs);
1174 std::vector<int> which = ordertest->whichVecs();
1178 std::vector<int> inlocked(0), insolver(0);
1179 for (
unsigned int i=0; i<which.size(); i++) {
1180 if (which[i] >= 0) {
1181 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
1182 insolver.push_back(which[i]);
1186 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
1187 inlocked.push_back(which[i] + curNumLocked);
1191 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (
unsigned int)sol.numVecs,std::logic_error,
"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
1194 if (insolver.size() > 0) {
1196 int lclnum = insolver.size();
1197 std::vector<int> tosol(lclnum);
1198 for (
int i=0; i<lclnum; i++) tosol[i] = i;
1199 Teuchos::RCP<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
1200 MVT::SetBlock(*v,tosol,*sol.Evecs);
1202 std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
1203 for (
unsigned int i=0; i<insolver.size(); i++) {
1204 vals[i] = fromsolver[insolver[i]].realpart;
1209 if (inlocked.size() > 0) {
1210 int solnum = insolver.size();
1212 int lclnum = inlocked.size();
1213 std::vector<int> tosol(lclnum);
1214 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1215 Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1216 MVT::SetBlock(*v,tosol,*sol.Evecs);
1218 for (
unsigned int i=0; i<inlocked.size(); i++) {
1219 vals[i+solnum] = lockvals[inlocked[i]];
1225 std::vector<int> order(sol.numVecs);
1226 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1228 for (
int i=0; i<sol.numVecs; i++) {
1229 sol.Evals[i].realpart = vals[i];
1230 sol.Evals[i].imagpart = MT::zero();
1233 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1237 sol.index.resize(sol.numVecs,0);
1242 bd_solver->currentStatus(printer_->stream(
FinalSummary));
1245#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1247 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails ) );
1251 problem_->setSolution(sol);
1252 printer_->stream(
Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1255 numIters_ = bd_solver->getNumIters();
1257 if (sol.numVecs < nev) {