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;
514 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
515 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"How To Choose Shift\"; valid options are \"Largest Converged\", \"Adjusted Ritz Values\", \"Ritz Values\".");
518 considerClusters_ = pl.get(
"Consider Clusters",
true);
521 useMultipleShifts_ = pl.get(
"Use Multiple Shifts",
true);
527 ortho_ = pl.get(
"Orthogonalization",
"SVQB");
528 TEUCHOS_TEST_FOR_EXCEPTION(ortho_ !=
"DGKS" && ortho_ !=
"SVQB" && ortho_ !=
"ICGS", std::invalid_argument,
529 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Orthogonalization\"; valid options are \"DGKS\", \"SVQB\", \"ICGS\".");
531 strtmp = pl.get(
"Saddle Solver Type",
"Projected Krylov");
533 if(strtmp ==
"Projected Krylov")
534 saddleSolType_ = PROJECTED_KRYLOV_SOLVER;
535 else if(strtmp ==
"Schur Complement")
536 saddleSolType_ = SCHUR_COMPLEMENT_SOLVER;
537 else if(strtmp ==
"Block Diagonal Preconditioned Minres")
538 saddleSolType_ = BD_PREC_MINRES;
539 else if(strtmp ==
"HSS Preconditioned Gmres")
540 saddleSolType_ = HSS_PREC_GMRES;
542 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
543 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Saddle Solver Type\"; valid options are \"Projected Krylov\", \"Schur Complement\", and \"Block Diagonal Preconditioned Minres\".");
545 projectAllVecs_ = pl.get(
"Project All Vectors",
true);
546 projectLockedVecs_ = pl.get(
"Project Locked Vectors",
true);
547 computeAllRes_ = pl.get(
"Compute All Residuals",
true);
548 useRHSR_ = pl.get(
"Use Residual as RHS",
false);
549 alpha_ = pl.get(
"HSS: alpha", 1.0);
551 TEUCHOS_TEST_FOR_EXCEPTION(projectLockedVecs_ && ! projectAllVecs_, std::invalid_argument,
552 "Anasazi::TraceMinBaseSolMgr: If you want to project out the locked vectors, you should really project out ALL the vectors of X.");
555 maxKrylovIter_ = pl.get(
"Maximum Krylov Iterations", 200);
556 TEUCHOS_TEST_FOR_EXCEPTION(maxKrylovIter_ < 1, std::invalid_argument,
557 "Anasazi::TraceMinBaseSolMgr: \"Maximum Krylov Iterations\" must be greater than 0.");
560 which_ = pl.get(
"Which",
"SM");
561 TEUCHOS_TEST_FOR_EXCEPTION(which_ !=
"SM" && which_ !=
"LM", std::invalid_argument,
562 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Which\"; valid options are \"SM\" and \"LM\".");
566 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getOperator() == Teuchos::null && whenToShift_ != NEVER_SHIFT, std::invalid_argument,
567 "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\".");
569 #ifdef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP 572 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER && useMultipleShifts_, std::invalid_argument,
573 "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.");
580 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER, std::invalid_argument,
581 "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)");
593 template<
class ScalarType,
class MV,
class OP>
599 const int nev = problem_->getNEV();
602 RCP<Teuchos::FancyOStream>
603 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(
Debug)));
604 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
605 *out <<
"Entering Anasazi::TraceMinBaseSolMgr::solve()\n";
617 RCP<const OP> swapHelper = problem_->getOperator();
618 problem_->setOperator(problem_->getM());
619 problem_->setM(swapHelper);
620 problem_->setProblem();
627 RCP<StatusTest<ScalarType,MV,OP> > convtest;
628 if (globalTest_ == Teuchos::null) {
632 convtest = rcp(
new StatusTestSpecTrans<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_,
true,problem_->getOperator()) );
635 convtest = globalTest_;
637 RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
640 RCP<StatusTest<ScalarType,MV,OP> > locktest;
642 if (lockingTest_ == Teuchos::null) {
646 locktest = rcp(
new StatusTestSpecTrans<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_,
true,problem_->getOperator()) );
649 locktest = lockingTest_;
653 Teuchos::Array<RCP<StatusTest<ScalarType,MV,OP> > > alltests;
654 alltests.push_back(ordertest);
655 if (locktest != Teuchos::null) alltests.push_back(locktest);
656 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
658 RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
661 RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
662 if ( printer_->isVerbosity(
Debug) ) {
671 RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
672 if (ortho_==
"SVQB") {
674 }
else if (ortho_==
"DGKS") {
676 }
else if (ortho_==
"ICGS") {
679 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): Invalid orthogonalization type.");
684 Teuchos::ParameterList plist;
685 setParameters(plist);
689 RCP<TraceMinBase<ScalarType,MV,OP> > tm_solver
690 = createSolver(sorter,outputtest,ortho,plist);
692 RCP< const MV > probauxvecs = problem_->getAuxVecs();
693 if (probauxvecs != Teuchos::null) {
694 tm_solver->setAuxVecs( Teuchos::tuple< RCP<const MV> >(probauxvecs) );
701 int curNumLocked = 0;
708 if (maxLocked_ > 0) {
709 lockvecs =
MVT::Clone(*problem_->getInitVec(),maxLocked_);
711 std::vector<MagnitudeType> lockvals;
756 const ScalarType ONE = SCT::one();
757 const ScalarType ZERO = SCT::zero();
762 problem_->setSolution(sol);
768 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 769 Teuchos::TimeMonitor slvtimer(*_timerSolve);
775 tm_solver->iterate();
781 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
782 throw AnasaziError(
"Anasazi::TraceMinBaseSolMgr::solve(): User-specified debug status test returned Passed.");
789 else if (ordertest->getStatus() ==
Passed ) {
800 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
802 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 803 Teuchos::TimeMonitor lcktimer(*_timerLocking);
809 const int curdim = state.
curDim;
813 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
814 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() non-positive.");
815 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
816 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
818 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
819 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: locking not deactivated.");
822 std::vector<int> tmp_vector_int;
823 if (curNumLocked + locktest->howMany() > maxLocked_) {
825 for(
int i=0; i<maxLocked_-curNumLocked; i++)
826 tmp_vector_int.push_back(locktest->whichVecs()[i]);
830 tmp_vector_int = locktest->whichVecs();
833 const std::vector<int> lockind(tmp_vector_int);
834 const int numNewLocked = lockind.size();
838 const int numUnlocked = curdim-numNewLocked;
839 tmp_vector_int.resize(curdim);
840 for (
int i=0; i<curdim; i++) tmp_vector_int[i] = i;
841 const std::vector<int> curind(tmp_vector_int);
842 tmp_vector_int.resize(numUnlocked);
843 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
844 const std::vector<int> unlockind(tmp_vector_int);
845 tmp_vector_int.clear();
849 if (printer_->isVerbosity(
Debug)) {
850 printer_->print(
Debug,
"Locking vectors: ");
851 for (
unsigned int i=0; i<lockind.size(); i++) {printer_->stream(
Debug) <<
" " << lockind[i];}
852 printer_->print(
Debug,
"\n");
853 printer_->print(
Debug,
"Not locking vectors: ");
854 for (
unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(
Debug) <<
" " << unlockind[i];}
855 printer_->print(
Debug,
"\n");
859 std::vector<Value<ScalarType> > allvals = tm_solver->getRitzValues();
860 for(
unsigned int i=0; i<allvals.size(); i++)
861 printer_->stream(
Debug) <<
"Ritz value[" << i <<
"] = " << allvals[i].realpart << std::endl;
862 for (
int i=0; i<numNewLocked; i++) {
863 lockvals.push_back(allvals[lockind[i]].realpart);
867 RCP<const MV> newLocked =
MVT::CloneView(*tm_solver->getRitzVectors(),lockind);
868 std::vector<int> indlock(numNewLocked);
869 for (
int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
873 ortho->normalizeMat(*tempMV);
886 for(
unsigned int aliciaInd=0; aliciaInd<lockvals.size(); aliciaInd++)
888 lockvals[aliciaInd] = 0.0;
891 ordertest->setAuxVals(lockvals);
895 curNumLocked += numNewLocked;
897 if(ordertest->getStatus() ==
Passed)
break;
899 std::vector<int> curlockind(curNumLocked);
900 for (
int i=0; i<curNumLocked; i++) curlockind[i] = i;
903 Teuchos::Array< RCP<const MV> > aux;
904 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
905 aux.push_back(curlocked);
906 tm_solver->setAuxVecs(aux);
909 printer_->stream(
Debug) <<
"curNumLocked: " << curNumLocked << std::endl;
910 printer_->stream(
Debug) <<
"maxLocked: " << maxLocked_ << std::endl;
911 if (curNumLocked == maxLocked_) {
913 combotest->removeTest(locktest);
914 locktest = Teuchos::null;
915 printer_->stream(
Debug) <<
"Removed locking test\n";
918 int newdim = numRestartBlocks_*blockSize_;
920 if(newdim <= numUnlocked)
924 std::vector<int> desiredSubscripts(newdim);
925 for(
int i=0; i<newdim; i++)
927 desiredSubscripts[i] = unlockind[i];
928 printer_->stream(
Debug) <<
"H desiredSubscripts[" << i <<
"] = " << desiredSubscripts[i] << std::endl;
930 newstate.
V =
MVT::CloneView(*tm_solver->getRitzVectors(),desiredSubscripts);
935 std::vector<int> desiredSubscripts(newdim);
936 for(
int i=0; i<newdim; i++)
938 desiredSubscripts[i] = unlockind[i];
939 printer_->stream(
Debug) <<
"desiredSubscripts[" << i <<
"] = " << desiredSubscripts[i] << std::endl;
942 copyPartOfState(state, newstate, desiredSubscripts);
950 int nrandom = newdim-numUnlocked;
952 RCP<const MV> helperMV;
953 RCP<MV> totalV =
MVT::Clone(*tm_solver->getRitzVectors(),newdim);
956 tmp_vector_int.resize(numUnlocked);
957 for(
int i=0; i<numUnlocked; i++) tmp_vector_int[i] = i;
962 helperMV =
MVT::CloneView(*tm_solver->getRitzVectors(),unlockind);
968 tmp_vector_int.resize(nrandom);
969 for(
int i=0; i<nrandom; i++) tmp_vector_int[i] = i+numUnlocked;
979 RCP< std::vector<ScalarType> > helperRS = rcp(
new std::vector<ScalarType>(blockSize_) );
980 for(
unsigned int i=0; i<(
unsigned int)blockSize_; i++)
982 if(i < unlockind.size() && unlockind[i] < blockSize_)
983 (*helperRS)[i] = (*state.
ritzShifts)[unlockind[i]];
985 (*helperRS)[i] = ZERO;
992 for(
size_t i=0; i<lockvals.size(); i++)
998 newstate.
NEV = state.
NEV - numNewLocked;
999 tm_solver->initialize(newstate);
1006 else if ( needToRestart(tm_solver) ) {
1008 if(performRestart(numRestarts, tm_solver) ==
false)
1018 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): Invalid return from tm_solver::iterate().");
1023 <<
"Anasazi::TraceMinBaseSolMgr::solve() caught unexpected exception from Anasazi::TraceMinBase::iterate() at iteration " << tm_solver->getNumIters() << std::endl
1024 << err.what() << std::endl
1025 <<
"Anasazi::TraceMinBaseSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1030 sol.
numVecs = ordertest->howMany();
1035 std::vector<MagnitudeType> vals(sol.
numVecs);
1038 std::vector<int> which = ordertest->whichVecs();
1042 std::vector<int> inlocked(0), insolver(0);
1043 for (
unsigned int i=0; i<which.size(); i++) {
1044 if (which[i] >= 0) {
1045 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): positive indexing mistake from ordertest.");
1046 insolver.push_back(which[i]);
1050 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): negative indexing mistake from ordertest.");
1051 inlocked.push_back(which[i] + curNumLocked);
1055 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (
unsigned int)sol.
numVecs,std::logic_error,
"Anasazi::TraceMinBaseSolMgr::solve(): indexing mistake.");
1058 if (insolver.size() > 0) {
1060 int lclnum = insolver.size();
1061 std::vector<int> tosol(lclnum);
1062 for (
int i=0; i<lclnum; i++) tosol[i] = i;
1063 RCP<const MV> v =
MVT::CloneView(*tm_solver->getRitzVectors(),insolver);
1066 std::vector<Value<ScalarType> > fromsolver = tm_solver->getRitzValues();
1067 for (
unsigned int i=0; i<insolver.size(); i++) {
1068 vals[i] = fromsolver[insolver[i]].realpart;
1073 if (inlocked.size() > 0) {
1074 int solnum = insolver.size();
1076 int lclnum = inlocked.size();
1077 std::vector<int> tosol(lclnum);
1078 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1082 for (
unsigned int i=0; i<inlocked.size(); i++) {
1083 vals[i+solnum] = lockvals[inlocked[i]];
1091 for(
size_t i=0; i<vals.size(); i++)
1092 vals[i] = ONE/vals[i];
1097 std::vector<int> order(sol.
numVecs);
1098 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.
numVecs);
1100 for (
int i=0; i<sol.
numVecs; i++) {
1101 sol.
Evals[i].realpart = vals[i];
1102 sol.
Evals[i].imagpart = MT::zero();
1114 tm_solver->currentStatus(printer_->stream(
FinalSummary));
1119 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1121 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails ) );
1125 problem_->setSolution(sol);
1126 printer_->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1129 numIters_ = tm_solver->getNumIters();
1138 template <
class ScalarType,
class MV,
class OP>
1143 globalTest_ = global;
1146 template <
class ScalarType,
class MV,
class OP>
1147 const RCP< StatusTest<ScalarType,MV,OP> > &
1153 template <
class ScalarType,
class MV,
class OP>
1161 template <
class ScalarType,
class MV,
class OP>
1162 const RCP< StatusTest<ScalarType,MV,OP> > &
1168 template <
class ScalarType,
class MV,
class OP>
1173 lockingTest_ = locking;
1176 template <
class ScalarType,
class MV,
class OP>
1177 const RCP< StatusTest<ScalarType,MV,OP> > &
1180 return lockingTest_;
1183 template <
class ScalarType,
class MV,
class OP>
1186 const ScalarType ONE = Teuchos::ScalarTraits<MagnitudeType>::one();
1187 const ScalarType ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
1189 newState.
curDim = indToCopy.size();
1190 std::vector<int> fullIndices(oldState.
curDim);
1191 for(
int i=0; i<oldState.
curDim; i++) fullIndices[i] = i;
1199 std::vector<int> oldIndices;
1200 std::vector<int> newIndices;
1201 for(
int i=0; i<newState.
curDim; i++)
1203 if(indToCopy[i] < blockSize_)
1204 oldIndices.push_back(indToCopy[i]);
1206 newIndices.push_back(indToCopy[i]);
1209 int olddim = oldIndices.size();
1210 int newdim = newIndices.size();
1217 newState.
X = newState.
V;
1219 if(problem_->getOperator() != Teuchos::null)
1222 newState.
KX = newState.
KV;
1226 newState.
KV = Teuchos::null;
1227 newState.
KX = Teuchos::null;
1230 if(problem_->getM() != Teuchos::null)
1233 newState.
MX = newState.
MopV;
1237 newState.
MopV = Teuchos::null;
1238 newState.
MX = Teuchos::null;
1241 else if(newdim == 0)
1243 std::vector<int> blockind(blockSize_);
1244 for(
int i=0; i<blockSize_; i++)
1254 if(problem_->getM() != Teuchos::null)
1261 newState.
MopV = Teuchos::null;
1262 newState.
MX = Teuchos::null;
1268 std::vector<int> oldPart(olddim);
1269 for(
int i=0; i<olddim; i++) oldPart[i] = i;
1270 std::vector<int> newPart(newdim);
1271 for(
int i=0; i<newdim; i++) newPart[i] = olddim+i;
1277 RCP<const MV> viewHelper;
1280 Teuchos::SerialDenseMatrix<int,ScalarType> newRV(oldState.
curDim,newdim);
1281 for(
int r=0; r<oldState.
curDim; r++)
1283 for(
int c=0; c<newdim; c++)
1284 newRV(r,c) = (*oldState.
RV)(r,newIndices[c]);
1302 if(problem_->getM() != Teuchos::null)
1311 newState.
MopV = newState.
V;
1314 std::vector<int> blockVec(blockSize_);
1315 for(
int i=0; i<blockSize_; i++) blockVec[i] = i;
1321 if(blockSize_-oldIndices.size() > 0)
1324 newPart.resize(blockSize_-oldIndices.size());
1331 std::vector<ScalarType> scalarVec(blockSize_-oldIndices.size());
1332 for(
unsigned int i=0; i<(
unsigned int)blockSize_-oldIndices.size(); i++) scalarVec[i] = (*oldState.
T)[newPart[i]];
1344 newState.
R = oldState.
R;
1351 RCP< std::vector<ScalarType> > helperT = rcp(
new std::vector<ScalarType>(newState.
curDim) );
1352 for(
int i=0; i<newState.
curDim; i++) (*helperT)[i] = (*oldState.
T)[indToCopy[i]];
1353 newState.
T = helperT;
1356 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newKK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.
curDim,newState.
curDim) );
1357 for(
int i=0; i<newState.
curDim; i++)
1358 (*newKK)(i,i) = (*newState.
T)[i];
1359 newState.
KK = newKK;
1362 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newRV = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.
curDim,newState.
curDim) );
1363 for(
int i=0; i<newState.
curDim; i++)
1364 (*newRV)(i,i) = ONE;
1365 newState.
RV = newRV;
1368 RCP< std::vector<ScalarType> > helperRS = rcp(
new std::vector<ScalarType>(blockSize_) );
1369 for(
int i=0; i<blockSize_; i++)
1371 if(indToCopy[i] < blockSize_)
1372 (*helperRS)[i] = (*oldState.
ritzShifts)[indToCopy[i]];
1374 (*helperRS)[i] = ZERO;
1380 template <
class ScalarType,
class MV,
class OP>
1383 pl.set(
"Block Size", blockSize_);
1384 pl.set(
"Num Blocks", numBlocks_);
1385 pl.set(
"Num Restart Blocks", numRestartBlocks_);
1386 pl.set(
"When To Shift", whenToShift_);
1387 pl.set(
"Trace Threshold", traceThresh_);
1388 pl.set(
"Shift Tolerance", shiftTol_);
1389 pl.set(
"Relative Shift Tolerance", relShiftTol_);
1390 pl.set(
"Shift Norm", shiftNorm_);
1391 pl.set(
"How To Choose Shift", howToShift_);
1392 pl.set(
"Consider Clusters", considerClusters_);
1393 pl.set(
"Use Multiple Shifts", useMultipleShifts_);
1394 pl.set(
"Saddle Solver Type", saddleSolType_);
1395 pl.set(
"Project All Vectors", projectAllVecs_);
1396 pl.set(
"Project Locked Vectors", projectLockedVecs_);
1397 pl.set(
"Compute All Residuals", computeAllRes_);
1398 pl.set(
"Use Residual as RHS", useRHSR_);
1399 pl.set(
"Use Harmonic Ritz Values", useHarmonic_);
1400 pl.set(
"Maximum Krylov Iterations", maxKrylovIter_);
1401 pl.set(
"HSS: alpha", alpha_);
1405 template <
class ScalarType,
class MV,
class OP>
1409 os <<
"========================================\n";
1410 os <<
"========= TraceMin parameters ==========\n";
1411 os <<
"========================================\n";
1412 os <<
"=========== Block parameters ===========\n";
1413 os <<
"Block Size: " << blockSize_ << std::endl;
1414 os <<
"Num Blocks: " << numBlocks_ << std::endl;
1415 os <<
"Num Restart Blocks: " << numRestartBlocks_ << std::endl;
1416 os <<
"======== Convergence parameters ========\n";
1417 os <<
"Convergence Tolerance: " << convTol_ << std::endl;
1418 os <<
"Relative Convergence Tolerance: " << relConvTol_ << std::endl;
1419 os <<
"========== Locking parameters ==========\n";
1420 os <<
"Use Locking: " << useLocking_ << std::endl;
1421 os <<
"Locking Tolerance: " << lockTol_ << std::endl;
1422 os <<
"Relative Locking Tolerance: " << relLockTol_ << std::endl;
1423 os <<
"Max Locked: " << maxLocked_ << std::endl;
1424 os <<
"Locking Quorum: " << lockQuorum_ << std::endl;
1425 os <<
"========== Shifting parameters =========\n";
1426 os <<
"When To Shift: ";
1427 if(whenToShift_ == NEVER_SHIFT) os <<
"Never\n";
1428 else if(whenToShift_ == SHIFT_WHEN_TRACE_LEVELS) os <<
"After Trace Levels\n";
1429 else if(whenToShift_ == SHIFT_WHEN_RESID_SMALL) os <<
"Residual Becomes Small\n";
1430 else if(whenToShift_ == ALWAYS_SHIFT) os <<
"Always\n";
1431 os <<
"Consider Clusters: " << considerClusters_ << std::endl;
1432 os <<
"Trace Threshohld: " << traceThresh_ << std::endl;
1433 os <<
"Shift Tolerance: " << shiftTol_ << std::endl;
1434 os <<
"Relative Shift Tolerance: " << relShiftTol_ << std::endl;
1435 os <<
"How To Choose Shift: ";
1436 if(howToShift_ == LARGEST_CONVERGED_SHIFT) os <<
"Largest Converged\n";
1437 else if(howToShift_ == ADJUSTED_RITZ_SHIFT) os <<
"Adjusted Ritz Values\n";
1438 else if(howToShift_ == RITZ_VALUES_SHIFT) os <<
"Ritz Values\n";
1439 os <<
"Use Multiple Shifts: " << useMultipleShifts_ << std::endl;
1440 os <<
"=========== Other parameters ===========\n";
1441 os <<
"Orthogonalization: " << ortho_ << std::endl;
1442 os <<
"Saddle Solver Type: ";
1443 if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER) os <<
"Projected Krylov\n";
1444 else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER) os <<
"Schur Complement\n";
1445 os <<
"Project All Vectors: " << projectAllVecs_ << std::endl;
1446 os <<
"Project Locked Vectors: " << projectLockedVecs_ << std::endl;
1447 os <<
"Compute All Residuals: " << computeAllRes_ << std::endl;
1448 os <<
"========================================\n\n\n";
Pure virtual base class which describes the basic interface for a solver manager. ...
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
const RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
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.
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
const RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
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...
static void Assign(const MV &A, MV &mv)
mv := A
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.
Teuchos::Array< RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A status test for testing the norm of the eigenvectors residuals.
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
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.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
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.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
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.
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...
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.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Basic implementation of the Anasazi::OrthoManager class.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Class which provides internal utilities for the Anasazi solvers.
const RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.