29 #ifndef ANASAZI_TraceMinBase_SOLMGR_HPP 30 #define ANASAZI_TraceMinBase_SOLMGR_HPP 47 #include "AnasaziStatusTestSpecTrans.hpp" 54 #include "Teuchos_TimeMonitor.hpp" 56 # include <Teuchos_FancyOStream.hpp> 101 template<
class ScalarType,
class MV,
class OP>
107 typedef Teuchos::ScalarTraits<ScalarType> SCT;
108 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
109 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
159 Teuchos::ParameterList &pl );
186 return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
240 RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
245 int blockSize_, numBlocks_, numRestartBlocks_;
248 RCP<BasicOutputManager<ScalarType> > printer_;
251 MagnitudeType convTol_;
256 MagnitudeType lockTol_;
257 int maxLocked_, lockQuorum_;
258 bool useLocking_, relLockTol_;
262 enum WhenToShiftType whenToShift_;
263 MagnitudeType traceThresh_, shiftTol_;
264 enum HowToShiftType howToShift_;
265 bool useMultipleShifts_, relShiftTol_, considerClusters_;
266 std::string shiftNorm_;
270 std::string ortho_, which_;
271 enum SaddleSolType saddleSolType_;
272 bool projectAllVecs_, projectLockedVecs_, computeAllRes_, useRHSR_, useHarmonic_, noSort_;
273 MagnitudeType alpha_;
276 RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
279 RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
280 RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
281 RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
286 void setParameters(Teuchos::ParameterList &pl)
const;
288 void printParameters(std::ostream &os)
const;
290 virtual RCP< TraceMinBase<ScalarType,MV,OP> > createSolver(
291 const RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
294 Teuchos::ParameterList &plist
306 template<
class ScalarType,
class MV,
class OP>
309 Teuchos::ParameterList &pl ) :
311 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
312 , _timerSolve(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBaseSolMgr::solve()")),
313 _timerRestarting(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBaseSolMgr restarting")),
314 _timerLocking(Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBaseSolMgr locking"))
317 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
318 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
319 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument,
"Problem not symmetric.");
320 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
342 std::string fntemplate =
"";
343 bool allProcs =
false;
344 if (pl.isParameter(
"Output on all processors")) {
345 if (Teuchos::isParameterType<bool>(pl,
"Output on all processors")) {
346 allProcs = pl.get(
"Output on all processors",allProcs);
348 allProcs = ( Teuchos::getParameter<int>(pl,
"Output on all processors") != 0 );
351 fntemplate = pl.get(
"Output filename template",fntemplate);
356 MPI_Initialized(&mpiStarted);
357 if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
362 if (fntemplate !=
"") {
363 std::ostringstream MyPIDstr;
367 while ( (pos = fntemplate.find(
"%d",start)) != -1 ) {
368 fntemplate.replace(pos,2,MyPIDstr.str());
373 if (fntemplate !=
"") {
374 osp = rcp(
new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
376 osp = Teuchos::rcpFromRef(std::cout);
377 std::cout <<
"Anasazi::TraceMinBaseSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
381 osp = Teuchos::rcpFromRef(std::cout);
385 if (pl.isParameter(
"Verbosity")) {
386 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
387 verbosity = pl.get(
"Verbosity", verbosity);
389 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
405 convTol_ = pl.get(
"Convergence Tolerance",MT::prec());
406 TEUCHOS_TEST_FOR_EXCEPTION(convTol_ < 0, std::invalid_argument,
407 "Anasazi::TraceMinBaseSolMgr: \"Convergence Tolerance\" must be nonnegative.");
409 relConvTol_ = pl.get(
"Relative Convergence Tolerance",
true);
410 strtmp = pl.get(
"Convergence Norm",std::string(
"2"));
412 convNorm_ = RES_2NORM;
414 else if (strtmp ==
"M") {
415 convNorm_ = RES_ORTH;
418 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
419 "Anasazi::TraceMinBaseSolMgr: Invalid Convergence Norm.");
424 useLocking_ = pl.get(
"Use Locking",
true);
425 relLockTol_ = pl.get(
"Relative Locking Tolerance",
true);
426 lockTol_ = pl.get(
"Locking Tolerance",convTol_/10);
428 TEUCHOS_TEST_FOR_EXCEPTION(relConvTol_ != relLockTol_, std::invalid_argument,
429 "Anasazi::TraceMinBaseSolMgr: \"Relative Convergence Tolerance\" and \"Relative Locking Tolerance\" have different values. If you set one, you should always set the other.");
431 TEUCHOS_TEST_FOR_EXCEPTION(lockTol_ < 0, std::invalid_argument,
432 "Anasazi::TraceMinBaseSolMgr: \"Locking Tolerance\" must be nonnegative.");
434 strtmp = pl.get(
"Locking Norm",std::string(
"2"));
436 lockNorm_ = RES_2NORM;
438 else if (strtmp ==
"M") {
439 lockNorm_ = RES_ORTH;
442 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
443 "Anasazi::TraceMinBaseSolMgr: Invalid Locking Norm.");
448 maxLocked_ = pl.get(
"Max Locked",problem_->getNEV());
449 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ <= 0, std::invalid_argument,
450 "Anasazi::TraceMinBaseSolMgr: \"Max Locked\" must be strictly positive.");
457 lockQuorum_ = pl.get(
"Locking Quorum",1);
458 TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0, std::invalid_argument,
459 "Anasazi::TraceMinBaseSolMgr: \"Locking Quorum\" must be strictly positive.");
466 strtmp = pl.get(
"When To Shift",
"Always");
468 if(strtmp ==
"Never")
469 whenToShift_ = NEVER_SHIFT;
470 else if(strtmp ==
"After Trace Levels")
471 whenToShift_ = SHIFT_WHEN_TRACE_LEVELS;
472 else if(strtmp ==
"Residual Becomes Small")
473 whenToShift_ = SHIFT_WHEN_RESID_SMALL;
474 else if(strtmp ==
"Always")
475 whenToShift_ = ALWAYS_SHIFT;
477 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
478 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"When To Shift\"; valid options are \"Never\", \"After Trace Levels\", \"Residual Becomes Small\", \"Always\".");
482 traceThresh_ = pl.get(
"Trace Threshold", 0.02);
484 TEUCHOS_TEST_FOR_EXCEPTION(traceThresh_ < 0, std::invalid_argument,
485 "Anasazi::TraceMinBaseSolMgr: \"Trace Threshold\" must be nonnegative.");
488 shiftTol_ = pl.get(
"Shift Tolerance", 0.1);
490 TEUCHOS_TEST_FOR_EXCEPTION(shiftTol_ < 0, std::invalid_argument,
491 "Anasazi::TraceMinBaseSolMgr: \"Shift Tolerance\" must be nonnegative.");
494 relShiftTol_ = pl.get(
"Relative Shift Tolerance",
true);
497 shiftNorm_ = pl.get(
"Shift Norm",
"2");
499 TEUCHOS_TEST_FOR_EXCEPTION(shiftNorm_ !=
"2" && shiftNorm_ !=
"M", std::invalid_argument,
500 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Shift Norm\"; valid options are \"2\" and \"M\".");
502 noSort_ = pl.get(
"No Sorting",
false);
505 strtmp = pl.get(
"How To Choose Shift",
"Adjusted Ritz Values");
507 if(strtmp ==
"Largest Converged")
508 howToShift_ = LARGEST_CONVERGED_SHIFT;
509 else if(strtmp ==
"Adjusted Ritz Values")
510 howToShift_ = ADJUSTED_RITZ_SHIFT;
511 else if(strtmp ==
"Ritz Values")
512 howToShift_ = RITZ_VALUES_SHIFT;
513 else if(strtmp ==
"Experimental Shift")
514 howToShift_ = EXPERIMENTAL_SHIFT;
516 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
517 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"How To Choose Shift\"; valid options are \"Largest Converged\", \"Adjusted Ritz Values\", \"Ritz Values\".");
520 considerClusters_ = pl.get(
"Consider Clusters",
true);
523 useMultipleShifts_ = pl.get(
"Use Multiple Shifts",
true);
529 ortho_ = pl.get(
"Orthogonalization",
"SVQB");
530 TEUCHOS_TEST_FOR_EXCEPTION(ortho_ !=
"DGKS" && ortho_ !=
"SVQB" && ortho_ !=
"ICGS", std::invalid_argument,
531 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Orthogonalization\"; valid options are \"DGKS\", \"SVQB\", \"ICGS\".");
533 strtmp = pl.get(
"Saddle Solver Type",
"Projected Krylov");
535 if(strtmp ==
"Projected Krylov")
536 saddleSolType_ = PROJECTED_KRYLOV_SOLVER;
537 else if(strtmp ==
"Schur Complement")
538 saddleSolType_ = SCHUR_COMPLEMENT_SOLVER;
539 else if(strtmp ==
"Block Diagonal Preconditioned Minres")
540 saddleSolType_ = BD_PREC_MINRES;
541 else if(strtmp ==
"HSS Preconditioned Gmres")
542 saddleSolType_ = HSS_PREC_GMRES;
544 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
545 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Saddle Solver Type\"; valid options are \"Projected Krylov\", \"Schur Complement\", and \"Block Diagonal Preconditioned Minres\".");
547 projectAllVecs_ = pl.get(
"Project All Vectors",
true);
548 projectLockedVecs_ = pl.get(
"Project Locked Vectors",
true);
549 computeAllRes_ = pl.get(
"Compute All Residuals",
true);
550 useRHSR_ = pl.get(
"Use Residual as RHS",
false);
551 alpha_ = pl.get(
"HSS: alpha", 1.0);
553 TEUCHOS_TEST_FOR_EXCEPTION(projectLockedVecs_ && ! projectAllVecs_, std::invalid_argument,
554 "Anasazi::TraceMinBaseSolMgr: If you want to project out the locked vectors, you should really project out ALL the vectors of X.");
557 maxKrylovIter_ = pl.get(
"Maximum Krylov Iterations", 200);
558 TEUCHOS_TEST_FOR_EXCEPTION(maxKrylovIter_ < 1, std::invalid_argument,
559 "Anasazi::TraceMinBaseSolMgr: \"Maximum Krylov Iterations\" must be greater than 0.");
562 which_ = pl.get(
"Which",
"SM");
563 TEUCHOS_TEST_FOR_EXCEPTION(which_ !=
"SM" && which_ !=
"LM", std::invalid_argument,
564 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Which\"; valid options are \"SM\" and \"LM\".");
568 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getOperator() == Teuchos::null && whenToShift_ != NEVER_SHIFT, std::invalid_argument,
569 "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\".");
571 #ifdef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP 574 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER && useMultipleShifts_, std::invalid_argument,
575 "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.");
582 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER, std::invalid_argument,
583 "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)");
595 template<
class ScalarType,
class MV,
class OP>
601 const int nev = problem_->getNEV();
604 RCP<Teuchos::FancyOStream>
605 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(
Debug)));
606 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
607 *out <<
"Entering Anasazi::TraceMinBaseSolMgr::solve()\n";
619 RCP<const OP> swapHelper = problem_->getOperator();
620 problem_->setOperator(problem_->getM());
621 problem_->setM(swapHelper);
622 problem_->setProblem();
629 RCP<StatusTest<ScalarType,MV,OP> > convtest;
630 if (globalTest_ == Teuchos::null) {
634 convtest = rcp(
new StatusTestSpecTrans<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_,
true,problem_->getOperator()) );
637 convtest = globalTest_;
639 RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
642 RCP<StatusTest<ScalarType,MV,OP> > locktest;
644 if (lockingTest_ == Teuchos::null) {
648 locktest = rcp(
new StatusTestSpecTrans<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_,
true,problem_->getOperator()) );
651 locktest = lockingTest_;
655 Teuchos::Array<RCP<StatusTest<ScalarType,MV,OP> > > alltests;
656 alltests.push_back(ordertest);
657 if (locktest != Teuchos::null) alltests.push_back(locktest);
658 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
660 RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
663 RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
664 if ( printer_->isVerbosity(
Debug) ) {
673 RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
674 if (ortho_==
"SVQB") {
676 }
else if (ortho_==
"DGKS") {
678 }
else if (ortho_==
"ICGS") {
681 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): Invalid orthogonalization type.");
686 Teuchos::ParameterList plist;
687 setParameters(plist);
691 RCP<TraceMinBase<ScalarType,MV,OP> > tm_solver
692 = createSolver(sorter,outputtest,ortho,plist);
694 RCP< const MV > probauxvecs = problem_->getAuxVecs();
695 if (probauxvecs != Teuchos::null) {
696 tm_solver->setAuxVecs( Teuchos::tuple< RCP<const MV> >(probauxvecs) );
703 int curNumLocked = 0;
710 if (maxLocked_ > 0) {
711 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
713 std::vector<MagnitudeType> lockvals;
758 const ScalarType ONE = SCT::one();
759 const ScalarType ZERO = SCT::zero();
764 problem_->setSolution(sol);
770 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 771 Teuchos::TimeMonitor slvtimer(*_timerSolve);
777 tm_solver->iterate();
783 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
784 throw AnasaziError(
"Anasazi::TraceMinBaseSolMgr::solve(): User-specified debug status test returned Passed.");
791 else if (ordertest->getStatus() ==
Passed ) {
802 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
804 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 805 Teuchos::TimeMonitor lcktimer(*_timerLocking);
811 const int curdim = state.
curDim;
815 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
816 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() non-positive.");
817 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
818 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
820 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
821 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: locking not deactivated.");
824 std::vector<int> tmp_vector_int;
825 if (curNumLocked + locktest->howMany() > maxLocked_) {
827 for(
int i=0; i<maxLocked_-curNumLocked; i++)
828 tmp_vector_int.push_back(locktest->whichVecs()[i]);
832 tmp_vector_int = locktest->whichVecs();
835 const std::vector<int> lockind(tmp_vector_int);
836 const int numNewLocked = lockind.size();
840 const int numUnlocked = curdim-numNewLocked;
841 tmp_vector_int.resize(curdim);
842 for (
int i=0; i<curdim; i++) tmp_vector_int[i] = i;
843 const std::vector<int> curind(tmp_vector_int);
844 tmp_vector_int.resize(numUnlocked);
845 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
846 const std::vector<int> unlockind(tmp_vector_int);
847 tmp_vector_int.clear();
851 if (printer_->isVerbosity(
Debug)) {
852 printer_->print(
Debug,
"Locking vectors: ");
853 for (
unsigned int i=0; i<lockind.size(); i++) {printer_->stream(
Debug) <<
" " << lockind[i];}
854 printer_->print(
Debug,
"\n");
855 printer_->print(
Debug,
"Not locking vectors: ");
856 for (
unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(
Debug) <<
" " << unlockind[i];}
857 printer_->print(
Debug,
"\n");
861 std::vector<Value<ScalarType> > allvals = tm_solver->getRitzValues();
862 for(
unsigned int i=0; i<allvals.size(); i++)
863 printer_->stream(
Debug) <<
"Ritz value[" << i <<
"] = " << allvals[i].realpart << std::endl;
864 for (
int i=0; i<numNewLocked; i++) {
865 lockvals.push_back(allvals[lockind[i]].realpart);
869 RCP<const MV> newLocked = MVT::CloneView(*tm_solver->getRitzVectors(),lockind);
870 std::vector<int> indlock(numNewLocked);
871 for (
int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
874 RCP<MV> tempMV = MVT::CloneCopy(*newLocked);
875 ortho->normalizeMat(*tempMV);
876 MVT::SetBlock(*tempMV,indlock,*lockvecs);
880 MVT::SetBlock(*newLocked,indlock,*lockvecs);
888 for(
unsigned int aliciaInd=0; aliciaInd<lockvals.size(); aliciaInd++)
890 lockvals[aliciaInd] = 0.0;
893 ordertest->setAuxVals(lockvals);
897 curNumLocked += numNewLocked;
899 if(ordertest->getStatus() ==
Passed)
break;
901 std::vector<int> curlockind(curNumLocked);
902 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
903 RCP<const MV> curlocked = MVT::CloneView(*lockvecs,curlockind);
905 Teuchos::Array< RCP<const MV> > aux;
906 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
907 aux.push_back(curlocked);
908 tm_solver->setAuxVecs(aux);
911 printer_->stream(
Debug) <<
"curNumLocked: " << curNumLocked << std::endl;
912 printer_->stream(
Debug) <<
"maxLocked: " << maxLocked_ << std::endl;
913 if (curNumLocked == maxLocked_) {
915 combotest->removeTest(locktest);
916 locktest = Teuchos::null;
917 printer_->stream(
Debug) <<
"Removed locking test\n";
920 int newdim = numRestartBlocks_*blockSize_;
922 if(newdim <= numUnlocked)
926 std::vector<int> desiredSubscripts(newdim);
927 for(
int i=0; i<newdim; i++)
929 desiredSubscripts[i] = unlockind[i];
930 printer_->stream(
Debug) <<
"H desiredSubscripts[" << i <<
"] = " << desiredSubscripts[i] << std::endl;
932 newstate.
V = MVT::CloneView(*tm_solver->getRitzVectors(),desiredSubscripts);
937 std::vector<int> desiredSubscripts(newdim);
938 for(
int i=0; i<newdim; i++)
940 desiredSubscripts[i] = unlockind[i];
941 printer_->stream(
Debug) <<
"desiredSubscripts[" << i <<
"] = " << desiredSubscripts[i] << std::endl;
944 copyPartOfState(state, newstate, desiredSubscripts);
952 int nrandom = newdim-numUnlocked;
954 RCP<const MV> helperMV;
955 RCP<MV> totalV = MVT::Clone(*tm_solver->getRitzVectors(),newdim);
958 tmp_vector_int.resize(numUnlocked);
959 for(
int i=0; i<numUnlocked; i++) tmp_vector_int[i] = i;
960 RCP<MV> oldV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
964 helperMV = MVT::CloneView(*tm_solver->getRitzVectors(),unlockind);
966 helperMV = MVT::CloneView(*state.
V,unlockind);
967 MVT::Assign(*helperMV,*oldV);
970 tmp_vector_int.resize(nrandom);
971 for(
int i=0; i<nrandom; i++) tmp_vector_int[i] = i+numUnlocked;
972 RCP<MV> newV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
975 MVT::MvRandom(*newV);
981 RCP< std::vector<ScalarType> > helperRS = rcp(
new std::vector<ScalarType>(blockSize_) );
982 for(
unsigned int i=0; i<(
unsigned int)blockSize_; i++)
984 if(i < unlockind.size() && unlockind[i] < blockSize_)
985 (*helperRS)[i] = (*state.
ritzShifts)[unlockind[i]];
987 (*helperRS)[i] = ZERO;
994 for(
size_t i=0; i<lockvals.size(); i++)
1000 newstate.
NEV = state.
NEV - numNewLocked;
1001 tm_solver->initialize(newstate);
1008 else if ( needToRestart(tm_solver) ) {
1010 if(performRestart(numRestarts, tm_solver) ==
false)
1020 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): Invalid return from tm_solver::iterate().");
1025 <<
"Anasazi::TraceMinBaseSolMgr::solve() caught unexpected exception from Anasazi::TraceMinBase::iterate() at iteration " << tm_solver->getNumIters() << std::endl
1026 << err.what() << std::endl
1027 <<
"Anasazi::TraceMinBaseSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1032 sol.
numVecs = ordertest->howMany();
1034 sol.
Evecs = MVT::Clone(*problem_->getInitVec(),sol.
numVecs);
1037 std::vector<MagnitudeType> vals(sol.
numVecs);
1040 std::vector<int> which = ordertest->whichVecs();
1044 std::vector<int> inlocked(0), insolver(0);
1045 for (
unsigned int i=0; i<which.size(); i++) {
1046 if (which[i] >= 0) {
1047 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): positive indexing mistake from ordertest.");
1048 insolver.push_back(which[i]);
1052 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): negative indexing mistake from ordertest.");
1053 inlocked.push_back(which[i] + curNumLocked);
1057 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (
unsigned int)sol.
numVecs,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): indexing mistake.");
1060 if (insolver.size() > 0) {
1062 int lclnum = insolver.size();
1063 std::vector<int> tosol(lclnum);
1064 for (
int i=0; i<lclnum; i++) tosol[i] = i;
1065 RCP<const MV> v = MVT::CloneView(*tm_solver->getRitzVectors(),insolver);
1066 MVT::SetBlock(*v,tosol,*sol.
Evecs);
1068 std::vector<Value<ScalarType> > fromsolver = tm_solver->getRitzValues();
1069 for (
unsigned int i=0; i<insolver.size(); i++) {
1070 vals[i] = fromsolver[insolver[i]].realpart;
1075 if (inlocked.size() > 0) {
1076 int solnum = insolver.size();
1078 int lclnum = inlocked.size();
1079 std::vector<int> tosol(lclnum);
1080 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1081 RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1082 MVT::SetBlock(*v,tosol,*sol.
Evecs);
1084 for (
unsigned int i=0; i<inlocked.size(); i++) {
1085 vals[i+solnum] = lockvals[inlocked[i]];
1093 for(
size_t i=0; i<vals.size(); i++)
1094 vals[i] = ONE/vals[i];
1099 std::vector<int> order(sol.
numVecs);
1100 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.
numVecs);
1102 for (
int i=0; i<sol.
numVecs; i++) {
1103 sol.
Evals[i].realpart = vals[i];
1104 sol.
Evals[i].imagpart = MT::zero();
1116 tm_solver->currentStatus(printer_->stream(
FinalSummary));
1121 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1123 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails ) );
1127 problem_->setSolution(sol);
1128 printer_->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1131 numIters_ = tm_solver->getNumIters();
1140 template <
class ScalarType,
class MV,
class OP>
1145 globalTest_ = global;
1148 template <
class ScalarType,
class MV,
class OP>
1149 const RCP< StatusTest<ScalarType,MV,OP> > &
1155 template <
class ScalarType,
class MV,
class OP>
1163 template <
class ScalarType,
class MV,
class OP>
1164 const RCP< StatusTest<ScalarType,MV,OP> > &
1170 template <
class ScalarType,
class MV,
class OP>
1175 lockingTest_ = locking;
1178 template <
class ScalarType,
class MV,
class OP>
1179 const RCP< StatusTest<ScalarType,MV,OP> > &
1182 return lockingTest_;
1185 template <
class ScalarType,
class MV,
class OP>
1188 const ScalarType ONE = Teuchos::ScalarTraits<MagnitudeType>::one();
1189 const ScalarType ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
1191 newState.
curDim = indToCopy.size();
1192 std::vector<int> fullIndices(oldState.
curDim);
1193 for(
int i=0; i<oldState.
curDim; i++) fullIndices[i] = i;
1201 std::vector<int> oldIndices;
1202 std::vector<int> newIndices;
1203 for(
int i=0; i<newState.
curDim; i++)
1205 if(indToCopy[i] < blockSize_)
1206 oldIndices.push_back(indToCopy[i]);
1208 newIndices.push_back(indToCopy[i]);
1211 int olddim = oldIndices.size();
1212 int newdim = newIndices.size();
1217 newState.
V = MVT::CloneView(*oldState.
X, indToCopy);
1218 newState.
R = MVT::CloneView(*oldState.
R, indToCopy);
1219 newState.
X = newState.
V;
1221 if(problem_->getOperator() != Teuchos::null)
1223 newState.
KV = MVT::CloneView(*oldState.
KX, indToCopy);
1224 newState.
KX = newState.
KV;
1228 newState.
KV = Teuchos::null;
1229 newState.
KX = Teuchos::null;
1232 if(problem_->getM() != Teuchos::null)
1234 newState.
MopV = MVT::CloneView(*oldState.
MX, indToCopy);
1235 newState.
MX = newState.
MopV;
1239 newState.
MopV = Teuchos::null;
1240 newState.
MX = Teuchos::null;
1243 else if(newdim == 0)
1245 std::vector<int> blockind(blockSize_);
1246 for(
int i=0; i<blockSize_; i++)
1250 newState.
V = MVT::CloneView(*oldState.
X, blockind);
1251 newState.
KV = MVT::CloneView(*oldState.
KX, blockind);
1252 newState.
R = MVT::CloneView(*oldState.
R, blockind);
1253 newState.
X = MVT::CloneView(*newState.
V, blockind);
1254 newState.
KX = MVT::CloneView(*newState.
KV, blockind);
1256 if(problem_->getM() != Teuchos::null)
1258 newState.
MopV = MVT::CloneView(*oldState.
MX, blockind);
1259 newState.
MX = MVT::CloneView(*newState.
MopV, blockind);
1263 newState.
MopV = Teuchos::null;
1264 newState.
MX = Teuchos::null;
1270 std::vector<int> oldPart(olddim);
1271 for(
int i=0; i<olddim; i++) oldPart[i] = i;
1272 std::vector<int> newPart(newdim);
1273 for(
int i=0; i<newdim; i++) newPart[i] = olddim+i;
1276 RCP<MV> helper = MVT::Clone(*oldState.
V,newState.
curDim);
1277 RCP<MV> oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1278 RCP<MV> newHelper = MVT::CloneViewNonConst(*helper,newPart);
1279 RCP<const MV> viewHelper;
1282 Teuchos::SerialDenseMatrix<int,ScalarType> newRV(oldState.
curDim,newdim);
1283 for(
int r=0; r<oldState.
curDim; r++)
1285 for(
int c=0; c<newdim; c++)
1286 newRV(r,c) = (*oldState.
RV)(r,newIndices[c]);
1290 viewHelper = MVT::CloneView(*oldState.
V,fullIndices);
1291 MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1292 viewHelper = MVT::CloneView(*oldState.
X,oldIndices);
1293 MVT::Assign(*viewHelper,*oldHelper);
1294 newState.
V = MVT::CloneCopy(*helper);
1297 viewHelper = MVT::CloneView(*oldState.
KV,fullIndices);
1298 MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1299 viewHelper = MVT::CloneView(*oldState.
KX,oldIndices);
1300 MVT::Assign(*viewHelper,*oldHelper);
1301 newState.
KV = MVT::CloneCopy(*helper);
1304 if(problem_->getM() != Teuchos::null)
1306 viewHelper = MVT::CloneView(*oldState.
MopV,fullIndices);
1307 MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1308 viewHelper = MVT::CloneView(*oldState.
MX,oldIndices);
1309 MVT::Assign(*viewHelper,*oldHelper);
1310 newState.
MopV = MVT::CloneCopy(*helper);
1313 newState.
MopV = newState.
V;
1316 std::vector<int> blockVec(blockSize_);
1317 for(
int i=0; i<blockSize_; i++) blockVec[i] = i;
1318 newState.
X = MVT::CloneView(*newState.
V,blockVec);
1319 newState.
KX = MVT::CloneView(*newState.
KV,blockVec);
1320 newState.
MX = MVT::CloneView(*newState.
MopV,blockVec);
1323 if(blockSize_-oldIndices.size() > 0)
1326 newPart.resize(blockSize_-oldIndices.size());
1327 helper = MVT::Clone(*oldState.
V,blockSize_);
1328 oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1329 newHelper = MVT::CloneViewNonConst(*helper,newPart);
1331 RCP<MV> scaledMV = MVT::CloneCopy(*newState.
MX,newPart);
1332 RCP<const MV> localKV = MVT::CloneView(*newState.
KX,newPart);
1333 std::vector<ScalarType> scalarVec(blockSize_-oldIndices.size());
1334 for(
unsigned int i=0; i<(
unsigned int)blockSize_-oldIndices.size(); i++) scalarVec[i] = (*oldState.
T)[newPart[i]];
1335 MVT::MvScale(*scaledMV,scalarVec);
1337 helper = MVT::Clone(*oldState.
V,blockSize_);
1338 oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1339 newHelper = MVT::CloneViewNonConst(*helper,newPart);
1340 MVT::MvAddMv(ONE,*localKV,-ONE,*scaledMV,*newHelper);
1341 viewHelper = MVT::CloneView(*oldState.
R,oldIndices);
1342 MVT::Assign(*viewHelper,*oldHelper);
1343 newState.
R = MVT::CloneCopy(*helper);
1346 newState.
R = oldState.
R;
1353 RCP< std::vector<ScalarType> > helperT = rcp(
new std::vector<ScalarType>(newState.
curDim) );
1354 for(
int i=0; i<newState.
curDim; i++) (*helperT)[i] = (*oldState.
T)[indToCopy[i]];
1355 newState.
T = helperT;
1358 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newKK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.
curDim,newState.
curDim) );
1359 for(
int i=0; i<newState.
curDim; i++)
1360 (*newKK)(i,i) = (*newState.
T)[i];
1361 newState.
KK = newKK;
1364 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newRV = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.
curDim,newState.
curDim) );
1365 for(
int i=0; i<newState.
curDim; i++)
1366 (*newRV)(i,i) = ONE;
1367 newState.
RV = newRV;
1370 RCP< std::vector<ScalarType> > helperRS = rcp(
new std::vector<ScalarType>(blockSize_) );
1371 for(
int i=0; i<blockSize_; i++)
1373 if(indToCopy[i] < blockSize_)
1374 (*helperRS)[i] = (*oldState.
ritzShifts)[indToCopy[i]];
1376 (*helperRS)[i] = ZERO;
1382 template <
class ScalarType,
class MV,
class OP>
1383 void TraceMinBaseSolMgr<ScalarType,MV,OP>::setParameters(Teuchos::ParameterList &pl)
const 1385 pl.set(
"Block Size", blockSize_);
1386 pl.set(
"Num Blocks", numBlocks_);
1387 pl.set(
"Num Restart Blocks", numRestartBlocks_);
1388 pl.set(
"When To Shift", whenToShift_);
1389 pl.set(
"Trace Threshold", traceThresh_);
1390 pl.set(
"Shift Tolerance", shiftTol_);
1391 pl.set(
"Relative Shift Tolerance", relShiftTol_);
1392 pl.set(
"Shift Norm", shiftNorm_);
1393 pl.set(
"How To Choose Shift", howToShift_);
1394 pl.set(
"Consider Clusters", considerClusters_);
1395 pl.set(
"Use Multiple Shifts", useMultipleShifts_);
1396 pl.set(
"Saddle Solver Type", saddleSolType_);
1397 pl.set(
"Project All Vectors", projectAllVecs_);
1398 pl.set(
"Project Locked Vectors", projectLockedVecs_);
1399 pl.set(
"Compute All Residuals", computeAllRes_);
1400 pl.set(
"Use Residual as RHS", useRHSR_);
1401 pl.set(
"Use Harmonic Ritz Values", useHarmonic_);
1402 pl.set(
"Maximum Krylov Iterations", maxKrylovIter_);
1403 pl.set(
"HSS: alpha", alpha_);
1407 template <
class ScalarType,
class MV,
class OP>
1408 void TraceMinBaseSolMgr<ScalarType,MV,OP>::printParameters(std::ostream &os)
const 1411 os <<
"========================================\n";
1412 os <<
"========= TraceMin parameters ==========\n";
1413 os <<
"========================================\n";
1414 os <<
"=========== Block parameters ===========\n";
1415 os <<
"Block Size: " << blockSize_ << std::endl;
1416 os <<
"Num Blocks: " << numBlocks_ << std::endl;
1417 os <<
"Num Restart Blocks: " << numRestartBlocks_ << std::endl;
1418 os <<
"======== Convergence parameters ========\n";
1419 os <<
"Convergence Tolerance: " << convTol_ << std::endl;
1420 os <<
"Relative Convergence Tolerance: " << relConvTol_ << std::endl;
1421 os <<
"========== Locking parameters ==========\n";
1422 os <<
"Use Locking: " << useLocking_ << std::endl;
1423 os <<
"Locking Tolerance: " << lockTol_ << std::endl;
1424 os <<
"Relative Locking Tolerance: " << relLockTol_ << std::endl;
1425 os <<
"Max Locked: " << maxLocked_ << std::endl;
1426 os <<
"Locking Quorum: " << lockQuorum_ << std::endl;
1427 os <<
"========== Shifting parameters =========\n";
1428 os <<
"When To Shift: ";
1429 if(whenToShift_ == NEVER_SHIFT) os <<
"Never\n";
1430 else if(whenToShift_ == SHIFT_WHEN_TRACE_LEVELS) os <<
"After Trace Levels\n";
1431 else if(whenToShift_ == SHIFT_WHEN_RESID_SMALL) os <<
"Residual Becomes Small\n";
1432 else if(whenToShift_ == ALWAYS_SHIFT) os <<
"Always\n";
1433 os <<
"Consider Clusters: " << considerClusters_ << std::endl;
1434 os <<
"Trace Threshohld: " << traceThresh_ << std::endl;
1435 os <<
"Shift Tolerance: " << shiftTol_ << std::endl;
1436 os <<
"Relative Shift Tolerance: " << relShiftTol_ << std::endl;
1437 os <<
"How To Choose Shift: ";
1438 if(howToShift_ == LARGEST_CONVERGED_SHIFT) os <<
"Largest Converged\n";
1439 else if(howToShift_ == ADJUSTED_RITZ_SHIFT) os <<
"Adjusted Ritz Values\n";
1440 else if(howToShift_ == RITZ_VALUES_SHIFT) os <<
"Ritz Values\n";
1441 os <<
"Use Multiple Shifts: " << useMultipleShifts_ << std::endl;
1442 os <<
"=========== Other parameters ===========\n";
1443 os <<
"Orthogonalization: " << ortho_ << std::endl;
1444 os <<
"Saddle Solver Type: ";
1445 if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER) os <<
"Projected Krylov\n";
1446 else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER) os <<
"Schur Complement\n";
1447 os <<
"Project All Vectors: " << projectAllVecs_ << std::endl;
1448 os <<
"Project Locked Vectors: " << projectLockedVecs_ << std::endl;
1449 os <<
"Compute All Residuals: " << computeAllRes_ << std::endl;
1450 os <<
"========================================\n\n\n";
Pure virtual base class which describes the basic interface for a solver manager. ...
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Abstract base class for trace minimization eigensolvers.
int NEV
Number of unconverged eigenvalues.
TraceMinBaseSolMgr(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for TraceMinBaseSolMgr.
const RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
RCP< const MV > X
The current eigenvectors.
A special StatusTest for printing other status tests.
This class defines the interface required by an eigensolver and status test class to compute solution...
Structure to contain pointers to TraceMinBase state variables.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
RCP< const MV > V
The current basis.
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
Status test for forming logical combinations of other status tests.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
void setGlobalStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
Basic implementation of the Anasazi::SortManager class.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
An exception class parent to all Anasazi exceptions.
bool isOrtho
Whether V has been projected and orthonormalized already.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Anasazi's templated, static class providing utilities for the solvers.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
int numVecs
The number of computed eigenpairs.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
Anasazi's basic output manager for sending information of select verbosity levels to the appropriate ...
Abstract base class which defines the interface required by an eigensolver and status test class to c...
ReturnType
Enumerated type used to pass back information from a solver manager.
const RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
A status test for testing the norm of the eigenvectors residuals.
Traits class which defines basic operations on multivectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
int curDim
The current dimension of the solver.
RCP< const MV > KV
The image of V under K.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
The Anasazi::TraceMinBaseSolMgr provides an abstract base class for the TraceMin series of solver man...
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > RV
The current Ritz vectors.
ScalarType largestSafeShift
Largest safe shift.
RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
Teuchos::Array< RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
virtual ~TraceMinBaseSolMgr()
Destructor.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
RCP< const MV > KX
The image of the current eigenvectors under K.
RCP< const MV > R
The current residual vectors.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
void setLockingStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
Status test for forming logical combinations of other status tests.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
int getNumIters() const
Get the iteration count for the most recent call to solve().
const RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
This is an abstract base class for the trace minimization eigensolvers.
void setDebugStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Common interface of stopping criteria for Anasazi's solvers.
A status test for testing the norm of the eigenvectors residuals.
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::OrthoManager class.
Class which provides internal utilities for the Anasazi solvers.