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 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef ANASAZI_LOBPCG_SOLMGR_HPP
43 #define ANASAZI_LOBPCG_SOLMGR_HPP
44 
49 #include "AnasaziConfigDefs.hpp"
50 #include "AnasaziTypes.hpp"
51 
52 #include "AnasaziEigenproblem.hpp"
53 #include "AnasaziSolverManager.hpp"
54 #include "AnasaziSolverUtils.hpp"
55 
56 #include "AnasaziLOBPCG.hpp"
57 #include "AnasaziBasicSort.hpp"
65 #include "AnasaziOutputManager.hpp"
67 
68 #include "Teuchos_BLAS.hpp"
69 #include "Teuchos_TimeMonitor.hpp"
70 #include "Teuchos_FancyOStream.hpp"
71 
122 
123 namespace Anasazi {
124 
172 template<class ScalarType, class MV, class OP>
173 class LOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
174 
175  private:
179  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
181 
182  public:
183 
185 
186 
216 
218  virtual ~LOBPCGSolMgr() {};
220 
222 
223 
226  return *problem_;
227  }
228 
230  int getNumIters() const {
231  return numIters_;
232  }
233 
241  return Teuchos::tuple(_timerSolve, _timerLocking);
242  }
243 
244 
246 
248 
249 
276  ReturnType solve();
277 
280 
283 
286 
289 
292 
295 
297 
298  private:
300 
301  std::string whch_, ortho_;
302 
303  MagnitudeType convtol_, locktol_;
304  int maxIters_, numIters_;
305  bool useLocking_;
306  bool relconvtol_, rellocktol_;
307  int blockSize_;
308  bool fullOrtho_;
309  int maxLocked_;
310  int verbosity_;
311  int lockQuorum_;
312  bool recover_;
314  enum ResType convNorm_, lockNorm_;
315  int osProc_;
316 
317  Teuchos::RCP<Teuchos::Time> _timerSolve, _timerLocking;
318 
320 
324 };
325 
326 
327 // Constructor
328 template<class ScalarType, class MV, class OP>
331  Teuchos::ParameterList &pl ) :
332  problem_(problem),
333  whch_("SR"),
334  ortho_("SVQB"),
335  convtol_(MT::prec()),
336  maxIters_(100),
337  numIters_(0),
338  useLocking_(false),
339  relconvtol_(true),
340  rellocktol_(true),
341  blockSize_(0),
342  fullOrtho_(true),
343  maxLocked_(0),
344  verbosity_(Anasazi::Errors),
345  lockQuorum_(1),
346  recover_(true),
347  osProc_(0)
348 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
349  , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr::solve()")),
350  _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr locking"))
351 #endif
352 {
353  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
354  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
355  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
356  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
357 
358 
359  std::string strtmp;
360 
361  // which values to solve for
362  whch_ = pl.get("Which",whch_);
363  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
364  std::invalid_argument, "Anasazi::LOBPCGSolMgr: Invalid sorting string.");
365 
366  // which orthogonalization to use
367  ortho_ = pl.get("Orthogonalization",ortho_);
368  if (ortho_ != "DGKS" && ortho_ != "SVQB") {
369  ortho_ = "SVQB";
370  }
371 
372  // convergence tolerance
373  convtol_ = pl.get("Convergence Tolerance",convtol_);
374  relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
375  strtmp = pl.get("Convergence Norm",std::string("2"));
376  if (strtmp == "2") {
377  convNorm_ = RES_2NORM;
378  }
379  else if (strtmp == "M") {
380  convNorm_ = RES_ORTH;
381  }
382  else {
383  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
384  "Anasazi::LOBPCGSolMgr: Invalid Convergence Norm.");
385  }
386 
387 
388  // locking tolerance
389  useLocking_ = pl.get("Use Locking",useLocking_);
390  rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
391  // default: should be less than convtol_
392  locktol_ = convtol_/10;
393  locktol_ = pl.get("Locking Tolerance",locktol_);
394  strtmp = pl.get("Locking Norm",std::string("2"));
395  if (strtmp == "2") {
396  lockNorm_ = RES_2NORM;
397  }
398  else if (strtmp == "M") {
399  lockNorm_ = RES_ORTH;
400  }
401  else {
402  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
403  "Anasazi::LOBPCGSolMgr: Invalid Locking Norm.");
404  }
405 
406  // maximum number of iterations
407  maxIters_ = pl.get("Maximum Iterations",maxIters_);
408 
409  // block size: default is nev()
410  blockSize_ = pl.get("Block Size",problem_->getNEV());
411  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
412  "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
413 
414  // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
415  if (useLocking_) {
416  maxLocked_ = pl.get("Max Locked",problem_->getNEV());
417  }
418  else {
419  maxLocked_ = 0;
420  }
421  if (maxLocked_ == 0) {
422  useLocking_ = false;
423  }
424  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
425  "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
426  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
427  std::invalid_argument,
428  "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
429 
430  if (useLocking_) {
431  lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
432  TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
433  std::invalid_argument,
434  "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
435  }
436 
437  // full orthogonalization: default true
438  fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
439 
440  // Create a formatted output stream to print to.
441  // See if user requests output processor.
442  osProc_ = pl.get("Output Processor", osProc_);
443 
444  // If not passed in by user, it will be chosen based upon operator type.
445  if (pl.isParameter("Output Stream")) {
446  osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
447  }
448  else {
449  osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc_);
450  }
451 
452  // verbosity level
453  if (pl.isParameter("Verbosity")) {
454  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
455  verbosity_ = pl.get("Verbosity", verbosity_);
456  } else {
457  verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
458  }
459  }
460 
461  // recover from LOBPCGRitzFailure
462  recover_ = pl.get("Recover",recover_);
463 
464  // get (optionally) an initial state
465  if (pl.isParameter("Init")) {
466  state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
467  }
468 }
469 
470 
471 // solve()
472 template<class ScalarType, class MV, class OP>
475 
476  typedef SolverUtils<ScalarType,MV,OP> msutils;
477 
478  const int nev = problem_->getNEV();
479 
480 
481 
483  // Sort manager
485 
487  // Output manager
489 
491  // Status tests
492  //
493  // maximum number of iterations: optional test
495  if (maxIters_ > 0) {
496  maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
497  }
498  // convergence
500  if (globalTest_ == Teuchos::null) {
501  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
502  }
503  else {
504  convtest = globalTest_;
505  }
507  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
508  // locking
510  if (useLocking_) {
511  if (lockingTest_ == Teuchos::null) {
512  locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
513  }
514  else {
515  locktest = lockingTest_;
516  }
517  }
518  // for a non-short-circuited OR test, the order doesn't matter
520  alltests.push_back(ordertest);
521  if (locktest != Teuchos::null) alltests.push_back(locktest);
522  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
523  if (maxtest != Teuchos::null) alltests.push_back(maxtest);
526  // printing StatusTest
528  if ( printer->isVerbosity(Debug) ) {
529  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
530  }
531  else {
532  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
533  }
534 
536  // Orthomanager
538  if (ortho_=="SVQB") {
539  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
540  } else if (ortho_=="DGKS") {
541  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
542  } else {
543  TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
544  }
545 
547  // Parameter list
549  plist.set("Block Size",blockSize_);
550  plist.set("Full Ortho",fullOrtho_);
551 
553  // LOBPCG solver
555  = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
556  // set any auxiliary vectors defined in the problem
557  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
558  if (probauxvecs != Teuchos::null) {
559  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
560  }
561 
563  // Storage
564  //
565  // lockvecs will contain eigenvectors that have been determined "locked" by the status test
566  int curNumLocked = 0;
567  Teuchos::RCP<MV> lockvecs;
568  if (useLocking_) {
569  lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
570  }
571  std::vector<MagnitudeType> lockvals;
572  // workMV will be used as work space for LOBPCGRitzFailure recovery and locking
573  // it will be partitioned in these cases as follows:
574  // for LOBPCGRitzFailure recovery:
575  // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
576  // total size: 2*3*blocksize
577  // for locking
578  // workMV = [X P MX MP], with MX,MP needing storage only if hasM==true
579  // total size: 2*blocksize or 4*blocksize
580  Teuchos::RCP<MV> workMV;
581  if (fullOrtho_ == false && recover_ == true) {
582  workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
583  }
584  else if (useLocking_) {
585  if (problem_->getM() != Teuchos::null) {
586  workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
587  }
588  else {
589  workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
590  }
591  }
592 
593  // initialize the solution to nothing in case we throw an exception
595  sol.numVecs = 0;
596  problem_->setSolution(sol);
597 
598  // initialize the solver if the user specified a state
599  if (state_ != Teuchos::null) {
600  lobpcg_solver->initialize(*state_);
601  }
602 
603  // enter solve() iterations
604  {
605 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
606  Teuchos::TimeMonitor slvtimer(*_timerSolve);
607 #endif
608 
609  // tell the lobpcg_solver to iterate
610  while (1) {
611  try {
612  lobpcg_solver->iterate();
613 
615  //
616  // check user-specified debug test; if it passed, return an exception
617  //
619  if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
620  throw AnasaziError("Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
621  }
623  //
624  // check convergence first
625  //
627  else if (ordertest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
628  // we have convergence or not
629  // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
630  // ordertest->howMany() will tell us how many
631  break;
632  }
634  //
635  // check locking if we didn't converge
636  //
638  else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
639 
640 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
641  Teuchos::TimeMonitor lcktimer(*_timerLocking);
642 #endif
643 
644  // remove the locked vectors,values from lobpcg_solver: put them in newvecs, newvals
645  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
646  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
647  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
648  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
649  TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
650  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
651  // get the indices
652  int numnew = locktest->howMany();
653  std::vector<int> indnew = locktest->whichVecs();
654 
655  // don't lock more than maxLocked_; we didn't allocate enough space.
656  if (curNumLocked + numnew > maxLocked_) {
657  numnew = maxLocked_ - curNumLocked;
658  indnew.resize(numnew);
659  }
660 
661  // the call below to lobpcg_solver->setAuxVecs() will reset the solver to unitialized with hasP() == false
662  // store the hasP() state for use below
663  bool hadP = lobpcg_solver->hasP();
664 
665  {
666  // debug printing
667  printer->print(Debug,"Locking vectors: ");
668  for (unsigned int i=0; i<indnew.size(); i++) {printer->stream(Debug) << " " << indnew[i];}
669  printer->print(Debug,"\n");
670  }
671  std::vector<MagnitudeType> newvals(numnew);
672  Teuchos::RCP<const MV> newvecs;
673  {
674  // work in a local scope, to hide the variabes needed for extracting this info
675  // get the vectors
676  newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
677  // get the values
678  std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
679  for (int i=0; i<numnew; i++) {
680  newvals[i] = allvals[indnew[i]].realpart;
681  }
682  }
683  // put newvecs into lockvecs
684  {
685  std::vector<int> indlock(numnew);
686  for (int i=0; i<numnew; i++) indlock[i] = curNumLocked+i;
687  MVT::SetBlock(*newvecs,indlock,*lockvecs);
688  newvecs = Teuchos::null;
689  }
690  // put newvals into lockvals
691  lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
692  curNumLocked += numnew;
693  // add locked vecs as aux vecs, along with aux vecs from problem
694  {
695  std::vector<int> indlock(curNumLocked);
696  for (int i=0; i<curNumLocked; i++) indlock[i] = i;
697  Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
698  if (probauxvecs != Teuchos::null) {
699  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs,curlocked) );
700  }
701  else {
702  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(curlocked) );
703  }
704  }
705  // add locked vals to ordertest
706  ordertest->setAuxVals(lockvals);
707  // fill out the empty state in the solver
708  {
709  LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
710  Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
711  //
712  // workMV will be partitioned as follows: workMV = [X P MX MP],
713  //
714  // make a copy of the current X,MX state
715  std::vector<int> bsind(blockSize_);
716  for (int i=0; i<blockSize_; i++) bsind[i] = i;
717  newstateX = MVT::CloneViewNonConst(*workMV,bsind);
718  MVT::SetBlock(*state.X,bsind,*newstateX);
719 
720  if (state.MX != Teuchos::null) {
721  std::vector<int> block3(blockSize_);
722  for (int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
723  newstateMX = MVT::CloneViewNonConst(*workMV,block3);
724  MVT::SetBlock(*state.MX,bsind,*newstateMX);
725  }
726  //
727  // get select part, set to random, apply M
728  {
729  Teuchos::RCP<MV> newX = MVT::CloneViewNonConst(*newstateX,indnew);
730  MVT::MvRandom(*newX);
731 
732  if (newstateMX != Teuchos::null) {
733  Teuchos::RCP<MV> newMX = MVT::CloneViewNonConst(*newstateMX,indnew);
734  OPT::Apply(*problem_->getM(),*newX,*newMX);
735  }
736  }
737 
738  Teuchos::Array<Teuchos::RCP<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
740  // ortho X against the aux vectors
741  ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
742 
743  if (hadP) {
744  //
745  // get P and optionally MP, orthogonalize against X and auxiliary vectors
746  std::vector<int> block2(blockSize_);
747  for (int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
748  newstateP = MVT::CloneViewNonConst(*workMV,block2);
749  MVT::SetBlock(*state.P,bsind,*newstateP);
750 
751  if (state.MP != Teuchos::null) {
752  std::vector<int> block4(blockSize_);
753  for (int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
754  newstateMP = MVT::CloneViewNonConst(*workMV,block4);
755  MVT::SetBlock(*state.MP,bsind,*newstateMP);
756  }
757 
758  if (fullOrtho_) {
759  // ortho P against the new aux vectors and new X
760  curauxvecs.push_back(newstateX);
761  ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
762  }
763  else {
764  // ortho P against the new aux vectors
765  ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
766  }
767  }
768  // set the new state
770  newstate.X = newstateX;
771  newstate.MX = newstateMX;
772  newstate.P = newstateP;
773  newstate.MP = newstateMP;
774  lobpcg_solver->initialize(newstate);
775  }
776 
777  if (curNumLocked == maxLocked_) {
778  // disable locking now; remove locking test from combo test
779  combotest->removeTest(locktest);
780  }
781  }
782  else {
783  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
784  }
785  }
787  //
788  // check Ritz Failure
789  //
791  catch (const LOBPCGRitzFailure &re) {
792  if (fullOrtho_==true || recover_==false) {
793  // if we are already using full orthogonalization, there isn't much we can do here.
794  // the most recent information in the status tests is still valid, and can be used to extract/return the
795  // eigenpairs that have converged.
796  printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
797  << "Will not try to recover." << std::endl;
798  break; // while(1)
799  }
800  printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
801  << "Full orthogonalization is off; will try to recover." << std::endl;
802  // get the current "basis" from the solver, orthonormalize it, do a rayleigh-ritz, and restart with the ritz vectors
803  // if there aren't enough, break and quit with what we have
804  //
805  // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
806  LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
807  Teuchos::RCP<MV> restart, Krestart, Mrestart;
808  int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
809  bool hasM = problem_->getM() != Teuchos::null;
810  {
811  std::vector<int> recind(localsize);
812  for (int i=0; i<localsize; i++) recind[i] = i;
813  restart = MVT::CloneViewNonConst(*workMV,recind);
814  }
815  {
816  std::vector<int> recind(localsize);
817  for (int i=0; i<localsize; i++) recind[i] = localsize+i;
818  Krestart = MVT::CloneViewNonConst(*workMV,recind);
819  }
820  if (hasM) {
821  Mrestart = Krestart;
822  }
823  else {
824  Mrestart = restart;
825  }
826  //
827  // set restart = [X H P] and Mrestart = M*[X H P]
828  //
829  // put X into [0 , blockSize)
830  {
831  std::vector<int> blk1(blockSize_);
832  for (int i=0; i < blockSize_; i++) blk1[i] = i;
833  MVT::SetBlock(*curstate.X,blk1,*restart);
834 
835  // put MX into [0 , blockSize)
836  if (hasM) {
837  MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
838  }
839  }
840  //
841  // put H into [blockSize_ , 2*blockSize)
842  {
843  std::vector<int> blk2(blockSize_);
844  for (int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
845  MVT::SetBlock(*curstate.H,blk2,*restart);
846 
847  // put MH into [blockSize_ , 2*blockSize)
848  if (hasM) {
849  MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
850  }
851  }
852  // optionally, put P into [2*blockSize,3*blockSize)
853  if (localsize == 3*blockSize_) {
854  std::vector<int> blk3(blockSize_);
855  for (int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
856  MVT::SetBlock(*curstate.P,blk3,*restart);
857 
858  // put MP into [2*blockSize,3*blockSize)
859  if (hasM) {
860  MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
861  }
862  }
863  // project against auxvecs and locked vecs, and orthonormalize the basis
866  {
867  if (curNumLocked > 0) {
868  std::vector<int> indlock(curNumLocked);
869  for (int i=0; i<curNumLocked; i++) indlock[i] = i;
870  Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
871  Q.push_back(curlocked);
872  }
873  if (probauxvecs != Teuchos::null) {
874  Q.push_back(probauxvecs);
875  }
876  }
877  int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
878  if (rank < blockSize_) {
879  // quit
880  printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
881  << "Recovery failed." << std::endl;
882  break;
883  }
884  // reduce multivec size if necessary
885  if (rank < localsize) {
886  localsize = rank;
887  std::vector<int> redind(localsize);
888  for (int i=0; i<localsize; i++) redind[i] = i;
889  // grab the first part of restart,Krestart
890  restart = MVT::CloneViewNonConst(*restart,redind);
891  Krestart = MVT::CloneViewNonConst(*Krestart,redind);
892  if (hasM) {
893  Mrestart = Krestart;
894  }
895  else {
896  Mrestart = restart;
897  }
898  }
899  Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
900  std::vector<MagnitudeType> theta(localsize);
901  // project the matrices
902  //
903  // MM = restart^H M restart
904  MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
905  //
906  // compute Krestart = K*restart
907  OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
908  //
909  // KK = restart^H K restart
910  MVT::MvTransMv(1.0,*restart,*Krestart,KK);
911  rank = localsize;
912  msutils::directSolver(localsize,KK,Teuchos::rcpFromRef(MM),S,theta,rank,1);
913  if (rank < blockSize_) {
914  printer->stream(Errors) << "Error! Recovered basis of rank " << rank << " produced only " << rank << "ritz vectors.\n"
915  << "Block size is " << blockSize_ << ".\n"
916  << "Recovery failed." << std::endl;
917  break;
918  }
919  theta.resize(rank);
920  //
921  // sort the ritz values using the sort manager
922  {
924  std::vector<int> order(rank);
925  // sort
926  sorter->sort( theta, Teuchos::rcpFromRef(order),rank ); // don't catch exception
927  // Sort the primitive ritz vectors
929  msutils::permuteVectors(order,curS);
930  }
931  //
933  //
934  // compute the ritz vectors: store them in Krestart
936  Teuchos::RCP<MV> newX;
937  {
938  std::vector<int> bsind(blockSize_);
939  for (int i=0; i<blockSize_; i++) bsind[i] = i;
940  newX = MVT::CloneViewNonConst(*Krestart,bsind);
941  }
942  MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
943  // send X and theta into the solver
944  newstate.X = newX;
945  theta.resize(blockSize_);
946  newstate.T = Teuchos::rcpFromRef(theta);
947  // initialize
948  lobpcg_solver->initialize(newstate);
949  }
950  catch (const AnasaziError &err) {
951  printer->stream(Errors)
952  << "Anasazi::LOBPCGSolMgr::solve() caught unexpected exception from Anasazi::LOBPCG::iterate() at iteration " << lobpcg_solver->getNumIters() << std::endl
953  << err.what() << std::endl
954  << "Anasazi::LOBPCGSolMgr::solve() returning Unconverged with no solutions." << std::endl;
955  return Unconverged;
956  }
957  // don't catch any other exceptions
958  }
959 
960  sol.numVecs = ordertest->howMany();
961  if (sol.numVecs > 0) {
962  sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
963  sol.Espace = sol.Evecs;
964  sol.Evals.resize(sol.numVecs);
965  std::vector<MagnitudeType> vals(sol.numVecs);
966 
967  // copy them into the solution
968  std::vector<int> which = ordertest->whichVecs();
969  // indices between [0,blockSize) refer to vectors/values in the solver
970  // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
971  // everything has already been ordered by the solver; we just have to partition the two references
972  std::vector<int> inlocked(0), insolver(0);
973  for (unsigned int i=0; i<which.size(); i++) {
974  if (which[i] >= 0) {
975  TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): positive indexing mistake from ordertest.");
976  insolver.push_back(which[i]);
977  }
978  else {
979  // sanity check
980  TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): negative indexing mistake from ordertest.");
981  inlocked.push_back(which[i] + curNumLocked);
982  }
983  }
984 
985  TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
986 
987  // set the vecs,vals in the solution
988  if (insolver.size() > 0) {
989  // set vecs
990  int lclnum = insolver.size();
991  std::vector<int> tosol(lclnum);
992  for (int i=0; i<lclnum; i++) tosol[i] = i;
993  Teuchos::RCP<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
994  MVT::SetBlock(*v,tosol,*sol.Evecs);
995  // set vals
996  std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
997  for (unsigned int i=0; i<insolver.size(); i++) {
998  vals[i] = fromsolver[insolver[i]].realpart;
999  }
1000  }
1001 
1002  // get the vecs,vals from locked storage
1003  if (inlocked.size() > 0) {
1004  int solnum = insolver.size();
1005  // set vecs
1006  int lclnum = inlocked.size();
1007  std::vector<int> tosol(lclnum);
1008  for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1009  Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1010  MVT::SetBlock(*v,tosol,*sol.Evecs);
1011  // set vals
1012  for (unsigned int i=0; i<inlocked.size(); i++) {
1013  vals[i+solnum] = lockvals[inlocked[i]];
1014  }
1015  }
1016 
1017  // sort the eigenvalues and permute the eigenvectors appropriately
1018  {
1019  std::vector<int> order(sol.numVecs);
1020  sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
1021  // store the values in the Eigensolution
1022  for (int i=0; i<sol.numVecs; i++) {
1023  sol.Evals[i].realpart = vals[i];
1024  sol.Evals[i].imagpart = MT::zero();
1025  }
1026  // now permute the eigenvectors according to order
1027  msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1028  }
1029 
1030  // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1031  sol.index.resize(sol.numVecs,0);
1032  }
1033  }
1034 
1035  // print final summary
1036  lobpcg_solver->currentStatus(printer->stream(FinalSummary));
1037 
1038  // print timing information
1039 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1040  if ( printer->isVerbosity( TimingDetails ) ) {
1041  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
1042  }
1043 #endif
1044 
1045  problem_->setSolution(sol);
1046  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1047 
1048  // get the number of iterations performed in this call to solve.
1049  numIters_ = lobpcg_solver->getNumIters();
1050 
1051  if (sol.numVecs < nev) {
1052  return Unconverged; // return from LOBPCGSolMgr::solve()
1053  }
1054  return Converged; // return from LOBPCGSolMgr::solve()
1055 }
1056 
1057 
1058 template <class ScalarType, class MV, class OP>
1059 void
1061  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
1062 {
1063  globalTest_ = global;
1064 }
1065 
1066 template <class ScalarType, class MV, class OP>
1069 {
1070  return globalTest_;
1071 }
1072 
1073 template <class ScalarType, class MV, class OP>
1074 void
1077 {
1078  debugTest_ = debug;
1079 }
1080 
1081 template <class ScalarType, class MV, class OP>
1084 {
1085  return debugTest_;
1086 }
1087 
1088 template <class ScalarType, class MV, class OP>
1089 void
1091  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
1092 {
1093  lockingTest_ = locking;
1094 }
1095 
1096 template <class ScalarType, class MV, class OP>
1099 {
1100  return lockingTest_;
1101 }
1102 
1103 } // end Anasazi namespace
1104 
1105 #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)
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&#39;s iterate() routine until ...
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.
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.
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.