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:
107  typedef Teuchos::ScalarTraits<ScalarType> SCT;
108  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
109  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
110 
111  public:
112 
114 
115 
159  Teuchos::ParameterList &pl );
160 
162  virtual ~TraceMinBaseSolMgr() {};
164 
166 
167 
170  return *problem_;
171  }
172 
174  int getNumIters() const {
175  return numIters_;
176  }
177 
185  Teuchos::Array<RCP<Teuchos::Time> > getTimers() const {
186  return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
187  }
188 
190 
192 
193 
217  ReturnType solve();
218 
220  void setGlobalStatusTest(const RCP< StatusTest<ScalarType,MV,OP> > &global);
221 
223  const RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
224 
226  void setLockingStatusTest(const RCP< StatusTest<ScalarType,MV,OP> > &locking);
227 
229  const RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
230 
232  void setDebugStatusTest(const RCP< StatusTest<ScalarType,MV,OP> > &debug);
233 
235  const RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
236 
238 
239  protected:
240  RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
241 
242  int numIters_;
243 
244  // Block variables
245  int blockSize_, numBlocks_, numRestartBlocks_;
246 
247  // Output variables
248  RCP<BasicOutputManager<ScalarType> > printer_;
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_;
281  RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
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(
291  const RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
292  const RCP<StatusTest<ScalarType,MV,OP> > &outputtest,
293  const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
294  Teuchos::ParameterList &plist
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
514  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
515  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"How To Choose Shift\"; valid options are \"Largest Converged\", \"Adjusted Ritz Values\", \"Ritz Values\".");
516 
517  // Consider clusters - if all eigenvalues are in one cluster, it's not expecially safe to shift
518  considerClusters_ = pl.get("Consider Clusters", true);
519 
520  // Use multiple shifts
521  useMultipleShifts_ = pl.get("Use Multiple Shifts", true);
522 
524  // Other parameters
525 
526  // which orthogonalization to use
527  ortho_ = pl.get("Orthogonalization", "SVQB");
528  TEUCHOS_TEST_FOR_EXCEPTION(ortho_ != "DGKS" && ortho_ != "SVQB" && ortho_ != "ICGS", std::invalid_argument,
529  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Orthogonalization\"; valid options are \"DGKS\", \"SVQB\", \"ICGS\".");
530 
531  strtmp = pl.get("Saddle Solver Type", "Projected Krylov");
532 
533  if(strtmp == "Projected Krylov")
534  saddleSolType_ = PROJECTED_KRYLOV_SOLVER;
535  else if(strtmp == "Schur Complement")
536  saddleSolType_ = SCHUR_COMPLEMENT_SOLVER;
537  else if(strtmp == "Block Diagonal Preconditioned Minres")
538  saddleSolType_ = BD_PREC_MINRES;
539  else if(strtmp == "HSS Preconditioned Gmres")
540  saddleSolType_ = HSS_PREC_GMRES;
541  else
542  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
543  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Saddle Solver Type\"; valid options are \"Projected Krylov\", \"Schur Complement\", and \"Block Diagonal Preconditioned Minres\".");
544 
545  projectAllVecs_ = pl.get("Project All Vectors", true);
546  projectLockedVecs_ = pl.get("Project Locked Vectors", true);
547  computeAllRes_ = pl.get("Compute All Residuals", true);
548  useRHSR_ = pl.get("Use Residual as RHS", false);
549  alpha_ = pl.get("HSS: alpha", 1.0);
550 
551  TEUCHOS_TEST_FOR_EXCEPTION(projectLockedVecs_ && ! projectAllVecs_, std::invalid_argument,
552  "Anasazi::TraceMinBaseSolMgr: If you want to project out the locked vectors, you should really project out ALL the vectors of X.");
553 
554  // Maximum number of inner iterations
555  maxKrylovIter_ = pl.get("Maximum Krylov Iterations", 200);
556  TEUCHOS_TEST_FOR_EXCEPTION(maxKrylovIter_ < 1, std::invalid_argument,
557  "Anasazi::TraceMinBaseSolMgr: \"Maximum Krylov Iterations\" must be greater than 0.");
558 
559  // Which eigenvalues we want to get
560  which_ = pl.get("Which", "SM");
561  TEUCHOS_TEST_FOR_EXCEPTION(which_ != "SM" && which_ != "LM", std::invalid_argument,
562  "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Which\"; valid options are \"SM\" and \"LM\".");
563 
564  // Test whether we are shifting without an operator K
565  // This is a really bad idea
566  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getOperator() == Teuchos::null && whenToShift_ != NEVER_SHIFT, std::invalid_argument,
567  "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\".");
568 
569 #ifdef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
570  // Test whether we are using a projected preconditioner with multiple Ritz shifts
571  // We can't currently do this for reasons that are complicated and are explained in the user manual
572  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER && useMultipleShifts_, std::invalid_argument,
573  "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.");
574 #else
575  // Test whether we are using a projected preconditioner without Belos.
576  // P Prec P should be positive definite if Prec is positive-definite,
577  // but it tends not to be in practice, presumably due to machine arithmetic
578  // As a result, we have to use pseudo-block gmres for now.
579  // Make sure it's available.
580  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER, std::invalid_argument,
581  "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)");
582 
583 
584 #endif
585 
586 
587 }
588 
589 
592 // solve()
593 template<class ScalarType, class MV, class OP>
594 ReturnType
596 {
597  typedef SolverUtils<ScalarType,MV,OP> msutils;
598 
599  const int nev = problem_->getNEV();
600 
601 #ifdef TEUCHOS_DEBUG
602  RCP<Teuchos::FancyOStream>
603  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
604  out->setShowAllFrontMatter(false).setShowProcRank(true);
605  *out << "Entering Anasazi::TraceMinBaseSolMgr::solve()\n";
606 #endif
607 
609  // Sort manager
610  RCP<BasicSort<MagnitudeType> > sorter = rcp( new BasicSort<MagnitudeType>("SM") );
611 
613  // Handle the spectral transformation if necessary
614  // TODO: Make sure we undo this before returning...
615  if(which_ == "LM")
616  {
617  RCP<const OP> swapHelper = problem_->getOperator();
618  problem_->setOperator(problem_->getM());
619  problem_->setM(swapHelper);
620  problem_->setProblem();
621  }
622 
624  // Status tests
625  //
626  // convergence
627  RCP<StatusTest<ScalarType,MV,OP> > convtest;
628  if (globalTest_ == Teuchos::null) {
629  if(which_ == "SM")
630  convtest = rcp( new StatusTestResNorm<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_) );
631  else
632  convtest = rcp( new StatusTestSpecTrans<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_,true,problem_->getOperator()) );
633  }
634  else {
635  convtest = globalTest_;
636  }
637  RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
638  = rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
639  // locking
640  RCP<StatusTest<ScalarType,MV,OP> > locktest;
641  if (useLocking_) {
642  if (lockingTest_ == Teuchos::null) {
643  if(which_ == "SM")
644  locktest = rcp( new StatusTestResNorm<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_) );
645  else
646  locktest = rcp( new StatusTestSpecTrans<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_,true,problem_->getOperator()) );
647  }
648  else {
649  locktest = lockingTest_;
650  }
651  }
652  // for a non-short-circuited OR test, the order doesn't matter
653  Teuchos::Array<RCP<StatusTest<ScalarType,MV,OP> > > alltests;
654  alltests.push_back(ordertest);
655  if (locktest != Teuchos::null) alltests.push_back(locktest);
656  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
657 
658  RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
660  // printing StatusTest
661  RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
662  if ( printer_->isVerbosity(Debug) ) {
663  outputtest = rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
664  }
665  else {
666  outputtest = rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
667  }
668 
670  // Orthomanager
671  RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
672  if (ortho_=="SVQB") {
673  ortho = rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
674  } else if (ortho_=="DGKS") {
675  ortho = rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
676  } else if (ortho_=="ICGS") {
677  ortho = rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
678  } else {
679  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): Invalid orthogonalization type.");
680  }
681 
683  // Parameter list
684  Teuchos::ParameterList plist;
685  setParameters(plist);
686 
688  // TraceMinBase solver
689  RCP<TraceMinBase<ScalarType,MV,OP> > tm_solver
690  = createSolver(sorter,outputtest,ortho,plist);
691  // set any auxiliary vectors defined in the problem
692  RCP< const MV > probauxvecs = problem_->getAuxVecs();
693  if (probauxvecs != Teuchos::null) {
694  tm_solver->setAuxVecs( Teuchos::tuple< RCP<const MV> >(probauxvecs) );
695  }
696 
698  // Storage
699  //
700  // lockvecs will contain eigenvectors that have been determined "locked" by the status test
701  int curNumLocked = 0;
702  RCP<MV> lockvecs;
703  // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking.
704  // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
705  // we will produce numnew random vectors, which will go into the space with the new basis.
706  // we will also need numnew storage for the image of these random vectors under A and M;
707  // columns [curlocked+1,curlocked+numnew] will be used for this storage
708  if (maxLocked_ > 0) {
709  lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
710  }
711  std::vector<MagnitudeType> lockvals;
712  //
713  // Restarting occurs under two scenarios: when the basis is full and after locking.
714  //
715  // For the former, a new basis of size blockSize*numRestartBlocks is generated using the current basis
716  // and the most significant primitive Ritz vectors (projected eigenvectors).
717  // [S,L] = eig(KK)
718  // S = [Sr St] // some for "r"estarting, some are "t"runcated
719  // newV = V*Sr
720  // KK_new = newV'*K*newV = Sr'*V'*K*V*Sr = Sr'*KK*Sr
721  // Therefore, the only multivector operation needed is for the generation of newV.
722  //
723  // * If the multiplication is explicit, it requires a workspace of blockSize*numRestartBlocks vectors.
724  // This space must be specifically allocated for that task, as we don't have any space of that size.
725  // It (workMV) will be allocated at the beginning of solve()
726  // * Optionally, the multiplication can be performed implicitly, via a Householder QR factorization of
727  // Sr. This can be done in situ, using the basis multivector contained in the solver. This requires
728  // that we cast away the const on the multivector returned from getState(). Workspace for this approach
729  // is a single vector. the solver's internal storage must be preserved (X,MX,KX,R), requiring us to
730  // allocate this vector.
731  //
732  // For the latter (restarting after locking), the new basis is the same size as existing basis. If numnew
733  // vectors are locked, they are deflated from the current basis and replaced with randomly generated
734  // vectors.
735  // [S,L] = eig(KK)
736  // S = [Sl Su] // partitioned: "l"ocked and "u"nlocked
737  // newL = V*Sl = X(locked)
738  // defV = V*Su
739  // augV = rand(numnew) // orthogonal to oldL,newL,defV,auxvecs
740  // newV = [defV augV]
741  // Kknew = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
742  // [augV'*K*defV augV'*K*augV]
743  // locked = [oldL newL]
744  // Clearly, this operation is more complicated than the previous.
745  // Here is a list of the significant computations that need to be performed:
746  // - newL will be put into space in lockvecs, but will be copied from getState().X at the end
747  // - defV,augV will be stored in workspace the size of the current basis.
748  // - M*augV and K*augV are needed; they will be stored in lockvecs. As a result, newL will
749  // not be put into lockvecs until the end.
750  //
751  // Therefore, we must allocate workMV when ((maxRestarts_ > 0) || (useLocking_ == true)) && inSituRestart == false
752  // It will be allocated to size (numBlocks-1)*blockSize
753  //
754 
755  // some consts and utils
756  const ScalarType ONE = SCT::one();
757  const ScalarType ZERO = SCT::zero();
758 
759  // go ahead and initialize the solution to nothing in case we throw an exception
761  sol.numVecs = 0;
762  problem_->setSolution(sol);
763 
764  int numRestarts = 0;
765 
766  // enter solve() iterations
767  {
768 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
769  Teuchos::TimeMonitor slvtimer(*_timerSolve);
770 #endif
771 
772  // tell tm_solver to iterate
773  while (1) {
774  try {
775  tm_solver->iterate();
776 
778  //
779  //
781  if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
782  throw AnasaziError("Anasazi::TraceMinBaseSolMgr::solve(): User-specified debug status test returned Passed.");
783  }
785  //
786  // check convergence next
787  //
789  else if (ordertest->getStatus() == Passed ) {
790  // we have convergence
791  // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
792  // ordertest->howMany() will tell us how many
793  break;
794  }
796  //
797  // check locking if we didn't converge or restart
798  //
800  else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
801 
802 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
803  Teuchos::TimeMonitor lcktimer(*_timerLocking);
804 #endif
805 
806  //
807  // get current state
808  TraceMinBaseState<ScalarType,MV> state = tm_solver->getState();
809  const int curdim = state.curDim;
810 
811  //
812  // get number,indices of vectors to be locked
813  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
814  "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() non-positive.");
815  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
816  "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
817  // we should have room to lock vectors; otherwise, locking should have been deactivated
818  TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
819  "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: locking not deactivated.");
820  //
821  // don't lock more than maxLocked_; we didn't allocate enough space.
822  std::vector<int> tmp_vector_int;
823  if (curNumLocked + locktest->howMany() > maxLocked_) {
824  // just use the first of them
825  for(int i=0; i<maxLocked_-curNumLocked; i++)
826  tmp_vector_int.push_back(locktest->whichVecs()[i]);
827 // tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
828  }
829  else {
830  tmp_vector_int = locktest->whichVecs();
831  }
832 
833  const std::vector<int> lockind(tmp_vector_int);
834  const int numNewLocked = lockind.size();
835  //
836  // generate indices of vectors left unlocked
837  // curind = [0,...,curdim-1] = UNION( lockind, unlockind )
838  const int numUnlocked = curdim-numNewLocked;
839  tmp_vector_int.resize(curdim);
840  for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
841  const std::vector<int> curind(tmp_vector_int); // curind = [0 ... curdim-1]
842  tmp_vector_int.resize(numUnlocked);
843  std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
844  const std::vector<int> unlockind(tmp_vector_int); // unlockind = [0 ... curdim-1] - lockind
845  tmp_vector_int.clear();
846 
847  //
848  // debug printing
849  if (printer_->isVerbosity(Debug)) {
850  printer_->print(Debug,"Locking vectors: ");
851  for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
852  printer_->print(Debug,"\n");
853  printer_->print(Debug,"Not locking vectors: ");
854  for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
855  printer_->print(Debug,"\n");
856  }
857 
858  // Copy eigenvalues we want to lock into lockvals
859  std::vector<Value<ScalarType> > allvals = tm_solver->getRitzValues();
860  for(unsigned int i=0; i<allvals.size(); i++)
861  printer_->stream(Debug) << "Ritz value[" << i << "] = " << allvals[i].realpart << std::endl;
862  for (int i=0; i<numNewLocked; i++) {
863  lockvals.push_back(allvals[lockind[i]].realpart);
864  }
865 
866  // Copy vectors we want to lock into lockvecs
867  RCP<const MV> newLocked = MVT::CloneView(*tm_solver->getRitzVectors(),lockind);
868  std::vector<int> indlock(numNewLocked);
869  for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
870  if(useHarmonic_)
871  {
872  RCP<MV> tempMV = MVT::CloneCopy(*newLocked);
873  ortho->normalizeMat(*tempMV);
874  MVT::SetBlock(*tempMV,indlock,*lockvecs);
875  }
876  else
877  {
878  MVT::SetBlock(*newLocked,indlock,*lockvecs);
879  }
880 
881  // Tell the StatusTestWithOrdering that things have been locked
882  // This is VERY important
883  // If this set of lines is removed, the code does not terminate correctly
884  if(noSort_)
885  {
886  for(unsigned int aliciaInd=0; aliciaInd<lockvals.size(); aliciaInd++)
887  {
888  lockvals[aliciaInd] = 0.0;
889  }
890  }
891  ordertest->setAuxVals(lockvals);
892 
893  // Set the auxiliary vectors so that we remain orthogonal to the ones we locked
894  // Remember to include any aux vecs provided by the user
895  curNumLocked += numNewLocked;
896 
897  if(ordertest->getStatus() == Passed) break;
898 
899  std::vector<int> curlockind(curNumLocked);
900  for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
901  RCP<const MV> curlocked = MVT::CloneView(*lockvecs,curlockind);
902 
903  Teuchos::Array< RCP<const MV> > aux;
904  if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
905  aux.push_back(curlocked);
906  tm_solver->setAuxVecs(aux);
907 
908  // Disable locking if we have locked the maximum number of things
909  printer_->stream(Debug) << "curNumLocked: " << curNumLocked << std::endl;
910  printer_->stream(Debug) << "maxLocked: " << maxLocked_ << std::endl;
911  if (curNumLocked == maxLocked_) {
912  // disabled locking now
913  combotest->removeTest(locktest);
914  locktest = Teuchos::null;
915  printer_->stream(Debug) << "Removed locking test\n";
916  }
917 
918  int newdim = numRestartBlocks_*blockSize_;
920  if(newdim <= numUnlocked)
921  {
922  if(useHarmonic_)
923  {
924  std::vector<int> desiredSubscripts(newdim);
925  for(int i=0; i<newdim; i++)
926  {
927  desiredSubscripts[i] = unlockind[i];
928  printer_->stream(Debug) << "H desiredSubscripts[" << i << "] = " << desiredSubscripts[i] << std::endl;
929  }
930  newstate.V = MVT::CloneView(*tm_solver->getRitzVectors(),desiredSubscripts);
931  newstate.curDim = newdim;
932  }
933  else
934  {
935  std::vector<int> desiredSubscripts(newdim);
936  for(int i=0; i<newdim; i++)
937  {
938  desiredSubscripts[i] = unlockind[i];
939  printer_->stream(Debug) << "desiredSubscripts[" << i << "] = " << desiredSubscripts[i] << std::endl;
940  }
941 
942  copyPartOfState(state, newstate, desiredSubscripts);
943  }
944  }
945  else
946  {
947  // TODO: Come back to this and make it more efficient
948 
949  // Replace the converged eigenvectors with random ones
950  int nrandom = newdim-numUnlocked;
951 
952  RCP<const MV> helperMV;
953  RCP<MV> totalV = MVT::Clone(*tm_solver->getRitzVectors(),newdim);
954 
955  // Holds old things that we're keeping
956  tmp_vector_int.resize(numUnlocked);
957  for(int i=0; i<numUnlocked; i++) tmp_vector_int[i] = i;
958  RCP<MV> oldV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
959 
960  // Copy over the old things
961  if(useHarmonic_)
962  helperMV = MVT::CloneView(*tm_solver->getRitzVectors(),unlockind);
963  else
964  helperMV = MVT::CloneView(*state.V,unlockind);
965  MVT::Assign(*helperMV,*oldV);
966 
967  // Holds random vectors we're generating
968  tmp_vector_int.resize(nrandom);
969  for(int i=0; i<nrandom; i++) tmp_vector_int[i] = i+numUnlocked;
970  RCP<MV> newV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
971 
972  // Create random things
973  MVT::MvRandom(*newV);
974 
975  newstate.V = totalV;
976  newstate.curDim = newdim;
977 
978  // Copy Ritz shifts
979  RCP< std::vector<ScalarType> > helperRS = rcp( new std::vector<ScalarType>(blockSize_) );
980  for(unsigned int i=0; i<(unsigned int)blockSize_; i++)
981  {
982  if(i < unlockind.size() && unlockind[i] < blockSize_)
983  (*helperRS)[i] = (*state.ritzShifts)[unlockind[i]];
984  else
985  (*helperRS)[i] = ZERO;
986  }
987  newstate.ritzShifts = helperRS;
988  }
989 
990  // Determine the largest safe shift
991  newstate.largestSafeShift = std::abs(lockvals[0]);
992  for(size_t i=0; i<lockvals.size(); i++)
993  newstate.largestSafeShift = std::max(std::abs(lockvals[i]), newstate.largestSafeShift);
994 
995  // Prepare new state, removing converged vectors
996  // TODO: Init will perform some unnecessary calculations; look into it
997  // TODO: The residual norms should be part of the state
998  newstate.NEV = state.NEV - numNewLocked;
999  tm_solver->initialize(newstate);
1000  } // end of locking
1002  //
1003  // check for restarting before locking: if we need to lock, it will happen after the restart
1004  //
1006  else if ( needToRestart(tm_solver) ) {
1007  // if performRestart returns false, we exceeded the maximum number of restarts
1008  if(performRestart(numRestarts, tm_solver) == false)
1009  break;
1010  } // end of restarting
1012  //
1013  // we returned from iterate(), but none of our status tests Passed.
1014  // something is wrong, and it is probably our fault.
1015  //
1017  else {
1018  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): Invalid return from tm_solver::iterate().");
1019  }
1020  }
1021  catch (const AnasaziError &err) {
1022  printer_->stream(Errors)
1023  << "Anasazi::TraceMinBaseSolMgr::solve() caught unexpected exception from Anasazi::TraceMinBase::iterate() at iteration " << tm_solver->getNumIters() << std::endl
1024  << err.what() << std::endl
1025  << "Anasazi::TraceMinBaseSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1026  return Unconverged;
1027  }
1028  }
1029 
1030  sol.numVecs = ordertest->howMany();
1031  if (sol.numVecs > 0) {
1032  sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
1033  sol.Espace = sol.Evecs;
1034  sol.Evals.resize(sol.numVecs);
1035  std::vector<MagnitudeType> vals(sol.numVecs);
1036 
1037  // copy them into the solution
1038  std::vector<int> which = ordertest->whichVecs();
1039  // indices between [0,blockSize) refer to vectors/values in the solver
1040  // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
1041  // everything has already been ordered by the solver; we just have to partition the two references
1042  std::vector<int> inlocked(0), insolver(0);
1043  for (unsigned int i=0; i<which.size(); i++) {
1044  if (which[i] >= 0) {
1045  TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): positive indexing mistake from ordertest.");
1046  insolver.push_back(which[i]);
1047  }
1048  else {
1049  // sanity check
1050  TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): negative indexing mistake from ordertest.");
1051  inlocked.push_back(which[i] + curNumLocked);
1052  }
1053  }
1054 
1055  TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): indexing mistake.");
1056 
1057  // set the vecs,vals in the solution
1058  if (insolver.size() > 0) {
1059  // set vecs
1060  int lclnum = insolver.size();
1061  std::vector<int> tosol(lclnum);
1062  for (int i=0; i<lclnum; i++) tosol[i] = i;
1063  RCP<const MV> v = MVT::CloneView(*tm_solver->getRitzVectors(),insolver);
1064  MVT::SetBlock(*v,tosol,*sol.Evecs);
1065  // set vals
1066  std::vector<Value<ScalarType> > fromsolver = tm_solver->getRitzValues();
1067  for (unsigned int i=0; i<insolver.size(); i++) {
1068  vals[i] = fromsolver[insolver[i]].realpart;
1069  }
1070  }
1071 
1072  // get the vecs,vals from locked storage
1073  if (inlocked.size() > 0) {
1074  int solnum = insolver.size();
1075  // set vecs
1076  int lclnum = inlocked.size();
1077  std::vector<int> tosol(lclnum);
1078  for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1079  RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1080  MVT::SetBlock(*v,tosol,*sol.Evecs);
1081  // set vals
1082  for (unsigned int i=0; i<inlocked.size(); i++) {
1083  vals[i+solnum] = lockvals[inlocked[i]];
1084  }
1085  }
1086 
1087  // undo the spectral transformation if necessary
1088  // if we really passed the solver Bx = \lambda A x, invert the eigenvalues
1089  if(which_ == "LM")
1090  {
1091  for(size_t i=0; i<vals.size(); i++)
1092  vals[i] = ONE/vals[i];
1093  }
1094 
1095  // sort the eigenvalues and permute the eigenvectors appropriately
1096  {
1097  std::vector<int> order(sol.numVecs);
1098  sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1099  // store the values in the Eigensolution
1100  for (int i=0; i<sol.numVecs; i++) {
1101  sol.Evals[i].realpart = vals[i];
1102  sol.Evals[i].imagpart = MT::zero();
1103  }
1104  // now permute the eigenvectors according to order
1105  msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1106  }
1107 
1108  // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1109  sol.index.resize(sol.numVecs,0);
1110  }
1111  }
1112 
1113  // print final summary
1114  tm_solver->currentStatus(printer_->stream(FinalSummary));
1115 
1116  printParameters(printer_->stream(FinalSummary));
1117 
1118  // print timing information
1119 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1120  if ( printer_->isVerbosity( TimingDetails ) ) {
1121  Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
1122  }
1123 #endif
1124 
1125  problem_->setSolution(sol);
1126  printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1127 
1128  // get the number of iterations taken for this call to solve().
1129  numIters_ = tm_solver->getNumIters();
1130 
1131  if (sol.numVecs < nev) {
1132  return Unconverged; // return from TraceMinBaseSolMgr::solve()
1133  }
1134  return Converged; // return from TraceMinBaseSolMgr::solve()
1135 }
1136 
1137 
1138 template <class ScalarType, class MV, class OP>
1139 void
1141  const RCP< StatusTest<ScalarType,MV,OP> > &global)
1142 {
1143  globalTest_ = global;
1144 }
1145 
1146 template <class ScalarType, class MV, class OP>
1147 const RCP< StatusTest<ScalarType,MV,OP> > &
1149 {
1150  return globalTest_;
1151 }
1152 
1153 template <class ScalarType, class MV, class OP>
1154 void
1156  const RCP< StatusTest<ScalarType,MV,OP> > &debug)
1157 {
1158  debugTest_ = debug;
1159 }
1160 
1161 template <class ScalarType, class MV, class OP>
1162 const RCP< StatusTest<ScalarType,MV,OP> > &
1164 {
1165  return debugTest_;
1166 }
1167 
1168 template <class ScalarType, class MV, class OP>
1169 void
1171  const RCP< StatusTest<ScalarType,MV,OP> > &locking)
1172 {
1173  lockingTest_ = locking;
1174 }
1175 
1176 template <class ScalarType, class MV, class OP>
1177 const RCP< StatusTest<ScalarType,MV,OP> > &
1179 {
1180  return lockingTest_;
1181 }
1182 
1183 template <class ScalarType, class MV, class OP>
1184 void TraceMinBaseSolMgr<ScalarType,MV,OP>::copyPartOfState(const TraceMinBaseState<ScalarType,MV>& oldState, TraceMinBaseState<ScalarType,MV>& newState, const std::vector<int> indToCopy) const
1185 {
1186  const ScalarType ONE = Teuchos::ScalarTraits<MagnitudeType>::one();
1187  const ScalarType ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
1188 
1189  newState.curDim = indToCopy.size();
1190  std::vector<int> fullIndices(oldState.curDim);
1191  for(int i=0; i<oldState.curDim; i++) fullIndices[i] = i;
1192 
1193  // Initialize with X.
1194  // Note that we didn't compute enough vectors of X, but we can very easily using the Ritz vectors.
1195  // That's why they're part of the state.
1196  // Note that there will ALWAYS be enough vectors
1197 
1198  // Helpful vectors for computing views and whatnot
1199  std::vector<int> oldIndices;
1200  std::vector<int> newIndices;
1201  for(int i=0; i<newState.curDim; i++)
1202  {
1203  if(indToCopy[i] < blockSize_)
1204  oldIndices.push_back(indToCopy[i]);
1205  else
1206  newIndices.push_back(indToCopy[i]);
1207  }
1208 
1209  int olddim = oldIndices.size();
1210  int newdim = newIndices.size();
1211 
1212  // If there are no new vectors being copied
1213  if(computeAllRes_)
1214  {
1215  newState.V = MVT::CloneView(*oldState.X, indToCopy);
1216  newState.R = MVT::CloneView(*oldState.R, indToCopy);
1217  newState.X = newState.V;
1218 
1219  if(problem_->getOperator() != Teuchos::null)
1220  {
1221  newState.KV = MVT::CloneView(*oldState.KX, indToCopy);
1222  newState.KX = newState.KV;
1223  }
1224  else
1225  {
1226  newState.KV = Teuchos::null;
1227  newState.KX = Teuchos::null;
1228  }
1229 
1230  if(problem_->getM() != Teuchos::null)
1231  {
1232  newState.MopV = MVT::CloneView(*oldState.MX, indToCopy);
1233  newState.MX = newState.MopV;
1234  }
1235  else
1236  {
1237  newState.MopV = Teuchos::null;
1238  newState.MX = Teuchos::null;
1239  }
1240  }
1241  else if(newdim == 0)
1242  {
1243  std::vector<int> blockind(blockSize_);
1244  for(int i=0; i<blockSize_; i++)
1245  blockind[i] = i;
1246 
1247  // Initialize with X
1248  newState.V = MVT::CloneView(*oldState.X, blockind);
1249  newState.KV = MVT::CloneView(*oldState.KX, blockind);
1250  newState.R = MVT::CloneView(*oldState.R, blockind);
1251  newState.X = MVT::CloneView(*newState.V, blockind);
1252  newState.KX = MVT::CloneView(*newState.KV, blockind);
1253 
1254  if(problem_->getM() != Teuchos::null)
1255  {
1256  newState.MopV = MVT::CloneView(*oldState.MX, blockind);
1257  newState.MX = MVT::CloneView(*newState.MopV, blockind);
1258  }
1259  else
1260  {
1261  newState.MopV = Teuchos::null;
1262  newState.MX = Teuchos::null;
1263  }
1264  }
1265  else
1266  {
1267  // More helpful vectors
1268  std::vector<int> oldPart(olddim);
1269  for(int i=0; i<olddim; i++) oldPart[i] = i;
1270  std::vector<int> newPart(newdim);
1271  for(int i=0; i<newdim; i++) newPart[i] = olddim+i;
1272 
1273  // Helpful multivectors for views and whatnot
1274  RCP<MV> helper = MVT::Clone(*oldState.V,newState.curDim);
1275  RCP<MV> oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1276  RCP<MV> newHelper = MVT::CloneViewNonConst(*helper,newPart);
1277  RCP<const MV> viewHelper;
1278 
1279  // Get the parts of the Ritz vectors we are interested in.
1280  Teuchos::SerialDenseMatrix<int,ScalarType> newRV(oldState.curDim,newdim);
1281  for(int r=0; r<oldState.curDim; r++)
1282  {
1283  for(int c=0; c<newdim; c++)
1284  newRV(r,c) = (*oldState.RV)(r,newIndices[c]);
1285  }
1286 
1287  // We're going to compute X as V*RitzVecs
1288  viewHelper = MVT::CloneView(*oldState.V,fullIndices);
1289  MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1290  viewHelper = MVT::CloneView(*oldState.X,oldIndices);
1291  MVT::Assign(*viewHelper,*oldHelper);
1292  newState.V = MVT::CloneCopy(*helper);
1293 
1294  // Also compute KX as KV*RitzVecs
1295  viewHelper = MVT::CloneView(*oldState.KV,fullIndices);
1296  MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1297  viewHelper = MVT::CloneView(*oldState.KX,oldIndices);
1298  MVT::Assign(*viewHelper,*oldHelper);
1299  newState.KV = MVT::CloneCopy(*helper);
1300 
1301  // Do the same with MX if necessary
1302  if(problem_->getM() != Teuchos::null)
1303  {
1304  viewHelper = MVT::CloneView(*oldState.MopV,fullIndices);
1305  MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1306  viewHelper = MVT::CloneView(*oldState.MX,oldIndices);
1307  MVT::Assign(*viewHelper,*oldHelper);
1308  newState.MopV = MVT::CloneCopy(*helper);
1309  }
1310  else
1311  newState.MopV = newState.V;
1312 
1313  // Get X, MX, KX
1314  std::vector<int> blockVec(blockSize_);
1315  for(int i=0; i<blockSize_; i++) blockVec[i] = i;
1316  newState.X = MVT::CloneView(*newState.V,blockVec);
1317  newState.KX = MVT::CloneView(*newState.KV,blockVec);
1318  newState.MX = MVT::CloneView(*newState.MopV,blockVec);
1319 
1320  // Update the residuals
1321  if(blockSize_-oldIndices.size() > 0)
1322  {
1323  // There are vectors we have not computed the residual for yet
1324  newPart.resize(blockSize_-oldIndices.size());
1325  helper = MVT::Clone(*oldState.V,blockSize_);
1326  oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1327  newHelper = MVT::CloneViewNonConst(*helper,newPart);
1328 
1329  RCP<MV> scaledMV = MVT::CloneCopy(*newState.MX,newPart);
1330  RCP<const MV> localKV = MVT::CloneView(*newState.KX,newPart);
1331  std::vector<ScalarType> scalarVec(blockSize_-oldIndices.size());
1332  for(unsigned int i=0; i<(unsigned int)blockSize_-oldIndices.size(); i++) scalarVec[i] = (*oldState.T)[newPart[i]];
1333  MVT::MvScale(*scaledMV,scalarVec);
1334 
1335  helper = MVT::Clone(*oldState.V,blockSize_);
1336  oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1337  newHelper = MVT::CloneViewNonConst(*helper,newPart);
1338  MVT::MvAddMv(ONE,*localKV,-ONE,*scaledMV,*newHelper);
1339  viewHelper = MVT::CloneView(*oldState.R,oldIndices);
1340  MVT::Assign(*viewHelper,*oldHelper);
1341  newState.R = MVT::CloneCopy(*helper);
1342  }
1343  else
1344  newState.R = oldState.R;
1345  }
1346 
1347  // Since we are setting V:=X, V is orthonormal
1348  newState.isOrtho = true;
1349 
1350  // Get the first eigenvalues
1351  RCP< std::vector<ScalarType> > helperT = rcp( new std::vector<ScalarType>(newState.curDim) );
1352  for(int i=0; i<newState.curDim; i++) (*helperT)[i] = (*oldState.T)[indToCopy[i]];
1353  newState.T = helperT;
1354 
1355  // X'KX is diag(T)
1356  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.curDim,newState.curDim) );
1357  for(int i=0; i<newState.curDim; i++)
1358  (*newKK)(i,i) = (*newState.T)[i];
1359  newState.KK = newKK;
1360 
1361  // The associated Ritz vectors are I
1362  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newRV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.curDim,newState.curDim) );
1363  for(int i=0; i<newState.curDim; i++)
1364  (*newRV)(i,i) = ONE;
1365  newState.RV = newRV;
1366 
1367  // Get the Ritz shifts
1368  RCP< std::vector<ScalarType> > helperRS = rcp( new std::vector<ScalarType>(blockSize_) );
1369  for(int i=0; i<blockSize_; i++)
1370  {
1371  if(indToCopy[i] < blockSize_)
1372  (*helperRS)[i] = (*oldState.ritzShifts)[indToCopy[i]];
1373  else
1374  (*helperRS)[i] = ZERO;
1375  }
1376  newState.ritzShifts = helperRS;
1377 }
1378 
1379 
1380 template <class ScalarType, class MV, class OP>
1381 void TraceMinBaseSolMgr<ScalarType,MV,OP>::setParameters(Teuchos::ParameterList &pl) const
1382 {
1383  pl.set("Block Size", blockSize_);
1384  pl.set("Num Blocks", numBlocks_);
1385  pl.set("Num Restart Blocks", numRestartBlocks_);
1386  pl.set("When To Shift", whenToShift_);
1387  pl.set("Trace Threshold", traceThresh_);
1388  pl.set("Shift Tolerance", shiftTol_);
1389  pl.set("Relative Shift Tolerance", relShiftTol_);
1390  pl.set("Shift Norm", shiftNorm_);
1391  pl.set("How To Choose Shift", howToShift_);
1392  pl.set("Consider Clusters", considerClusters_);
1393  pl.set("Use Multiple Shifts", useMultipleShifts_);
1394  pl.set("Saddle Solver Type", saddleSolType_);
1395  pl.set("Project All Vectors", projectAllVecs_);
1396  pl.set("Project Locked Vectors", projectLockedVecs_);
1397  pl.set("Compute All Residuals", computeAllRes_);
1398  pl.set("Use Residual as RHS", useRHSR_);
1399  pl.set("Use Harmonic Ritz Values", useHarmonic_);
1400  pl.set("Maximum Krylov Iterations", maxKrylovIter_);
1401  pl.set("HSS: alpha", alpha_);
1402 }
1403 
1404 
1405 template <class ScalarType, class MV, class OP>
1406 void TraceMinBaseSolMgr<ScalarType,MV,OP>::printParameters(std::ostream &os) const
1407 {
1408  os << "\n\n\n";
1409  os << "========================================\n";
1410  os << "========= TraceMin parameters ==========\n";
1411  os << "========================================\n";
1412  os << "=========== Block parameters ===========\n";
1413  os << "Block Size: " << blockSize_ << std::endl;
1414  os << "Num Blocks: " << numBlocks_ << std::endl;
1415  os << "Num Restart Blocks: " << numRestartBlocks_ << std::endl;
1416  os << "======== Convergence parameters ========\n";
1417  os << "Convergence Tolerance: " << convTol_ << std::endl;
1418  os << "Relative Convergence Tolerance: " << relConvTol_ << std::endl;
1419  os << "========== Locking parameters ==========\n";
1420  os << "Use Locking: " << useLocking_ << std::endl;
1421  os << "Locking Tolerance: " << lockTol_ << std::endl;
1422  os << "Relative Locking Tolerance: " << relLockTol_ << std::endl;
1423  os << "Max Locked: " << maxLocked_ << std::endl;
1424  os << "Locking Quorum: " << lockQuorum_ << std::endl;
1425  os << "========== Shifting parameters =========\n";
1426  os << "When To Shift: ";
1427  if(whenToShift_ == NEVER_SHIFT) os << "Never\n";
1428  else if(whenToShift_ == SHIFT_WHEN_TRACE_LEVELS) os << "After Trace Levels\n";
1429  else if(whenToShift_ == SHIFT_WHEN_RESID_SMALL) os << "Residual Becomes Small\n";
1430  else if(whenToShift_ == ALWAYS_SHIFT) os << "Always\n";
1431  os << "Consider Clusters: " << considerClusters_ << std::endl;
1432  os << "Trace Threshohld: " << traceThresh_ << std::endl;
1433  os << "Shift Tolerance: " << shiftTol_ << std::endl;
1434  os << "Relative Shift Tolerance: " << relShiftTol_ << std::endl;
1435  os << "How To Choose Shift: ";
1436  if(howToShift_ == LARGEST_CONVERGED_SHIFT) os << "Largest Converged\n";
1437  else if(howToShift_ == ADJUSTED_RITZ_SHIFT) os << "Adjusted Ritz Values\n";
1438  else if(howToShift_ == RITZ_VALUES_SHIFT) os << "Ritz Values\n";
1439  os << "Use Multiple Shifts: " << useMultipleShifts_ << std::endl;
1440  os << "=========== Other parameters ===========\n";
1441  os << "Orthogonalization: " << ortho_ << std::endl;
1442  os << "Saddle Solver Type: ";
1443  if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER) os << "Projected Krylov\n";
1444  else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER) os << "Schur Complement\n";
1445  os << "Project All Vectors: " << projectAllVecs_ << std::endl;
1446  os << "Project Locked Vectors: " << projectLockedVecs_ << std::endl;
1447  os << "Compute All Residuals: " << computeAllRes_ << std::endl;
1448  os << "========================================\n\n\n";
1449 }
1450 
1451 
1452 }} // end Anasazi namespace
1453 
1454 #endif /* ANASAZI_TraceMinBase_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 .
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
const RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
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.
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
const 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...
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.
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
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...
static void Assign(const MV &A, MV &mv)
mv := A
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.
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.
Teuchos::Array< RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A status test for testing the norm of the eigenvectors residuals.
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
int curDim
The current dimension of the solver.
RCP< const MV > KV
The image of V under K.
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...
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.
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.
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...
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.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
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.
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...
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.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Basic implementation of the Anasazi::OrthoManager class.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Class which provides internal utilities for the Anasazi solvers.
const RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.