29 #ifndef ANASAZI_LOBPCG_SOLMGR_HPP 30 #define ANASAZI_LOBPCG_SOLMGR_HPP 156 template<
class ScalarType,
class MV,
class OP>
222 return Teuchos::tuple(_timerSolve, _timerLocking);
282 std::string whch_, ortho_;
284 MagnitudeType convtol_, locktol_;
285 int maxIters_, numIters_;
287 bool relconvtol_, rellocktol_;
295 enum ResType convNorm_, lockNorm_;
306 template<
class ScalarType,
class MV,
class OP>
313 convtol_(
MT::prec()),
325 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
326 , _timerSolve(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCGSolMgr::solve()")),
327 _timerLocking(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: LOBPCGSolMgr locking"))
333 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
339 whch_ = pl.
get(
"Which",whch_);
341 std::invalid_argument,
"Anasazi::LOBPCGSolMgr: Invalid sorting string.");
344 ortho_ = pl.
get(
"Orthogonalization",ortho_);
345 if (ortho_ !=
"DGKS" && ortho_ !=
"SVQB") {
350 convtol_ = pl.
get(
"Convergence Tolerance",convtol_);
351 relconvtol_ = pl.
get(
"Relative Convergence Tolerance",relconvtol_);
352 strtmp = pl.
get(
"Convergence Norm",std::string(
"2"));
354 convNorm_ = RES_2NORM;
356 else if (strtmp ==
"M") {
357 convNorm_ = RES_ORTH;
361 "Anasazi::LOBPCGSolMgr: Invalid Convergence Norm.");
366 useLocking_ = pl.
get(
"Use Locking",useLocking_);
367 rellocktol_ = pl.
get(
"Relative Locking Tolerance",rellocktol_);
369 locktol_ = convtol_/10;
370 locktol_ = pl.
get(
"Locking Tolerance",locktol_);
371 strtmp = pl.
get(
"Locking Norm",std::string(
"2"));
373 lockNorm_ = RES_2NORM;
375 else if (strtmp ==
"M") {
376 lockNorm_ = RES_ORTH;
380 "Anasazi::LOBPCGSolMgr: Invalid Locking Norm.");
384 maxIters_ = pl.
get(
"Maximum Iterations",maxIters_);
387 blockSize_ = pl.
get(
"Block Size",problem_->getNEV());
389 "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
393 maxLocked_ = pl.
get(
"Max Locked",problem_->getNEV());
398 if (maxLocked_ == 0) {
402 "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
404 std::invalid_argument,
405 "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
408 lockQuorum_ = pl.
get(
"Locking Quorum",lockQuorum_);
410 std::invalid_argument,
411 "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
415 fullOrtho_ = pl.
get(
"Full Ortho",fullOrtho_);
419 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
420 verbosity_ = pl.
get(
"Verbosity", verbosity_);
422 verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
427 recover_ = pl.
get(
"Recover",recover_);
431 state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,
"Init");
437 template<
class ScalarType,
class MV,
class OP>
443 const int nev = problem_->getNEV();
465 if (globalTest_ == Teuchos::null) {
469 convtest = globalTest_;
476 if (lockingTest_ == Teuchos::null) {
480 locktest = lockingTest_;
486 if (locktest != Teuchos::null) alltests.
push_back(locktest);
487 if (debugTest_ != Teuchos::null) alltests.
push_back(debugTest_);
488 if (maxtest != Teuchos::null) alltests.
push_back(maxtest);
493 if ( printer->isVerbosity(
Debug) ) {
503 if (ortho_==
"SVQB") {
505 }
else if (ortho_==
"DGKS") {
508 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!=
"SVQB"&&ortho_!=
"DGKS",std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
514 plist.
set(
"Block Size",blockSize_);
515 plist.
set(
"Full Ortho",fullOrtho_);
523 if (probauxvecs != Teuchos::null) {
531 int curNumLocked = 0;
534 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
536 std::vector<MagnitudeType> lockvals;
546 if (fullOrtho_ ==
false && recover_ ==
true) {
547 workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
549 else if (useLocking_) {
550 if (problem_->getM() != Teuchos::null) {
551 workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
554 workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
561 problem_->setSolution(sol);
564 if (state_ != Teuchos::null) {
565 lobpcg_solver->initialize(*state_);
570 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 577 lobpcg_solver->iterate();
584 if (debugTest_ != Teuchos::null && debugTest_->getStatus() ==
Passed) {
585 throw AnasaziError(
"Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
592 else if (ordertest->getStatus() ==
Passed || (maxtest != Teuchos::null && maxtest->getStatus() ==
Passed) ) {
603 else if (locktest != Teuchos::null && locktest->getStatus() ==
Passed) {
605 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 611 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
613 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
615 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
617 int numnew = locktest->howMany();
618 std::vector<int> indnew = locktest->whichVecs();
621 if (curNumLocked + numnew > maxLocked_) {
622 numnew = maxLocked_ - curNumLocked;
623 indnew.resize(numnew);
628 bool hadP = lobpcg_solver->hasP();
632 printer->print(
Debug,
"Locking vectors: ");
633 for (
unsigned int i=0; i<indnew.size(); i++) {printer->stream(
Debug) <<
" " << indnew[i];}
634 printer->print(
Debug,
"\n");
636 std::vector<MagnitudeType> newvals(numnew);
641 newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
643 std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
644 for (
int i=0; i<numnew; i++) {
645 newvals[i] = allvals[indnew[i]].realpart;
650 std::vector<int> indlock(numnew);
651 for (
int i=0; i<numnew; i++) indlock[i] = curNumLocked+i;
652 MVT::SetBlock(*newvecs,indlock,*lockvecs);
653 newvecs = Teuchos::null;
656 lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
657 curNumLocked += numnew;
660 std::vector<int> indlock(curNumLocked);
661 for (
int i=0; i<curNumLocked; i++) indlock[i] = i;
663 if (probauxvecs != Teuchos::null) {
671 ordertest->setAuxVals(lockvals);
680 std::vector<int> bsind(blockSize_);
681 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
682 newstateX = MVT::CloneViewNonConst(*workMV,bsind);
683 MVT::SetBlock(*state.
X,bsind,*newstateX);
685 if (state.
MX != Teuchos::null) {
686 std::vector<int> block3(blockSize_);
687 for (
int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
688 newstateMX = MVT::CloneViewNonConst(*workMV,block3);
689 MVT::SetBlock(*state.
MX,bsind,*newstateMX);
695 MVT::MvRandom(*newX);
697 if (newstateMX != Teuchos::null) {
699 OPT::Apply(*problem_->getM(),*newX,*newMX);
706 ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
711 std::vector<int> block2(blockSize_);
712 for (
int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
713 newstateP = MVT::CloneViewNonConst(*workMV,block2);
714 MVT::SetBlock(*state.
P,bsind,*newstateP);
716 if (state.
MP != Teuchos::null) {
717 std::vector<int> block4(blockSize_);
718 for (
int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
719 newstateMP = MVT::CloneViewNonConst(*workMV,block4);
720 MVT::SetBlock(*state.
MP,bsind,*newstateMP);
726 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
730 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
735 newstate.
X = newstateX;
736 newstate.
MX = newstateMX;
737 newstate.
P = newstateP;
738 newstate.
MP = newstateMP;
739 lobpcg_solver->initialize(newstate);
742 if (curNumLocked == maxLocked_) {
744 combotest->removeTest(locktest);
748 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
757 if (fullOrtho_==
true || recover_==
false) {
761 printer->stream(
Warnings) <<
"Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
762 <<
"Will not try to recover." << std::endl;
765 printer->stream(
Warnings) <<
"Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
766 <<
"Full orthogonalization is off; will try to recover." << std::endl;
773 int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
774 bool hasM = problem_->getM() != Teuchos::null;
776 std::vector<int> recind(localsize);
777 for (
int i=0; i<localsize; i++) recind[i] = i;
778 restart = MVT::CloneViewNonConst(*workMV,recind);
781 std::vector<int> recind(localsize);
782 for (
int i=0; i<localsize; i++) recind[i] = localsize+i;
783 Krestart = MVT::CloneViewNonConst(*workMV,recind);
796 std::vector<int> blk1(blockSize_);
797 for (
int i=0; i < blockSize_; i++) blk1[i] = i;
798 MVT::SetBlock(*curstate.
X,blk1,*restart);
802 MVT::SetBlock(*curstate.
MX,blk1,*Mrestart);
808 std::vector<int> blk2(blockSize_);
809 for (
int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
810 MVT::SetBlock(*curstate.
H,blk2,*restart);
814 MVT::SetBlock(*curstate.
MH,blk2,*Mrestart);
818 if (localsize == 3*blockSize_) {
819 std::vector<int> blk3(blockSize_);
820 for (
int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
821 MVT::SetBlock(*curstate.
P,blk3,*restart);
825 MVT::SetBlock(*curstate.
MP,blk3,*Mrestart);
832 if (curNumLocked > 0) {
833 std::vector<int> indlock(curNumLocked);
834 for (
int i=0; i<curNumLocked; i++) indlock[i] = i;
838 if (probauxvecs != Teuchos::null) {
842 int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
843 if (rank < blockSize_) {
845 printer->stream(
Errors) <<
"Error! Recovered basis only rank " << rank <<
". Block size is " << blockSize_ <<
".\n" 846 <<
"Recovery failed." << std::endl;
850 if (rank < localsize) {
852 std::vector<int> redind(localsize);
853 for (
int i=0; i<localsize; i++) redind[i] = i;
855 restart = MVT::CloneViewNonConst(*restart,redind);
856 Krestart = MVT::CloneViewNonConst(*Krestart,redind);
865 std::vector<MagnitudeType> theta(localsize);
869 MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
872 OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
875 MVT::MvTransMv(1.0,*restart,*Krestart,KK);
877 msutils::directSolver(localsize,KK,Teuchos::rcpFromRef(MM),S,theta,rank,1);
878 if (rank < blockSize_) {
879 printer->stream(
Errors) <<
"Error! Recovered basis of rank " << rank <<
" produced only " << rank <<
"ritz vectors.\n" 880 <<
"Block size is " << blockSize_ <<
".\n" 881 <<
"Recovery failed." << std::endl;
889 std::vector<int> order(rank);
891 sorter->sort( theta, Teuchos::rcpFromRef(order),rank );
894 msutils::permuteVectors(order,curS);
903 std::vector<int> bsind(blockSize_);
904 for (
int i=0; i<blockSize_; i++) bsind[i] = i;
905 newX = MVT::CloneViewNonConst(*Krestart,bsind);
907 MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
910 theta.resize(blockSize_);
911 newstate.
T = Teuchos::rcpFromRef(theta);
913 lobpcg_solver->initialize(newstate);
917 <<
"Anasazi::LOBPCGSolMgr::solve() caught unexpected exception from Anasazi::LOBPCG::iterate() at iteration " << lobpcg_solver->getNumIters() << std::endl
918 << err.what() << std::endl
919 <<
"Anasazi::LOBPCGSolMgr::solve() returning Unconverged with no solutions." << std::endl;
925 sol.
numVecs = ordertest->howMany();
927 sol.
Evecs = MVT::Clone(*problem_->getInitVec(),sol.
numVecs);
930 std::vector<MagnitudeType> vals(sol.
numVecs);
933 std::vector<int> which = ordertest->whichVecs();
937 std::vector<int> inlocked(0), insolver(0);
938 for (
unsigned int i=0; i<which.size(); i++) {
940 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): positive indexing mistake from ordertest.");
941 insolver.push_back(which[i]);
945 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,
"Anasazi::LOBPCGSolMgr::solve(): negative indexing mistake from ordertest.");
946 inlocked.push_back(which[i] + curNumLocked);
953 if (insolver.size() > 0) {
955 int lclnum = insolver.size();
956 std::vector<int> tosol(lclnum);
957 for (
int i=0; i<lclnum; i++) tosol[i] = i;
959 MVT::SetBlock(*v,tosol,*sol.
Evecs);
961 std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
962 for (
unsigned int i=0; i<insolver.size(); i++) {
963 vals[i] = fromsolver[insolver[i]].realpart;
968 if (inlocked.size() > 0) {
969 int solnum = insolver.size();
971 int lclnum = inlocked.size();
972 std::vector<int> tosol(lclnum);
973 for (
int i=0; i<lclnum; i++) tosol[i] = solnum + i;
975 MVT::SetBlock(*v,tosol,*sol.
Evecs);
977 for (
unsigned int i=0; i<inlocked.size(); i++) {
978 vals[i+solnum] = lockvals[inlocked[i]];
984 std::vector<int> order(sol.
numVecs);
985 sorter->sort( vals, Teuchos::rcpFromRef(order), sol.
numVecs);
987 for (
int i=0; i<sol.
numVecs; i++) {
988 sol.
Evals[i].realpart = vals[i];
989 sol.
Evals[i].imagpart = MT::zero();
1001 lobpcg_solver->currentStatus(printer->stream(
FinalSummary));
1004 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1010 problem_->setSolution(sol);
1011 printer->stream(
Debug) <<
"Returning " << sol.
numVecs <<
" eigenpairs to eigenproblem." << std::endl;
1014 numIters_ = lobpcg_solver->getNumIters();
1023 template <
class ScalarType,
class MV,
class OP>
1028 globalTest_ = global;
1031 template <
class ScalarType,
class MV,
class OP>
1038 template <
class ScalarType,
class MV,
class OP>
1046 template <
class ScalarType,
class MV,
class OP>
1053 template <
class ScalarType,
class MV,
class OP>
1058 lockingTest_ = locking;
1061 template <
class ScalarType,
class MV,
class OP>
1065 return lockingTest_;
Pure virtual base class which describes the basic interface for a solver manager. ...
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
A special StatusTest for printing other status tests.
This class defines the interface required by an eigensolver and status test class to compute solution...
LOBPCGSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for LOBPCGSolMgr.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration...
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Status test for forming logical combinations of other status tests.
Teuchos::RCP< const MultiVector > H
The current preconditioned residual vectors.
int getNumIters() const
Get the iteration count for the most recent call to solve().
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Teuchos::RCP< const MultiVector > P
The current search direction.
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.
Implementation of the locally-optimal block preconditioned conjugate gradient (LOBPCG) method...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
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.
int numVecs
The number of computed eigenpairs.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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...
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
ReturnType
Enumerated type used to pass back information from a solver manager.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
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...
Teuchos::RCP< const MultiVector > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
void setLockingStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
Teuchos::RCP< const MultiVector > MH
The image of the current preconditioned residual vectors under M, or Teuchos::null if M was not speci...
A status test for testing the number of iterations.
Status test for testing the number of iterations.
void push_back(const value_type &x)
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
LOBPCGRitzFailure is thrown when the LOBPCG solver is unable to continue a call to LOBPCG::iterate() ...
Status test for forming logical combinations of other status tests.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
bool isParameter(const std::string &name) const
virtual ~LOBPCGSolMgr()
Destructor.
Teuchos::RCP< const MultiVector > MP
The image of the current search direction under M, or Teuchos::null if M was not specified.
Common interface of stopping criteria for Anasazi's solvers.
A status test for testing the norm of the eigenvectors residuals.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
Basic implementation of the Anasazi::OrthoManager class.
User interface for the LOBPCG eigensolver.
Structure to contain pointers to Anasazi state variables.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
Teuchos::RCP< const MultiVector > X
The current eigenvectors.
Class which provides internal utilities for the Anasazi solvers.