Anasazi  Version of the Day
AnasaziTraceMinBaseSolMgr.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_TraceMinBase_SOLMGR_HPP
30 #define ANASAZI_TraceMinBase_SOLMGR_HPP
31 
38 #include "AnasaziBasicSort.hpp"
39 #include "AnasaziConfigDefs.hpp"
40 #include "AnasaziEigenproblem.hpp"
42 #include "AnasaziSolverManager.hpp"
43 #include "AnasaziSolverUtils.hpp"
47 #include "AnasaziStatusTestSpecTrans.hpp"
50 #include "AnasaziTraceMinBase.hpp"
51 #include "AnasaziTraceMinTypes.hpp"
52 #include "AnasaziTypes.hpp"
53 
54 #include "Teuchos_TimeMonitor.hpp"
55 #ifdef TEUCHOS_DEBUG
56 # include <Teuchos_FancyOStream.hpp>
57 #endif
58 #ifdef HAVE_MPI
59 #include <mpi.h>
60 #endif
61 
62 using Teuchos::RCP;
63 using Teuchos::rcp;
64 
65 namespace Anasazi {
66 namespace Experimental {
67 
68 
101 template<class ScalarType, class MV, class OP>
102 class TraceMinBaseSolMgr : public SolverManager<ScalarType,MV,OP> {
103 
104  private:
108  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
110 
111  public:
112 
114 
115 
160 
162  virtual ~TraceMinBaseSolMgr() {};
164 
166 
167 
170  return *problem_;
171  }
172 
174  int getNumIters() const {
175  return numIters_;
176  }
177 
186  return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
187  }
188 
190 
192 
193 
217  ReturnType solve();
218 
221 
224 
226  void setLockingStatusTest(const RCP< StatusTest<ScalarType,MV,OP> > &locking);
227 
230 
233 
236 
238 
239  protected:
241 
242  int numIters_;
243 
244  // Block variables
245  int blockSize_, numBlocks_, numRestartBlocks_;
246 
247  // Output variables
249 
250  // Convergence variables
251  MagnitudeType convTol_;
252  bool relConvTol_;
253  enum ResType convNorm_;
254 
255  // Locking variables
256  MagnitudeType lockTol_;
257  int maxLocked_, lockQuorum_;
258  bool useLocking_, relLockTol_;
259  enum ResType lockNorm_;
260 
261  // Shifting variables
262  enum WhenToShiftType whenToShift_;
263  MagnitudeType traceThresh_, shiftTol_;
264  enum HowToShiftType howToShift_;
265  bool useMultipleShifts_, relShiftTol_, considerClusters_;
266  std::string shiftNorm_;
267 
268  // Other variables
269  int maxKrylovIter_;
270  std::string ortho_, which_;
271  enum SaddleSolType saddleSolType_;
272  bool projectAllVecs_, projectLockedVecs_, computeAllRes_, useRHSR_, useHarmonic_, noSort_;
273  MagnitudeType alpha_;
274 
275  // Timers
276  RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
277 
278  // Status tests
279  RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
280  RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
282 
283  // TraceMin specific functions
284  void copyPartOfState(const TraceMinBaseState<ScalarType,MV>& oldState, TraceMinBaseState<ScalarType,MV>& newState, const std::vector<int> indToCopy) const;
285 
286  void setParameters(Teuchos::ParameterList &pl) const;
287 
288  void printParameters(std::ostream &os) const;
289 
290  virtual RCP< TraceMinBase<ScalarType,MV,OP> > createSolver(
292  const RCP<StatusTest<ScalarType,MV,OP> > &outputtest,
293  const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
295  ) =0;
296 
297  virtual bool needToRestart(const RCP< TraceMinBase<ScalarType,MV,OP> > solver) =0;
298 
299  virtual bool performRestart(int &numRestarts, RCP< TraceMinBase<ScalarType,MV,OP> > solver) =0;
300 };
301 
302 
305 // Constructor
306 template<class ScalarType, class MV, class OP>
308  const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
309  Teuchos::ParameterList &pl ) :
310  problem_(problem)
311 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
312  , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBaseSolMgr::solve()")),
313  _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBaseSolMgr restarting")),
314  _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBaseSolMgr locking"))
315 #endif
316 {
317  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
318  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
319  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
320  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
321 
322  std::string strtmp;
323 
325  // Block parameters
326 
327  // TODO: the default is different for TraceMin and TraceMin-Davidson
328  // block size: default is nev()
329 // blockSize_ = pl.get("Block Size",problem_->getNEV());
330 // TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
331 // "Anasazi::TraceMinBaseSolMgr: \"Block Size\" must be strictly positive.");
332 
333  // TODO: add Num Blocks as a parameter to both child classes, since they have different default values
334 // numBlocks_ = pl.get("Num Blocks",5);
335 // TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ < 2, std::invalid_argument,
336 // "Anasazi::TraceMinBaseSolMgr: \"Num Blocks\" must be >= 2.");
337 
339  // Output parameters
340 
341  // output stream
342  std::string fntemplate = "";
343  bool allProcs = false;
344  if (pl.isParameter("Output on all processors")) {
345  if (Teuchos::isParameterType<bool>(pl,"Output on all processors")) {
346  allProcs = pl.get("Output on all processors",allProcs);
347  } else {
348  allProcs = ( Teuchos::getParameter<int>(pl,"Output on all processors") != 0 );
349  }
350  }
351  fntemplate = pl.get("Output filename template",fntemplate);
352  int MyPID;
353 # ifdef HAVE_MPI
354  // Initialize MPI
355  int mpiStarted = 0;
356  MPI_Initialized(&mpiStarted);
357  if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
358  else MyPID=0;
359 # else
360  MyPID = 0;
361 # endif
362  if (fntemplate != "") {
363  std::ostringstream MyPIDstr;
364  MyPIDstr << MyPID;
365  // replace %d in fntemplate with MyPID
366  int pos, start=0;
367  while ( (pos = fntemplate.find("%d",start)) != -1 ) {
368  fntemplate.replace(pos,2,MyPIDstr.str());
369  start = pos+2;
370  }
371  }
372  RCP<ostream> osp;
373  if (fntemplate != "") {
374  osp = rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
375  if (!*osp) {
376  osp = Teuchos::rcpFromRef(std::cout);
377  std::cout << "Anasazi::TraceMinBaseSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
378  }
379  }
380  else {
381  osp = Teuchos::rcpFromRef(std::cout);
382  }
383  // Output manager
384  int verbosity = Anasazi::Errors;
385  if (pl.isParameter("Verbosity")) {
386  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
387  verbosity = pl.get("Verbosity", verbosity);
388  } else {
389  verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
390  }
391  }
392  if (allProcs) {
393  // print on all procs
394  printer_ = rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
395  }
396  else {
397  // print only on proc 0
398  printer_ = rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
399  }
400 
401  // TODO: Add restart parameters to TraceMin-Davidson
402 
404  // Convergence parameters
405  convTol_ = pl.get("Convergence Tolerance",MT::prec());
406  TEUCHOS_TEST_FOR_EXCEPTION(convTol_ < 0, std::invalid_argument,
407  "Anasazi::TraceMinBaseSolMgr: \"Convergence Tolerance\" must be nonnegative.");
408 
409  relConvTol_ = pl.get("Relative Convergence Tolerance",true);
410  strtmp = pl.get("Convergence Norm",std::string("2"));
411  if (strtmp == "2") {
412  convNorm_ = RES_2NORM;
413  }
414  else if (strtmp == "M") {
415  convNorm_ = RES_ORTH;
416  }
417  else {
418  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
419  "Anasazi::TraceMinBaseSolMgr: Invalid Convergence Norm.");
420  }
421 
423  // Locking parameters
424  useLocking_ = pl.get("Use Locking",true);
425  relLockTol_ = pl.get("Relative Locking Tolerance",true);
426  lockTol_ = pl.get("Locking Tolerance",convTol_/10);
427 
428  TEUCHOS_TEST_FOR_EXCEPTION(relConvTol_ != relLockTol_, std::invalid_argument,
429  "Anasazi::TraceMinBaseSolMgr: \"Relative Convergence Tolerance\" and \"Relative Locking Tolerance\" have different values. If you set one, you should always set the other.");
430 
431  TEUCHOS_TEST_FOR_EXCEPTION(lockTol_ < 0, std::invalid_argument,
432  "Anasazi::TraceMinBaseSolMgr: \"Locking Tolerance\" must be nonnegative.");
433 
434  strtmp = pl.get("Locking Norm",std::string("2"));
435  if (strtmp == "2") {
436  lockNorm_ = RES_2NORM;
437  }
438  else if (strtmp == "M") {
439  lockNorm_ = RES_ORTH;
440  }
441  else {
442  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
443  "Anasazi::TraceMinBaseSolMgr: Invalid Locking Norm.");
444  }
445 
446  // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
447  if (useLocking_) {
448  maxLocked_ = pl.get("Max Locked",problem_->getNEV());
449  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ <= 0, std::invalid_argument,
450  "Anasazi::TraceMinBaseSolMgr: \"Max Locked\" must be strictly positive.");
451  }
452  else {
453  maxLocked_ = 0;
454  }
455 
456  if (useLocking_) {
457  lockQuorum_ = pl.get("Locking Quorum",1);
458  TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0, std::invalid_argument,
459  "Anasazi::TraceMinBaseSolMgr: \"Locking Quorum\" must be strictly positive.");
460  }
461 
463  // Ritz shift parameters
464 
465  // When to shift - what triggers a shift?
466  strtmp = pl.get("When To Shift", "Always");
467 
468  if(strtmp == "Never")
469  whenToShift_ = NEVER_SHIFT;
470  else if(strtmp == "After Trace Levels")
471  whenToShift_ = SHIFT_WHEN_TRACE_LEVELS;
472  else if(strtmp == "Residual Becomes Small")
473  whenToShift_ = SHIFT_WHEN_RESID_SMALL;
474  else if(strtmp == "Always")
475  whenToShift_ = ALWAYS_SHIFT;
476  else
477  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
478  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"When To Shift\"; valid options are \"Never\", \"After Trace Levels\", \"Residual Becomes Small\", \"Always\".");
479 
480 
481  // How small does the % change in trace have to get before shifting?
482  traceThresh_ = pl.get("Trace Threshold", 0.02);
483 
484  TEUCHOS_TEST_FOR_EXCEPTION(traceThresh_ < 0, std::invalid_argument,
485  "Anasazi::TraceMinBaseSolMgr: \"Trace Threshold\" must be nonnegative.");
486 
487  // Shift threshold - if the residual of an eigenpair is less than this, then shift
488  shiftTol_ = pl.get("Shift Tolerance", 0.1);
489 
490  TEUCHOS_TEST_FOR_EXCEPTION(shiftTol_ < 0, std::invalid_argument,
491  "Anasazi::TraceMinBaseSolMgr: \"Shift Tolerance\" must be nonnegative.");
492 
493  // Use relative convergence tolerance - scale by eigenvalue?
494  relShiftTol_ = pl.get("Relative Shift Tolerance", true);
495 
496  // Which norm to use in determining whether to shift
497  shiftNorm_ = pl.get("Shift Norm", "2");
498 
499  TEUCHOS_TEST_FOR_EXCEPTION(shiftNorm_ != "2" && shiftNorm_ != "M", std::invalid_argument,
500  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Shift Norm\"; valid options are \"2\" and \"M\".");
501 
502  noSort_ = pl.get("No Sorting", false);
503 
504  // How to choose shift
505  strtmp = pl.get("How To Choose Shift", "Adjusted Ritz Values");
506 
507  if(strtmp == "Largest Converged")
508  howToShift_ = LARGEST_CONVERGED_SHIFT;
509  else if(strtmp == "Adjusted Ritz Values")
510  howToShift_ = ADJUSTED_RITZ_SHIFT;
511  else if(strtmp == "Ritz Values")
512  howToShift_ = RITZ_VALUES_SHIFT;
513  else if(strtmp == "Experimental Shift")
514  howToShift_ = EXPERIMENTAL_SHIFT;
515  else
516  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
517  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"How To Choose Shift\"; valid options are \"Largest Converged\", \"Adjusted Ritz Values\", \"Ritz Values\".");
518 
519  // Consider clusters - if all eigenvalues are in one cluster, it's not expecially safe to shift
520  considerClusters_ = pl.get("Consider Clusters", true);
521 
522  // Use multiple shifts
523  useMultipleShifts_ = pl.get("Use Multiple Shifts", true);
524 
526  // Other parameters
527 
528  // which orthogonalization to use
529  ortho_ = pl.get("Orthogonalization", "SVQB");
530  TEUCHOS_TEST_FOR_EXCEPTION(ortho_ != "DGKS" && ortho_ != "SVQB" && ortho_ != "ICGS", std::invalid_argument,
531  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Orthogonalization\"; valid options are \"DGKS\", \"SVQB\", \"ICGS\".");
532 
533  strtmp = pl.get("Saddle Solver Type", "Projected Krylov");
534 
535  if(strtmp == "Projected Krylov")
536  saddleSolType_ = PROJECTED_KRYLOV_SOLVER;
537  else if(strtmp == "Schur Complement")
538  saddleSolType_ = SCHUR_COMPLEMENT_SOLVER;
539  else if(strtmp == "Block Diagonal Preconditioned Minres")
540  saddleSolType_ = BD_PREC_MINRES;
541  else if(strtmp == "HSS Preconditioned Gmres")
542  saddleSolType_ = HSS_PREC_GMRES;
543  else
544  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
545  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Saddle Solver Type\"; valid options are \"Projected Krylov\", \"Schur Complement\", and \"Block Diagonal Preconditioned Minres\".");
546 
547  projectAllVecs_ = pl.get("Project All Vectors", true);
548  projectLockedVecs_ = pl.get("Project Locked Vectors", true);
549  computeAllRes_ = pl.get("Compute All Residuals", true);
550  useRHSR_ = pl.get("Use Residual as RHS", false);
551  alpha_ = pl.get("HSS: alpha", 1.0);
552 
553  TEUCHOS_TEST_FOR_EXCEPTION(projectLockedVecs_ && ! projectAllVecs_, std::invalid_argument,
554  "Anasazi::TraceMinBaseSolMgr: If you want to project out the locked vectors, you should really project out ALL the vectors of X.");
555 
556  // Maximum number of inner iterations
557  maxKrylovIter_ = pl.get("Maximum Krylov Iterations", 200);
558  TEUCHOS_TEST_FOR_EXCEPTION(maxKrylovIter_ < 1, std::invalid_argument,
559  "Anasazi::TraceMinBaseSolMgr: \"Maximum Krylov Iterations\" must be greater than 0.");
560 
561  // Which eigenvalues we want to get
562  which_ = pl.get("Which", "SM");
563  TEUCHOS_TEST_FOR_EXCEPTION(which_ != "SM" && which_ != "LM", std::invalid_argument,
564  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Which\"; valid options are \"SM\" and \"LM\".");
565 
566  // Test whether we are shifting without an operator K
567  // This is a really bad idea
568  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getOperator() == Teuchos::null && whenToShift_ != NEVER_SHIFT, std::invalid_argument,
569  "Anasazi::TraceMinBaseSolMgr: It is an exceptionally bad idea to use Ritz shifts when finding the largest eigenpairs of a standard eigenvalue problem. If we don't use Ritz shifts, it may take extra iterations to converge, but we NEVER have to solve a single linear system. Using Ritz shifts forces us to solve systems of the form (I + sigma A)x=f, and it probably doesn't benefit us enough to outweigh the extra cost. We may add support for this feature in the future, but for now, please set \"When To Shift\" to \"Never\".");
570 
571 #ifdef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
572  // Test whether we are using a projected preconditioner with multiple Ritz shifts
573  // We can't currently do this for reasons that are complicated and are explained in the user manual
574  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER && useMultipleShifts_, std::invalid_argument,
575  "Anasazi::TraceMinBaseSolMgr: When you use the projected Krylov solver with preconditioning, the preconditioner must be projected as well. In theory, if the preconditioner is SPD, the projected preconditioner will be SPSD, but in practice, it can have small negative eigenvalues, presumably due to machine arithmetic. This means we can't use TraceMin's built-in MINRES, and we are forced to use Belos for now. When you use multiple Ritz shifts, you are essentially using a different operator to solve each linear system. Belos can't handle this right now, but we're working on a solution. For now, please set \"Use Multiple Shifts\" to false.");
576 #else
577  // Test whether we are using a projected preconditioner without Belos.
578  // P Prec P should be positive definite if Prec is positive-definite,
579  // but it tends not to be in practice, presumably due to machine arithmetic
580  // As a result, we have to use pseudo-block gmres for now.
581  // Make sure it's available.
582  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER, std::invalid_argument,
583  "Anasazi::TraceMinBaseSolMgr: When you use the projected Krylov solver with preconditioning, the preconditioner must be projected as well. In theory, if the preconditioner is SPD, the projected preconditioner will be SPSD, but in practice, it can have small negative eigenvalues, presumably due to machine arithmetic. This means we can't use TraceMin's built-in MINRES, and we are forced to use Belos for now. You didn't install Belos. You have three options to correct this problem:\n1. Reinstall Trilinos with Belos enabled\n2. Don't use a preconditioner\n3. Choose a different method for solving the saddle-point problem (Recommended)");
584 
585 
586 #endif
587 
588 
589 }
590 
591 
594 // solve()
595 template<class ScalarType, class MV, class OP>
596 ReturnType
598 {
599  typedef SolverUtils<ScalarType,MV,OP> msutils;
600 
601  const int nev = problem_->getNEV();
602 
603 #ifdef TEUCHOS_DEBUG
605  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
606  out->setShowAllFrontMatter(false).setShowProcRank(true);
607  *out << "Entering Anasazi::TraceMinBaseSolMgr::solve()\n";
608 #endif
609 
611  // Sort manager
613 
615  // Handle the spectral transformation if necessary
616  // TODO: Make sure we undo this before returning...
617  if(which_ == "LM")
618  {
619  RCP<const OP> swapHelper = problem_->getOperator();
620  problem_->setOperator(problem_->getM());
621  problem_->setM(swapHelper);
622  problem_->setProblem();
623  }
624 
626  // Status tests
627  //
628  // convergence
630  if (globalTest_ == Teuchos::null) {
631  if(which_ == "SM")
632  convtest = rcp( new StatusTestResNorm<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_) );
633  else
634  convtest = rcp( new StatusTestSpecTrans<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_,true,problem_->getOperator()) );
635  }
636  else {
637  convtest = globalTest_;
638  }
640  = rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
641  // locking
643  if (useLocking_) {
644  if (lockingTest_ == Teuchos::null) {
645  if(which_ == "SM")
646  locktest = rcp( new StatusTestResNorm<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_) );
647  else
648  locktest = rcp( new StatusTestSpecTrans<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_,true,problem_->getOperator()) );
649  }
650  else {
651  locktest = lockingTest_;
652  }
653  }
654  // for a non-short-circuited OR test, the order doesn't matter
656  alltests.push_back(ordertest);
657  if (locktest != Teuchos::null) alltests.push_back(locktest);
658  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
659 
662  // printing StatusTest
664  if ( printer_->isVerbosity(Debug) ) {
665  outputtest = rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
666  }
667  else {
668  outputtest = rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
669  }
670 
672  // Orthomanager
674  if (ortho_=="SVQB") {
675  ortho = rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
676  } else if (ortho_=="DGKS") {
677  ortho = rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
678  } else if (ortho_=="ICGS") {
679  ortho = rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
680  } else {
681  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): Invalid orthogonalization type.");
682  }
683 
685  // Parameter list
687  setParameters(plist);
688 
690  // TraceMinBase solver
692  = createSolver(sorter,outputtest,ortho,plist);
693  // set any auxiliary vectors defined in the problem
694  RCP< const MV > probauxvecs = problem_->getAuxVecs();
695  if (probauxvecs != Teuchos::null) {
696  tm_solver->setAuxVecs( Teuchos::tuple< RCP<const MV> >(probauxvecs) );
697  }
698 
700  // Storage
701  //
702  // lockvecs will contain eigenvectors that have been determined "locked" by the status test
703  int curNumLocked = 0;
704  RCP<MV> lockvecs;
705  // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking.
706  // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
707  // we will produce numnew random vectors, which will go into the space with the new basis.
708  // we will also need numnew storage for the image of these random vectors under A and M;
709  // columns [curlocked+1,curlocked+numnew] will be used for this storage
710  if (maxLocked_ > 0) {
711  lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
712  }
713  std::vector<MagnitudeType> lockvals;
714  //
715  // Restarting occurs under two scenarios: when the basis is full and after locking.
716  //
717  // For the former, a new basis of size blockSize*numRestartBlocks is generated using the current basis
718  // and the most significant primitive Ritz vectors (projected eigenvectors).
719  // [S,L] = eig(KK)
720  // S = [Sr St] // some for "r"estarting, some are "t"runcated
721  // newV = V*Sr
722  // KK_new = newV'*K*newV = Sr'*V'*K*V*Sr = Sr'*KK*Sr
723  // Therefore, the only multivector operation needed is for the generation of newV.
724  //
725  // * If the multiplication is explicit, it requires a workspace of blockSize*numRestartBlocks vectors.
726  // This space must be specifically allocated for that task, as we don't have any space of that size.
727  // It (workMV) will be allocated at the beginning of solve()
728  // * Optionally, the multiplication can be performed implicitly, via a Householder QR factorization of
729  // Sr. This can be done in situ, using the basis multivector contained in the solver. This requires
730  // that we cast away the const on the multivector returned from getState(). Workspace for this approach
731  // is a single vector. the solver's internal storage must be preserved (X,MX,KX,R), requiring us to
732  // allocate this vector.
733  //
734  // For the latter (restarting after locking), the new basis is the same size as existing basis. If numnew
735  // vectors are locked, they are deflated from the current basis and replaced with randomly generated
736  // vectors.
737  // [S,L] = eig(KK)
738  // S = [Sl Su] // partitioned: "l"ocked and "u"nlocked
739  // newL = V*Sl = X(locked)
740  // defV = V*Su
741  // augV = rand(numnew) // orthogonal to oldL,newL,defV,auxvecs
742  // newV = [defV augV]
743  // Kknew = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
744  // [augV'*K*defV augV'*K*augV]
745  // locked = [oldL newL]
746  // Clearly, this operation is more complicated than the previous.
747  // Here is a list of the significant computations that need to be performed:
748  // - newL will be put into space in lockvecs, but will be copied from getState().X at the end
749  // - defV,augV will be stored in workspace the size of the current basis.
750  // - M*augV and K*augV are needed; they will be stored in lockvecs. As a result, newL will
751  // not be put into lockvecs until the end.
752  //
753  // Therefore, we must allocate workMV when ((maxRestarts_ > 0) || (useLocking_ == true)) && inSituRestart == false
754  // It will be allocated to size (numBlocks-1)*blockSize
755  //
756 
757  // some consts and utils
758  const ScalarType ONE = SCT::one();
759  const ScalarType ZERO = SCT::zero();
760 
761  // go ahead and initialize the solution to nothing in case we throw an exception
763  sol.numVecs = 0;
764  problem_->setSolution(sol);
765 
766  int numRestarts = 0;
767 
768  // enter solve() iterations
769  {
770 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
771  Teuchos::TimeMonitor slvtimer(*_timerSolve);
772 #endif
773 
774  // tell tm_solver to iterate
775  while (1) {
776  try {
777  tm_solver->iterate();
778 
780  //
781  //
783  if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
784  throw AnasaziError("Anasazi::TraceMinBaseSolMgr::solve(): User-specified debug status test returned Passed.");
785  }
787  //
788  // check convergence next
789  //
791  else if (ordertest->getStatus() == Passed ) {
792  // we have convergence
793  // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
794  // ordertest->howMany() will tell us how many
795  break;
796  }
798  //
799  // check locking if we didn't converge or restart
800  //
802  else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
803 
804 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
805  Teuchos::TimeMonitor lcktimer(*_timerLocking);
806 #endif
807 
808  //
809  // get current state
810  TraceMinBaseState<ScalarType,MV> state = tm_solver->getState();
811  const int curdim = state.curDim;
812 
813  //
814  // get number,indices of vectors to be locked
815  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
816  "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() non-positive.");
817  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
818  "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
819  // we should have room to lock vectors; otherwise, locking should have been deactivated
820  TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
821  "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: locking not deactivated.");
822  //
823  // don't lock more than maxLocked_; we didn't allocate enough space.
824  std::vector<int> tmp_vector_int;
825  if (curNumLocked + locktest->howMany() > maxLocked_) {
826  // just use the first of them
827  for(int i=0; i<maxLocked_-curNumLocked; i++)
828  tmp_vector_int.push_back(locktest->whichVecs()[i]);
829 // tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
830  }
831  else {
832  tmp_vector_int = locktest->whichVecs();
833  }
834 
835  const std::vector<int> lockind(tmp_vector_int);
836  const int numNewLocked = lockind.size();
837  //
838  // generate indices of vectors left unlocked
839  // curind = [0,...,curdim-1] = UNION( lockind, unlockind )
840  const int numUnlocked = curdim-numNewLocked;
841  tmp_vector_int.resize(curdim);
842  for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
843  const std::vector<int> curind(tmp_vector_int); // curind = [0 ... curdim-1]
844  tmp_vector_int.resize(numUnlocked);
845  std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
846  const std::vector<int> unlockind(tmp_vector_int); // unlockind = [0 ... curdim-1] - lockind
847  tmp_vector_int.clear();
848 
849  //
850  // debug printing
851  if (printer_->isVerbosity(Debug)) {
852  printer_->print(Debug,"Locking vectors: ");
853  for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
854  printer_->print(Debug,"\n");
855  printer_->print(Debug,"Not locking vectors: ");
856  for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
857  printer_->print(Debug,"\n");
858  }
859 
860  // Copy eigenvalues we want to lock into lockvals
861  std::vector<Value<ScalarType> > allvals = tm_solver->getRitzValues();
862  for(unsigned int i=0; i<allvals.size(); i++)
863  printer_->stream(Debug) << "Ritz value[" << i << "] = " << allvals[i].realpart << std::endl;
864  for (int i=0; i<numNewLocked; i++) {
865  lockvals.push_back(allvals[lockind[i]].realpart);
866  }
867 
868  // Copy vectors we want to lock into lockvecs
869  RCP<const MV> newLocked = MVT::CloneView(*tm_solver->getRitzVectors(),lockind);
870  std::vector<int> indlock(numNewLocked);
871  for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
872  if(useHarmonic_)
873  {
874  RCP<MV> tempMV = MVT::CloneCopy(*newLocked);
875  ortho->normalizeMat(*tempMV);
876  MVT::SetBlock(*tempMV,indlock,*lockvecs);
877  }
878  else
879  {
880  MVT::SetBlock(*newLocked,indlock,*lockvecs);
881  }
882 
883  // Tell the StatusTestWithOrdering that things have been locked
884  // This is VERY important
885  // If this set of lines is removed, the code does not terminate correctly
886  if(noSort_)
887  {
888  for(unsigned int aliciaInd=0; aliciaInd<lockvals.size(); aliciaInd++)
889  {
890  lockvals[aliciaInd] = 0.0;
891  }
892  }
893  ordertest->setAuxVals(lockvals);
894 
895  // Set the auxiliary vectors so that we remain orthogonal to the ones we locked
896  // Remember to include any aux vecs provided by the user
897  curNumLocked += numNewLocked;
898 
899  if(ordertest->getStatus() == Passed) break;
900 
901  std::vector<int> curlockind(curNumLocked);
902  for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
903  RCP<const MV> curlocked = MVT::CloneView(*lockvecs,curlockind);
904 
906  if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
907  aux.push_back(curlocked);
908  tm_solver->setAuxVecs(aux);
909 
910  // Disable locking if we have locked the maximum number of things
911  printer_->stream(Debug) << "curNumLocked: " << curNumLocked << std::endl;
912  printer_->stream(Debug) << "maxLocked: " << maxLocked_ << std::endl;
913  if (curNumLocked == maxLocked_) {
914  // disabled locking now
915  combotest->removeTest(locktest);
916  locktest = Teuchos::null;
917  printer_->stream(Debug) << "Removed locking test\n";
918  }
919 
920  int newdim = numRestartBlocks_*blockSize_;
922  if(newdim <= numUnlocked)
923  {
924  if(useHarmonic_)
925  {
926  std::vector<int> desiredSubscripts(newdim);
927  for(int i=0; i<newdim; i++)
928  {
929  desiredSubscripts[i] = unlockind[i];
930  printer_->stream(Debug) << "H desiredSubscripts[" << i << "] = " << desiredSubscripts[i] << std::endl;
931  }
932  newstate.V = MVT::CloneView(*tm_solver->getRitzVectors(),desiredSubscripts);
933  newstate.curDim = newdim;
934  }
935  else
936  {
937  std::vector<int> desiredSubscripts(newdim);
938  for(int i=0; i<newdim; i++)
939  {
940  desiredSubscripts[i] = unlockind[i];
941  printer_->stream(Debug) << "desiredSubscripts[" << i << "] = " << desiredSubscripts[i] << std::endl;
942  }
943 
944  copyPartOfState(state, newstate, desiredSubscripts);
945  }
946  }
947  else
948  {
949  // TODO: Come back to this and make it more efficient
950 
951  // Replace the converged eigenvectors with random ones
952  int nrandom = newdim-numUnlocked;
953 
954  RCP<const MV> helperMV;
955  RCP<MV> totalV = MVT::Clone(*tm_solver->getRitzVectors(),newdim);
956 
957  // Holds old things that we're keeping
958  tmp_vector_int.resize(numUnlocked);
959  for(int i=0; i<numUnlocked; i++) tmp_vector_int[i] = i;
960  RCP<MV> oldV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
961 
962  // Copy over the old things
963  if(useHarmonic_)
964  helperMV = MVT::CloneView(*tm_solver->getRitzVectors(),unlockind);
965  else
966  helperMV = MVT::CloneView(*state.V,unlockind);
967  MVT::Assign(*helperMV,*oldV);
968 
969  // Holds random vectors we're generating
970  tmp_vector_int.resize(nrandom);
971  for(int i=0; i<nrandom; i++) tmp_vector_int[i] = i+numUnlocked;
972  RCP<MV> newV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
973 
974  // Create random things
975  MVT::MvRandom(*newV);
976 
977  newstate.V = totalV;
978  newstate.curDim = newdim;
979 
980  // Copy Ritz shifts
981  RCP< std::vector<ScalarType> > helperRS = rcp( new std::vector<ScalarType>(blockSize_) );
982  for(unsigned int i=0; i<(unsigned int)blockSize_; i++)
983  {
984  if(i < unlockind.size() && unlockind[i] < blockSize_)
985  (*helperRS)[i] = (*state.ritzShifts)[unlockind[i]];
986  else
987  (*helperRS)[i] = ZERO;
988  }
989  newstate.ritzShifts = helperRS;
990  }
991 
992  // Determine the largest safe shift
993  newstate.largestSafeShift = std::abs(lockvals[0]);
994  for(size_t i=0; i<lockvals.size(); i++)
995  newstate.largestSafeShift = std::max(std::abs(lockvals[i]), newstate.largestSafeShift);
996 
997  // Prepare new state, removing converged vectors
998  // TODO: Init will perform some unnecessary calculations; look into it
999  // TODO: The residual norms should be part of the state
1000  newstate.NEV = state.NEV - numNewLocked;
1001  tm_solver->initialize(newstate);
1002  } // end of locking
1004  //
1005  // check for restarting before locking: if we need to lock, it will happen after the restart
1006  //
1008  else if ( needToRestart(tm_solver) ) {
1009  // if performRestart returns false, we exceeded the maximum number of restarts
1010  if(performRestart(numRestarts, tm_solver) == false)
1011  break;
1012  } // end of restarting
1014  //
1015  // we returned from iterate(), but none of our status tests Passed.
1016  // something is wrong, and it is probably our fault.
1017  //
1019  else {
1020  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): Invalid return from tm_solver::iterate().");
1021  }
1022  }
1023  catch (const AnasaziError &err) {
1024  printer_->stream(Errors)
1025  << "Anasazi::TraceMinBaseSolMgr::solve() caught unexpected exception from Anasazi::TraceMinBase::iterate() at iteration " << tm_solver->getNumIters() << std::endl
1026  << err.what() << std::endl
1027  << "Anasazi::TraceMinBaseSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1028  return Unconverged;
1029  }
1030  }
1031 
1032  sol.numVecs = ordertest->howMany();
1033  if (sol.numVecs > 0) {
1034  sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
1035  sol.Espace = sol.Evecs;
1036  sol.Evals.resize(sol.numVecs);
1037  std::vector<MagnitudeType> vals(sol.numVecs);
1038 
1039  // copy them into the solution
1040  std::vector<int> which = ordertest->whichVecs();
1041  // indices between [0,blockSize) refer to vectors/values in the solver
1042  // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
1043  // everything has already been ordered by the solver; we just have to partition the two references
1044  std::vector<int> inlocked(0), insolver(0);
1045  for (unsigned int i=0; i<which.size(); i++) {
1046  if (which[i] >= 0) {
1047  TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): positive indexing mistake from ordertest.");
1048  insolver.push_back(which[i]);
1049  }
1050  else {
1051  // sanity check
1052  TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): negative indexing mistake from ordertest.");
1053  inlocked.push_back(which[i] + curNumLocked);
1054  }
1055  }
1056 
1057  TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): indexing mistake.");
1058 
1059  // set the vecs,vals in the solution
1060  if (insolver.size() > 0) {
1061  // set vecs
1062  int lclnum = insolver.size();
1063  std::vector<int> tosol(lclnum);
1064  for (int i=0; i<lclnum; i++) tosol[i] = i;
1065  RCP<const MV> v = MVT::CloneView(*tm_solver->getRitzVectors(),insolver);
1066  MVT::SetBlock(*v,tosol,*sol.Evecs);
1067  // set vals
1068  std::vector<Value<ScalarType> > fromsolver = tm_solver->getRitzValues();
1069  for (unsigned int i=0; i<insolver.size(); i++) {
1070  vals[i] = fromsolver[insolver[i]].realpart;
1071  }
1072  }
1073 
1074  // get the vecs,vals from locked storage
1075  if (inlocked.size() > 0) {
1076  int solnum = insolver.size();
1077  // set vecs
1078  int lclnum = inlocked.size();
1079  std::vector<int> tosol(lclnum);
1080  for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1081  RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1082  MVT::SetBlock(*v,tosol,*sol.Evecs);
1083  // set vals
1084  for (unsigned int i=0; i<inlocked.size(); i++) {
1085  vals[i+solnum] = lockvals[inlocked[i]];
1086  }
1087  }
1088 
1089  // undo the spectral transformation if necessary
1090  // if we really passed the solver Bx = \lambda A x, invert the eigenvalues
1091  if(which_ == "LM")
1092  {
1093  for(size_t i=0; i<vals.size(); i++)
1094  vals[i] = ONE/vals[i];
1095  }
1096 
1097  // sort the eigenvalues and permute the eigenvectors appropriately
1098  {
1099  std::vector<int> order(sol.numVecs);
1100  sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1101  // store the values in the Eigensolution
1102  for (int i=0; i<sol.numVecs; i++) {
1103  sol.Evals[i].realpart = vals[i];
1104  sol.Evals[i].imagpart = MT::zero();
1105  }
1106  // now permute the eigenvectors according to order
1107  msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1108  }
1109 
1110  // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1111  sol.index.resize(sol.numVecs,0);
1112  }
1113  }
1114 
1115  // print final summary
1116  tm_solver->currentStatus(printer_->stream(FinalSummary));
1117 
1118  printParameters(printer_->stream(FinalSummary));
1119 
1120  // print timing information
1121 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1122  if ( printer_->isVerbosity( TimingDetails ) ) {
1123  Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
1124  }
1125 #endif
1126 
1127  problem_->setSolution(sol);
1128  printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1129 
1130  // get the number of iterations taken for this call to solve().
1131  numIters_ = tm_solver->getNumIters();
1132 
1133  if (sol.numVecs < nev) {
1134  return Unconverged; // return from TraceMinBaseSolMgr::solve()
1135  }
1136  return Converged; // return from TraceMinBaseSolMgr::solve()
1137 }
1138 
1139 
1140 template <class ScalarType, class MV, class OP>
1141 void
1143  const RCP< StatusTest<ScalarType,MV,OP> > &global)
1144 {
1145  globalTest_ = global;
1146 }
1147 
1148 template <class ScalarType, class MV, class OP>
1151 {
1152  return globalTest_;
1153 }
1154 
1155 template <class ScalarType, class MV, class OP>
1156 void
1158  const RCP< StatusTest<ScalarType,MV,OP> > &debug)
1159 {
1160  debugTest_ = debug;
1161 }
1162 
1163 template <class ScalarType, class MV, class OP>
1166 {
1167  return debugTest_;
1168 }
1169 
1170 template <class ScalarType, class MV, class OP>
1171 void
1173  const RCP< StatusTest<ScalarType,MV,OP> > &locking)
1174 {
1175  lockingTest_ = locking;
1176 }
1177 
1178 template <class ScalarType, class MV, class OP>
1181 {
1182  return lockingTest_;
1183 }
1184 
1185 template <class ScalarType, class MV, class OP>
1186 void TraceMinBaseSolMgr<ScalarType,MV,OP>::copyPartOfState(const TraceMinBaseState<ScalarType,MV>& oldState, TraceMinBaseState<ScalarType,MV>& newState, const std::vector<int> indToCopy) const
1187 {
1188  const ScalarType ONE = Teuchos::ScalarTraits<MagnitudeType>::one();
1189  const ScalarType ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
1190 
1191  newState.curDim = indToCopy.size();
1192  std::vector<int> fullIndices(oldState.curDim);
1193  for(int i=0; i<oldState.curDim; i++) fullIndices[i] = i;
1194 
1195  // Initialize with X.
1196  // Note that we didn't compute enough vectors of X, but we can very easily using the Ritz vectors.
1197  // That's why they're part of the state.
1198  // Note that there will ALWAYS be enough vectors
1199 
1200  // Helpful vectors for computing views and whatnot
1201  std::vector<int> oldIndices;
1202  std::vector<int> newIndices;
1203  for(int i=0; i<newState.curDim; i++)
1204  {
1205  if(indToCopy[i] < blockSize_)
1206  oldIndices.push_back(indToCopy[i]);
1207  else
1208  newIndices.push_back(indToCopy[i]);
1209  }
1210 
1211  int olddim = oldIndices.size();
1212  int newdim = newIndices.size();
1213 
1214  // If there are no new vectors being copied
1215  if(computeAllRes_)
1216  {
1217  newState.V = MVT::CloneView(*oldState.X, indToCopy);
1218  newState.R = MVT::CloneView(*oldState.R, indToCopy);
1219  newState.X = newState.V;
1220 
1221  if(problem_->getOperator() != Teuchos::null)
1222  {
1223  newState.KV = MVT::CloneView(*oldState.KX, indToCopy);
1224  newState.KX = newState.KV;
1225  }
1226  else
1227  {
1228  newState.KV = Teuchos::null;
1229  newState.KX = Teuchos::null;
1230  }
1231 
1232  if(problem_->getM() != Teuchos::null)
1233  {
1234  newState.MopV = MVT::CloneView(*oldState.MX, indToCopy);
1235  newState.MX = newState.MopV;
1236  }
1237  else
1238  {
1239  newState.MopV = Teuchos::null;
1240  newState.MX = Teuchos::null;
1241  }
1242  }
1243  else if(newdim == 0)
1244  {
1245  std::vector<int> blockind(blockSize_);
1246  for(int i=0; i<blockSize_; i++)
1247  blockind[i] = i;
1248 
1249  // Initialize with X
1250  newState.V = MVT::CloneView(*oldState.X, blockind);
1251  newState.KV = MVT::CloneView(*oldState.KX, blockind);
1252  newState.R = MVT::CloneView(*oldState.R, blockind);
1253  newState.X = MVT::CloneView(*newState.V, blockind);
1254  newState.KX = MVT::CloneView(*newState.KV, blockind);
1255 
1256  if(problem_->getM() != Teuchos::null)
1257  {
1258  newState.MopV = MVT::CloneView(*oldState.MX, blockind);
1259  newState.MX = MVT::CloneView(*newState.MopV, blockind);
1260  }
1261  else
1262  {
1263  newState.MopV = Teuchos::null;
1264  newState.MX = Teuchos::null;
1265  }
1266  }
1267  else
1268  {
1269  // More helpful vectors
1270  std::vector<int> oldPart(olddim);
1271  for(int i=0; i<olddim; i++) oldPart[i] = i;
1272  std::vector<int> newPart(newdim);
1273  for(int i=0; i<newdim; i++) newPart[i] = olddim+i;
1274 
1275  // Helpful multivectors for views and whatnot
1276  RCP<MV> helper = MVT::Clone(*oldState.V,newState.curDim);
1277  RCP<MV> oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1278  RCP<MV> newHelper = MVT::CloneViewNonConst(*helper,newPart);
1279  RCP<const MV> viewHelper;
1280 
1281  // Get the parts of the Ritz vectors we are interested in.
1282  Teuchos::SerialDenseMatrix<int,ScalarType> newRV(oldState.curDim,newdim);
1283  for(int r=0; r<oldState.curDim; r++)
1284  {
1285  for(int c=0; c<newdim; c++)
1286  newRV(r,c) = (*oldState.RV)(r,newIndices[c]);
1287  }
1288 
1289  // We're going to compute X as V*RitzVecs
1290  viewHelper = MVT::CloneView(*oldState.V,fullIndices);
1291  MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1292  viewHelper = MVT::CloneView(*oldState.X,oldIndices);
1293  MVT::Assign(*viewHelper,*oldHelper);
1294  newState.V = MVT::CloneCopy(*helper);
1295 
1296  // Also compute KX as KV*RitzVecs
1297  viewHelper = MVT::CloneView(*oldState.KV,fullIndices);
1298  MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1299  viewHelper = MVT::CloneView(*oldState.KX,oldIndices);
1300  MVT::Assign(*viewHelper,*oldHelper);
1301  newState.KV = MVT::CloneCopy(*helper);
1302 
1303  // Do the same with MX if necessary
1304  if(problem_->getM() != Teuchos::null)
1305  {
1306  viewHelper = MVT::CloneView(*oldState.MopV,fullIndices);
1307  MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1308  viewHelper = MVT::CloneView(*oldState.MX,oldIndices);
1309  MVT::Assign(*viewHelper,*oldHelper);
1310  newState.MopV = MVT::CloneCopy(*helper);
1311  }
1312  else
1313  newState.MopV = newState.V;
1314 
1315  // Get X, MX, KX
1316  std::vector<int> blockVec(blockSize_);
1317  for(int i=0; i<blockSize_; i++) blockVec[i] = i;
1318  newState.X = MVT::CloneView(*newState.V,blockVec);
1319  newState.KX = MVT::CloneView(*newState.KV,blockVec);
1320  newState.MX = MVT::CloneView(*newState.MopV,blockVec);
1321 
1322  // Update the residuals
1323  if(blockSize_-oldIndices.size() > 0)
1324  {
1325  // There are vectors we have not computed the residual for yet
1326  newPart.resize(blockSize_-oldIndices.size());
1327  helper = MVT::Clone(*oldState.V,blockSize_);
1328  oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1329  newHelper = MVT::CloneViewNonConst(*helper,newPart);
1330 
1331  RCP<MV> scaledMV = MVT::CloneCopy(*newState.MX,newPart);
1332  RCP<const MV> localKV = MVT::CloneView(*newState.KX,newPart);
1333  std::vector<ScalarType> scalarVec(blockSize_-oldIndices.size());
1334  for(unsigned int i=0; i<(unsigned int)blockSize_-oldIndices.size(); i++) scalarVec[i] = (*oldState.T)[newPart[i]];
1335  MVT::MvScale(*scaledMV,scalarVec);
1336 
1337  helper = MVT::Clone(*oldState.V,blockSize_);
1338  oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1339  newHelper = MVT::CloneViewNonConst(*helper,newPart);
1340  MVT::MvAddMv(ONE,*localKV,-ONE,*scaledMV,*newHelper);
1341  viewHelper = MVT::CloneView(*oldState.R,oldIndices);
1342  MVT::Assign(*viewHelper,*oldHelper);
1343  newState.R = MVT::CloneCopy(*helper);
1344  }
1345  else
1346  newState.R = oldState.R;
1347  }
1348 
1349  // Since we are setting V:=X, V is orthonormal
1350  newState.isOrtho = true;
1351 
1352  // Get the first eigenvalues
1353  RCP< std::vector<ScalarType> > helperT = rcp( new std::vector<ScalarType>(newState.curDim) );
1354  for(int i=0; i<newState.curDim; i++) (*helperT)[i] = (*oldState.T)[indToCopy[i]];
1355  newState.T = helperT;
1356 
1357  // X'KX is diag(T)
1359  for(int i=0; i<newState.curDim; i++)
1360  (*newKK)(i,i) = (*newState.T)[i];
1361  newState.KK = newKK;
1362 
1363  // The associated Ritz vectors are I
1365  for(int i=0; i<newState.curDim; i++)
1366  (*newRV)(i,i) = ONE;
1367  newState.RV = newRV;
1368 
1369  // Get the Ritz shifts
1370  RCP< std::vector<ScalarType> > helperRS = rcp( new std::vector<ScalarType>(blockSize_) );
1371  for(int i=0; i<blockSize_; i++)
1372  {
1373  if(indToCopy[i] < blockSize_)
1374  (*helperRS)[i] = (*oldState.ritzShifts)[indToCopy[i]];
1375  else
1376  (*helperRS)[i] = ZERO;
1377  }
1378  newState.ritzShifts = helperRS;
1379 }
1380 
1381 
1382 template <class ScalarType, class MV, class OP>
1383 void TraceMinBaseSolMgr<ScalarType,MV,OP>::setParameters(Teuchos::ParameterList &pl) const
1384 {
1385  pl.set("Block Size", blockSize_);
1386  pl.set("Num Blocks", numBlocks_);
1387  pl.set("Num Restart Blocks", numRestartBlocks_);
1388  pl.set("When To Shift", whenToShift_);
1389  pl.set("Trace Threshold", traceThresh_);
1390  pl.set("Shift Tolerance", shiftTol_);
1391  pl.set("Relative Shift Tolerance", relShiftTol_);
1392  pl.set("Shift Norm", shiftNorm_);
1393  pl.set("How To Choose Shift", howToShift_);
1394  pl.set("Consider Clusters", considerClusters_);
1395  pl.set("Use Multiple Shifts", useMultipleShifts_);
1396  pl.set("Saddle Solver Type", saddleSolType_);
1397  pl.set("Project All Vectors", projectAllVecs_);
1398  pl.set("Project Locked Vectors", projectLockedVecs_);
1399  pl.set("Compute All Residuals", computeAllRes_);
1400  pl.set("Use Residual as RHS", useRHSR_);
1401  pl.set("Use Harmonic Ritz Values", useHarmonic_);
1402  pl.set("Maximum Krylov Iterations", maxKrylovIter_);
1403  pl.set("HSS: alpha", alpha_);
1404 }
1405 
1406 
1407 template <class ScalarType, class MV, class OP>
1408 void TraceMinBaseSolMgr<ScalarType,MV,OP>::printParameters(std::ostream &os) const
1409 {
1410  os << "\n\n\n";
1411  os << "========================================\n";
1412  os << "========= TraceMin parameters ==========\n";
1413  os << "========================================\n";
1414  os << "=========== Block parameters ===========\n";
1415  os << "Block Size: " << blockSize_ << std::endl;
1416  os << "Num Blocks: " << numBlocks_ << std::endl;
1417  os << "Num Restart Blocks: " << numRestartBlocks_ << std::endl;
1418  os << "======== Convergence parameters ========\n";
1419  os << "Convergence Tolerance: " << convTol_ << std::endl;
1420  os << "Relative Convergence Tolerance: " << relConvTol_ << std::endl;
1421  os << "========== Locking parameters ==========\n";
1422  os << "Use Locking: " << useLocking_ << std::endl;
1423  os << "Locking Tolerance: " << lockTol_ << std::endl;
1424  os << "Relative Locking Tolerance: " << relLockTol_ << std::endl;
1425  os << "Max Locked: " << maxLocked_ << std::endl;
1426  os << "Locking Quorum: " << lockQuorum_ << std::endl;
1427  os << "========== Shifting parameters =========\n";
1428  os << "When To Shift: ";
1429  if(whenToShift_ == NEVER_SHIFT) os << "Never\n";
1430  else if(whenToShift_ == SHIFT_WHEN_TRACE_LEVELS) os << "After Trace Levels\n";
1431  else if(whenToShift_ == SHIFT_WHEN_RESID_SMALL) os << "Residual Becomes Small\n";
1432  else if(whenToShift_ == ALWAYS_SHIFT) os << "Always\n";
1433  os << "Consider Clusters: " << considerClusters_ << std::endl;
1434  os << "Trace Threshohld: " << traceThresh_ << std::endl;
1435  os << "Shift Tolerance: " << shiftTol_ << std::endl;
1436  os << "Relative Shift Tolerance: " << relShiftTol_ << std::endl;
1437  os << "How To Choose Shift: ";
1438  if(howToShift_ == LARGEST_CONVERGED_SHIFT) os << "Largest Converged\n";
1439  else if(howToShift_ == ADJUSTED_RITZ_SHIFT) os << "Adjusted Ritz Values\n";
1440  else if(howToShift_ == RITZ_VALUES_SHIFT) os << "Ritz Values\n";
1441  os << "Use Multiple Shifts: " << useMultipleShifts_ << std::endl;
1442  os << "=========== Other parameters ===========\n";
1443  os << "Orthogonalization: " << ortho_ << std::endl;
1444  os << "Saddle Solver Type: ";
1445  if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER) os << "Projected Krylov\n";
1446  else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER) os << "Schur Complement\n";
1447  os << "Project All Vectors: " << projectAllVecs_ << std::endl;
1448  os << "Project Locked Vectors: " << projectLockedVecs_ << std::endl;
1449  os << "Compute All Residuals: " << computeAllRes_ << std::endl;
1450  os << "========================================\n\n\n";
1451 }
1452 
1453 
1454 }} // end Anasazi namespace
1455 
1456 #endif /* ANASAZI_TraceMinBase_SOLMGR_HPP */
Pure virtual base class which describes the basic interface for a solver manager. ...
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Abstract base class for trace minimization eigensolvers.
int NEV
Number of unconverged eigenvalues.
TraceMinBaseSolMgr(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for TraceMinBaseSolMgr.
const RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
RCP< const MV > X
The current eigenvectors.
A special StatusTest for printing other status tests.
This class defines the interface required by an eigensolver and status test class to compute solution...
Structure to contain pointers to TraceMinBase state variables.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
RCP< const MV > V
The current basis.
basic_FancyOStream & setShowProcRank(const bool showProcRank)
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)
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
Status test for forming logical combinations of other status tests.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
void setGlobalStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
Basic implementation of the Anasazi::SortManager class.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
An exception class parent to all Anasazi exceptions.
bool isOrtho
Whether V has been projected and orthonormalized already.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Anasazi&#39;s templated, static class providing utilities for the solvers.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver&#39;s iterate() routine until ...
int numVecs
The number of computed eigenpairs.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Anasazi&#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...
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)
static magnitudeType prec()
ReturnType
Enumerated type used to pass back information from a solver manager.
const RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
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...
int curDim
The current dimension of the solver.
RCP< const MV > KV
The image of V under K.
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
The Anasazi::TraceMinBaseSolMgr provides an abstract base class for the TraceMin series of solver man...
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > RV
The current Ritz vectors.
ScalarType largestSafeShift
Largest safe shift.
RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
void push_back(const value_type &x)
Teuchos::Array< RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
RCP< const MV > KX
The image of the current eigenvectors under K.
RCP< const MV > R
The current residual vectors.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
void setLockingStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
Status test for forming logical combinations of other status tests.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
int getNumIters() const
Get the iteration count for the most recent call to solve().
bool isParameter(const std::string &name) const
const RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
This is an abstract base class for the trace minimization eigensolvers.
void setDebugStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
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.
Class which provides internal utilities for the Anasazi solvers.