114 typedef MultiVecTraits<ScalarType,MV> MVT;
115 typedef OperatorTraits<ScalarType,MV,OP> OPT;
116 typedef Teuchos::ScalarTraits<ScalarType> SCT;
117 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
118 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
171 Teuchos::ParameterList &pl );
198 return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
252 RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
257 int blockSize_, numBlocks_, numRestartBlocks_;
260 RCP<OutputManager<ScalarType> > printer_;
263 MagnitudeType convTol_;
268 MagnitudeType lockTol_;
269 int maxLocked_, lockQuorum_;
270 bool useLocking_, relLockTol_;
274 enum WhenToShiftType whenToShift_;
275 MagnitudeType traceThresh_, shiftTol_;
276 enum HowToShiftType howToShift_;
277 bool useMultipleShifts_, relShiftTol_, considerClusters_;
278 std::string shiftNorm_;
282 std::string ortho_, which_;
283 enum SaddleSolType saddleSolType_;
284 bool projectAllVecs_, projectLockedVecs_, computeAllRes_, useRHSR_, useHarmonic_, noSort_;
285 MagnitudeType alpha_;
288 RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
291 RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
292 RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
293 RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
296 void copyPartOfState(
const TraceMinBaseState<ScalarType,MV>& oldState, TraceMinBaseState<ScalarType,MV>& newState,
const std::vector<int> indToCopy)
const;
298 void setParameters(Teuchos::ParameterList &pl)
const;
300 void printParameters(std::ostream &os)
const;
302 virtual RCP< TraceMinBase<ScalarType,MV,OP> > createSolver(
303 const RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
304 const RCP<StatusTest<ScalarType,MV,OP> > &outputtest,
305 const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
306 Teuchos::ParameterList &plist
309 virtual bool needToRestart(
const RCP< TraceMinBase<ScalarType,MV,OP> > solver) =0;
311 virtual bool performRestart(
int &numRestarts, RCP< TraceMinBase<ScalarType,MV,OP> > solver) =0;
320 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
321 Teuchos::ParameterList &pl ) :
323#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
324 , _timerSolve(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBaseSolMgr::solve()")),
325 _timerRestarting(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBaseSolMgr restarting")),
326 _timerLocking(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBaseSolMgr locking"))
329 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
330 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
331 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument,
"Problem not symmetric.");
332 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
355 int osProc = pl.get(
"Output Processor", 0);
358 Teuchos::RCP<Teuchos::FancyOStream> osp;
360 if (pl.isParameter(
"Output Stream")) {
361 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
368 if (pl.isParameter(
"Verbosity")) {
369 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
370 verbosity = pl.get(
"Verbosity", verbosity);
372 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
375 printer_ = rcp(
new OutputManager<ScalarType>(verbosity,osp) );
381 convTol_ = pl.get(
"Convergence Tolerance",MT::prec());
382 TEUCHOS_TEST_FOR_EXCEPTION(convTol_ < 0, std::invalid_argument,
383 "Anasazi::TraceMinBaseSolMgr: \"Convergence Tolerance\" must be nonnegative.");
385 relConvTol_ = pl.get(
"Relative Convergence Tolerance",
true);
386 strtmp = pl.get(
"Convergence Norm",std::string(
"2"));
388 convNorm_ = RES_2NORM;
390 else if (strtmp ==
"M") {
391 convNorm_ = RES_ORTH;
394 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
395 "Anasazi::TraceMinBaseSolMgr: Invalid Convergence Norm.");
400 useLocking_ = pl.get(
"Use Locking",
true);
401 relLockTol_ = pl.get(
"Relative Locking Tolerance",
true);
402 lockTol_ = pl.get(
"Locking Tolerance",convTol_/10);
404 TEUCHOS_TEST_FOR_EXCEPTION(relConvTol_ != relLockTol_, std::invalid_argument,
405 "Anasazi::TraceMinBaseSolMgr: \"Relative Convergence Tolerance\" and \"Relative Locking Tolerance\" have different values. If you set one, you should always set the other.");
407 TEUCHOS_TEST_FOR_EXCEPTION(lockTol_ < 0, std::invalid_argument,
408 "Anasazi::TraceMinBaseSolMgr: \"Locking Tolerance\" must be nonnegative.");
410 strtmp = pl.get(
"Locking Norm",std::string(
"2"));
412 lockNorm_ = RES_2NORM;
414 else if (strtmp ==
"M") {
415 lockNorm_ = RES_ORTH;
418 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
419 "Anasazi::TraceMinBaseSolMgr: Invalid Locking Norm.");
424 maxLocked_ = pl.get(
"Max Locked",problem_->getNEV());
425 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ <= 0, std::invalid_argument,
426 "Anasazi::TraceMinBaseSolMgr: \"Max Locked\" must be strictly positive.");
433 lockQuorum_ = pl.get(
"Locking Quorum",1);
434 TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0, std::invalid_argument,
435 "Anasazi::TraceMinBaseSolMgr: \"Locking Quorum\" must be strictly positive.");
442 strtmp = pl.get(
"When To Shift",
"Always");
444 if(strtmp ==
"Never")
445 whenToShift_ = NEVER_SHIFT;
446 else if(strtmp ==
"After Trace Levels")
447 whenToShift_ = SHIFT_WHEN_TRACE_LEVELS;
448 else if(strtmp ==
"Residual Becomes Small")
449 whenToShift_ = SHIFT_WHEN_RESID_SMALL;
450 else if(strtmp ==
"Always")
451 whenToShift_ = ALWAYS_SHIFT;
453 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
454 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"When To Shift\"; valid options are \"Never\", \"After Trace Levels\", \"Residual Becomes Small\", \"Always\".");
458 traceThresh_ = pl.get(
"Trace Threshold", 0.02);
460 TEUCHOS_TEST_FOR_EXCEPTION(traceThresh_ < 0, std::invalid_argument,
461 "Anasazi::TraceMinBaseSolMgr: \"Trace Threshold\" must be nonnegative.");
464 shiftTol_ = pl.get(
"Shift Tolerance", 0.1);
466 TEUCHOS_TEST_FOR_EXCEPTION(shiftTol_ < 0, std::invalid_argument,
467 "Anasazi::TraceMinBaseSolMgr: \"Shift Tolerance\" must be nonnegative.");
470 relShiftTol_ = pl.get(
"Relative Shift Tolerance",
true);
473 shiftNorm_ = pl.get(
"Shift Norm",
"2");
475 TEUCHOS_TEST_FOR_EXCEPTION(shiftNorm_ !=
"2" && shiftNorm_ !=
"M", std::invalid_argument,
476 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Shift Norm\"; valid options are \"2\" and \"M\".");
478 noSort_ = pl.get(
"No Sorting",
false);
481 strtmp = pl.get(
"How To Choose Shift",
"Adjusted Ritz Values");
483 if(strtmp ==
"Largest Converged")
484 howToShift_ = LARGEST_CONVERGED_SHIFT;
485 else if(strtmp ==
"Adjusted Ritz Values")
486 howToShift_ = ADJUSTED_RITZ_SHIFT;
487 else if(strtmp ==
"Ritz Values")
488 howToShift_ = RITZ_VALUES_SHIFT;
489 else if(strtmp ==
"Experimental Shift")
490 howToShift_ = EXPERIMENTAL_SHIFT;
492 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
493 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"How To Choose Shift\"; valid options are \"Largest Converged\", \"Adjusted Ritz Values\", \"Ritz Values\".");
496 considerClusters_ = pl.get(
"Consider Clusters",
true);
499 useMultipleShifts_ = pl.get(
"Use Multiple Shifts",
true);
505 ortho_ = pl.get(
"Orthogonalization",
"SVQB");
506 TEUCHOS_TEST_FOR_EXCEPTION(ortho_ !=
"DGKS" && ortho_ !=
"SVQB" && ortho_ !=
"ICGS", std::invalid_argument,
507 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Orthogonalization\"; valid options are \"DGKS\", \"SVQB\", \"ICGS\".");
509 strtmp = pl.get(
"Saddle Solver Type",
"Projected Krylov");
511 if(strtmp ==
"Projected Krylov")
512 saddleSolType_ = PROJECTED_KRYLOV_SOLVER;
513 else if(strtmp ==
"Schur Complement")
514 saddleSolType_ = SCHUR_COMPLEMENT_SOLVER;
515 else if(strtmp ==
"Block Diagonal Preconditioned Minres")
516 saddleSolType_ = BD_PREC_MINRES;
517 else if(strtmp ==
"HSS Preconditioned Gmres")
518 saddleSolType_ = HSS_PREC_GMRES;
520 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
521 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Saddle Solver Type\"; valid options are \"Projected Krylov\", \"Schur Complement\", and \"Block Diagonal Preconditioned Minres\".");
523 projectAllVecs_ = pl.get(
"Project All Vectors",
true);
524 projectLockedVecs_ = pl.get(
"Project Locked Vectors",
true);
525 computeAllRes_ = pl.get(
"Compute All Residuals",
true);
526 useRHSR_ = pl.get(
"Use Residual as RHS",
false);
527 alpha_ = pl.get(
"HSS: alpha", 1.0);
529 TEUCHOS_TEST_FOR_EXCEPTION(projectLockedVecs_ && ! projectAllVecs_, std::invalid_argument,
530 "Anasazi::TraceMinBaseSolMgr: If you want to project out the locked vectors, you should really project out ALL the vectors of X.");
533 maxKrylovIter_ = pl.get(
"Maximum Krylov Iterations", 200);
534 TEUCHOS_TEST_FOR_EXCEPTION(maxKrylovIter_ < 1, std::invalid_argument,
535 "Anasazi::TraceMinBaseSolMgr: \"Maximum Krylov Iterations\" must be greater than 0.");
538 which_ = pl.get(
"Which",
"SM");
539 TEUCHOS_TEST_FOR_EXCEPTION(which_ !=
"SM" && which_ !=
"LM", std::invalid_argument,
540 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Which\"; valid options are \"SM\" and \"LM\".");
544 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getOperator() == Teuchos::null && whenToShift_ != NEVER_SHIFT, std::invalid_argument,
545 "Anasazi::TraceMinBaseSolMgr: It is an exceptionally bad idea to use Ritz shifts when finding the largest eigenpairs of a standard eigenvalue problem. If we don't use Ritz shifts, it may take extra iterations to converge, but we NEVER have to solve a single linear system. Using Ritz shifts forces us to solve systems of the form (I + sigma A)x=f, and it probably doesn't benefit us enough to outweigh the extra cost. We may add support for this feature in the future, but for now, please set \"When To Shift\" to \"Never\".");
547#ifdef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
550 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER && useMultipleShifts_, std::invalid_argument,
551 "Anasazi::TraceMinBaseSolMgr: When you use the projected Krylov solver with preconditioning, the preconditioner must be projected as well. In theory, if the preconditioner is SPD, the projected preconditioner will be SPSD, but in practice, it can have small negative eigenvalues, presumably due to machine arithmetic. This means we can't use TraceMin's built-in MINRES, and we are forced to use Belos for now. When you use multiple Ritz shifts, you are essentially using a different operator to solve each linear system. Belos can't handle this right now, but we're working on a solution. For now, please set \"Use Multiple Shifts\" to false.");
558 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER, std::invalid_argument,
559 "Anasazi::TraceMinBaseSolMgr: When you use the projected Krylov solver with preconditioning, the preconditioner must be projected as well. In theory, if the preconditioner is SPD, the projected preconditioner will be SPSD, but in practice, it can have small negative eigenvalues, presumably due to machine arithmetic. This means we can't use TraceMin's built-in MINRES, and we are forced to use Belos for now. You didn't install Belos. You have three options to correct this problem:\n1. Reinstall Trilinos with Belos enabled\n2. Don't use a preconditioner\n3. Choose a different method for solving the saddle-point problem (Recommended)");
575 typedef SolverUtils<ScalarType,MV,OP> msutils;
577 const int nev = problem_->getNEV();
580 RCP<Teuchos::FancyOStream>
581 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(
Debug)));
582 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
583 *out <<
"Entering Anasazi::TraceMinBaseSolMgr::solve()\n";
588 RCP<BasicSort<MagnitudeType> > sorter = rcp(
new BasicSort<MagnitudeType>(
"SM") );
595 RCP<const OP> swapHelper = problem_->getOperator();
596 problem_->setOperator(problem_->getM());
597 problem_->setM(swapHelper);
598 problem_->setProblem();
605 RCP<StatusTest<ScalarType,MV,OP> > convtest;
606 if (globalTest_ == Teuchos::null) {
608 convtest = rcp(
new StatusTestResNorm<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_) );
610 convtest = rcp(
new StatusTestSpecTrans<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_,
true,problem_->getOperator()) );
613 convtest = globalTest_;
615 RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
616 = rcp(
new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
618 RCP<StatusTest<ScalarType,MV,OP> > locktest;
620 if (lockingTest_ == Teuchos::null) {
622 locktest = rcp(
new StatusTestResNorm<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_) );
624 locktest = rcp(
new StatusTestSpecTrans<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_,
true,problem_->getOperator()) );
627 locktest = lockingTest_;
631 Teuchos::Array<RCP<StatusTest<ScalarType,MV,OP> > > alltests;
632 alltests.push_back(ordertest);
633 if (locktest != Teuchos::null) alltests.push_back(locktest);
634 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
636 RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
639 RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
640 if ( printer_->isVerbosity(
Debug) ) {
641 outputtest = rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,
Passed+
Failed+
Undefined ) );
644 outputtest = rcp(
new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,
Passed ) );
649 RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
650 if (ortho_==
"SVQB") {
651 ortho = rcp(
new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
652 }
else if (ortho_==
"DGKS") {
653 ortho = rcp(
new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
654 }
else if (ortho_==
"ICGS") {
655 ortho = rcp(
new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
657 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): Invalid orthogonalization type.");
662 Teuchos::ParameterList plist;
663 setParameters(plist);
667 RCP<TraceMinBase<ScalarType,MV,OP> > tm_solver
668 = createSolver(sorter,outputtest,ortho,plist);
670 RCP< const MV > probauxvecs = problem_->getAuxVecs();
671 if (probauxvecs != Teuchos::null) {
672 tm_solver->setAuxVecs( Teuchos::tuple< RCP<const MV> >(probauxvecs) );
679 int curNumLocked = 0;
686 if (maxLocked_ > 0) {
687 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
689 std::vector<MagnitudeType> lockvals;
734 const ScalarType ONE = SCT::one();
735 const ScalarType ZERO = SCT::zero();
738 Eigensolution<ScalarType,MV> sol;
740 problem_->setSolution(sol);
746#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
747 Teuchos::TimeMonitor slvtimer(*_timerSolve);
753 tm_solver->iterate();
759 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
760 throw AnasaziError(
"Anasazi::TraceMinBaseSolMgr::solve(): User-specified debug status test returned Passed.");
767 else if (ordertest->getStatus() ==
Passed ) {
778 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
780#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
781 Teuchos::TimeMonitor lcktimer(*_timerLocking);
786 TraceMinBaseState<ScalarType,MV> state = tm_solver->getState();
787 const int curdim = state.curDim;
791 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
792 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() non-positive.");
793 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (
int)locktest->whichVecs().size(),std::logic_error,
794 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
796 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
797 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: locking not deactivated.");
800 std::vector<int> tmp_vector_int;
801 if (curNumLocked + locktest->howMany() > maxLocked_) {
803 for(
int i=0; i<maxLocked_-curNumLocked; i++)
804 tmp_vector_int.push_back(locktest->whichVecs()[i]);
808 tmp_vector_int = locktest->whichVecs();
811 const std::vector<int> lockind(tmp_vector_int);
812 const int numNewLocked = lockind.size();
816 const int numUnlocked = curdim-numNewLocked;
817 tmp_vector_int.resize(curdim);
818 for (
int i=0; i<curdim; i++) tmp_vector_int[i] = i;
819 const std::vector<int> curind(tmp_vector_int);
820 tmp_vector_int.resize(numUnlocked);
821 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
822 const std::vector<int> unlockind(tmp_vector_int);
823 tmp_vector_int.clear();
827 if (printer_->isVerbosity(
Debug)) {
828 printer_->print(
Debug,
"Locking vectors: ");
829 for (
unsigned int i=0; i<lockind.size(); i++) {printer_->stream(
Debug) <<
" " << lockind[i];}
830 printer_->print(
Debug,
"\n");
831 printer_->print(
Debug,
"Not locking vectors: ");
832 for (
unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(
Debug) <<
" " << unlockind[i];}
833 printer_->print(
Debug,
"\n");
837 std::vector<Value<ScalarType> > allvals = tm_solver->getRitzValues();
838 for(
unsigned int i=0; i<allvals.size(); i++)
839 printer_->stream(
Debug) <<
"Ritz value[" << i <<
"] = " << allvals[i].realpart << std::endl;
840 for (
int i=0; i<numNewLocked; i++) {
841 lockvals.push_back(allvals[lockind[i]].realpart);
845 RCP<const MV> newLocked = MVT::CloneView(*tm_solver->getRitzVectors(),lockind);
846 std::vector<int> indlock(numNewLocked);
847 for (
int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
850 RCP<MV> tempMV = MVT::CloneCopy(*newLocked);
851 ortho->normalizeMat(*tempMV);
852 MVT::SetBlock(*tempMV,indlock,*lockvecs);
856 MVT::SetBlock(*newLocked,indlock,*lockvecs);
864 for(
unsigned int aliciaInd=0; aliciaInd<lockvals.size(); aliciaInd++)
866 lockvals[aliciaInd] = 0.0;
869 ordertest->setAuxVals(lockvals);
873 curNumLocked += numNewLocked;
875 if(ordertest->getStatus() ==
Passed)
break;
877 std::vector<int> curlockind(curNumLocked);
878 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
879 RCP<const MV> curlocked = MVT::CloneView(*lockvecs,curlockind);
881 Teuchos::Array< RCP<const MV> > aux;
882 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
883 aux.push_back(curlocked);
884 tm_solver->setAuxVecs(aux);
887 printer_->stream(
Debug) <<
"curNumLocked: " << curNumLocked << std::endl;
888 printer_->stream(
Debug) <<
"maxLocked: " << maxLocked_ << std::endl;
889 if (curNumLocked == maxLocked_) {
891 combotest->removeTest(locktest);
892 locktest = Teuchos::null;
893 printer_->stream(
Debug) <<
"Removed locking test\n";
896 int newdim = numRestartBlocks_*blockSize_;
897 TraceMinBaseState<ScalarType,MV> newstate;
898 if(newdim <= numUnlocked)
902 std::vector<int> desiredSubscripts(newdim);
903 for(
int i=0; i<newdim; i++)
905 desiredSubscripts[i] = unlockind[i];
906 printer_->stream(
Debug) <<
"H desiredSubscripts[" << i <<
"] = " << desiredSubscripts[i] << std::endl;
908 newstate.V = MVT::CloneView(*tm_solver->getRitzVectors(),desiredSubscripts);
909 newstate.curDim = newdim;
913 std::vector<int> desiredSubscripts(newdim);
914 for(
int i=0; i<newdim; i++)
916 desiredSubscripts[i] = unlockind[i];
917 printer_->stream(
Debug) <<
"desiredSubscripts[" << i <<
"] = " << desiredSubscripts[i] << std::endl;
920 copyPartOfState(state, newstate, desiredSubscripts);
928 int nrandom = newdim-numUnlocked;
930 RCP<const MV> helperMV;
931 RCP<MV> totalV = MVT::Clone(*tm_solver->getRitzVectors(),newdim);
934 tmp_vector_int.resize(numUnlocked);
935 for(
int i=0; i<numUnlocked; i++) tmp_vector_int[i] = i;
936 RCP<MV> oldV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
940 helperMV = MVT::CloneView(*tm_solver->getRitzVectors(),unlockind);
942 helperMV = MVT::CloneView(*state.V,unlockind);
943 MVT::Assign(*helperMV,*oldV);
946 tmp_vector_int.resize(nrandom);
947 for(
int i=0; i<nrandom; i++) tmp_vector_int[i] = i+numUnlocked;
948 RCP<MV> newV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
951 MVT::MvRandom(*newV);
954 newstate.curDim = newdim;
957 RCP< std::vector<ScalarType> > helperRS = rcp(
new std::vector<ScalarType>(blockSize_) );
958 for(
unsigned int i=0; i<(
unsigned int)blockSize_; i++)
960 if(i < unlockind.size() && unlockind[i] < blockSize_)
961 (*helperRS)[i] = (*state.ritzShifts)[unlockind[i]];
963 (*helperRS)[i] = ZERO;
965 newstate.ritzShifts = helperRS;
969 newstate.largestSafeShift = std::abs(lockvals[0]);
970 for(
size_t i=0; i<lockvals.size(); i++)
971 newstate.largestSafeShift = std::max(std::abs(lockvals[i]), newstate.largestSafeShift);
976 newstate.NEV = state.NEV - numNewLocked;
977 tm_solver->initialize(newstate);
984 else if ( needToRestart(tm_solver) ) {
986 if(performRestart(numRestarts, tm_solver) ==
false)
996 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): Invalid return from tm_solver::iterate().");
1001 <<
"Anasazi::TraceMinBaseSolMgr::solve() caught unexpected exception from Anasazi::TraceMinBase::iterate() at iteration " << tm_solver->getNumIters() << std::endl
1002 << err.what() << std::endl
1003 <<
"Anasazi::TraceMinBaseSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1008 sol.numVecs = ordertest->howMany();
1009 if (sol.numVecs > 0) {
1010 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
1011 sol.Espace = sol.Evecs;
1012 sol.Evals.resize(sol.numVecs);
1013 std::vector<MagnitudeType> vals(sol.numVecs);
1016 std::vector<int> which = ordertest->whichVecs();
1020 std::vector<int> inlocked(0), insolver(0);
1021 for (
unsigned int i=0; i<which.size(); i++) {
1022 if (which[i] >= 0) {
1023 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): positive indexing mistake from ordertest.");
1024 insolver.push_back(which[i]);
1028 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): negative indexing mistake from ordertest.");
1029 inlocked.push_back(which[i] + curNumLocked);
1033 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (
unsigned int)sol.numVecs,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): indexing mistake.");
1036 if (insolver.size() > 0) {
1038 int lclnum = insolver.size();
1039 std::vector<int> tosol(lclnum);
1040 for (
int i=0; i<lclnum; i++) tosol[i] = i;
1041 RCP<const MV> v = MVT::CloneView(*tm_solver->getRitzVectors(),insolver);
1042 MVT::SetBlock(*v,tosol,*sol.Evecs);
1044 std::vector<Value<ScalarType> > fromsolver = tm_solver->getRitzValues();
1045 for (
unsigned int i=0; i<insolver.size(); i++) {
1046 vals[i] = fromsolver[insolver[i]].realpart;
1051 if (inlocked.size() > 0) {
1052 int solnum = insolver.size();
1054 int lclnum = inlocked.size();
1055 std::vector<int> tosol(lclnum);
1056 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1057 RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1058 MVT::SetBlock(*v,tosol,*sol.Evecs);
1060 for (
unsigned int i=0; i<inlocked.size(); i++) {
1061 vals[i+solnum] = lockvals[inlocked[i]];
1069 for(
size_t i=0; i<vals.size(); i++)
1070 vals[i] = ONE/vals[i];
1075 std::vector<int> order(sol.numVecs);
1076 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1078 for (
int i=0; i<sol.numVecs; i++) {
1079 sol.Evals[i].realpart = vals[i];
1080 sol.Evals[i].imagpart = MT::zero();
1083 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1087 sol.index.resize(sol.numVecs,0);
1092 tm_solver->currentStatus(printer_->stream(
FinalSummary));
1097#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1099 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails ) );
1103 problem_->setSolution(sol);
1104 printer_->stream(
Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1107 numIters_ = tm_solver->getNumIters();
1109 if (sol.numVecs < nev) {