Anasazi  Version of the Day
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, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_LOBPCG_SOLMGR_HPP
30 #define ANASAZI_LOBPCG_SOLMGR_HPP
31 
36 #include "AnasaziConfigDefs.hpp"
37 #include "AnasaziTypes.hpp"
38 
39 #include "AnasaziEigenproblem.hpp"
40 #include "AnasaziSolverManager.hpp"
41 #include "AnasaziSolverUtils.hpp"
42 
43 #include "AnasaziLOBPCG.hpp"
44 #include "AnasaziBasicSort.hpp"
53 #include "Teuchos_BLAS.hpp"
54 #include "Teuchos_TimeMonitor.hpp"
55 
106 
107 namespace Anasazi {
108 
156 template<class ScalarType, class MV, class OP>
157 class LOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
158 
159  private:
162  typedef Teuchos::ScalarTraits<ScalarType> SCT;
163  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
164  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
165 
166  public:
167 
169 
170 
195  LOBPCGSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
196  Teuchos::ParameterList &pl );
197 
199  virtual ~LOBPCGSolMgr() {};
201 
203 
204 
207  return *problem_;
208  }
209 
211  int getNumIters() const {
212  return numIters_;
213  }
214 
221  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
222  return Teuchos::tuple(_timerSolve, _timerLocking);
223  }
224 
225 
227 
229 
230 
257  ReturnType solve();
258 
260  void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
261 
263  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
264 
266  void setLockingStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking);
267 
269  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
270 
272  void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
273 
275  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
276 
278 
279  private:
280  Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
281 
282  std::string whch_, ortho_;
283 
284  MagnitudeType convtol_, locktol_;
285  int maxIters_, numIters_;
286  bool useLocking_;
287  bool relconvtol_, rellocktol_;
288  int blockSize_;
289  bool fullOrtho_;
290  int maxLocked_;
291  int verbosity_;
292  int lockQuorum_;
293  bool recover_;
294  Teuchos::RCP<LOBPCGState<ScalarType,MV> > state_;
295  enum ResType convNorm_, lockNorm_;
296 
297  Teuchos::RCP<Teuchos::Time> _timerSolve, _timerLocking;
298 
299  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
300  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
301  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
302 };
303 
304 
305 // Constructor
306 template<class ScalarType, class MV, class OP>
308  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
309  Teuchos::ParameterList &pl ) :
310  problem_(problem),
311  whch_("SR"),
312  ortho_("SVQB"),
313  convtol_(MT::prec()),
314  maxIters_(100),
315  numIters_(0),
316  useLocking_(false),
317  relconvtol_(true),
318  rellocktol_(true),
319  blockSize_(0),
320  fullOrtho_(true),
321  maxLocked_(0),
322  verbosity_(Anasazi::Errors),
323  lockQuorum_(1),
324  recover_(true)
325 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
326  , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr::solve()")),
327  _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr locking"))
328 #endif
329 {
330  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
331  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
332  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
333  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
334 
335 
336  std::string strtmp;
337 
338  // which values to solve for
339  whch_ = pl.get("Which",whch_);
340  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
341  std::invalid_argument, "Anasazi::LOBPCGSolMgr: Invalid sorting string.");
342 
343  // which orthogonalization to use
344  ortho_ = pl.get("Orthogonalization",ortho_);
345  if (ortho_ != "DGKS" && ortho_ != "SVQB") {
346  ortho_ = "SVQB";
347  }
348 
349  // convergence tolerance
350  convtol_ = pl.get("Convergence Tolerance",convtol_);
351  relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
352  strtmp = pl.get("Convergence Norm",std::string("2"));
353  if (strtmp == "2") {
354  convNorm_ = RES_2NORM;
355  }
356  else if (strtmp == "M") {
357  convNorm_ = RES_ORTH;
358  }
359  else {
360  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
361  "Anasazi::LOBPCGSolMgr: Invalid Convergence Norm.");
362  }
363 
364 
365  // locking tolerance
366  useLocking_ = pl.get("Use Locking",useLocking_);
367  rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
368  // default: should be less than convtol_
369  locktol_ = convtol_/10;
370  locktol_ = pl.get("Locking Tolerance",locktol_);
371  strtmp = pl.get("Locking Norm",std::string("2"));
372  if (strtmp == "2") {
373  lockNorm_ = RES_2NORM;
374  }
375  else if (strtmp == "M") {
376  lockNorm_ = RES_ORTH;
377  }
378  else {
379  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
380  "Anasazi::LOBPCGSolMgr: Invalid Locking Norm.");
381  }
382 
383  // maximum number of iterations
384  maxIters_ = pl.get("Maximum Iterations",maxIters_);
385 
386  // block size: default is nev()
387  blockSize_ = pl.get("Block Size",problem_->getNEV());
388  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
389  "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
390 
391  // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
392  if (useLocking_) {
393  maxLocked_ = pl.get("Max Locked",problem_->getNEV());
394  }
395  else {
396  maxLocked_ = 0;
397  }
398  if (maxLocked_ == 0) {
399  useLocking_ = false;
400  }
401  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
402  "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
403  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
404  std::invalid_argument,
405  "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
406 
407  if (useLocking_) {
408  lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
409  TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
410  std::invalid_argument,
411  "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
412  }
413 
414  // full orthogonalization: default true
415  fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
416 
417  // verbosity level
418  if (pl.isParameter("Verbosity")) {
419  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
420  verbosity_ = pl.get("Verbosity", verbosity_);
421  } else {
422  verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
423  }
424  }
425 
426  // recover from LOBPCGRitzFailure
427  recover_ = pl.get("Recover",recover_);
428 
429  // get (optionally) an initial state
430  if (pl.isParameter("Init")) {
431  state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
432  }
433 }
434 
435 
436 // solve()
437 template<class ScalarType, class MV, class OP>
440 
441  typedef SolverUtils<ScalarType,MV,OP> msutils;
442 
443  const int nev = problem_->getNEV();
444 
445 
446 
448  // Sort manager
449  Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
450 
452  // Output manager
453  Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity_) );
454 
456  // Status tests
457  //
458  // maximum number of iterations: optional test
459  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxtest;
460  if (maxIters_ > 0) {
461  maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
462  }
463  // convergence
464  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
465  if (globalTest_ == Teuchos::null) {
466  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
467  }
468  else {
469  convtest = globalTest_;
470  }
471  Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
472  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
473  // locking
474  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
475  if (useLocking_) {
476  if (lockingTest_ == Teuchos::null) {
477  locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
478  }
479  else {
480  locktest = lockingTest_;
481  }
482  }
483  // for a non-short-circuited OR test, the order doesn't matter
484  Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
485  alltests.push_back(ordertest);
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);
489  Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
491  // printing StatusTest
492  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
493  if ( printer->isVerbosity(Debug) ) {
494  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
495  }
496  else {
497  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
498  }
499 
501  // Orthomanager
502  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
503  if (ortho_=="SVQB") {
504  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
505  } else if (ortho_=="DGKS") {
506  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
507  } else {
508  TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
509  }
510 
512  // Parameter list
513  Teuchos::ParameterList plist;
514  plist.set("Block Size",blockSize_);
515  plist.set("Full Ortho",fullOrtho_);
516 
518  // LOBPCG solver
519  Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
520  = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
521  // set any auxiliary vectors defined in the problem
522  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
523  if (probauxvecs != Teuchos::null) {
524  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
525  }
526 
528  // Storage
529  //
530  // lockvecs will contain eigenvectors that have been determined "locked" by the status test
531  int curNumLocked = 0;
532  Teuchos::RCP<MV> lockvecs;
533  if (useLocking_) {
534  lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
535  }
536  std::vector<MagnitudeType> lockvals;
537  // workMV will be used as work space for LOBPCGRitzFailure recovery and locking
538  // it will be partitioned in these cases as follows:
539  // for LOBPCGRitzFailure recovery:
540  // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
541  // total size: 2*3*blocksize
542  // for locking
543  // workMV = [X P MX MP], with MX,MP needing storage only if hasM==true
544  // total size: 2*blocksize or 4*blocksize
545  Teuchos::RCP<MV> workMV;
546  if (fullOrtho_ == false && recover_ == true) {
547  workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
548  }
549  else if (useLocking_) {
550  if (problem_->getM() != Teuchos::null) {
551  workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
552  }
553  else {
554  workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
555  }
556  }
557 
558  // initialize the solution to nothing in case we throw an exception
560  sol.numVecs = 0;
561  problem_->setSolution(sol);
562 
563  // initialize the solver if the user specified a state
564  if (state_ != Teuchos::null) {
565  lobpcg_solver->initialize(*state_);
566  }
567 
568  // enter solve() iterations
569  {
570 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
571  Teuchos::TimeMonitor slvtimer(*_timerSolve);
572 #endif
573 
574  // tell the lobpcg_solver to iterate
575  while (1) {
576  try {
577  lobpcg_solver->iterate();
578 
580  //
581  // check user-specified debug test; if it passed, return an exception
582  //
584  if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
585  throw AnasaziError("Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
586  }
588  //
589  // check convergence first
590  //
592  else if (ordertest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
593  // we have convergence or not
594  // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
595  // ordertest->howMany() will tell us how many
596  break;
597  }
599  //
600  // check locking if we didn't converge
601  //
603  else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
604 
605 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
606  Teuchos::TimeMonitor lcktimer(*_timerLocking);
607 #endif
608 
609  // remove the locked vectors,values from lobpcg_solver: put them in newvecs, newvals
610  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
611  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
612  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
613  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
614  TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
615  "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
616  // get the indices
617  int numnew = locktest->howMany();
618  std::vector<int> indnew = locktest->whichVecs();
619 
620  // don't lock more than maxLocked_; we didn't allocate enough space.
621  if (curNumLocked + numnew > maxLocked_) {
622  numnew = maxLocked_ - curNumLocked;
623  indnew.resize(numnew);
624  }
625 
626  // the call below to lobpcg_solver->setAuxVecs() will reset the solver to unitialized with hasP() == false
627  // store the hasP() state for use below
628  bool hadP = lobpcg_solver->hasP();
629 
630  {
631  // debug printing
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");
635  }
636  std::vector<MagnitudeType> newvals(numnew);
637  Teuchos::RCP<const MV> newvecs;
638  {
639  // work in a local scope, to hide the variabes needed for extracting this info
640  // get the vectors
641  newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
642  // get the values
643  std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
644  for (int i=0; i<numnew; i++) {
645  newvals[i] = allvals[indnew[i]].realpart;
646  }
647  }
648  // put newvecs into lockvecs
649  {
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;
654  }
655  // put newvals into lockvals
656  lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
657  curNumLocked += numnew;
658  // add locked vecs as aux vecs, along with aux vecs from problem
659  {
660  std::vector<int> indlock(curNumLocked);
661  for (int i=0; i<curNumLocked; i++) indlock[i] = i;
662  Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
663  if (probauxvecs != Teuchos::null) {
664  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs,curlocked) );
665  }
666  else {
667  lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(curlocked) );
668  }
669  }
670  // add locked vals to ordertest
671  ordertest->setAuxVals(lockvals);
672  // fill out the empty state in the solver
673  {
674  LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
675  Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
676  //
677  // workMV will be partitioned as follows: workMV = [X P MX MP],
678  //
679  // make a copy of the current X,MX state
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);
684 
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);
690  }
691  //
692  // get select part, set to random, apply M
693  {
694  Teuchos::RCP<MV> newX = MVT::CloneViewNonConst(*newstateX,indnew);
695  MVT::MvRandom(*newX);
696 
697  if (newstateMX != Teuchos::null) {
698  Teuchos::RCP<MV> newMX = MVT::CloneViewNonConst(*newstateMX,indnew);
699  OPT::Apply(*problem_->getM(),*newX,*newMX);
700  }
701  }
702 
703  Teuchos::Array<Teuchos::RCP<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
704  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
705  // ortho X against the aux vectors
706  ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
707 
708  if (hadP) {
709  //
710  // get P and optionally MP, orthogonalize against X and auxiliary vectors
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);
715 
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);
721  }
722 
723  if (fullOrtho_) {
724  // ortho P against the new aux vectors and new X
725  curauxvecs.push_back(newstateX);
726  ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
727  }
728  else {
729  // ortho P against the new aux vectors
730  ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
731  }
732  }
733  // set the new state
735  newstate.X = newstateX;
736  newstate.MX = newstateMX;
737  newstate.P = newstateP;
738  newstate.MP = newstateMP;
739  lobpcg_solver->initialize(newstate);
740  }
741 
742  if (curNumLocked == maxLocked_) {
743  // disable locking now; remove locking test from combo test
744  combotest->removeTest(locktest);
745  }
746  }
747  else {
748  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
749  }
750  }
752  //
753  // check Ritz Failure
754  //
756  catch (const LOBPCGRitzFailure &re) {
757  if (fullOrtho_==true || recover_==false) {
758  // if we are already using full orthogonalization, there isn't much we can do here.
759  // the most recent information in the status tests is still valid, and can be used to extract/return the
760  // eigenpairs that have converged.
761  printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
762  << "Will not try to recover." << std::endl;
763  break; // while(1)
764  }
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;
767  // get the current "basis" from the solver, orthonormalize it, do a rayleigh-ritz, and restart with the ritz vectors
768  // if there aren't enough, break and quit with what we have
769  //
770  // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
771  LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
772  Teuchos::RCP<MV> restart, Krestart, Mrestart;
773  int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
774  bool hasM = problem_->getM() != Teuchos::null;
775  {
776  std::vector<int> recind(localsize);
777  for (int i=0; i<localsize; i++) recind[i] = i;
778  restart = MVT::CloneViewNonConst(*workMV,recind);
779  }
780  {
781  std::vector<int> recind(localsize);
782  for (int i=0; i<localsize; i++) recind[i] = localsize+i;
783  Krestart = MVT::CloneViewNonConst(*workMV,recind);
784  }
785  if (hasM) {
786  Mrestart = Krestart;
787  }
788  else {
789  Mrestart = restart;
790  }
791  //
792  // set restart = [X H P] and Mrestart = M*[X H P]
793  //
794  // put X into [0 , blockSize)
795  {
796  std::vector<int> blk1(blockSize_);
797  for (int i=0; i < blockSize_; i++) blk1[i] = i;
798  MVT::SetBlock(*curstate.X,blk1,*restart);
799 
800  // put MX into [0 , blockSize)
801  if (hasM) {
802  MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
803  }
804  }
805  //
806  // put H into [blockSize_ , 2*blockSize)
807  {
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);
811 
812  // put MH into [blockSize_ , 2*blockSize)
813  if (hasM) {
814  MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
815  }
816  }
817  // optionally, put P into [2*blockSize,3*blockSize)
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);
822 
823  // put MP into [2*blockSize,3*blockSize)
824  if (hasM) {
825  MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
826  }
827  }
828  // project against auxvecs and locked vecs, and orthonormalize the basis
829  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
830  Teuchos::Array<Teuchos::RCP<const MV> > Q;
831  {
832  if (curNumLocked > 0) {
833  std::vector<int> indlock(curNumLocked);
834  for (int i=0; i<curNumLocked; i++) indlock[i] = i;
835  Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
836  Q.push_back(curlocked);
837  }
838  if (probauxvecs != Teuchos::null) {
839  Q.push_back(probauxvecs);
840  }
841  }
842  int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
843  if (rank < blockSize_) {
844  // quit
845  printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
846  << "Recovery failed." << std::endl;
847  break;
848  }
849  // reduce multivec size if necessary
850  if (rank < localsize) {
851  localsize = rank;
852  std::vector<int> redind(localsize);
853  for (int i=0; i<localsize; i++) redind[i] = i;
854  // grab the first part of restart,Krestart
855  restart = MVT::CloneViewNonConst(*restart,redind);
856  Krestart = MVT::CloneViewNonConst(*Krestart,redind);
857  if (hasM) {
858  Mrestart = Krestart;
859  }
860  else {
861  Mrestart = restart;
862  }
863  }
864  Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
865  std::vector<MagnitudeType> theta(localsize);
866  // project the matrices
867  //
868  // MM = restart^H M restart
869  MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
870  //
871  // compute Krestart = K*restart
872  OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
873  //
874  // KK = restart^H K restart
875  MVT::MvTransMv(1.0,*restart,*Krestart,KK);
876  rank = localsize;
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;
882  break;
883  }
884  theta.resize(rank);
885  //
886  // sort the ritz values using the sort manager
887  {
888  Teuchos::BLAS<int,ScalarType> blas;
889  std::vector<int> order(rank);
890  // sort
891  sorter->sort( theta, Teuchos::rcpFromRef(order),rank ); // don't catch exception
892  // Sort the primitive ritz vectors
893  Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,rank,rank);
894  msutils::permuteVectors(order,curS);
895  }
896  //
897  Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,localsize,blockSize_);
898  //
899  // compute the ritz vectors: store them in Krestart
901  Teuchos::RCP<MV> newX;
902  {
903  std::vector<int> bsind(blockSize_);
904  for (int i=0; i<blockSize_; i++) bsind[i] = i;
905  newX = MVT::CloneViewNonConst(*Krestart,bsind);
906  }
907  MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
908  // send X and theta into the solver
909  newstate.X = newX;
910  theta.resize(blockSize_);
911  newstate.T = Teuchos::rcpFromRef(theta);
912  // initialize
913  lobpcg_solver->initialize(newstate);
914  }
915  catch (const AnasaziError &err) {
916  printer->stream(Errors)
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;
920  return Unconverged;
921  }
922  // don't catch any other exceptions
923  }
924 
925  sol.numVecs = ordertest->howMany();
926  if (sol.numVecs > 0) {
927  sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
928  sol.Espace = sol.Evecs;
929  sol.Evals.resize(sol.numVecs);
930  std::vector<MagnitudeType> vals(sol.numVecs);
931 
932  // copy them into the solution
933  std::vector<int> which = ordertest->whichVecs();
934  // indices between [0,blockSize) refer to vectors/values in the solver
935  // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
936  // everything has already been ordered by the solver; we just have to partition the two references
937  std::vector<int> inlocked(0), insolver(0);
938  for (unsigned int i=0; i<which.size(); i++) {
939  if (which[i] >= 0) {
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]);
942  }
943  else {
944  // sanity check
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);
947  }
948  }
949 
950  TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
951 
952  // set the vecs,vals in the solution
953  if (insolver.size() > 0) {
954  // set vecs
955  int lclnum = insolver.size();
956  std::vector<int> tosol(lclnum);
957  for (int i=0; i<lclnum; i++) tosol[i] = i;
958  Teuchos::RCP<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
959  MVT::SetBlock(*v,tosol,*sol.Evecs);
960  // set vals
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;
964  }
965  }
966 
967  // get the vecs,vals from locked storage
968  if (inlocked.size() > 0) {
969  int solnum = insolver.size();
970  // set vecs
971  int lclnum = inlocked.size();
972  std::vector<int> tosol(lclnum);
973  for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
974  Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
975  MVT::SetBlock(*v,tosol,*sol.Evecs);
976  // set vals
977  for (unsigned int i=0; i<inlocked.size(); i++) {
978  vals[i+solnum] = lockvals[inlocked[i]];
979  }
980  }
981 
982  // sort the eigenvalues and permute the eigenvectors appropriately
983  {
984  std::vector<int> order(sol.numVecs);
985  sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
986  // store the values in the Eigensolution
987  for (int i=0; i<sol.numVecs; i++) {
988  sol.Evals[i].realpart = vals[i];
989  sol.Evals[i].imagpart = MT::zero();
990  }
991  // now permute the eigenvectors according to order
992  msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
993  }
994 
995  // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
996  sol.index.resize(sol.numVecs,0);
997  }
998  }
999 
1000  // print final summary
1001  lobpcg_solver->currentStatus(printer->stream(FinalSummary));
1002 
1003  // print timing information
1004 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1005  if ( printer->isVerbosity( TimingDetails ) ) {
1006  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
1007  }
1008 #endif
1009 
1010  problem_->setSolution(sol);
1011  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1012 
1013  // get the number of iterations performed in this call to solve.
1014  numIters_ = lobpcg_solver->getNumIters();
1015 
1016  if (sol.numVecs < nev) {
1017  return Unconverged; // return from LOBPCGSolMgr::solve()
1018  }
1019  return Converged; // return from LOBPCGSolMgr::solve()
1020 }
1021 
1022 
1023 template <class ScalarType, class MV, class OP>
1024 void
1026  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
1027 {
1028  globalTest_ = global;
1029 }
1030 
1031 template <class ScalarType, class MV, class OP>
1032 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1034 {
1035  return globalTest_;
1036 }
1037 
1038 template <class ScalarType, class MV, class OP>
1039 void
1041  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
1042 {
1043  debugTest_ = debug;
1044 }
1045 
1046 template <class ScalarType, class MV, class OP>
1047 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1049 {
1050  return debugTest_;
1051 }
1052 
1053 template <class ScalarType, class MV, class OP>
1054 void
1056  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
1057 {
1058  lockingTest_ = locking;
1059 }
1060 
1061 template <class ScalarType, class MV, class OP>
1062 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1064 {
1065  return lockingTest_;
1066 }
1067 
1068 } // end Anasazi namespace
1069 
1070 #endif /* ANASAZI_LOBPCG_SOLMGR_HPP */
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 .
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...
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration...
Virtual base class which defines basic traits for the operator type.
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.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
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&#39;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...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
Anasazi&#39;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.
A status test for testing the norm of the eigenvectors residuals.
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
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...
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.
void setLockingStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
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...
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.
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.
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.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Class which provides internal utilities for the Anasazi solvers.