Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AnasaziLOBPCGSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Anasazi: Block Eigensolvers Package
4 //
5 // Copyright 2004 NTESS and the Anasazi contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ANASAZI_LOBPCG_SOLMGR_HPP
11 #define ANASAZI_LOBPCG_SOLMGR_HPP
12 
17 #include "AnasaziConfigDefs.hpp"
18 #include "AnasaziTypes.hpp"
19 
20 #include "AnasaziEigenproblem.hpp"
21 #include "AnasaziSolverManager.hpp"
22 #include "AnasaziSolverUtils.hpp"
23 
24 #include "AnasaziLOBPCG.hpp"
25 #include "AnasaziBasicSort.hpp"
34 #include "AnasaziOutputManager.hpp"
36 
37 #include "Teuchos_BLAS.hpp"
38 #include "Teuchos_TimeMonitor.hpp"
39 #include "Teuchos_FancyOStream.hpp"
40 
91 
92 namespace Anasazi {
93 
141 template<class ScalarType, class MV, class OP>
142 class LOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
143 
144  private:
148  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
150 
151  public:
152 
154 
155 
185 
187  virtual ~LOBPCGSolMgr() {};
189 
191 
192 
195  return *problem_;
196  }
197 
199  int getNumIters() const {
200  return numIters_;
201  }
202 
210  return Teuchos::tuple(_timerSolve, _timerLocking);
211  }
212 
213 
215 
217 
218 
245  ReturnType solve();
246 
249 
252 
255 
258 
261 
264 
266 
267  private:
269 
270  std::string whch_, ortho_;
271 
272  MagnitudeType convtol_, locktol_;
273  int maxIters_, numIters_;
274  bool useLocking_;
275  bool relconvtol_, rellocktol_;
276  int blockSize_;
277  bool fullOrtho_;
278  int maxLocked_;
279  int verbosity_;
280  int lockQuorum_;
281  bool recover_;
283  enum ResType convNorm_, lockNorm_;
284  int osProc_;
285 
286  Teuchos::RCP<Teuchos::Time> _timerSolve, _timerLocking;
287 
289 
293 };
294 
295 
296 // Constructor
297 template<class ScalarType, class MV, class OP>
300  Teuchos::ParameterList &pl ) :
301  problem_(problem),
302  whch_("SR"),
303  ortho_("SVQB"),
304  convtol_(MT::prec()),
305  maxIters_(100),
306  numIters_(0),
307  useLocking_(false),
308  relconvtol_(true),
309  rellocktol_(true),
310  blockSize_(0),
311  fullOrtho_(true),
312  maxLocked_(0),
313  verbosity_(Anasazi::Errors),
314  lockQuorum_(1),
315  recover_(true),
316  osProc_(0)
317 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
318  , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr::solve()")),
319  _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr locking"))
320 #endif
321 {
322  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
323  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
324  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
325  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
326 
327 
328  std::string strtmp;
329 
330  // which values to solve for
331  whch_ = pl.get("Which",whch_);
332  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
333  std::invalid_argument, "Anasazi::LOBPCGSolMgr: Invalid sorting string.");
334 
335  // which orthogonalization to use
336  ortho_ = pl.get("Orthogonalization",ortho_);
337  if (ortho_ != "DGKS" && ortho_ != "SVQB" && ortho_ != "ICGS") {
338  ortho_ = "SVQB";
339  }
340 
341  // convergence tolerance
342  convtol_ = pl.get("Convergence Tolerance",convtol_);
343  relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
344  strtmp = pl.get("Convergence Norm",std::string("2"));
345  if (strtmp == "2") {
346  convNorm_ = RES_2NORM;
347  }
348  else if (strtmp == "M") {
349  convNorm_ = RES_ORTH;
350  }
351  else {
352  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
353  "Anasazi::LOBPCGSolMgr: Invalid Convergence Norm.");
354  }
355 
356 
357  // locking tolerance
358  useLocking_ = pl.get("Use Locking",useLocking_);
359  rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
360  // default: should be less than convtol_
361  locktol_ = convtol_/10;
362  locktol_ = pl.get("Locking Tolerance",locktol_);
363  strtmp = pl.get("Locking Norm",std::string("2"));
364  if (strtmp == "2") {
365  lockNorm_ = RES_2NORM;
366  }
367  else if (strtmp == "M") {
368  lockNorm_ = RES_ORTH;
369  }
370  else {
371  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
372  "Anasazi::LOBPCGSolMgr: Invalid Locking Norm.");
373  }
374 
375  // maximum number of iterations
376  maxIters_ = pl.get("Maximum Iterations",maxIters_);
377 
378  // block size: default is nev()
379  blockSize_ = pl.get("Block Size",problem_->getNEV());
380  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
381  "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
382 
383  // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
384  if (useLocking_) {
385  maxLocked_ = pl.get("Max Locked",problem_->getNEV());
386  }
387  else {
388  maxLocked_ = 0;
389  }
390  if (maxLocked_ == 0) {
391  useLocking_ = false;
392  }
393  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
394  "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
395  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
396  std::invalid_argument,
397  "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
398 
399  if (useLocking_) {
400  lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
401  TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
402  std::invalid_argument,
403  "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
404  }
405 
406  // full orthogonalization: default true
407  fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
408 
409  // Create a formatted output stream to print to.
410  // See if user requests output processor.
411  osProc_ = pl.get("Output Processor", osProc_);
412 
413  // If not passed in by user, it will be chosen based upon operator type.
414  if (pl.isParameter("Output Stream")) {
415  osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
416  }
417  else {
418  osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc_);
419  }
420 
421  // verbosity level
422  if (pl.isParameter("Verbosity")) {
423  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
424  verbosity_ = pl.get("Verbosity", verbosity_);
425  } else {
426  verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
427  }
428  }
429 
430  // recover from LOBPCGRitzFailure
431  recover_ = pl.get("Recover",recover_);
432 
433  // get (optionally) an initial state
434  if (pl.isParameter("Init")) {
435  state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
436  }
437 }
438 
439 
440 // solve()
441 template<class ScalarType, class MV, class OP>
444 
445  typedef SolverUtils<ScalarType,MV,OP> msutils;
446 
447  const int nev = problem_->getNEV();
448 
449 
450 
452  // Sort manager
454 
456  // Output manager
458 
460  // Status tests
461  //
462  // maximum number of iterations: optional test
464  if (maxIters_ > 0) {
465  maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
466  }
467  // convergence
469  if (globalTest_ == Teuchos::null) {
470  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
471  }
472  else {
473  convtest = globalTest_;
474  }
476  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
477  // locking
479  if (useLocking_) {
480  if (lockingTest_ == Teuchos::null) {
481  locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
482  }
483  else {
484  locktest = lockingTest_;
485  }
486  }
487  // for a non-short-circuited OR test, the order doesn't matter
489  alltests.push_back(ordertest);
490  if (locktest != Teuchos::null) alltests.push_back(locktest);
491  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
492  if (maxtest != Teuchos::null) alltests.push_back(maxtest);
495  // printing StatusTest
497  if ( printer->isVerbosity(Debug) ) {
498  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
499  }
500  else {
501  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
502  }
503 
505  // Orthomanager
507  if (ortho_=="SVQB") {
508  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
509  } else if (ortho_=="DGKS") {
510  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
511  } else if (ortho_=="ICGS") {
512  ortho = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
513  } else {
514  TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS"&&ortho_!="ICGS",std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
515  }
516 
518  // Parameter list
520  plist.set("Block Size",blockSize_);
521  plist.set("Full Ortho",fullOrtho_);
522 
524  // LOBPCG solver
526  = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
527  // set any auxiliary vectors defined in the problem
528  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
529  if (probauxvecs != Teuchos::null) {
530  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
531  }
532 
534  // Storage
535  //
536  // lockvecs will contain eigenvectors that have been determined "locked" by the status test
537  int curNumLocked = 0;
538  Teuchos::RCP<MV> lockvecs;
539  if (useLocking_) {
540  lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
541  }
542  std::vector<MagnitudeType> lockvals;
543  // workMV will be used as work space for LOBPCGRitzFailure recovery and locking
544  // it will be partitioned in these cases as follows:
545  // for LOBPCGRitzFailure recovery:
546  // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
547  // total size: 2*3*blocksize
548  // for locking
549  // workMV = [X P MX MP], with MX,MP needing storage only if hasM==true
550  // total size: 2*blocksize or 4*blocksize
551  Teuchos::RCP<MV> workMV;
552  if (fullOrtho_ == false && recover_ == true) {
553  workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
554  }
555  else if (useLocking_) {
556  if (problem_->getM() != Teuchos::null) {
557  workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
558  }
559  else {
560  workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
561  }
562  }
563 
564  // initialize the solution to nothing in case we throw an exception
566  sol.numVecs = 0;
567  problem_->setSolution(sol);
568 
569  // initialize the solver if the user specified a state
570  if (state_ != Teuchos::null) {
571  lobpcg_solver->initialize(*state_);
572  }
573 
574  // enter solve() iterations
575  {
576 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
577  Teuchos::TimeMonitor slvtimer(*_timerSolve);
578 #endif
579 
580  // tell the lobpcg_solver to iterate
581  while (1) {
582  try {
583  lobpcg_solver->iterate();
584 
586  //
587  // check user-specified debug test; if it passed, return an exception
588  //
590  if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
591  throw AnasaziError("Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
592  }
594  //
595  // check convergence first
596  //
598  else if (ordertest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
599  // we have convergence or not
600  // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
601  // ordertest->howMany() will tell us how many
602  break;
603  }
605  //
606  // check locking if we didn't converge
607  //
609  else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
610 
611 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
612  Teuchos::TimeMonitor lcktimer(*_timerLocking);
613 #endif
614 
615  // remove the locked vectors,values from lobpcg_solver: put them in newvecs, newvals
616  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
617  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
618  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
619  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
620  TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
621  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
622  // get the indices
623  int numnew = locktest->howMany();
624  std::vector<int> indnew = locktest->whichVecs();
625 
626  // don't lock more than maxLocked_; we didn't allocate enough space.
627  if (curNumLocked + numnew > maxLocked_) {
628  numnew = maxLocked_ - curNumLocked;
629  indnew.resize(numnew);
630  }
631 
632  // the call below to lobpcg_solver->setAuxVecs() will reset the solver to unitialized with hasP() == false
633  // store the hasP() state for use below
634  bool hadP = lobpcg_solver->hasP();
635 
636  {
637  // debug printing
638  printer->print(Debug,"Locking vectors: ");
639  for (unsigned int i=0; i<indnew.size(); i++) {printer->stream(Debug) << " " << indnew[i];}
640  printer->print(Debug,"\n");
641  }
642  std::vector<MagnitudeType> newvals(numnew);
643  Teuchos::RCP<const MV> newvecs;
644  {
645  // work in a local scope, to hide the variabes needed for extracting this info
646  // get the vectors
647  newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
648  // get the values
649  std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
650  for (int i=0; i<numnew; i++) {
651  newvals[i] = allvals[indnew[i]].realpart;
652  }
653  }
654  // put newvecs into lockvecs
655  {
656  std::vector<int> indlock(numnew);
657  for (int i=0; i<numnew; i++) indlock[i] = curNumLocked+i;
658  MVT::SetBlock(*newvecs,indlock,*lockvecs);
659  newvecs = Teuchos::null;
660  }
661  // put newvals into lockvals
662  lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
663  curNumLocked += numnew;
664  // add locked vecs as aux vecs, along with aux vecs from problem
665  {
666  std::vector<int> indlock(curNumLocked);
667  for (int i=0; i<curNumLocked; i++) indlock[i] = i;
668  Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
669  if (probauxvecs != Teuchos::null) {
670  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs,curlocked) );
671  }
672  else {
673  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(curlocked) );
674  }
675  }
676  // add locked vals to ordertest
677  ordertest->setAuxVals(lockvals);
678  // fill out the empty state in the solver
679  {
680  LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
681  Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
682  //
683  // workMV will be partitioned as follows: workMV = [X P MX MP],
684  //
685  // make a copy of the current X,MX state
686  std::vector<int> bsind(blockSize_);
687  for (int i=0; i<blockSize_; i++) bsind[i] = i;
688  newstateX = MVT::CloneViewNonConst(*workMV,bsind);
689  MVT::SetBlock(*state.X,bsind,*newstateX);
690 
691  if (state.MX != Teuchos::null) {
692  std::vector<int> block3(blockSize_);
693  for (int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
694  newstateMX = MVT::CloneViewNonConst(*workMV,block3);
695  MVT::SetBlock(*state.MX,bsind,*newstateMX);
696  }
697  //
698  // get select part, set to random, apply M
699  {
700  Teuchos::RCP<MV> newX = MVT::CloneViewNonConst(*newstateX,indnew);
701  MVT::MvRandom(*newX);
702 
703  if (newstateMX != Teuchos::null) {
704  Teuchos::RCP<MV> newMX = MVT::CloneViewNonConst(*newstateMX,indnew);
705  OPT::Apply(*problem_->getM(),*newX,*newMX);
706  }
707  }
708 
709  Teuchos::Array<Teuchos::RCP<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
711  // ortho X against the aux vectors
712  ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
713 
714  if (hadP) {
715  //
716  // get P and optionally MP, orthogonalize against X and auxiliary vectors
717  std::vector<int> block2(blockSize_);
718  for (int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
719  newstateP = MVT::CloneViewNonConst(*workMV,block2);
720  MVT::SetBlock(*state.P,bsind,*newstateP);
721 
722  if (state.MP != Teuchos::null) {
723  std::vector<int> block4(blockSize_);
724  for (int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
725  newstateMP = MVT::CloneViewNonConst(*workMV,block4);
726  MVT::SetBlock(*state.MP,bsind,*newstateMP);
727  }
728 
729  if (fullOrtho_) {
730  // ortho P against the new aux vectors and new X
731  curauxvecs.push_back(newstateX);
732  ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
733  }
734  else {
735  // ortho P against the new aux vectors
736  ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
737  }
738  }
739  // set the new state
741  newstate.X = newstateX;
742  newstate.MX = newstateMX;
743  newstate.P = newstateP;
744  newstate.MP = newstateMP;
745  lobpcg_solver->initialize(newstate);
746  }
747 
748  if (curNumLocked == maxLocked_) {
749  // disable locking now; remove locking test from combo test
750  combotest->removeTest(locktest);
751  }
752  }
753  else {
754  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
755  }
756  }
758  //
759  // check Ritz Failure
760  //
762  catch (const LOBPCGRitzFailure &re) {
763  if (fullOrtho_==true || recover_==false) {
764  // if we are already using full orthogonalization, there isn't much we can do here.
765  // the most recent information in the status tests is still valid, and can be used to extract/return the
766  // eigenpairs that have converged.
767  printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
768  << "Will not try to recover." << std::endl;
769  break; // while(1)
770  }
771  printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
772  << "Full orthogonalization is off; will try to recover." << std::endl;
773  // get the current "basis" from the solver, orthonormalize it, do a rayleigh-ritz, and restart with the ritz vectors
774  // if there aren't enough, break and quit with what we have
775  //
776  // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
777  LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
778  Teuchos::RCP<MV> restart, Krestart, Mrestart;
779  int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
780  bool hasM = problem_->getM() != Teuchos::null;
781  {
782  std::vector<int> recind(localsize);
783  for (int i=0; i<localsize; i++) recind[i] = i;
784  restart = MVT::CloneViewNonConst(*workMV,recind);
785  }
786  {
787  std::vector<int> recind(localsize);
788  for (int i=0; i<localsize; i++) recind[i] = localsize+i;
789  Krestart = MVT::CloneViewNonConst(*workMV,recind);
790  }
791  if (hasM) {
792  Mrestart = Krestart;
793  }
794  else {
795  Mrestart = restart;
796  }
797  //
798  // set restart = [X H P] and Mrestart = M*[X H P]
799  //
800  // put X into [0 , blockSize)
801  {
802  std::vector<int> blk1(blockSize_);
803  for (int i=0; i < blockSize_; i++) blk1[i] = i;
804  MVT::SetBlock(*curstate.X,blk1,*restart);
805 
806  // put MX into [0 , blockSize)
807  if (hasM) {
808  MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
809  }
810  }
811  //
812  // put H into [blockSize_ , 2*blockSize)
813  {
814  std::vector<int> blk2(blockSize_);
815  for (int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
816  MVT::SetBlock(*curstate.H,blk2,*restart);
817 
818  // put MH into [blockSize_ , 2*blockSize)
819  if (hasM) {
820  MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
821  }
822  }
823  // optionally, put P into [2*blockSize,3*blockSize)
824  if (localsize == 3*blockSize_) {
825  std::vector<int> blk3(blockSize_);
826  for (int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
827  MVT::SetBlock(*curstate.P,blk3,*restart);
828 
829  // put MP into [2*blockSize,3*blockSize)
830  if (hasM) {
831  MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
832  }
833  }
834  // project against auxvecs and locked vecs, and orthonormalize the basis
837  {
838  if (curNumLocked > 0) {
839  std::vector<int> indlock(curNumLocked);
840  for (int i=0; i<curNumLocked; i++) indlock[i] = i;
841  Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
842  Q.push_back(curlocked);
843  }
844  if (probauxvecs != Teuchos::null) {
845  Q.push_back(probauxvecs);
846  }
847  }
848  int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
849  if (rank < blockSize_) {
850  // quit
851  printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
852  << "Recovery failed." << std::endl;
853  break;
854  }
855  // reduce multivec size if necessary
856  if (rank < localsize) {
857  localsize = rank;
858  std::vector<int> redind(localsize);
859  for (int i=0; i<localsize; i++) redind[i] = i;
860  // grab the first part of restart,Krestart
861  restart = MVT::CloneViewNonConst(*restart,redind);
862  Krestart = MVT::CloneViewNonConst(*Krestart,redind);
863  if (hasM) {
864  Mrestart = Krestart;
865  }
866  else {
867  Mrestart = restart;
868  }
869  }
870  Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
871  std::vector<MagnitudeType> theta(localsize);
872  // project the matrices
873  //
874  // MM = restart^H M restart
875  MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
876  //
877  // compute Krestart = K*restart
878  OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
879  //
880  // KK = restart^H K restart
881  MVT::MvTransMv(1.0,*restart,*Krestart,KK);
882  rank = localsize;
883  msutils::directSolver(localsize,KK,Teuchos::rcpFromRef(MM),S,theta,rank,1);
884  if (rank < blockSize_) {
885  printer->stream(Errors) << "Error! Recovered basis of rank " << rank << " produced only " << rank << "ritz vectors.\n"
886  << "Block size is " << blockSize_ << ".\n"
887  << "Recovery failed." << std::endl;
888  break;
889  }
890  theta.resize(rank);
891  //
892  // sort the ritz values using the sort manager
893  {
895  std::vector<int> order(rank);
896  // sort
897  sorter->sort( theta, Teuchos::rcpFromRef(order),rank ); // don't catch exception
898  // Sort the primitive ritz vectors
900  msutils::permuteVectors(order,curS);
901  }
902  //
904  //
905  // compute the ritz vectors: store them in Krestart
907  Teuchos::RCP<MV> newX;
908  {
909  std::vector<int> bsind(blockSize_);
910  for (int i=0; i<blockSize_; i++) bsind[i] = i;
911  newX = MVT::CloneViewNonConst(*Krestart,bsind);
912  }
913  MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
914  // send X and theta into the solver
915  newstate.X = newX;
916  theta.resize(blockSize_);
917  newstate.T = Teuchos::rcpFromRef(theta);
918  // initialize
919  lobpcg_solver->initialize(newstate);
920  }
921  catch (const AnasaziError &err) {
922  printer->stream(Errors)
923  << "Anasazi::LOBPCGSolMgr::solve() caught unexpected exception from Anasazi::LOBPCG::iterate() at iteration " << lobpcg_solver->getNumIters() << std::endl
924  << err.what() << std::endl
925  << "Anasazi::LOBPCGSolMgr::solve() returning Unconverged with no solutions." << std::endl;
926  return Unconverged;
927  }
928  // don't catch any other exceptions
929  }
930 
931  sol.numVecs = ordertest->howMany();
932  if (sol.numVecs > 0) {
933  sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
934  sol.Espace = sol.Evecs;
935  sol.Evals.resize(sol.numVecs);
936  std::vector<MagnitudeType> vals(sol.numVecs);
937 
938  // copy them into the solution
939  std::vector<int> which = ordertest->whichVecs();
940  // indices between [0,blockSize) refer to vectors/values in the solver
941  // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
942  // everything has already been ordered by the solver; we just have to partition the two references
943  std::vector<int> inlocked(0), insolver(0);
944  for (unsigned int i=0; i<which.size(); i++) {
945  if (which[i] >= 0) {
946  TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): positive indexing mistake from ordertest.");
947  insolver.push_back(which[i]);
948  }
949  else {
950  // sanity check
951  TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): negative indexing mistake from ordertest.");
952  inlocked.push_back(which[i] + curNumLocked);
953  }
954  }
955 
956  TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
957 
958  // set the vecs,vals in the solution
959  if (insolver.size() > 0) {
960  // set vecs
961  int lclnum = insolver.size();
962  std::vector<int> tosol(lclnum);
963  for (int i=0; i<lclnum; i++) tosol[i] = i;
964  Teuchos::RCP<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
965  MVT::SetBlock(*v,tosol,*sol.Evecs);
966  // set vals
967  std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
968  for (unsigned int i=0; i<insolver.size(); i++) {
969  vals[i] = fromsolver[insolver[i]].realpart;
970  }
971  }
972 
973  // get the vecs,vals from locked storage
974  if (inlocked.size() > 0) {
975  int solnum = insolver.size();
976  // set vecs
977  int lclnum = inlocked.size();
978  std::vector<int> tosol(lclnum);
979  for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
980  Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
981  MVT::SetBlock(*v,tosol,*sol.Evecs);
982  // set vals
983  for (unsigned int i=0; i<inlocked.size(); i++) {
984  vals[i+solnum] = lockvals[inlocked[i]];
985  }
986  }
987 
988  // sort the eigenvalues and permute the eigenvectors appropriately
989  {
990  std::vector<int> order(sol.numVecs);
991  sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
992  // store the values in the Eigensolution
993  for (int i=0; i<sol.numVecs; i++) {
994  sol.Evals[i].realpart = vals[i];
995  sol.Evals[i].imagpart = MT::zero();
996  }
997  // now permute the eigenvectors according to order
998  msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
999  }
1000 
1001  // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1002  sol.index.resize(sol.numVecs,0);
1003  }
1004  }
1005 
1006  // print final summary
1007  lobpcg_solver->currentStatus(printer->stream(FinalSummary));
1008 
1009  // print timing information
1010 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1011  if ( printer->isVerbosity( TimingDetails ) ) {
1012  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
1013  }
1014 #endif
1015 
1016  problem_->setSolution(sol);
1017  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1018 
1019  // get the number of iterations performed in this call to solve.
1020  numIters_ = lobpcg_solver->getNumIters();
1021 
1022  if (sol.numVecs < nev) {
1023  return Unconverged; // return from LOBPCGSolMgr::solve()
1024  }
1025  return Converged; // return from LOBPCGSolMgr::solve()
1026 }
1027 
1028 
1029 template <class ScalarType, class MV, class OP>
1030 void
1032  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
1033 {
1034  globalTest_ = global;
1035 }
1036 
1037 template <class ScalarType, class MV, class OP>
1040 {
1041  return globalTest_;
1042 }
1043 
1044 template <class ScalarType, class MV, class OP>
1045 void
1048 {
1049  debugTest_ = debug;
1050 }
1051 
1052 template <class ScalarType, class MV, class OP>
1055 {
1056  return debugTest_;
1057 }
1058 
1059 template <class ScalarType, class MV, class OP>
1060 void
1062  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
1063 {
1064  lockingTest_ = locking;
1065 }
1066 
1067 template <class ScalarType, class MV, class OP>
1070 {
1071  return lockingTest_;
1072 }
1073 
1074 } // end Anasazi namespace
1075 
1076 #endif /* ANASAZI_LOBPCG_SOLMGR_HPP */
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.
const Teuchos::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...
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)
#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&#39;s iterate() routine until ...
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
Status test for forming logical combinations of other status tests.
Teuchos::RCP< const MultiVector > H
The current preconditioned residual vectors.
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.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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.
Output managers remove the need for the eigensolver to know any information about the required output...
Implementation of the locally-optimal block preconditioned conjugate gradient (LOBPCG) method...
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&#39;s templated, static class providing utilities for the solvers.
bool isParameter(const std::string &name) const
int numVecs
The number of computed eigenpairs.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
Abstract class definition for Anasazi Output Managers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
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.
Output managers remove the need for the eigensolver to know any information about the required output...
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 &quot;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.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
Abstract class definition for Anasazi output stream.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
virtual ~LOBPCGSolMgr()
Destructor.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< const MultiVector > MP
The image of the current search direction under M, or Teuchos::null if M was not specified.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Common interface of stopping criteria for Anasazi&#39;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.
User interface for the LOBPCG eigensolver.
Structure to contain pointers to Anasazi state variables.
Teuchos::RCP< const MultiVector > X
The current eigenvectors.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Class which provides internal utilities for the Anasazi solvers.