Anasazi  Version of the Day
AnasaziTraceMinBase.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 // TODO: Modify code so R is not computed unless needed
30 
35 #ifndef ANASAZI_TRACEMIN_BASE_HPP
36 #define ANASAZI_TRACEMIN_BASE_HPP
37 
38 #include "AnasaziBasicSort.hpp"
39 #include "AnasaziConfigDefs.hpp"
40 #include "AnasaziEigensolver.hpp"
46 #include "AnasaziSolverUtils.hpp"
48 #include "AnasaziTraceMinTypes.hpp"
49 #include "AnasaziTypes.hpp"
50 
51 #include "Teuchos_ParameterList.hpp"
52 #include "Teuchos_ScalarTraits.hpp"
53 #include "Teuchos_SerialDenseMatrix.hpp"
54 #include "Teuchos_SerialDenseSolver.hpp"
55 #include "Teuchos_TimeMonitor.hpp"
56 
57 using Teuchos::RCP;
58 using Teuchos::rcp;
59 
60 namespace Anasazi {
71 namespace Experimental {
72 
74 
75 
80  template <class ScalarType, class MV>
83  int curDim;
88  RCP<const MV> V;
90  RCP<const MV> KV;
92  RCP<const MV> MopV;
94  RCP<const MV> X;
96  RCP<const MV> KX;
98  RCP<const MV> MX;
100  RCP<const MV> R;
102  RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
108  RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > KK;
110  RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > RV;
112  bool isOrtho;
114  int NEV;
116  ScalarType largestSafeShift;
118  RCP< const std::vector<ScalarType> > ritzShifts;
119  TraceMinBaseState() : curDim(0), V(Teuchos::null), KV(Teuchos::null), MopV(Teuchos::null),
120  X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null), R(Teuchos::null),
121  T(Teuchos::null), KK(Teuchos::null), RV(Teuchos::null), isOrtho(false),
122  NEV(0), largestSafeShift(Teuchos::ScalarTraits<ScalarType>::zero()),
123  ritzShifts(Teuchos::null) {}
124  };
125 
127 
129 
130 
143  class TraceMinBaseInitFailure : public AnasaziError {public:
144  TraceMinBaseInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
145  {}};
146 
150  class TraceMinBaseOrthoFailure : public AnasaziError {public:
151  TraceMinBaseOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
152  {}};
153 
155 
168  template <class ScalarType, class MV, class OP>
169  class TraceMinBase : public Eigensolver<ScalarType,MV,OP> {
170  public:
172 
173 
201  TraceMinBase( const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
202  const RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
203  const RCP<OutputManager<ScalarType> > &printer,
204  const RCP<StatusTest<ScalarType,MV,OP> > &tester,
205  const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
206  Teuchos::ParameterList &params
207  );
208 
210  virtual ~TraceMinBase();
212 
213 
215 
216 
240  void iterate();
241 
242  void harmonicIterate();
243 
266  void initialize(TraceMinBaseState<ScalarType,MV>& newstate);
267 
268  void harmonicInitialize(TraceMinBaseState<ScalarType,MV> newstate);
269 
273  void initialize();
274 
290  bool isInitialized() const;
291 
300  TraceMinBaseState<ScalarType,MV> getState() const;
301 
303 
304 
306 
307 
309  int getNumIters() const;
310 
312  void resetNumIters();
313 
321  RCP<const MV> getRitzVectors();
322 
328  std::vector<Value<ScalarType> > getRitzValues();
329 
330 
339  std::vector<int> getRitzIndex();
340 
341 
347  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
348 
349 
355  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
356 
357 
365  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
366 
372  int getCurSubspaceDim() const;
373 
375  int getMaxSubspaceDim() const;
376 
378 
379 
381 
382 
384  void setStatusTest(RCP<StatusTest<ScalarType,MV,OP> > test);
385 
387  RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
388 
390  const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
391 
401  void setBlockSize(int blockSize);
402 
404  int getBlockSize() const;
405 
423  void setAuxVecs(const Teuchos::Array<RCP<const MV> > &auxvecs);
424 
426  Teuchos::Array<RCP<const MV> > getAuxVecs() const;
427 
429 
431 
432 
439  void setSize(int blockSize, int numBlocks);
440 
442 
443 
445  void currentStatus(std::ostream &os);
446 
448 
449  protected:
450  //
451  // Convenience typedefs
452  //
456  typedef Teuchos::ScalarTraits<ScalarType> SCT;
457  typedef typename SCT::magnitudeType MagnitudeType;
458  typedef TraceMinRitzOp<ScalarType,MV,OP> tracemin_ritz_op_type;
459  typedef SaddleContainer<ScalarType,MV> saddle_container_type;
460  typedef SaddleOperator<ScalarType,MV,tracemin_ritz_op_type> saddle_op_type;
461  const MagnitudeType ONE;
462  const MagnitudeType ZERO;
463  const MagnitudeType NANVAL;
464  //
465  // Classes inputed through constructor that define the eigenproblem to be solved.
466  //
467  const RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
468  const RCP<SortManager<MagnitudeType> > sm_;
469  const RCP<OutputManager<ScalarType> > om_;
470  RCP<StatusTest<ScalarType,MV,OP> > tester_;
471  const RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
472  //
473  // Information obtained from the eigenproblem
474  //
475  RCP<const OP> Op_;
476  RCP<const OP> MOp_;
477  RCP<const OP> Prec_;
478  bool hasM_;
479  //
480  // Internal timers
481  // TODO: Fix the timers
482  //
483  RCP<Teuchos::Time> timerOp_, timerMOp_, timerSaddle_, timerSortEval_, timerDS_,
484  timerLocal_, timerCompRes_, timerOrtho_, timerInit_;
485  //
486  // Internal structs
487  // TODO: Fix the checks
488  //
489  struct CheckList {
490  bool checkV, checkX, checkMX,
491  checkKX, checkQ, checkKK;
492  CheckList() : checkV(false),checkX(false),
493  checkMX(false),checkKX(false),
494  checkQ(false),checkKK(false) {};
495  };
496  //
497  // Internal methods
498  //
499  std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
500 
501  //
502  // Counters
503  //
504  int count_ApplyOp_, count_ApplyM_;
505 
506  //
507  // Algorithmic parameters.
508  //
509  // blockSize_ is the solver block size; it controls the number of vectors added to the basis on each iteration.
510  int blockSize_;
511  // numBlocks_ is the size of the allocated space for the basis, in blocks.
512  int numBlocks_;
513 
514  //
515  // Current solver state
516  //
517  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
518  // is capable of running; _initialize is controlled by the initialize() member method
519  // For the implications of the state of initialized_, please see documentation for initialize()
520  bool initialized_;
521  //
522  // curDim_ reflects how much of the current basis is valid
523  // NOTE: 0 <= curDim_ <= blockSize_*numBlocks_
524  // this also tells us how many of the values in theta_ are valid Ritz values
525  int curDim_;
526  //
527  // State Multivecs
528  // V is the current basis
529  // X is the current set of eigenvectors
530  // R is the residual
531  RCP<MV> X_, KX_, MX_, KV_, MV_, R_, V_;
532  //
533  // Projected matrices
534  //
535  RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_, ritzVecs_;
536  //
537  // auxiliary vectors
538  Teuchos::Array<RCP<const MV> > auxVecs_, MauxVecs_;
539  int numAuxVecs_;
540  //
541  // Number of iterations that have been performed.
542  int iter_;
543  //
544  // Current eigenvalues, residual norms
545  std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
546  //
547  // are the residual norms current with the residual?
548  bool Rnorms_current_, R2norms_current_;
549 
550  // TraceMin specific variables
551  RCP<tracemin_ritz_op_type> ritzOp_; // Operator which incorporates Ritz shifts
552  enum SaddleSolType saddleSolType_; // Type of saddle point solver being used
553  bool previouslyLeveled_; // True if the trace already leveled off
554  MagnitudeType previousTrace_; // Trace of X'KX at the previous iteration
555  bool posSafeToShift_, negSafeToShift_; // Whether there are multiple clusters
556  MagnitudeType largestSafeShift_; // The largest shift considered to be safe - generally the biggest converged eigenvalue
557  int NEV_; // Number of eigenvalues we seek - used in computation of trace
558  std::vector<ScalarType> ritzShifts_; // Current Ritz shifts
559 
560  // This is only used if we're using the Schur complement solver
561  RCP<MV> Z_;
562 
563  // User provided TraceMin parameters
564  enum WhenToShiftType whenToShift_; // What triggers a Ritz shift
565  enum HowToShiftType howToShift_; // Strategy for choosing the Ritz shifts
566  bool useMultipleShifts_; // Whether to use one Ritz shift or many
567  bool considerClusters_; // Don't shift if there is only one cluster
568  bool projectAllVecs_; // Use all vectors in projected minres, or just 1
569  bool projectLockedVecs_; // Project locked vectors too in minres; does nothing if projectAllVecs = false
570  bool computeAllRes_; // Compute all residuals, or just blocksize ones - helps make Ritz shifts safer
571  bool useRHSR_; // Use -R as the RHS of projected minres rather than AX
572  bool useHarmonic_;
573  MagnitudeType traceThresh_;
574  MagnitudeType alpha_;
575 
576  // Variables that are only used if we're shifting when the residual gets small
577  // TODO: These are currently unused
578  std::string shiftNorm_; // Measure 2-norm or M-norm of residual
579  MagnitudeType shiftThresh_; // How small must the residual be?
580  bool useRelShiftThresh_; // Are we scaling the shift threshold by the eigenvalues?
581 
582  // TraceMin specific functions
583  // Returns the trace of KK = X'KX
584  ScalarType getTrace() const;
585  // Returns true if the change in trace is very small (or was very small at one point)
586  // TODO: Check whether I want to return true if the trace WAS small
587  bool traceLeveled();
588  // Get the residuals of each cluster of eigenvalues
589  // TODO: Figure out whether I want to use these for all eigenvalues or just the first
590  std::vector<ScalarType> getClusterResids();
591  // Computes the Ritz shifts, which is not a trivial task
592  // TODO: Make it safer for indefinite matrices maybe?
593  void computeRitzShifts(const std::vector<ScalarType>& clusterResids);
594  // Compute the tolerance for the inner solves
595  // TODO: Come back to this and make sure it does what I want it to
596  std::vector<ScalarType> computeTol();
597  // Solves a saddle point problem
598  void solveSaddlePointProblem(RCP<MV> Delta);
599  // Solves a saddle point problem by using unpreconditioned projected minres
600  void solveSaddleProj(RCP<MV> Delta) const;
601  // Solves a saddle point problem by using preconditioned projected...Krylov...something
602  // TODO: Fix this. Replace minres with gmres or something.
603  void solveSaddleProjPrec(RCP<MV> Delta) const;
604  // Solves a saddle point problem by explicitly forming the inexact Schur complement
605  void solveSaddleSchur (RCP<MV> Delta) const;
606  // Solves a saddle point problem with a block diagonal preconditioner
607  void solveSaddleBDPrec (RCP<MV> Delta) const;
608  // Solves a saddle point problem with a Hermitian/non-Hermitian splitting preconditioner
609  void solveSaddleHSSPrec (RCP<MV> Delta) const;
610  // Computes KK = X'KX
611  void computeKK();
612  // Computes the eigenpairs of KK
613  void computeRitzPairs();
614  // Computes X = V evecs
615  void computeX();
616  // Updates KX and MX without a matvec by MAGIC (or by multiplying KV and MV by evecs)
617  void updateKXMX();
618  // Updates the residual
619  void updateResidual();
620  // Adds Delta to the basis
621  // TraceMin and TraceMinDavidson each do this differently, which is why this is pure virtual
622  virtual void addToBasis(const RCP<const MV> Delta) =0;
623  virtual void harmonicAddToBasis(const RCP<const MV> Delta) =0;
624  };
625 
628  //
629  // Implementations
630  //
633 
634 
636  // Constructor
637  // TODO: Allow the users to supply their own linear solver
638  // TODO: Add additional checking for logic errors (like trying to use gmres with multiple shifts)
639  template <class ScalarType, class MV, class OP>
641  const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
642  const RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
643  const RCP<OutputManager<ScalarType> > &printer,
644  const RCP<StatusTest<ScalarType,MV,OP> > &tester,
645  const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
646  Teuchos::ParameterList &params
647  ) :
648  ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
649  ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
650  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
651  // problem, tools
652  problem_(problem),
653  sm_(sorter),
654  om_(printer),
655  tester_(tester),
656  orthman_(ortho),
657  // timers, counters
658 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
659  timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Operation Op*x")),
660  timerMOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Operation M*x")),
661  timerSaddle_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Solving saddle point problem")),
662  timerSortEval_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Sorting eigenvalues")),
663  timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Direct solve")),
664  timerLocal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Local update")),
665  timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Computing residuals")),
666  timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Orthogonalization")),
667  timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Initialization")),
668 #endif
669  count_ApplyOp_(0),
670  count_ApplyM_(0),
671  // internal data
672  blockSize_(0),
673  numBlocks_(0),
674  initialized_(false),
675  curDim_(0),
676  auxVecs_( Teuchos::Array<RCP<const MV> >(0) ),
677  MauxVecs_( Teuchos::Array<RCP<const MV> >(0) ),
678  numAuxVecs_(0),
679  iter_(0),
680  Rnorms_current_(false),
681  R2norms_current_(false),
682  // TraceMin variables
683  previouslyLeveled_(false),
684  previousTrace_(ZERO),
685  posSafeToShift_(false),
686  negSafeToShift_(false),
687  largestSafeShift_(ZERO),
688  Z_(Teuchos::null)
689  {
690  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
691  "Anasazi::TraceMinBase::constructor: user passed null problem pointer.");
692  TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
693  "Anasazi::TraceMinBase::constructor: user passed null sort manager pointer.");
694  TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
695  "Anasazi::TraceMinBase::constructor: user passed null output manager pointer.");
696  TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
697  "Anasazi::TraceMinBase::constructor: user passed null status test pointer.");
698  TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
699  "Anasazi::TraceMinBase::constructor: user passed null orthogonalization manager pointer.");
700  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
701  "Anasazi::TraceMinBase::constructor: problem is not hermitian.");
702 
703  // get the problem operators
704  Op_ = problem_->getOperator();
705  MOp_ = problem_->getM();
706  Prec_ = problem_->getPrec();
707  hasM_ = (MOp_ != Teuchos::null);
708 
709  // Set the saddle point solver parameters
710  saddleSolType_ = params.get("Saddle Solver Type", PROJECTED_KRYLOV_SOLVER);
711  TEUCHOS_TEST_FOR_EXCEPTION(saddleSolType_ != PROJECTED_KRYLOV_SOLVER && saddleSolType_ != SCHUR_COMPLEMENT_SOLVER && saddleSolType_ != BD_PREC_MINRES && saddleSolType_ != HSS_PREC_GMRES, std::invalid_argument,
712  "Anasazi::TraceMin::constructor: Invalid value for \"Saddle Solver Type\"; valid options are PROJECTED_KRYLOV_SOLVER, SCHUR_COMPLEMENT_SOLVER, and BD_PREC_MINRES.");
713 
714  // Set the Ritz shift parameters
715  whenToShift_ = params.get("When To Shift", ALWAYS_SHIFT);
716  TEUCHOS_TEST_FOR_EXCEPTION(whenToShift_ != NEVER_SHIFT && whenToShift_ != SHIFT_WHEN_TRACE_LEVELS && whenToShift_ != SHIFT_WHEN_RESID_SMALL && whenToShift_ != ALWAYS_SHIFT, std::invalid_argument,
717  "Anasazi::TraceMin::constructor: Invalid value for \"When To Shift\"; valid options are \"NEVER_SHIFT\", \"SHIFT_WHEN_TRACE_LEVELS\", \"SHIFT_WHEN_RESID_SMALL\", and \"ALWAYS_SHIFT\".");
718 
719  traceThresh_ = params.get("Trace Threshold", 2e-2);
720  TEUCHOS_TEST_FOR_EXCEPTION(traceThresh_ < 0, std::invalid_argument,
721  "Anasazi::TraceMin::constructor: Invalid value for \"Trace Threshold\"; Must be positive.");
722 
723  howToShift_ = params.get("How To Choose Shift", ADJUSTED_RITZ_SHIFT);
724  TEUCHOS_TEST_FOR_EXCEPTION(howToShift_ != LARGEST_CONVERGED_SHIFT && howToShift_ != ADJUSTED_RITZ_SHIFT && howToShift_ != RITZ_VALUES_SHIFT, std::invalid_argument,
725  "Anasazi::TraceMin::constructor: Invalid value for \"How To Choose Shift\"; valid options are \"LARGEST_CONVERGED_SHIFT\", \"ADJUSTED_RITZ_SHIFT\", \"RITZ_VALUES_SHIFT\".");
726 
727  useMultipleShifts_ = params.get("Use Multiple Shifts", true);
728 
729  shiftThresh_ = params.get("Shift Threshold", 1e-2);
730  useRelShiftThresh_ = params.get("Relative Shift Threshold", true);
731  shiftNorm_ = params.get("Shift Norm", "2");
732  TEUCHOS_TEST_FOR_EXCEPTION(shiftNorm_ != "2" && shiftNorm_ != "M", std::invalid_argument,
733  "Anasazi::TraceMin::constructor: Invalid value for \"Shift Norm\"; valid options are \"2\", \"M\".");
734 
735  considerClusters_ = params.get("Consider Clusters", true);
736 
737  projectAllVecs_ = params.get("Project All Vectors", true);
738  projectLockedVecs_ = params.get("Project Locked Vectors", true);
739  useRHSR_ = params.get("Use Residual as RHS", true);
740  useHarmonic_ = params.get("Use Harmonic Ritz Values", false);
741  computeAllRes_ = params.get("Compute All Residuals", true);
742 
743  // set the block size and allocate data
744  int bs = params.get("Block Size", problem_->getNEV());
745  int nb = params.get("Num Blocks", 1);
746  setSize(bs,nb);
747 
748  NEV_ = problem_->getNEV();
749 
750  // Create the Ritz shift operator
751  ritzOp_ = rcp (new tracemin_ritz_op_type (Op_, MOp_, Prec_));
752 
753  // Set the maximum number of inner iterations
754  const int innerMaxIts = params.get ("Maximum Krylov Iterations", 200);
755  ritzOp_->setMaxIts (innerMaxIts);
756 
757  alpha_ = params.get ("HSS: alpha", ONE);
758  }
759 
760 
762  // Destructor
763  template <class ScalarType, class MV, class OP>
765 
766 
768  // Set the block size
769  // This simply calls setSize(), modifying the block size while retaining the number of blocks.
770  template <class ScalarType, class MV, class OP>
772  {
773  TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::TraceMinBase::setSize(blocksize,numblocks): blocksize must be strictly positive.");
774  setSize(blockSize,numBlocks_);
775  }
776 
777 
779  // Return the current auxiliary vectors
780  template <class ScalarType, class MV, class OP>
781  Teuchos::Array<RCP<const MV> > TraceMinBase<ScalarType,MV,OP>::getAuxVecs() const {
782  return auxVecs_;
783  }
784 
785 
787  // return the current block size
788  template <class ScalarType, class MV, class OP>
790  return(blockSize_);
791  }
792 
793 
795  // return eigenproblem
796  template <class ScalarType, class MV, class OP>
798  return(*problem_);
799  }
800 
801 
803  // return max subspace dim
804  template <class ScalarType, class MV, class OP>
806  return blockSize_*numBlocks_;
807  }
808 
810  // return current subspace dim
811  template <class ScalarType, class MV, class OP>
813  if (!initialized_) return 0;
814  return curDim_;
815  }
816 
817 
819  // return ritz residual 2-norms
820  template <class ScalarType, class MV, class OP>
821  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
823  return getRes2Norms();
824  }
825 
826 
828  // return ritz index
829  template <class ScalarType, class MV, class OP>
831  std::vector<int> ret(curDim_,0);
832  return ret;
833  }
834 
835 
837  // return ritz values
838  template <class ScalarType, class MV, class OP>
839  std::vector<Value<ScalarType> > TraceMinBase<ScalarType,MV,OP>::getRitzValues() {
840  std::vector<Value<ScalarType> > ret(curDim_);
841  for (int i=0; i<curDim_; ++i) {
842  ret[i].realpart = theta_[i];
843  ret[i].imagpart = ZERO;
844  }
845  return ret;
846  }
847 
848 
850  // return pointer to ritz vectors
851  template <class ScalarType, class MV, class OP>
853  return X_;
854  }
855 
856 
858  // reset number of iterations
859  template <class ScalarType, class MV, class OP>
861  iter_=0;
862  }
863 
864 
866  // return number of iterations
867  template <class ScalarType, class MV, class OP>
869  return(iter_);
870  }
871 
872 
874  // return state pointers
875  template <class ScalarType, class MV, class OP>
878  state.curDim = curDim_;
879  state.V = V_;
880  state.X = X_;
881  if (Op_ != Teuchos::null) {
882  state.KV = KV_;
883  state.KX = KX_;
884  }
885  else {
886  state.KV = Teuchos::null;
887  state.KX = Teuchos::null;
888  }
889  if (hasM_) {
890  state.MopV = MV_;
891  state.MX = MX_;
892  }
893  else {
894  state.MopV = Teuchos::null;
895  state.MX = Teuchos::null;
896  }
897  state.R = R_;
898  state.KK = KK_;
899  state.RV = ritzVecs_;
900  if (curDim_ > 0) {
901  state.T = rcp(new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
902  } else {
903  state.T = rcp(new std::vector<MagnitudeType>(0));
904  }
905  state.ritzShifts = rcp(new std::vector<MagnitudeType>(ritzShifts_.begin(),ritzShifts_.begin()+blockSize_) );
906  state.isOrtho = true;
907  state.largestSafeShift = largestSafeShift_;
908 
909  return state;
910  }
911 
912 
914  // Perform TraceMinBase iterations until the StatusTest tells us to stop.
915  template <class ScalarType, class MV, class OP>
917  {
918  if(useHarmonic_)
919  {
920  harmonicIterate();
921  return;
922  }
923 
924  //
925  // Initialize solver state
926  if (initialized_ == false) {
927  initialize();
928  }
929 
930  // as a data member, this would be redundant and require synchronization with
931  // blockSize_ and numBlocks_; we'll use a constant here.
932  const int searchDim = blockSize_*numBlocks_;
933 
934  // Update obtained from solving saddle point problem
935  RCP<MV> Delta = MVT::Clone(*X_,blockSize_);
936 
938  // iterate until the status test tells us to stop.
939  // also break if our basis is full
940  while (tester_->checkStatus(this) != Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
941 
942  // Print information on current iteration
943  if (om_->isVerbosity(Debug)) {
944  currentStatus( om_->stream(Debug) );
945  }
946  else if (om_->isVerbosity(IterationDetails)) {
947  currentStatus( om_->stream(IterationDetails) );
948  }
949 
950  ++iter_;
951 
952  // Solve the saddle point problem
953  solveSaddlePointProblem(Delta);
954 
955  // Insert Delta at the end of V
956  addToBasis(Delta);
957 
958  if (om_->isVerbosity( Debug ) ) {
959  // Check almost everything here
960  CheckList chk;
961  chk.checkV = true;
962  om_->print( Debug, accuracyCheck(chk, ": after addToBasis(Delta)") );
963  }
964 
965  // Compute KK := V'KV
966  computeKK();
967 
968  if (om_->isVerbosity( Debug ) ) {
969  // Check almost everything here
970  CheckList chk;
971  chk.checkKK = true;
972  om_->print( Debug, accuracyCheck(chk, ": after computeKK()") );
973  }
974 
975  // Compute the Ritz pairs
976  computeRitzPairs();
977 
978  // Compute X := V RitzVecs
979  computeX();
980 
981  if (om_->isVerbosity( Debug ) ) {
982  // Check almost everything here
983  CheckList chk;
984  chk.checkX = true;
985  om_->print( Debug, accuracyCheck(chk, ": after computeX()") );
986  }
987 
988  // Compute KX := KV RitzVecs and MX := MV RitzVecs (if necessary)
989  updateKXMX();
990 
991  if (om_->isVerbosity( Debug ) ) {
992  // Check almost everything here
993  CheckList chk;
994  chk.checkKX = true;
995  chk.checkMX = true;
996  om_->print( Debug, accuracyCheck(chk, ": after updateKXMX()") );
997  }
998 
999  // Update the residual vectors
1000  updateResidual();
1001  } // end while (statusTest == false)
1002 
1003  } // end of iterate()
1004 
1005 
1006 
1008  // Perform TraceMinBase iterations until the StatusTest tells us to stop.
1009  template <class ScalarType, class MV, class OP>
1011  {
1012  //
1013  // Initialize solver state
1014  if (initialized_ == false) {
1015  initialize();
1016  }
1017 
1018  // as a data member, this would be redundant and require synchronization with
1019  // blockSize_ and numBlocks_; we'll use a constant here.
1020  const int searchDim = blockSize_*numBlocks_;
1021 
1022  // Update obtained from solving saddle point problem
1023  RCP<MV> Delta = MVT::Clone(*X_,blockSize_);
1024 
1026  // iterate until the status test tells us to stop.
1027  // also break if our basis is full
1028  while (tester_->checkStatus(this) != Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
1029 
1030  // Print information on current iteration
1031  if (om_->isVerbosity(Debug)) {
1032  currentStatus( om_->stream(Debug) );
1033  }
1034  else if (om_->isVerbosity(IterationDetails)) {
1035  currentStatus( om_->stream(IterationDetails) );
1036  }
1037 
1038  ++iter_;
1039 
1040  // Solve the saddle point problem
1041  solveSaddlePointProblem(Delta);
1042 
1043  // Insert Delta at the end of V
1044  harmonicAddToBasis(Delta);
1045 
1046  if (om_->isVerbosity( Debug ) ) {
1047  // Check almost everything here
1048  CheckList chk;
1049  chk.checkV = true;
1050  om_->print( Debug, accuracyCheck(chk, ": after addToBasis(Delta)") );
1051  }
1052 
1053  // Compute KK := V'KV
1054  computeKK();
1055 
1056  if (om_->isVerbosity( Debug ) ) {
1057  // Check almost everything here
1058  CheckList chk;
1059  chk.checkKK = true;
1060  om_->print( Debug, accuracyCheck(chk, ": after computeKK()") );
1061  }
1062 
1063  // Compute the Ritz pairs
1064  computeRitzPairs();
1065 
1066  // Compute X := V RitzVecs
1067  computeX();
1068 
1069  // Get norm of each vector in X
1070  int nvecs;
1071  if(computeAllRes_)
1072  nvecs = curDim_;
1073  else
1074  nvecs = blockSize_;
1075  std::vector<int> dimind(nvecs);
1076  for(int i=0; i<nvecs; i++)
1077  dimind[i] = i;
1078  RCP<MV> lclX = MVT::CloneViewNonConst(*X_,dimind);
1079  std::vector<ScalarType> normvec(nvecs);
1080  orthman_->normMat(*lclX,normvec);
1081 
1082  // Scale X
1083  for(int i=0; i<nvecs; i++)
1084  normvec[i] = ONE/normvec[i];
1085  MVT::MvScale(*lclX,normvec);
1086 
1087  // Scale eigenvalues
1088  for(int i=0; i<nvecs; i++)
1089  {
1090  theta_[i] = theta_[i] * normvec[i] * normvec[i];
1091  }
1092 
1093  if (om_->isVerbosity( Debug ) ) {
1094  // Check almost everything here
1095  CheckList chk;
1096  chk.checkX = true;
1097  om_->print( Debug, accuracyCheck(chk, ": after computeX()") );
1098  }
1099 
1100  // Compute KX := KV RitzVecs and MX := MV RitzVecs (if necessary)
1101  updateKXMX();
1102 
1103  // Scale KX and MX
1104  if(Op_ != Teuchos::null)
1105  {
1106  RCP<MV> lclKX = MVT::CloneViewNonConst(*KX_,dimind);
1107  MVT::MvScale(*lclKX,normvec);
1108  }
1109  if(hasM_)
1110  {
1111  RCP<MV> lclMX = MVT::CloneViewNonConst(*MX_,dimind);
1112  MVT::MvScale(*lclMX,normvec);
1113  }
1114 
1115  if (om_->isVerbosity( Debug ) ) {
1116  // Check almost everything here
1117  CheckList chk;
1118  chk.checkKX = true;
1119  chk.checkMX = true;
1120  om_->print( Debug, accuracyCheck(chk, ": after updateKXMX()") );
1121  }
1122 
1123  // Update the residual vectors
1124  updateResidual();
1125  } // end while (statusTest == false)
1126 
1127  } // end of harmonicIterate()
1128 
1129 
1131  // initialize the solver with default state
1132  template <class ScalarType, class MV, class OP>
1134  {
1136  initialize(empty);
1137  }
1138 
1139 
1141  /* Initialize the state of the solver
1142  *
1143  * POST-CONDITIONS:
1144  *
1145  * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
1146  * theta_ contains Ritz w.r.t. V_(1:curDim_)
1147  * X is Ritz vectors w.r.t. V_(1:curDim_)
1148  * KX = Op*X
1149  * MX = M*X if hasM_
1150  * R = KX - MX*diag(theta_)
1151  *
1152  */
1153  template <class ScalarType, class MV, class OP>
1155  {
1156  // TODO: Check for bad input
1157  // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
1158  // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
1159 
1160 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1161  Teuchos::TimeMonitor inittimer( *timerInit_ );
1162 #endif
1163 
1164  previouslyLeveled_ = false;
1165 
1166  if(useHarmonic_)
1167  {
1168  harmonicInitialize(newstate);
1169  return;
1170  }
1171 
1172  std::vector<int> bsind(blockSize_);
1173  for (int i=0; i<blockSize_; ++i) bsind[i] = i;
1174 
1175  // in TraceMinBase, V is primary
1176  // the order of dependence follows like so.
1177  // --init-> V,KK
1178  // --ritz analysis-> theta,X
1179  // --op apply-> KX,MX
1180  // --compute-> R
1181  //
1182  // if the user specifies all data for a level, we will accept it.
1183  // otherwise, we will generate the whole level, and all subsequent levels.
1184  //
1185  // the data members are ordered based on dependence, and the levels are
1186  // partitioned according to the amount of work required to produce the
1187  // items in a level.
1188  //
1189  // inconsistent multivectors widths and lengths will not be tolerated, and
1190  // will be treated with exceptions.
1191  //
1192  // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to
1193  // multivectors in the solver, the copy will not be affected.
1194 
1195  // set up V and KK: get them from newstate if user specified them
1196  // otherwise, set them manually
1197  RCP<MV> lclV, lclKV, lclMV;
1198  RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK, lclRV;
1199 
1200  // If V was supplied by the user...
1201  if (newstate.V != Teuchos::null) {
1202  om_->stream(Debug) << "Copying V from the user\n";
1203 
1204  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
1205  "Anasazi::TraceMinBase::initialize(newstate): Vector length of V not correct." );
1206  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
1207  "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1208  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
1209  "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
1210 
1211  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
1212  "Anasazi::TraceMinBase::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
1213 
1214  curDim_ = newstate.curDim;
1215  // pick an integral amount
1216  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1217 
1218  TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
1219  "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
1220 
1221  // put data in V
1222  std::vector<int> nevind(curDim_);
1223  for (int i=0; i<curDim_; ++i) nevind[i] = i;
1224  if (newstate.V != V_) {
1225  if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
1226  newstate.V = MVT::CloneView(*newstate.V,nevind);
1227  }
1228  MVT::SetBlock(*newstate.V,nevind,*V_);
1229  }
1230  lclV = MVT::CloneViewNonConst(*V_,nevind);
1231  } // end if user specified V
1232  else
1233  {
1234  // get vectors from problem or generate something, projectAndNormalize
1235  RCP<const MV> ivec = problem_->getInitVec();
1236  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
1237  "Anasazi::TraceMinBase::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
1238  // clear newstate so we won't use any data from it below
1239  newstate.X = Teuchos::null;
1240  newstate.MX = Teuchos::null;
1241  newstate.KX = Teuchos::null;
1242  newstate.R = Teuchos::null;
1243  newstate.T = Teuchos::null;
1244  newstate.RV = Teuchos::null;
1245  newstate.KK = Teuchos::null;
1246  newstate.KV = Teuchos::null;
1247  newstate.MopV = Teuchos::null;
1248  newstate.isOrtho = false;
1249 
1250  // If the user did not specify a current dimension, use that of the initial vectors they provided
1251  if(newstate.curDim > 0)
1252  curDim_ = newstate.curDim;
1253  else
1254  curDim_ = MVT::GetNumberVecs(*ivec);
1255 
1256  // pick the largest multiple of blockSize_
1257  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1258  if (curDim_ > blockSize_*numBlocks_) {
1259  // user specified too many vectors... truncate
1260  // this produces a full subspace, but that is okay
1261  curDim_ = blockSize_*numBlocks_;
1262  }
1263  bool userand = false;
1264  if (curDim_ == 0) {
1265  // we need at least blockSize_ vectors
1266  // use a random multivec: ignore everything from InitVec
1267  userand = true;
1268  curDim_ = blockSize_;
1269  }
1270 
1271  std::vector<int> nevind(curDim_);
1272  for (int i=0; i<curDim_; ++i) nevind[i] = i;
1273 
1274  // get a pointer into V
1275  // lclV has curDim vectors
1276  //
1277  // get pointer to first curDim vectors in V_
1278  lclV = MVT::CloneViewNonConst(*V_,nevind);
1279  if (userand)
1280  {
1281  // generate random vector data
1282  MVT::MvRandom(*lclV);
1283  }
1284  else
1285  {
1286  if(newstate.curDim > 0)
1287  {
1288  if(MVT::GetNumberVecs(*newstate.V) > curDim_) {
1289  RCP<const MV> helperMV = MVT::CloneView(*newstate.V,nevind);
1290  MVT::SetBlock(*helperMV,nevind,*lclV);
1291  }
1292  // assign ivec to first part of lclV
1293  MVT::SetBlock(*newstate.V,nevind,*lclV);
1294  }
1295  else
1296  {
1297  if(MVT::GetNumberVecs(*ivec) > curDim_) {
1298  ivec = MVT::CloneView(*ivec,nevind);
1299  }
1300  // assign ivec to first part of lclV
1301  MVT::SetBlock(*ivec,nevind,*lclV);
1302  }
1303  }
1304  } // end if user did not specify V
1305 
1306  // A vector of relevant indices
1307  std::vector<int> dimind(curDim_);
1308  for (int i=0; i<curDim_; ++i) dimind[i] = i;
1309 
1310  // Compute MV if necessary
1311  if(hasM_ && newstate.MopV == Teuchos::null)
1312  {
1313  om_->stream(Debug) << "Computing MV\n";
1314 
1315 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1316  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1317 #endif
1318  count_ApplyM_+= curDim_;
1319 
1320  newstate.isOrtho = false;
1321  // Get a pointer to the relevant parts of MV
1322  lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1323  OPT::Apply(*MOp_,*lclV,*lclMV);
1324  }
1325  // Copy MV if necessary
1326  else if(hasM_)
1327  {
1328  om_->stream(Debug) << "Copying MV\n";
1329 
1330  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MopV) != MVT::GetGlobalLength(*MV_), std::invalid_argument,
1331  "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1332 
1333  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MopV) < curDim_, std::invalid_argument,
1334  "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in MV not correct.");
1335 
1336  if(newstate.MopV != MV_) {
1337  if(curDim_ < MVT::GetNumberVecs(*newstate.MopV)) {
1338  newstate.MopV = MVT::CloneView(*newstate.MopV,dimind);
1339  }
1340  MVT::SetBlock(*newstate.MopV,dimind,*MV_);
1341  }
1342  lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1343  }
1344  // There is no M, so there is no MV
1345  else
1346  {
1347  om_->stream(Debug) << "There is no MV\n";
1348 
1349  lclMV = lclV;
1350  MV_ = V_;
1351  }
1352 
1353  // Project and normalize so that V'V = I and Q'V=0
1354  if(newstate.isOrtho == false)
1355  {
1356  om_->stream(Debug) << "Project and normalize\n";
1357 
1358 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1359  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1360 #endif
1361 
1362  // These things are now invalid
1363  newstate.X = Teuchos::null;
1364  newstate.MX = Teuchos::null;
1365  newstate.KX = Teuchos::null;
1366  newstate.R = Teuchos::null;
1367  newstate.T = Teuchos::null;
1368  newstate.RV = Teuchos::null;
1369  newstate.KK = Teuchos::null;
1370  newstate.KV = Teuchos::null;
1371 
1372  int rank;
1373  if(auxVecs_.size() > 0)
1374  {
1375  rank = orthman_->projectAndNormalizeMat(*lclV, auxVecs_,
1376  Teuchos::tuple(RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
1377  Teuchos::null, lclMV, MauxVecs_);
1378  }
1379  else
1380  {
1381  rank = orthman_->normalizeMat(*lclV,Teuchos::null,lclMV);
1382  }
1383 
1384  TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,TraceMinBaseInitFailure,
1385  "Anasazi::TraceMinBase::initialize(): Couldn't generate initial basis of full rank.");
1386  }
1387 
1388  // Compute KV
1389  if(Op_ != Teuchos::null && newstate.KV == Teuchos::null)
1390  {
1391  om_->stream(Debug) << "Computing KV\n";
1392 
1393 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1394  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1395 #endif
1396  count_ApplyOp_+= curDim_;
1397 
1398  // These things are now invalid
1399  newstate.X = Teuchos::null;
1400  newstate.MX = Teuchos::null;
1401  newstate.KX = Teuchos::null;
1402  newstate.R = Teuchos::null;
1403  newstate.T = Teuchos::null;
1404  newstate.RV = Teuchos::null;
1405  newstate.KK = Teuchos::null;
1406  newstate.KV = Teuchos::null;
1407 
1408  lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1409  OPT::Apply(*Op_,*lclV,*lclKV);
1410  }
1411  // Copy KV
1412  else if(Op_ != Teuchos::null)
1413  {
1414  om_->stream(Debug) << "Copying MV\n";
1415 
1416  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KV) != MVT::GetGlobalLength(*KV_), std::invalid_argument,
1417  "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1418 
1419  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KV) < curDim_, std::invalid_argument,
1420  "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in KV not correct.");
1421 
1422  if (newstate.KV != KV_) {
1423  if (curDim_ < MVT::GetNumberVecs(*newstate.KV)) {
1424  newstate.KV = MVT::CloneView(*newstate.KV,dimind);
1425  }
1426  MVT::SetBlock(*newstate.KV,dimind,*KV_);
1427  }
1428  lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1429  }
1430  else
1431  {
1432  om_->stream(Debug) << "There is no KV\n";
1433 
1434  lclKV = lclV;
1435  KV_ = V_;
1436  }
1437 
1438  // Compute KK
1439  if(newstate.KK == Teuchos::null)
1440  {
1441  om_->stream(Debug) << "Computing KK\n";
1442 
1443  // These things are now invalid
1444  newstate.X = Teuchos::null;
1445  newstate.MX = Teuchos::null;
1446  newstate.KX = Teuchos::null;
1447  newstate.R = Teuchos::null;
1448  newstate.T = Teuchos::null;
1449  newstate.RV = Teuchos::null;
1450 
1451  // Note: there used to be a bug here.
1452  // We can't just call computeKK because it only computes the new parts of KK
1453 
1454  // Get a pointer to the part of KK we're interested in
1455  lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1456 
1457  // KK := V'KV
1458  MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
1459  }
1460  // Copy KK
1461  else
1462  {
1463  om_->stream(Debug) << "Copying KK\n";
1464 
1465  // check size of KK
1466  TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
1467  "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
1468 
1469  // put data into KK_
1470  lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1471  if (newstate.KK != KK_) {
1472  if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
1473  newstate.KK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
1474  }
1475  lclKK->assign(*newstate.KK);
1476  }
1477  }
1478 
1479  // Compute Ritz pairs
1480  if(newstate.T == Teuchos::null || newstate.RV == Teuchos::null)
1481  {
1482  om_->stream(Debug) << "Computing Ritz pairs\n";
1483 
1484  // These things are now invalid
1485  newstate.X = Teuchos::null;
1486  newstate.MX = Teuchos::null;
1487  newstate.KX = Teuchos::null;
1488  newstate.R = Teuchos::null;
1489  newstate.T = Teuchos::null;
1490  newstate.RV = Teuchos::null;
1491 
1492  computeRitzPairs();
1493  }
1494  // Copy Ritz pairs
1495  else
1496  {
1497  om_->stream(Debug) << "Copying Ritz pairs\n";
1498 
1499  TEUCHOS_TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
1500  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
1501 
1502  TEUCHOS_TEST_FOR_EXCEPTION( newstate.RV->numRows() < curDim_ || newstate.RV->numCols() < curDim_, std::invalid_argument,
1503  "Anasazi::TraceMinBase::initialize(newstate): Ritz vectors in new state must be as large as specified state rank.");
1504 
1505  std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
1506 
1507  lclRV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
1508  if (newstate.RV != ritzVecs_) {
1509  if (newstate.RV->numRows() > curDim_ || newstate.RV->numCols() > curDim_) {
1510  newstate.RV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.RV,curDim_,curDim_) );
1511  }
1512  lclRV->assign(*newstate.RV);
1513  }
1514  }
1515 
1516  // Compute X
1517  if(newstate.X == Teuchos::null)
1518  {
1519  om_->stream(Debug) << "Computing X\n";
1520 
1521  // These things are now invalid
1522  newstate.MX = Teuchos::null;
1523  newstate.KX = Teuchos::null;
1524  newstate.R = Teuchos::null;
1525 
1526  computeX();
1527  }
1528  // Copy X
1529  else
1530  {
1531  om_->stream(Debug) << "Copying X\n";
1532 
1533  if(computeAllRes_ == false)
1534  {
1535  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1536  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with block size and length of V.");
1537 
1538  if(newstate.X != X_) {
1539  MVT::SetBlock(*newstate.X,bsind,*X_);
1540  }
1541  }
1542  else
1543  {
1544  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != curDim_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1545  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with current dimension and length of V.");
1546 
1547  if(newstate.X != X_) {
1548  MVT::SetBlock(*newstate.X,dimind,*X_);
1549  }
1550  }
1551  }
1552 
1553  // Compute KX and MX if necessary
1554  // TODO: These technically should be separate; it won't matter much in terms of running time though
1555  if((Op_ != Teuchos::null && newstate.KX == Teuchos::null) || (hasM_ && newstate.MX == Teuchos::null))
1556  {
1557  om_->stream(Debug) << "Computing KX and MX\n";
1558 
1559  // These things are now invalid
1560  newstate.R = Teuchos::null;
1561 
1562  updateKXMX();
1563  }
1564  // Copy KX and MX if necessary
1565  else
1566  {
1567  om_->stream(Debug) << "Copying KX and MX\n";
1568  if(Op_ != Teuchos::null)
1569  {
1570  if(computeAllRes_ == false)
1571  {
1572  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
1573  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with block size and length of V.");
1574 
1575  if(newstate.KX != KX_) {
1576  MVT::SetBlock(*newstate.KX,bsind,*KX_);
1577  }
1578  }
1579  else
1580  {
1581  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != curDim_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
1582  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with current dimension and length of V.");
1583 
1584  if (newstate.KX != KX_) {
1585  MVT::SetBlock(*newstate.KX,dimind,*KX_);
1586  }
1587  }
1588  }
1589 
1590  if(hasM_)
1591  {
1592  if(computeAllRes_ == false)
1593  {
1594  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
1595  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with block size and length of V.");
1596 
1597  if (newstate.MX != MX_) {
1598  MVT::SetBlock(*newstate.MX,bsind,*MX_);
1599  }
1600  }
1601  else
1602  {
1603  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != curDim_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
1604  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with current dimension and length of V.");
1605 
1606  if (newstate.MX != MX_) {
1607  MVT::SetBlock(*newstate.MX,dimind,*MX_);
1608  }
1609  }
1610  }
1611  }
1612 
1613  // Compute R
1614  if(newstate.R == Teuchos::null)
1615  {
1616  om_->stream(Debug) << "Computing R\n";
1617 
1618  updateResidual();
1619  }
1620  // Copy R
1621  else
1622  {
1623  om_->stream(Debug) << "Copying R\n";
1624 
1625  if(computeAllRes_ == false)
1626  {
1627  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1628  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
1629  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
1630  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least block size vectors." );
1631  if (newstate.R != R_) {
1632  MVT::SetBlock(*newstate.R,bsind,*R_);
1633  }
1634  }
1635  else
1636  {
1637  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1638  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
1639  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != curDim_,
1640  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least curDim vectors." );
1641  if (newstate.R != R_) {
1642  MVT::SetBlock(*newstate.R,dimind,*R_);
1643  }
1644  }
1645  }
1646 
1647  // R has been updated; mark the norms as out-of-date
1648  Rnorms_current_ = false;
1649  R2norms_current_ = false;
1650 
1651  // Set the largest safe shift
1652  largestSafeShift_ = newstate.largestSafeShift;
1653 
1654  // Copy over the Ritz shifts
1655  if(newstate.ritzShifts != Teuchos::null)
1656  {
1657  om_->stream(Debug) << "Copying Ritz shifts\n";
1658  std::copy(newstate.ritzShifts->begin(),newstate.ritzShifts->end(),ritzShifts_.begin());
1659  }
1660  else
1661  {
1662  om_->stream(Debug) << "Setting Ritz shifts to 0\n";
1663  for(size_t i=0; i<ritzShifts_.size(); i++)
1664  ritzShifts_[i] = ZERO;
1665  }
1666 
1667  for(size_t i=0; i<ritzShifts_.size(); i++)
1668  om_->stream(Debug) << "Ritz shifts[" << i << "] = " << ritzShifts_[i] << std::endl;
1669 
1670  // finally, we are initialized
1671  initialized_ = true;
1672 
1673  if (om_->isVerbosity( Debug ) ) {
1674  // Check almost everything here
1675  CheckList chk;
1676  chk.checkV = true;
1677  chk.checkX = true;
1678  chk.checkKX = true;
1679  chk.checkMX = true;
1680  chk.checkQ = true;
1681  chk.checkKK = true;
1682  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
1683  }
1684 
1685  // Print information on current status
1686  if (om_->isVerbosity(Debug)) {
1687  currentStatus( om_->stream(Debug) );
1688  }
1689  else if (om_->isVerbosity(IterationDetails)) {
1690  currentStatus( om_->stream(IterationDetails) );
1691  }
1692  }
1693 
1694 
1695 
1697  /* Initialize the state of the solver
1698  *
1699  * POST-CONDITIONS:
1700  *
1701  * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
1702  * theta_ contains Ritz w.r.t. V_(1:curDim_)
1703  * X is Ritz vectors w.r.t. V_(1:curDim_)
1704  * KX = Op*X
1705  * MX = M*X if hasM_
1706  * R = KX - MX*diag(theta_)
1707  *
1708  */
1709  template <class ScalarType, class MV, class OP>
1711  {
1712  // TODO: Check for bad input
1713  // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
1714  // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
1715 
1716  std::vector<int> bsind(blockSize_);
1717  for (int i=0; i<blockSize_; ++i) bsind[i] = i;
1718 
1719  // in TraceMinBase, V is primary
1720  // the order of dependence follows like so.
1721  // --init-> V,KK
1722  // --ritz analysis-> theta,X
1723  // --op apply-> KX,MX
1724  // --compute-> R
1725  //
1726  // if the user specifies all data for a level, we will accept it.
1727  // otherwise, we will generate the whole level, and all subsequent levels.
1728  //
1729  // the data members are ordered based on dependence, and the levels are
1730  // partitioned according to the amount of work required to produce the
1731  // items in a level.
1732  //
1733  // inconsistent multivectors widths and lengths will not be tolerated, and
1734  // will be treated with exceptions.
1735  //
1736  // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to
1737  // multivectors in the solver, the copy will not be affected.
1738 
1739  // set up V and KK: get them from newstate if user specified them
1740  // otherwise, set them manually
1741  RCP<MV> lclV, lclKV, lclMV;
1742  RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK, lclRV;
1743 
1744  // If V was supplied by the user...
1745  if (newstate.V != Teuchos::null) {
1746  om_->stream(Debug) << "Copying V from the user\n";
1747 
1748  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
1749  "Anasazi::TraceMinBase::initialize(newstate): Vector length of V not correct." );
1750  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
1751  "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1752  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
1753  "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
1754 
1755  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
1756  "Anasazi::TraceMinBase::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
1757 
1758  curDim_ = newstate.curDim;
1759  // pick an integral amount
1760  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1761 
1762  TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
1763  "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
1764 
1765  // put data in V
1766  std::vector<int> nevind(curDim_);
1767  for (int i=0; i<curDim_; ++i) nevind[i] = i;
1768  if (newstate.V != V_) {
1769  if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
1770  newstate.V = MVT::CloneView(*newstate.V,nevind);
1771  }
1772  MVT::SetBlock(*newstate.V,nevind,*V_);
1773  }
1774  lclV = MVT::CloneViewNonConst(*V_,nevind);
1775  } // end if user specified V
1776  else
1777  {
1778  // get vectors from problem or generate something, projectAndNormalize
1779  RCP<const MV> ivec = problem_->getInitVec();
1780  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
1781  "Anasazi::TraceMinBase::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
1782  // clear newstate so we won't use any data from it below
1783  newstate.X = Teuchos::null;
1784  newstate.MX = Teuchos::null;
1785  newstate.KX = Teuchos::null;
1786  newstate.R = Teuchos::null;
1787  newstate.T = Teuchos::null;
1788  newstate.RV = Teuchos::null;
1789  newstate.KK = Teuchos::null;
1790  newstate.KV = Teuchos::null;
1791  newstate.MopV = Teuchos::null;
1792  newstate.isOrtho = false;
1793 
1794  // If the user did not specify a current dimension, use that of the initial vectors they provided
1795  if(newstate.curDim > 0)
1796  curDim_ = newstate.curDim;
1797  else
1798  curDim_ = MVT::GetNumberVecs(*ivec);
1799 
1800  // pick the largest multiple of blockSize_
1801  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1802  if (curDim_ > blockSize_*numBlocks_) {
1803  // user specified too many vectors... truncate
1804  // this produces a full subspace, but that is okay
1805  curDim_ = blockSize_*numBlocks_;
1806  }
1807  bool userand = false;
1808  if (curDim_ == 0) {
1809  // we need at least blockSize_ vectors
1810  // use a random multivec: ignore everything from InitVec
1811  userand = true;
1812  curDim_ = blockSize_;
1813  }
1814 
1815  std::vector<int> nevind(curDim_);
1816  for (int i=0; i<curDim_; ++i) nevind[i] = i;
1817 
1818  // get a pointer into V
1819  // lclV has curDim vectors
1820  //
1821  // get pointer to first curDim vectors in V_
1822  lclV = MVT::CloneViewNonConst(*V_,nevind);
1823  if (userand)
1824  {
1825  // generate random vector data
1826  MVT::MvRandom(*lclV);
1827  }
1828  else
1829  {
1830  if(newstate.curDim > 0)
1831  {
1832  if(MVT::GetNumberVecs(*newstate.V) > curDim_) {
1833  RCP<const MV> helperMV = MVT::CloneView(*newstate.V,nevind);
1834  MVT::SetBlock(*helperMV,nevind,*lclV);
1835  }
1836  // assign ivec to first part of lclV
1837  MVT::SetBlock(*newstate.V,nevind,*lclV);
1838  }
1839  else
1840  {
1841  if(MVT::GetNumberVecs(*ivec) > curDim_) {
1842  ivec = MVT::CloneView(*ivec,nevind);
1843  }
1844  // assign ivec to first part of lclV
1845  MVT::SetBlock(*ivec,nevind,*lclV);
1846  }
1847  }
1848  } // end if user did not specify V
1849 
1850  // Nuke everything from orbit
1851  // This is a temporary measure due to a bug in the code that I have not found yet
1852  // It adds a minimal amount of work
1853  newstate.X = Teuchos::null;
1854  newstate.MX = Teuchos::null;
1855  newstate.KX = Teuchos::null;
1856  newstate.R = Teuchos::null;
1857  newstate.T = Teuchos::null;
1858  newstate.RV = Teuchos::null;
1859  newstate.KK = Teuchos::null;
1860  newstate.KV = Teuchos::null;
1861  newstate.MopV = Teuchos::null;
1862  newstate.isOrtho = false;
1863 
1864  // A vector of relevant indices
1865  std::vector<int> dimind(curDim_);
1866  for (int i=0; i<curDim_; ++i) dimind[i] = i;
1867 
1868  // Project the auxVecs out of V
1869  if(auxVecs_.size() > 0)
1870  orthman_->projectMat(*lclV,auxVecs_);
1871 
1872  // Compute KV
1873  if(Op_ != Teuchos::null && newstate.KV == Teuchos::null)
1874  {
1875  om_->stream(Debug) << "Computing KV\n";
1876 
1877 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1878  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1879 #endif
1880  count_ApplyOp_+= curDim_;
1881 
1882  // These things are now invalid
1883  newstate.X = Teuchos::null;
1884  newstate.MX = Teuchos::null;
1885  newstate.KX = Teuchos::null;
1886  newstate.R = Teuchos::null;
1887  newstate.T = Teuchos::null;
1888  newstate.RV = Teuchos::null;
1889  newstate.KK = Teuchos::null;
1890 
1891  lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1892  OPT::Apply(*Op_,*lclV,*lclKV);
1893  }
1894  // Copy KV
1895  else if(Op_ != Teuchos::null)
1896  {
1897  om_->stream(Debug) << "Copying KV\n";
1898 
1899  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KV) != MVT::GetGlobalLength(*KV_), std::invalid_argument,
1900  "Anasazi::TraceMinBase::initialize(newstate): Vector length of KV not correct." );
1901 
1902  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KV) < curDim_, std::invalid_argument,
1903  "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in KV not correct.");
1904 
1905  if (newstate.KV != KV_) {
1906  if (curDim_ < MVT::GetNumberVecs(*newstate.KV)) {
1907  newstate.KV = MVT::CloneView(*newstate.KV,dimind);
1908  }
1909  MVT::SetBlock(*newstate.KV,dimind,*KV_);
1910  }
1911  lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1912  }
1913  else
1914  {
1915  om_->stream(Debug) << "There is no KV\n";
1916 
1917  lclKV = lclV;
1918  KV_ = V_;
1919  }
1920 
1921 
1922 
1923  // Project and normalize so that V'V = I and Q'V=0
1924  if(newstate.isOrtho == false)
1925  {
1926  om_->stream(Debug) << "Project and normalize\n";
1927 
1928 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1929  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1930 #endif
1931 
1932  // These things are now invalid
1933  newstate.MopV = Teuchos::null;
1934  newstate.X = Teuchos::null;
1935  newstate.MX = Teuchos::null;
1936  newstate.KX = Teuchos::null;
1937  newstate.R = Teuchos::null;
1938  newstate.T = Teuchos::null;
1939  newstate.RV = Teuchos::null;
1940  newstate.KK = Teuchos::null;
1941 
1942  // Normalize lclKV
1943  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(curDim_,curDim_));
1944  int rank = orthman_->normalizeMat(*lclKV,gamma);
1945 
1946  // lclV = lclV/gamma
1947  Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
1948  SDsolver.setMatrix(gamma);
1949  SDsolver.invert();
1950  RCP<MV> tempMV = MVT::CloneCopy(*lclV);
1951  MVT::MvTimesMatAddMv(ONE,*tempMV,*gamma,ZERO,*lclV);
1952 
1953  TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,TraceMinBaseInitFailure,
1954  "Anasazi::TraceMinBase::initialize(): Couldn't generate initial basis of full rank.");
1955  }
1956 
1957  // Compute MV if necessary
1958  if(hasM_ && newstate.MopV == Teuchos::null)
1959  {
1960  om_->stream(Debug) << "Computing MV\n";
1961 
1962 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1963  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1964 #endif
1965  count_ApplyM_+= curDim_;
1966 
1967  // Get a pointer to the relevant parts of MV
1968  lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1969  OPT::Apply(*MOp_,*lclV,*lclMV);
1970  }
1971  // Copy MV if necessary
1972  else if(hasM_)
1973  {
1974  om_->stream(Debug) << "Copying MV\n";
1975 
1976  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MopV) != MVT::GetGlobalLength(*MV_), std::invalid_argument,
1977  "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1978 
1979  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MopV) < curDim_, std::invalid_argument,
1980  "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in MV not correct.");
1981 
1982  if(newstate.MopV != MV_) {
1983  if(curDim_ < MVT::GetNumberVecs(*newstate.MopV)) {
1984  newstate.MopV = MVT::CloneView(*newstate.MopV,dimind);
1985  }
1986  MVT::SetBlock(*newstate.MopV,dimind,*MV_);
1987  }
1988  lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1989  }
1990  // There is no M, so there is no MV
1991  else
1992  {
1993  om_->stream(Debug) << "There is no MV\n";
1994 
1995  lclMV = lclV;
1996  MV_ = V_;
1997  }
1998 
1999  // Compute KK
2000  if(newstate.KK == Teuchos::null)
2001  {
2002  om_->stream(Debug) << "Computing KK\n";
2003 
2004  // These things are now invalid
2005  newstate.X = Teuchos::null;
2006  newstate.MX = Teuchos::null;
2007  newstate.KX = Teuchos::null;
2008  newstate.R = Teuchos::null;
2009  newstate.T = Teuchos::null;
2010  newstate.RV = Teuchos::null;
2011 
2012  // Note: there used to be a bug here.
2013  // We can't just call computeKK because it only computes the new parts of KK
2014 
2015  // Get a pointer to the part of KK we're interested in
2016  lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
2017 
2018  // KK := V'KV
2019  MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
2020  }
2021  // Copy KK
2022  else
2023  {
2024  om_->stream(Debug) << "Copying KK\n";
2025 
2026  // check size of KK
2027  TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
2028  "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
2029 
2030  // put data into KK_
2031  lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
2032  if (newstate.KK != KK_) {
2033  if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
2034  newstate.KK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
2035  }
2036  lclKK->assign(*newstate.KK);
2037  }
2038  }
2039 
2040  // Compute Ritz pairs
2041  if(newstate.T == Teuchos::null || newstate.RV == Teuchos::null)
2042  {
2043  om_->stream(Debug) << "Computing Ritz pairs\n";
2044 
2045  // These things are now invalid
2046  newstate.X = Teuchos::null;
2047  newstate.MX = Teuchos::null;
2048  newstate.KX = Teuchos::null;
2049  newstate.R = Teuchos::null;
2050  newstate.T = Teuchos::null;
2051  newstate.RV = Teuchos::null;
2052 
2053  computeRitzPairs();
2054  }
2055  // Copy Ritz pairs
2056  else
2057  {
2058  om_->stream(Debug) << "Copying Ritz pairs\n";
2059 
2060  TEUCHOS_TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
2061  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
2062 
2063  TEUCHOS_TEST_FOR_EXCEPTION( newstate.RV->numRows() < curDim_ || newstate.RV->numCols() < curDim_, std::invalid_argument,
2064  "Anasazi::TraceMinBase::initialize(newstate): Ritz vectors in new state must be as large as specified state rank.");
2065 
2066  std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
2067 
2068  lclRV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
2069  if (newstate.RV != ritzVecs_) {
2070  if (newstate.RV->numRows() > curDim_ || newstate.RV->numCols() > curDim_) {
2071  newstate.RV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.RV,curDim_,curDim_) );
2072  }
2073  lclRV->assign(*newstate.RV);
2074  }
2075  }
2076 
2077  // Compute X
2078  if(newstate.X == Teuchos::null)
2079  {
2080  om_->stream(Debug) << "Computing X\n";
2081 
2082  // These things are now invalid
2083  newstate.MX = Teuchos::null;
2084  newstate.KX = Teuchos::null;
2085  newstate.R = Teuchos::null;
2086 
2087  computeX();
2088  }
2089  // Copy X
2090  else
2091  {
2092  om_->stream(Debug) << "Copying X\n";
2093 
2094  if(computeAllRes_ == false)
2095  {
2096  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
2097  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with block size and length of V.");
2098 
2099  if(newstate.X != X_) {
2100  MVT::SetBlock(*newstate.X,bsind,*X_);
2101  }
2102  }
2103  else
2104  {
2105  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != curDim_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
2106  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with current dimension and length of V.");
2107 
2108  if(newstate.X != X_) {
2109  MVT::SetBlock(*newstate.X,dimind,*X_);
2110  }
2111  }
2112  }
2113 
2114  // Compute KX and MX if necessary
2115  // TODO: These technically should be separate; it won't matter much in terms of running time though
2116  if((Op_ != Teuchos::null && newstate.KX == Teuchos::null) || (hasM_ && newstate.MX == Teuchos::null))
2117  {
2118  om_->stream(Debug) << "Computing KX and MX\n";
2119 
2120  // These things are now invalid
2121  newstate.R = Teuchos::null;
2122 
2123  updateKXMX();
2124  }
2125  // Copy KX and MX if necessary
2126  else
2127  {
2128  om_->stream(Debug) << "Copying KX and MX\n";
2129  if(Op_ != Teuchos::null)
2130  {
2131  if(computeAllRes_ == false)
2132  {
2133  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
2134  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with block size and length of V.");
2135 
2136  if(newstate.KX != KX_) {
2137  MVT::SetBlock(*newstate.KX,bsind,*KX_);
2138  }
2139  }
2140  else
2141  {
2142  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != curDim_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
2143  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with current dimension and length of V.");
2144 
2145  if (newstate.KX != KX_) {
2146  MVT::SetBlock(*newstate.KX,dimind,*KX_);
2147  }
2148  }
2149  }
2150 
2151  if(hasM_)
2152  {
2153  if(computeAllRes_ == false)
2154  {
2155  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
2156  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with block size and length of V.");
2157 
2158  if (newstate.MX != MX_) {
2159  MVT::SetBlock(*newstate.MX,bsind,*MX_);
2160  }
2161  }
2162  else
2163  {
2164  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != curDim_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
2165  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with current dimension and length of V.");
2166 
2167  if (newstate.MX != MX_) {
2168  MVT::SetBlock(*newstate.MX,dimind,*MX_);
2169  }
2170  }
2171  }
2172  }
2173 
2174  // Scale X so each vector is of length 1
2175  {
2176  // Get norm of each vector in X
2177  const int nvecs = computeAllRes_ ? curDim_ : blockSize_;
2178  Teuchos::Range1D dimind2 (0, nvecs-1);
2179  RCP<MV> lclX = MVT::CloneViewNonConst(*X_, dimind2);
2180  std::vector<ScalarType> normvec(nvecs);
2181  orthman_->normMat(*lclX,normvec);
2182 
2183  // Scale X, KX, and MX accordingly
2184  for (int i = 0; i < nvecs; ++i) {
2185  normvec[i] = ONE / normvec[i];
2186  }
2187  MVT::MvScale (*lclX, normvec);
2188  if (Op_ != Teuchos::null) {
2189  RCP<MV> lclKX = MVT::CloneViewNonConst (*KX_, dimind2);
2190  MVT::MvScale (*lclKX, normvec);
2191  }
2192  if (hasM_) {
2193  RCP<MV> lclMX = MVT::CloneViewNonConst (*MX_, dimind2);
2194  MVT::MvScale (*lclMX, normvec);
2195  }
2196 
2197  // Scale eigenvalues
2198  for (int i = 0; i < nvecs; ++i) {
2199  theta_[i] = theta_[i] * normvec[i] * normvec[i];
2200  }
2201  }
2202 
2203  // Compute R
2204  if(newstate.R == Teuchos::null)
2205  {
2206  om_->stream(Debug) << "Computing R\n";
2207 
2208  updateResidual();
2209  }
2210  // Copy R
2211  else
2212  {
2213  om_->stream(Debug) << "Copying R\n";
2214 
2215  if(computeAllRes_ == false)
2216  {
2217  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
2218  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
2219  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
2220  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least block size vectors." );
2221  if (newstate.R != R_) {
2222  MVT::SetBlock(*newstate.R,bsind,*R_);
2223  }
2224  }
2225  else
2226  {
2227  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
2228  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
2229  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != curDim_,
2230  std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least curDim vectors." );
2231  if (newstate.R != R_) {
2232  MVT::SetBlock(*newstate.R,dimind,*R_);
2233  }
2234  }
2235  }
2236 
2237  // R has been updated; mark the norms as out-of-date
2238  Rnorms_current_ = false;
2239  R2norms_current_ = false;
2240 
2241  // Set the largest safe shift
2242  largestSafeShift_ = newstate.largestSafeShift;
2243 
2244  // Copy over the Ritz shifts
2245  if(newstate.ritzShifts != Teuchos::null)
2246  {
2247  om_->stream(Debug) << "Copying Ritz shifts\n";
2248  std::copy(newstate.ritzShifts->begin(),newstate.ritzShifts->end(),ritzShifts_.begin());
2249  }
2250  else
2251  {
2252  om_->stream(Debug) << "Setting Ritz shifts to 0\n";
2253  for(size_t i=0; i<ritzShifts_.size(); i++)
2254  ritzShifts_[i] = ZERO;
2255  }
2256 
2257  for(size_t i=0; i<ritzShifts_.size(); i++)
2258  om_->stream(Debug) << "Ritz shifts[" << i << "] = " << ritzShifts_[i] << std::endl;
2259 
2260  // finally, we are initialized
2261  initialized_ = true;
2262 
2263  if (om_->isVerbosity( Debug ) ) {
2264  // Check almost everything here
2265  CheckList chk;
2266  chk.checkV = true;
2267  chk.checkX = true;
2268  chk.checkKX = true;
2269  chk.checkMX = true;
2270  chk.checkQ = true;
2271  chk.checkKK = true;
2272  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
2273  }
2274 
2275  // Print information on current status
2276  if (om_->isVerbosity(Debug)) {
2277  currentStatus( om_->stream(Debug) );
2278  }
2279  else if (om_->isVerbosity(IterationDetails)) {
2280  currentStatus( om_->stream(IterationDetails) );
2281  }
2282  }
2283 
2284 
2286  // Return initialized state
2287  template <class ScalarType, class MV, class OP>
2288  bool TraceMinBase<ScalarType,MV,OP>::isInitialized() const { return initialized_; }
2289 
2290 
2292  // Set the block size and make necessary adjustments.
2293  template <class ScalarType, class MV, class OP>
2294  void TraceMinBase<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
2295  {
2296  // This routine only allocates space; it doesn't not perform any computation
2297  // any change in size will invalidate the state of the solver.
2298 
2299  TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::TraceMinBase::setSize(blocksize,numblocks): blocksize must be strictly positive.");
2300 
2301  if (blockSize == blockSize_ && numBlocks == numBlocks_) {
2302  // do nothing
2303  return;
2304  }
2305 
2306  blockSize_ = blockSize;
2307  numBlocks_ = numBlocks;
2308 
2309  RCP<const MV> tmp;
2310  // grab some Multivector to Clone
2311  // in practice, getInitVec() should always provide this, but it is possible to use a
2312  // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
2313  // in case of that strange scenario, we will try to Clone from X_ first, then resort to getInitVec()
2314  if (X_ != Teuchos::null) { // this is equivalent to blockSize_ > 0
2315  tmp = X_;
2316  }
2317  else {
2318  tmp = problem_->getInitVec();
2319  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
2320  "Anasazi::TraceMinBase::setSize(): eigenproblem did not specify initial vectors to clone from.");
2321  }
2322 
2323  TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*static_cast<ptrdiff_t>(numBlocks) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
2324  "Anasazi::TraceMinBase::setSize(): max subspace dimension and auxilliary subspace too large. Potentially impossible orthogonality constraints.");
2325 
2326  // New subspace dimension
2327  int newsd = blockSize_*numBlocks_;
2328 
2330  // blockSize dependent
2331  //
2332  ritzShifts_.resize(blockSize_,ZERO);
2333  if(computeAllRes_ == false)
2334  {
2335  // grow/allocate vectors
2336  Rnorms_.resize(blockSize_,NANVAL);
2337  R2norms_.resize(blockSize_,NANVAL);
2338  //
2339  // clone multivectors off of tmp
2340  //
2341  // free current allocation first, to make room for new allocation
2342  X_ = Teuchos::null;
2343  KX_ = Teuchos::null;
2344  MX_ = Teuchos::null;
2345  R_ = Teuchos::null;
2346  V_ = Teuchos::null;
2347  KV_ = Teuchos::null;
2348  MV_ = Teuchos::null;
2349 
2350  om_->print(Debug," >> Allocating X_\n");
2351  X_ = MVT::Clone(*tmp,blockSize_);
2352  if(Op_ != Teuchos::null) {
2353  om_->print(Debug," >> Allocating KX_\n");
2354  KX_ = MVT::Clone(*tmp,blockSize_);
2355  }
2356  else {
2357  KX_ = X_;
2358  }
2359  if (hasM_) {
2360  om_->print(Debug," >> Allocating MX_\n");
2361  MX_ = MVT::Clone(*tmp,blockSize_);
2362  }
2363  else {
2364  MX_ = X_;
2365  }
2366  om_->print(Debug," >> Allocating R_\n");
2367  R_ = MVT::Clone(*tmp,blockSize_);
2368  }
2369  else
2370  {
2371  // grow/allocate vectors
2372  Rnorms_.resize(newsd,NANVAL);
2373  R2norms_.resize(newsd,NANVAL);
2374  //
2375  // clone multivectors off of tmp
2376  //
2377  // free current allocation first, to make room for new allocation
2378  X_ = Teuchos::null;
2379  KX_ = Teuchos::null;
2380  MX_ = Teuchos::null;
2381  R_ = Teuchos::null;
2382  V_ = Teuchos::null;
2383  KV_ = Teuchos::null;
2384  MV_ = Teuchos::null;
2385 
2386  om_->print(Debug," >> Allocating X_\n");
2387  X_ = MVT::Clone(*tmp,newsd);
2388  if(Op_ != Teuchos::null) {
2389  om_->print(Debug," >> Allocating KX_\n");
2390  KX_ = MVT::Clone(*tmp,newsd);
2391  }
2392  else {
2393  KX_ = X_;
2394  }
2395  if (hasM_) {
2396  om_->print(Debug," >> Allocating MX_\n");
2397  MX_ = MVT::Clone(*tmp,newsd);
2398  }
2399  else {
2400  MX_ = X_;
2401  }
2402  om_->print(Debug," >> Allocating R_\n");
2403  R_ = MVT::Clone(*tmp,newsd);
2404  }
2405 
2407  // blockSize*numBlocks dependent
2408  //
2409  theta_.resize(newsd,NANVAL);
2410  om_->print(Debug," >> Allocating V_\n");
2411  V_ = MVT::Clone(*tmp,newsd);
2412  KK_ = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
2413  ritzVecs_ = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
2414 
2415  if(Op_ != Teuchos::null) {
2416  om_->print(Debug," >> Allocating KV_\n");
2417  KV_ = MVT::Clone(*tmp,newsd);
2418  }
2419  else {
2420  KV_ = V_;
2421  }
2422  if (hasM_) {
2423  om_->print(Debug," >> Allocating MV_\n");
2424  MV_ = MVT::Clone(*tmp,newsd);
2425  }
2426  else {
2427  MV_ = V_;
2428  }
2429 
2430  om_->print(Debug," >> done allocating.\n");
2431 
2432  initialized_ = false;
2433  curDim_ = 0;
2434  }
2435 
2436 
2438  // Set the auxiliary vectors
2439  template <class ScalarType, class MV, class OP>
2440  void TraceMinBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<RCP<const MV> > &auxvecs) {
2441  typedef typename Teuchos::Array<RCP<const MV> >::iterator tarcpmv;
2442 
2443  // set new auxiliary vectors
2444  auxVecs_ = auxvecs;
2445 
2446  if(hasM_)
2447  MauxVecs_.resize(0);
2448  numAuxVecs_ = 0;
2449 
2450  for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
2451  numAuxVecs_ += MVT::GetNumberVecs(**i);
2452 
2453  if(hasM_)
2454  {
2455 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
2456  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
2457 #endif
2458  count_ApplyM_+= MVT::GetNumberVecs(**i);
2459 
2460  RCP<MV> helperMV = MVT::Clone(**i,MVT::GetNumberVecs(**i));
2461  OPT::Apply(*MOp_,**i,*helperMV);
2462  MauxVecs_.push_back(helperMV);
2463  }
2464  }
2465 
2466  // If the solver has been initialized, V is not necessarily orthogonal to new auxiliary vectors
2467  if (numAuxVecs_ > 0 && initialized_) {
2468  initialized_ = false;
2469  }
2470 
2471  if (om_->isVerbosity( Debug ) ) {
2472  // Check almost everything here
2473  CheckList chk;
2474  chk.checkQ = true;
2475  om_->print( Debug, accuracyCheck(chk, ": after setAuxVecs()") );
2476  }
2477  }
2478 
2479 
2481  // compute/return residual M-norms
2482  template <class ScalarType, class MV, class OP>
2483  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2485  if (Rnorms_current_ == false) {
2486  // Update the residual norms
2487  if(computeAllRes_)
2488  {
2489  std::vector<int> curind(curDim_);
2490  for(int i=0; i<curDim_; i++)
2491  curind[i] = i;
2492 
2493  RCP<const MV> locR = MVT::CloneView(*R_,curind);
2494  std::vector<ScalarType> locNorms(curDim_);
2495  orthman_->norm(*locR,locNorms);
2496 
2497  for(int i=0; i<curDim_; i++)
2498  Rnorms_[i] = locNorms[i];
2499  for(int i=curDim_+1; i<blockSize_*numBlocks_; i++)
2500  Rnorms_[i] = NANVAL;
2501 
2502  Rnorms_current_ = true;
2503  locNorms.resize(blockSize_);
2504  return locNorms;
2505  }
2506  else
2507  orthman_->norm(*R_,Rnorms_);
2508  Rnorms_current_ = true;
2509  }
2510  else if(computeAllRes_)
2511  {
2512  std::vector<ScalarType> locNorms(blockSize_);
2513  for(int i=0; i<blockSize_; i++)
2514  locNorms[i] = Rnorms_[i];
2515  return locNorms;
2516  }
2517 
2518  return Rnorms_;
2519  }
2520 
2521 
2523  // compute/return residual 2-norms
2524  template <class ScalarType, class MV, class OP>
2525  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2527  if (R2norms_current_ == false) {
2528  // Update the residual 2-norms
2529  if(computeAllRes_)
2530  {
2531  std::vector<int> curind(curDim_);
2532  for(int i=0; i<curDim_; i++)
2533  curind[i] = i;
2534 
2535  RCP<const MV> locR = MVT::CloneView(*R_,curind);
2536  std::vector<ScalarType> locNorms(curDim_);
2537  MVT::MvNorm(*locR,locNorms);
2538 
2539  for(int i=0; i<curDim_; i++)
2540  {
2541  R2norms_[i] = locNorms[i];
2542  }
2543  for(int i=curDim_+1; i<blockSize_*numBlocks_; i++)
2544  R2norms_[i] = NANVAL;
2545 
2546  R2norms_current_ = true;
2547  locNorms.resize(blockSize_);
2548  return locNorms;
2549  }
2550  else
2551  MVT::MvNorm(*R_,R2norms_);
2552  R2norms_current_ = true;
2553  }
2554  else if(computeAllRes_)
2555  {
2556  std::vector<ScalarType> locNorms(blockSize_);
2557  for(int i=0; i<blockSize_; i++)
2558  locNorms[i] = R2norms_[i];
2559  return locNorms;
2560  }
2561 
2562  return R2norms_;
2563  }
2564 
2566  // Set a new StatusTest for the solver.
2567  template <class ScalarType, class MV, class OP>
2569  TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
2570  "Anasazi::TraceMinBase::setStatusTest() was passed a null StatusTest.");
2571  tester_ = test;
2572  }
2573 
2575  // Get the current StatusTest used by the solver.
2576  template <class ScalarType, class MV, class OP>
2577  RCP<StatusTest<ScalarType,MV,OP> > TraceMinBase<ScalarType,MV,OP>::getStatusTest() const {
2578  return tester_;
2579  }
2580 
2582  // Print the current status of the solver
2583  template <class ScalarType, class MV, class OP>
2584  void
2586  {
2587  using std::endl;
2588 
2589  os.setf(std::ios::scientific, std::ios::floatfield);
2590  os.precision(6);
2591  os <<endl;
2592  os <<"================================================================================" << endl;
2593  os << endl;
2594  os <<" TraceMinBase Solver Status" << endl;
2595  os << endl;
2596  os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
2597  os <<"The number of iterations performed is " <<iter_<<endl;
2598  os <<"The block size is " << blockSize_<<endl;
2599  os <<"The number of blocks is " << numBlocks_<<endl;
2600  os <<"The current basis size is " << curDim_<<endl;
2601  os <<"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
2602  os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
2603  os <<"The number of operations M*x is "<<count_ApplyM_<<endl;
2604 
2605  os.setf(std::ios_base::right, std::ios_base::adjustfield);
2606 
2607  if (initialized_) {
2608  os << endl;
2609  os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
2610  os << std::setw(20) << "Eigenvalue"
2611  << std::setw(20) << "Residual(M)"
2612  << std::setw(20) << "Residual(2)"
2613  << endl;
2614  os <<"--------------------------------------------------------------------------------"<<endl;
2615  for (int i=0; i<blockSize_; ++i) {
2616  os << std::setw(20) << theta_[i];
2617  if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
2618  else os << std::setw(20) << "not current";
2619  if (R2norms_current_) os << std::setw(20) << R2norms_[i];
2620  else os << std::setw(20) << "not current";
2621  os << endl;
2622  }
2623  }
2624  os <<"================================================================================" << endl;
2625  os << endl;
2626  }
2627 
2629  template <class ScalarType, class MV, class OP>
2630  ScalarType TraceMinBase<ScalarType,MV,OP>::getTrace() const
2631  {
2632  ScalarType currentTrace = ZERO;
2633 
2634  for(int i=0; i < blockSize_; i++)
2635  currentTrace += theta_[i];
2636 
2637  return currentTrace;
2638  }
2639 
2641  template <class ScalarType, class MV, class OP>
2643  {
2644  ScalarType ratioOfChange = traceThresh_;
2645 
2646  if(previouslyLeveled_)
2647  {
2648  om_->stream(Debug) << "The trace already leveled, so we're not going to check it again\n";
2649  return true;
2650  }
2651 
2652  ScalarType currentTrace = getTrace();
2653 
2654  om_->stream(Debug) << "The current trace is " << currentTrace << std::endl;
2655 
2656  // Compute the ratio of the change
2657  // We seek the point where the trace has leveled off
2658  // It should be reasonably safe to shift at this point
2659  if(previousTrace_ != ZERO)
2660  {
2661  om_->stream(Debug) << "The previous trace was " << previousTrace_ << std::endl;
2662 
2663  ratioOfChange = std::abs(previousTrace_-currentTrace)/std::abs(previousTrace_);
2664  om_->stream(Debug) << "The ratio of change is " << ratioOfChange << std::endl;
2665  }
2666 
2667  previousTrace_ = currentTrace;
2668 
2669  if(ratioOfChange < traceThresh_)
2670  {
2671  previouslyLeveled_ = true;
2672  return true;
2673  }
2674  return false;
2675  }
2676 
2678  // Compute the residual of each CLUSTER of eigenvalues
2679  // This is important for selecting the Ritz shifts
2680  template <class ScalarType, class MV, class OP>
2681  std::vector<ScalarType> TraceMinBase<ScalarType,MV,OP>::getClusterResids()
2682  {
2683  int nvecs;
2684 
2685  if(computeAllRes_)
2686  nvecs = curDim_;
2687  else
2688  nvecs = blockSize_;
2689 
2690  getRes2Norms();
2691 
2692  std::vector<ScalarType> clusterResids(nvecs);
2693  std::vector<int> clusterIndices;
2694  for(int i=0; i < nvecs; i++)
2695  {
2696  // test for cluster
2697  if(clusterIndices.empty() || (theta_[i-1] + R2norms_[i-1] >= theta_[i] - R2norms_[i]))
2698  {
2699  // Add to cluster
2700  if(!clusterIndices.empty()) om_->stream(Debug) << theta_[i-1] << " is in a cluster with " << theta_[i] << " because " << theta_[i-1] + R2norms_[i-1] << " >= " << theta_[i] - R2norms_[i] << std::endl;
2701  clusterIndices.push_back(i);
2702  }
2703  // Cluster completed
2704  else
2705  {
2706  om_->stream(Debug) << theta_[i-1] << " is NOT in a cluster with " << theta_[i] << " because " << theta_[i-1] + R2norms_[i-1] << " < " << theta_[i] - R2norms_[i] << std::endl;
2707  ScalarType totalRes = ZERO;
2708  for(size_t j=0; j < clusterIndices.size(); j++)
2709  totalRes += (R2norms_[clusterIndices[j]]*R2norms_[clusterIndices[j]]);
2710 
2711  // If the smallest magnitude value of this sign is in a cluster with the
2712  // largest magnitude cluster of this sign, it is not safe for the smallest
2713  // eigenvalue to use a shift
2714  if(theta_[clusterIndices[0]] < 0 && theta_[i] < 0)
2715  negSafeToShift_ = true;
2716  else if(theta_[clusterIndices[0]] > 0 && theta_[i] > 0)
2717  posSafeToShift_ = true;
2718 
2719  for(size_t j=0; j < clusterIndices.size(); j++)
2720  clusterResids[clusterIndices[j]] = sqrt(totalRes);
2721 
2722  clusterIndices.clear();
2723  clusterIndices.push_back(i);
2724  }
2725  }
2726 
2727  // Handle last cluster
2728  ScalarType totalRes = ZERO;
2729  for(size_t j=0; j < clusterIndices.size(); j++)
2730  totalRes += R2norms_[clusterIndices[j]];
2731  for(size_t j=0; j < clusterIndices.size(); j++)
2732  clusterResids[clusterIndices[j]] = totalRes;
2733 
2734  return clusterResids;
2735  }
2736 
2738  // Compute the Ritz shifts based on the Ritz values and residuals
2739  // A note on shifting: if the matrix is indefinite, you NEED to use a large block size
2740  // TODO: resids[i] on its own is unsafe for the generalized EVP
2741  // See "A Parallel Implementation of the Trace Minimization Eigensolver"
2742  // by Eloy Romero and Jose E. Roman
2743  template <class ScalarType, class MV, class OP>
2744  void TraceMinBase<ScalarType,MV,OP>::computeRitzShifts(const std::vector<ScalarType>& clusterResids)
2745  {
2746  std::vector<ScalarType> thetaMag(theta_);
2747  bool traceHasLeveled = traceLeveled();
2748 
2749  // Compute the magnitude of the eigenvalues
2750  for(int i=0; i<blockSize_; i++)
2751  thetaMag[i] = std::abs(thetaMag[i]);
2752 
2753  // Test whether it is safe to shift
2754  // TODO: Add residual shift option
2755  // Note: If you choose single shift with an indefinite matrix, you're gonna have a bad time...
2756  // Note: This is not safe for indefinite matrices, and I don't even know that it CAN be made safe.
2757  // There just isn't any theoretical reason it should work.
2758  // TODO: It feels like this should be different for TraceMin-Base; we should be able to use all eigenvalues in the current subspace to determine whether we have a cluster
2759  if(whenToShift_ == ALWAYS_SHIFT || (whenToShift_ == SHIFT_WHEN_TRACE_LEVELS && traceHasLeveled))
2760  {
2761  // Set the shift to the largest safe shift
2762  if(howToShift_ == LARGEST_CONVERGED_SHIFT)
2763  {
2764  for(int i=0; i<blockSize_; i++)
2765  ritzShifts_[i] = largestSafeShift_;
2766  }
2767  // Set the shifts to the Ritz values
2768  else if(howToShift_ == RITZ_VALUES_SHIFT)
2769  {
2770  ritzShifts_[0] = theta_[0];
2771 
2772  // If we're using mulitple shifts, set them to EACH Ritz value.
2773  // Otherwise, only use the smallest
2774  if(useMultipleShifts_)
2775  {
2776  for(int i=1; i<blockSize_; i++)
2777  ritzShifts_[i] = theta_[i];
2778  }
2779  else
2780  {
2781  for(int i=1; i<blockSize_; i++)
2782  ritzShifts_[i] = ritzShifts_[0];
2783  }
2784  }
2785  // Use Dr. Sameh's original shifting strategy
2786  else if(howToShift_ == ADJUSTED_RITZ_SHIFT)
2787  {
2788  om_->stream(Debug) << "\nSeeking a shift for theta[0]=" << thetaMag[0] << std::endl;
2789 
2790  // This is my adjustment. If all eigenvalues are in a single cluster, it's probably a bad idea to shift the smallest one.
2791  // If all of your eigenvalues are in one cluster, it's either way to early to shift or your subspace is too small
2792  if((theta_[0] > 0 && posSafeToShift_) || (theta_[0] < 0 && negSafeToShift_) || considerClusters_ == false)
2793  {
2794  // Initialize with a conservative shift, either the biggest safe shift or the eigenvalue adjusted by its cluster's residual
2795  ritzShifts_[0] = std::max(largestSafeShift_,thetaMag[0]-clusterResids[0]);
2796 
2797  om_->stream(Debug) << "Initializing with a conservative shift, either the most positive converged eigenvalue ("
2798  << largestSafeShift_ << ") or the eigenvalue adjusted by the residual (" << thetaMag[0] << "-"
2799  << clusterResids[0] << ").\n";
2800 
2801  // If this eigenvalue is NOT in a cluster, do an aggressive shift
2802  if(R2norms_[0] == clusterResids[0])
2803  {
2804  ritzShifts_[0] = thetaMag[0];
2805  om_->stream(Debug) << "Since this eigenvalue is NOT in a cluster, we can use the eigenvalue itself as a shift: ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2806  }
2807  else
2808  om_->stream(Debug) << "This eigenvalue is in a cluster, so it would not be safe to use the eigenvalue itself as a shift\n";
2809  }
2810  else
2811  {
2812  if(largestSafeShift_ > std::abs(ritzShifts_[0]))
2813  {
2814  om_->stream(Debug) << "Initializing with a conservative shift...the most positive converged eigenvalue: " << largestSafeShift_ << std::endl;
2815  ritzShifts_[0] = largestSafeShift_;
2816  }
2817  else
2818  om_->stream(Debug) << "Using the previous value of ritzShifts[0]=" << ritzShifts_[0];
2819 
2820  }
2821 
2822  om_->stream(Debug) << "ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2823 
2824  if(useMultipleShifts_)
2825  {
2827  // Compute shifts for other eigenvalues
2828  for(int i=1; i < blockSize_; i++)
2829  {
2830  om_->stream(Debug) << "\nSeeking a shift for theta[" << i << "]=" << thetaMag[i] << std::endl;
2831 
2832  // If the previous shift was aggressive and we are not in a cluster, do an aggressive shift
2833  if(ritzShifts_[i-1] == thetaMag[i-1] && i < blockSize_-1 && thetaMag[i] < thetaMag[i+1] - clusterResids[i+1])
2834  {
2835  ritzShifts_[i] = thetaMag[i];
2836  om_->stream(Debug) << "Using an aggressive shift: ritzShifts_[" << i << "]=" << ritzShifts_[i] << std::endl;
2837  }
2838  else
2839  {
2840  if(ritzShifts_[0] > std::abs(ritzShifts_[i]))
2841  {
2842  om_->stream(Debug) << "It was unsafe to use the aggressive shift. Choose the shift used by theta[0]="
2843  << thetaMag[0] << ": ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2844 
2845  // Choose a conservative shift, that of the smallest positive eigenvalue
2846  ritzShifts_[i] = ritzShifts_[0];
2847  }
2848  else
2849  om_->stream(Debug) << "It was unsafe to use the aggressive shift. We will use the shift from the previous iteration: " << ritzShifts_[i] << std::endl;
2850 
2851  om_->stream(Debug) << "Check whether any less conservative shifts would work (such as the biggest eigenvalue outside of the cluster, namely theta[ell] < "
2852  << thetaMag[i] << "-" << clusterResids[i] << " (" << thetaMag[i] - clusterResids[i] << ")\n";
2853 
2854  // If possible, choose a less conservative shift, that of the biggest eigenvalue outside of the cluster
2855  for(int ell=0; ell < i; ell++)
2856  {
2857  if(thetaMag[ell] < thetaMag[i] - clusterResids[i])
2858  {
2859  ritzShifts_[i] = thetaMag[ell];
2860  om_->stream(Debug) << "ritzShifts_[" << i << "]=" << ritzShifts_[i] << " is valid\n";
2861  }
2862  else
2863  break;
2864  }
2865  } // end else
2866 
2867  om_->stream(Debug) << "ritzShifts[" << i << "]=" << ritzShifts_[i] << std::endl;
2868  } // end for
2869  } // end if(useMultipleShifts_)
2870  else
2871  {
2872  for(int i=1; i<blockSize_; i++)
2873  ritzShifts_[i] = ritzShifts_[0];
2874  }
2875  } // end if(howToShift_ == "Adjusted Ritz Values")
2876  } // end if(whenToShift_ == "Always" || (whenToShift_ == "After Trace Levels" && traceHasLeveled))
2877 
2878  // Set the correct sign
2879  for(int i=0; i<blockSize_; i++)
2880  {
2881  if(theta_[i] < 0)
2882  ritzShifts_[i] = -abs(ritzShifts_[i]);
2883  else
2884  ritzShifts_[i] = abs(ritzShifts_[i]);
2885  }
2886  }
2887 
2889  template <class ScalarType, class MV, class OP>
2890  std::vector<ScalarType> TraceMinBase<ScalarType,MV,OP>::computeTol()
2891  {
2892  ScalarType temp1, temp2;
2893  int nvecs = ritzShifts_.size();
2894  std::vector<ScalarType> tolerances;
2895 
2896  for(int i=0; i < nvecs; i++)
2897  {
2898  if(std::abs(theta_[0]) != std::abs(ritzShifts_[i]))
2899  temp1 = std::abs(theta_[i]-ritzShifts_[i])/std::abs(std::abs(theta_[0])-std::abs(ritzShifts_[i]));
2900  else
2901  temp1 = ZERO;
2902 
2903  temp2 = pow(2.0,-iter_);
2904 
2905  // TODO: The min and max tolerances should not be hard coded
2906  // Neither should the maximum number of iterations
2907  tolerances.push_back(std::max(std::min(temp1*temp1,temp2),1e-8));
2908  }
2909 
2910  if(nvecs > 1)
2911  tolerances[nvecs-1] = tolerances[nvecs-2];
2912 
2913  return tolerances;
2914  }
2915 
2917  template <class ScalarType, class MV, class OP>
2919  {
2920 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
2921  Teuchos::TimeMonitor lcltimer( *timerSaddle_ );
2922 #endif
2923 
2924  // This case can arise when looking for the largest eigenpairs
2925  if(Op_ == Teuchos::null)
2926  {
2927  // dense solver
2928  Teuchos::SerialDenseSolver<int,ScalarType> My_Solver;
2929 
2930  // Schur complement
2931  RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclS = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(blockSize_,blockSize_) );
2932 
2933  if(computeAllRes_)
2934  {
2935  // Get the valid indices of X
2936  std::vector<int> curind(blockSize_);
2937  for(int i=0; i<blockSize_; i++)
2938  curind[i] = i;
2939 
2940  // Get a view of MX
2941  RCP<const MV> lclMX = MVT::CloneView(*MX_, curind);
2942 
2943  // form S = X' M^2 X
2944  MVT::MvTransMv(ONE,*lclMX,*lclMX,*lclS);
2945 
2946  // compute the inverse of S
2947  My_Solver.setMatrix(lclS);
2948  My_Solver.invert();
2949 
2950  // Delta = X - MX/S
2951  RCP<const MV> lclX = MVT::CloneView(*X_, curind);
2952  MVT::Assign(*lclX,*Delta);
2953  MVT::MvTimesMatAddMv( -ONE, *lclMX, *lclS, ONE, *Delta);
2954  }
2955  else
2956  {
2957  // form S = X' M^2 X
2958  MVT::MvTransMv(ONE,*MX_,*MX_,*lclS);
2959 
2960  // compute the inverse of S
2961  My_Solver.setMatrix(lclS);
2962  My_Solver.invert();
2963 
2964  // Delta = X - MX/S
2965  MVT::Assign(*X_,*Delta);
2966  MVT::MvTimesMatAddMv( -ONE, *MX_, *lclS, ONE, *Delta);
2967  }
2968  }
2969  else
2970  {
2971  std::vector<int> order(curDim_);
2972  std::vector<ScalarType> tempvec(blockSize_);
2973 // RCP<BasicSort<MagnitudeType> > sorter = rcp( new BasicSort<MagnitudeType>("SR") );
2974 
2975  // Stores the residual of each CLUSTER of eigenvalues
2976  std::vector<ScalarType> clusterResids;
2977 
2978 /* // Sort the eigenvalues in ascending order for the Ritz shift selection
2979  sorter->sort(theta_, Teuchos::rcpFromRef(order), curDim_); // don't catch exception
2980 
2981  // Apply the same ordering to the residual norms
2982  getRes2Norms();
2983  for (int i=0; i<blockSize_; i++)
2984  tempvec[i] = R2norms_[order[i]];
2985  R2norms_ = tempvec;*/
2986 
2987  // Compute the residual of each CLUSTER of eigenvalues
2988  // This is important for selecting the Ritz shifts
2989  clusterResids = getClusterResids();
2990 
2991 /* // Sort the eigenvalues based on what the user wanted
2992  sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
2993 
2994  // Apply the same ordering to the residual norms and cluster residuals
2995  for (int i=0; i<blockSize_; i++)
2996  tempvec[i] = R2norms_[order[i]];
2997  R2norms_ = tempvec;
2998 
2999  for (int i=0; i<blockSize_; i++)
3000  tempvec[i] = clusterResids[order[i]];
3001  clusterResids = tempvec;*/
3002 
3003  // Compute the Ritz shifts
3004  computeRitzShifts(clusterResids);
3005 
3006  // Compute the tolerances for the inner solves
3007  std::vector<ScalarType> tolerances = computeTol();
3008 
3009  for(int i=0; i<blockSize_; i++)
3010  {
3011  om_->stream(IterationDetails) << "Choosing Ritz shifts...theta[" << i << "]="
3012  << theta_[i] << ", resids[" << i << "]=" << R2norms_[i] << ", clusterResids[" << i << "]=" << clusterResids[i]
3013  << ", ritzShifts[" << i << "]=" << ritzShifts_[i] << ", and tol[" << i << "]=" << tolerances[i] << std::endl;
3014  }
3015 
3016  // Set the Ritz shifts for the solver
3017  ritzOp_->setRitzShifts(ritzShifts_);
3018 
3019  // Set the inner stopping tolerance
3020  // This uses the Ritz values to determine when to stop
3021  ritzOp_->setInnerTol(tolerances);
3022 
3023  // Solve the saddle point problem
3024  if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER)
3025  {
3026  if(Prec_ != Teuchos::null)
3027  solveSaddleProjPrec(Delta);
3028  else
3029  solveSaddleProj(Delta);
3030  }
3031  else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER)
3032  {
3033  if(Z_ == Teuchos::null || MVT::GetNumberVecs(*Z_) != blockSize_)
3034  {
3035  // We do NOT want Z to be 0, because that could result in stagnation
3036  // I know it's tempting to take out the MvRandom, but seriously, don't do it.
3037  Z_ = MVT::Clone(*X_,blockSize_);
3038  MVT::MvRandom(*Z_);
3039  }
3040  solveSaddleSchur(Delta);
3041  }
3042  else if(saddleSolType_ == BD_PREC_MINRES)
3043  {
3044  solveSaddleBDPrec(Delta);
3045 // Delta->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
3046  }
3047  else if(saddleSolType_ == HSS_PREC_GMRES)
3048  {
3049  solveSaddleHSSPrec(Delta);
3050  }
3051  else
3052  std::cout << "Invalid saddle solver type\n";
3053  }
3054  }
3055 
3057  // Solve the saddle point problem using projected minres
3058  // TODO: We should be able to choose KX or -R as RHS.
3059  template <class ScalarType, class MV, class OP>
3060  void TraceMinBase<ScalarType,MV,OP>::solveSaddleProj (RCP<MV> Delta) const
3061  {
3062  RCP<TraceMinProjRitzOp<ScalarType,MV,OP> > projOp;
3063 
3064  if(computeAllRes_)
3065  {
3066  // Get the valid indices of X
3067  std::vector<int> curind(blockSize_);
3068  for(int i=0; i<blockSize_; i++)
3069  curind[i] = i;
3070 
3071  RCP<const MV> locMX = MVT::CloneView(*MX_, curind);
3072 
3073  // We should really project out the converged eigenvectors too
3074  if(projectAllVecs_)
3075  {
3076  if(projectLockedVecs_ && numAuxVecs_ > 0)
3077  projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3078  else
3079  projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3080  }
3081  else
3082  projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX) );
3083 
3084  // Remember, Delta0 must equal 0
3085  // This ensures B-orthogonality between Delta and X
3086  MVT::MvInit(*Delta);
3087 
3088  if(useRHSR_)
3089  {
3090  RCP<const MV> locR = MVT::CloneView(*R_, curind);
3091  projOp->ApplyInverse(*locR, *Delta);
3092  }
3093  else
3094  {
3095  RCP<const MV> locKX = MVT::CloneView(*KX_, curind);
3096  projOp->ApplyInverse(*locKX, *Delta);
3097  }
3098  }
3099  else
3100  {
3101  // We should really project out the converged eigenvectors too
3102  if(projectAllVecs_)
3103  {
3104  if(projectLockedVecs_ && numAuxVecs_ > 0)
3105  projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3106  else
3107  projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3108  }
3109  else
3110  projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_) );
3111 
3112  // Remember, Delta0 must equal 0
3113  // This ensures B-orthogonality between Delta and X
3114  MVT::MvInit(*Delta);
3115 
3116  if(useRHSR_) {
3117  projOp->ApplyInverse(*R_, *Delta);
3118  }
3119  else {
3120  projOp->ApplyInverse(*KX_, *Delta);
3121  }
3122  }
3123  }
3124 
3126  // TODO: Fix preconditioning
3127  template <class ScalarType, class MV, class OP>
3128  void TraceMinBase<ScalarType,MV,OP>::solveSaddleProjPrec (RCP<MV> Delta) const
3129  {
3130  // If we don't have Belos installed, we can't use TraceMinProjRitzOpWithPrec
3131  // Of course, this problem will be detected in the constructor and an exception will be thrown
3132  // This is only here to make sure the code compiles properly
3133 #ifdef HAVE_ANASAZI_BELOS
3134  RCP<TraceMinProjRitzOpWithPrec<ScalarType,MV,OP> > projOp;
3135 
3136  if(computeAllRes_)
3137  {
3138  int dimension;
3139  if(projectAllVecs_)
3140  dimension = curDim_;
3141  else
3142  dimension = blockSize_;
3143 
3144  // Get the valid indices of X
3145  std::vector<int> curind(dimension);
3146  for(int i=0; i<dimension; i++)
3147  curind[i] = i;
3148 
3149  RCP<const MV> locMX = MVT::CloneView(*MX_, curind);
3150 
3151  // We should really project out the converged eigenvectors too
3152  if(projectAllVecs_)
3153  {
3154  if(projectLockedVecs_ && numAuxVecs_ > 0)
3155  projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3156  else
3157  projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3158  }
3159  else
3160  projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX) );
3161 
3162  // Remember, Delta0 must equal 0
3163  // This ensures B-orthogonality between Delta and X
3164  MVT::MvInit(*Delta);
3165 
3166  std::vector<int> dimind(blockSize_);
3167  for(int i=0; i<blockSize_; i++)
3168  dimind[i] = i;
3169 
3170  if(useRHSR_)
3171  {
3172  RCP<const MV> locR = MVT::CloneView(*R_, dimind);
3173  projOp->ApplyInverse(*locR, *Delta);
3174  MVT::MvScale(*Delta, -ONE);
3175  }
3176  else
3177  {
3178  RCP<const MV> locKX = MVT::CloneView(*KX_, dimind);
3179  projOp->ApplyInverse(*locKX, *Delta);
3180  }
3181  }
3182  else
3183  {
3184  // We should really project out the converged eigenvectors too
3185  if(projectAllVecs_)
3186  {
3187  if(projectLockedVecs_ && numAuxVecs_ > 0)
3188  projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3189  else
3190  projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3191  }
3192  else
3193  projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_) );
3194 
3195  // Remember, Delta0 must equal 0
3196  // This ensures B-orthogonality between Delta and X
3197  MVT::MvInit(*Delta);
3198 
3199  if(useRHSR_)
3200  {
3201  projOp->ApplyInverse(*R_, *Delta);
3202  MVT::MvScale(*Delta,-ONE);
3203  }
3204  else
3205  projOp->ApplyInverse(*KX_, *Delta);
3206  }
3207 #endif
3208  }
3209 
3211  // TODO: We can hold the Schur complement constant in later iterations
3212  // TODO: Make sure we're using the preconditioner correctly
3213  template <class ScalarType, class MV, class OP>
3214  void TraceMinBase<ScalarType,MV,OP>::solveSaddleSchur (RCP<MV> Delta) const
3215  {
3216  // dense solver
3217  Teuchos::SerialDenseSolver<int,ScalarType> My_Solver;
3218 
3219  RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclL;
3220  RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclS = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(blockSize_,blockSize_) );
3221 
3222  if(computeAllRes_)
3223  {
3224  // Get the valid indices of X
3225  std::vector<int> curind(blockSize_);
3226  for(int i=0; i<blockSize_; i++)
3227  curind[i] = i;
3228 
3229  // Z = K \ MX
3230  // Why would I have wanted to set Z <- X? I want to leave Z's previous value alone...
3231  RCP<const MV> lclMX = MVT::CloneView(*MX_, curind);
3232 
3233 #ifdef USE_APPLY_INVERSE
3234  Op_->ApplyInverse(*lclMX,*Z_);
3235 #else
3236  ritzOp_->ApplyInverse(*lclMX,*Z_);
3237 #endif
3238 
3239  // form S = X' M Z
3240  MVT::MvTransMv(ONE,*Z_,*lclMX,*lclS);
3241 
3242  // solve S L = I
3243  My_Solver.setMatrix(lclS);
3244  My_Solver.invert();
3245  lclL = lclS;
3246 
3247  // Delta = X - Z L
3248  RCP<const MV> lclX = MVT::CloneView(*X_, curind);
3249  MVT::Assign(*lclX,*Delta);
3250  MVT::MvTimesMatAddMv( -ONE, *Z_, *lclL, ONE, *Delta);
3251  }
3252  else
3253  {
3254  // Z = K \ MX
3255 #ifdef USE_APPLY_INVERSE
3256  Op_->ApplyInverse(*MX_,*Z_);
3257 #else
3258  ritzOp_->ApplyInverse(*MX_,*Z_);
3259 #endif
3260 
3261  // form S = X' M Z
3262  MVT::MvTransMv(ONE,*Z_,*MX_,*lclS);
3263 
3264  // solve S L = I
3265  My_Solver.setMatrix(lclS);
3266  My_Solver.invert();
3267  lclL = lclS;
3268 
3269  // Delta = X - Z L
3270  MVT::Assign(*X_,*Delta);
3271  MVT::MvTimesMatAddMv( -ONE, *Z_, *lclL, ONE, *Delta);
3272  }
3273  }
3274 
3275 
3277  // TODO: We can hold the Schur complement constant in later iterations
3278  template <class ScalarType, class MV, class OP>
3279  void TraceMinBase<ScalarType,MV,OP>::solveSaddleBDPrec (RCP<MV> Delta) const
3280  {
3281  RCP<MV> locKX, locMX;
3282  if(computeAllRes_)
3283  {
3284  std::vector<int> curind(blockSize_);
3285  for(int i=0; i<blockSize_; i++)
3286  curind[i] = i;
3287  locKX = MVT::CloneViewNonConst(*KX_, curind);
3288  locMX = MVT::CloneViewNonConst(*MX_, curind);
3289  }
3290  else
3291  {
3292  locKX = KX_;
3293  locMX = MX_;
3294  }
3295 
3296  // Create the operator [A BX; X'B 0]
3297  RCP<saddle_op_type> sadOp = rcp(new saddle_op_type(ritzOp_,locMX));
3298 
3299  // Create the RHS [AX; 0]
3300  RCP<saddle_container_type> sadRHS = rcp(new saddle_container_type(locKX));
3301 
3302 // locKX->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
3303 
3304 // locMX->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
3305 
3306  // Create the solution vector [Delta; L]
3307  MVT::MvInit(*Delta);
3308  RCP<saddle_container_type> sadSol = rcp(new saddle_container_type(Delta));
3309 
3310  // Create a minres solver
3311  RCP<PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type > > sadSolver;
3312  if(Prec_ != Teuchos::null)
3313  {
3314  RCP<saddle_op_type> sadPrec = rcp(new saddle_op_type(ritzOp_->getPrec(),locMX,BD_PREC));
3315  sadSolver = rcp(new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp, sadPrec));
3316  }
3317  else {
3318  sadSolver = rcp(new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp));
3319  }
3320 
3321  // Set the tolerance for the minres solver
3322  std::vector<ScalarType> tol;
3323  ritzOp_->getInnerTol(tol);
3324  sadSolver->setTol(tol);
3325 
3326  // Set the maximum number of iterations
3327  sadSolver->setMaxIter(ritzOp_->getMaxIts());
3328 
3329  // Set the solution vector
3330  sadSolver->setSol(sadSol);
3331 
3332  // Set the RHS
3333  sadSolver->setRHS(sadRHS);
3334 
3335  // Solve the saddle point problem
3336  sadSolver->solve();
3337  }
3338 
3339 
3340 
3342  // TODO: We can hold the Schur complement constant in later iterations
3343  template <class ScalarType, class MV, class OP>
3344  void TraceMinBase<ScalarType,MV,OP>::solveSaddleHSSPrec (RCP<MV> Delta) const
3345  {
3346 #ifdef HAVE_ANASAZI_BELOS
3347  typedef Belos::LinearProblem<ScalarType,saddle_container_type,saddle_op_type> LP;
3348  typedef Belos::PseudoBlockGmresSolMgr<ScalarType,saddle_container_type,saddle_op_type> GmresSolMgr;
3349 
3350  RCP<MV> locKX, locMX;
3351  if(computeAllRes_)
3352  {
3353  std::vector<int> curind(blockSize_);
3354  for(int i=0; i<blockSize_; i++)
3355  curind[i] = i;
3356  locKX = MVT::CloneViewNonConst(*KX_, curind);
3357  locMX = MVT::CloneViewNonConst(*MX_, curind);
3358  }
3359  else
3360  {
3361  locKX = KX_;
3362  locMX = MX_;
3363  }
3364 
3365  // Create the operator [A BX; X'B 0]
3366  RCP<saddle_op_type> sadOp = rcp(new saddle_op_type(ritzOp_,locMX,NONSYM));
3367 
3368  // Create the RHS [AX; 0]
3369  RCP<saddle_container_type> sadRHS = rcp(new saddle_container_type(locKX));
3370 
3371  // Create the solution vector [Delta; L]
3372  MVT::MvInit(*Delta);
3373  RCP<saddle_container_type> sadSol = rcp(new saddle_container_type(Delta));
3374 
3375  // Create a parameter list for the gmres solver
3376  RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
3377 
3378  // Set the tolerance for the gmres solver
3379  std::vector<ScalarType> tol;
3380  ritzOp_->getInnerTol(tol);
3381  pl->set("Convergence Tolerance",tol[0]);
3382 
3383  // Set the maximum number of iterations
3384  // TODO: Come back to this
3385  pl->set("Maximum Iterations", ritzOp_->getMaxIts());
3386  pl->set("Num Blocks", ritzOp_->getMaxIts());
3387 
3388  // Set the block size
3389  // TODO: Come back to this
3390  // TODO: This breaks the code right now, presumably because of a MVT cloneview issue.
3391  pl->set("Block Size", blockSize_);
3392 
3393  // Set the verbosity of gmres
3394 // pl->set("Verbosity", Belos::IterationDetails + Belos::StatusTestDetails + Belos::Debug);
3395 // pl->set("Output Frequency", 1);
3396 
3397  // Create the linear problem
3398  RCP<LP> problem = rcp(new LP(sadOp,sadSol,sadRHS));
3399 
3400  // Set the preconditioner
3401  if(Prec_ != Teuchos::null)
3402  {
3403  RCP<saddle_op_type> sadPrec = rcp(new saddle_op_type(ritzOp_->getPrec(),locMX,HSS_PREC,alpha_));
3404  problem->setLeftPrec(sadPrec);
3405  }
3406 
3407  // Set the problem
3408  problem->setProblem();
3409 
3410  // Create a minres solver
3411  RCP<GmresSolMgr> sadSolver = rcp(new GmresSolMgr(problem,pl)) ;
3412 
3413  // Solve the saddle point problem
3414  sadSolver->solve();
3415 #else
3416  std::cout << "No Belos. This is bad\n";
3417 #endif
3418  }
3419 
3420 
3421 
3423  // Compute KK := V'KV
3424  // We only compute the NEW elements
3425  template <class ScalarType, class MV, class OP>
3427  {
3428  // Get the valid indices of V
3429  std::vector<int> curind(curDim_);
3430  for(int i=0; i<curDim_; i++)
3431  curind[i] = i;
3432 
3433  // Get a pointer to the valid parts of V
3434  RCP<const MV> lclV = MVT::CloneView(*V_,curind);
3435 
3436  // Get the valid indices of KV
3437  curind.resize(blockSize_);
3438  for(int i=0; i<blockSize_; i++)
3439  curind[i] = curDim_-blockSize_+i;
3440  RCP<const MV> lclKV = MVT::CloneView(*KV_,curind);
3441 
3442  // Get a pointer to the valid part of KK
3443  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK =
3444  rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_-blockSize_) );
3445 
3446  // KK := V'KV
3447  MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
3448 
3449  // We only constructed the upper triangular part of the matrix, but that's okay because KK is symmetric!
3450  for(int r=0; r<curDim_; r++)
3451  {
3452  for(int c=0; c<r; c++)
3453  {
3454  (*KK_)(r,c) = (*KK_)(c,r);
3455  }
3456  }
3457  }
3458 
3460  // Compute the eigenpairs of KK, i.e. the Ritz pairs
3461  template <class ScalarType, class MV, class OP>
3463  {
3464  // Get a pointer to the valid part of KK
3465  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK =
3466  rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
3467 
3468  // Get a pointer to the valid part of ritzVecs
3469  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclRV =
3470  rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3471 
3472  // Compute Ritz pairs from KK
3473  {
3474 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3475  Teuchos::TimeMonitor lcltimer( *timerDS_ );
3476 #endif
3477  int rank = curDim_;
3478  Utils::directSolver(curDim_, *lclKK, Teuchos::null, *lclRV, theta_, rank, 10);
3479  // we want all ritz values back
3480  // TODO: This probably should not be an ortho failure
3481  TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,TraceMinBaseOrthoFailure,
3482  "Anasazi::TraceMinBase::computeRitzPairs(): Failed to compute all eigenpairs of KK.");
3483  }
3484 
3485  // Sort ritz pairs
3486  {
3487 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3488  Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
3489 #endif
3490 
3491  std::vector<int> order(curDim_);
3492  //
3493  // sort the first curDim_ values in theta_
3494  if(useHarmonic_)
3495  {
3497  sm.sort(theta_, Teuchos::rcpFromRef(order), curDim_);
3498  }
3499  else
3500  {
3501  sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_); // don't catch exception
3502  }
3503  //
3504  // apply the same ordering to the primitive ritz vectors
3505  Utils::permuteVectors(order,*lclRV);
3506  }
3507  }
3508 
3510  // Compute X := V evecs
3511  template <class ScalarType, class MV, class OP>
3513  {
3514 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3515  Teuchos::TimeMonitor lcltimer( *timerLocal_ );
3516 #endif
3517 
3518  // Get the valid indices of V
3519  std::vector<int> curind(curDim_);
3520  for(int i=0; i<curDim_; i++)
3521  curind[i] = i;
3522 
3523  // Get a pointer to the valid parts of V
3524  RCP<const MV> lclV = MVT::CloneView(*V_,curind);
3525 
3526  if(computeAllRes_)
3527  {
3528  // Capture the relevant eigenvectors
3529  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3530  rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3531 
3532  // X <- lclV*S
3533  RCP<MV> lclX = MVT::CloneViewNonConst(*X_,curind);
3534  MVT::MvTimesMatAddMv( ONE, *lclV, *relevantEvecs, ZERO, *lclX );
3535  }
3536  else
3537  {
3538  // Capture the relevant eigenvectors
3539  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3540  rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,blockSize_) );
3541 
3542  // X <- lclV*S
3543  MVT::MvTimesMatAddMv( ONE, *lclV, *relevantEvecs, ZERO, *X_ );
3544  }
3545  }
3546 
3548  // Compute KX := KV evecs and (if necessary) MX := MV evecs
3549  template <class ScalarType, class MV, class OP>
3551  {
3552 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3553  Teuchos::TimeMonitor lcltimer( *timerLocal_ );
3554 #endif
3555 
3556  // Get the valid indices of V
3557  std::vector<int> curind(curDim_);
3558  for(int i=0; i<curDim_; i++)
3559  curind[i] = i;
3560 
3561  // Get pointers to the valid parts of V, KV, and MV (if necessary)
3562  RCP<const MV> lclV = MVT::CloneView(*V_,curind);
3563  RCP<const MV> lclKV = MVT::CloneView(*KV_,curind);
3564 
3565  if(computeAllRes_)
3566  {
3567  // Capture the relevant eigenvectors
3568  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3569  rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3570 
3571  // Update KX and MX
3572  RCP<MV> lclKX = MVT::CloneViewNonConst(*KX_,curind);
3573  MVT::MvTimesMatAddMv( ONE, *lclKV, *relevantEvecs, ZERO, *lclKX );
3574  if(hasM_)
3575  {
3576  RCP<const MV> lclMV = MVT::CloneView(*MV_,curind);
3577  RCP<MV> lclMX = MVT::CloneViewNonConst(*MX_,curind);
3578  MVT::MvTimesMatAddMv( ONE, *lclMV, *relevantEvecs, ZERO, *lclMX );
3579  }
3580  }
3581  else
3582  {
3583  // Capture the relevant eigenvectors
3584  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3585  rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,blockSize_) );
3586 
3587  // Update KX and MX
3588  MVT::MvTimesMatAddMv( ONE, *lclKV, *relevantEvecs, ZERO, *KX_ );
3589  if(hasM_)
3590  {
3591  RCP<const MV> lclMV = MVT::CloneView(*MV_,curind);
3592  MVT::MvTimesMatAddMv( ONE, *lclMV, *relevantEvecs, ZERO, *MX_ );
3593  }
3594  }
3595  }
3596 
3598  // Update the residual R := KX - MX*T
3599  template <class ScalarType, class MV, class OP>
3601 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3602  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
3603 #endif
3604 
3605  if(computeAllRes_)
3606  {
3607  // Get the valid indices of X
3608  std::vector<int> curind(curDim_);
3609  for(int i=0; i<curDim_; i++)
3610  curind[i] = i;
3611 
3612  // Holds MX*T
3613  RCP<MV> MXT = MVT::CloneCopy(*MX_,curind);
3614 
3615  // Holds the relevant part of theta
3616  std::vector<ScalarType> locTheta(curDim_);
3617  for(int i=0; i<curDim_; i++)
3618  locTheta[i] = theta_[i];
3619 
3620  // Compute MX*T
3621  MVT::MvScale(*MXT,locTheta);
3622 
3623  // form R <- KX - MX*T
3624  RCP<const MV> locKX = MVT::CloneView(*KX_,curind);
3625  RCP<MV> locR = MVT::CloneViewNonConst(*R_,curind);
3626  MVT::MvAddMv(ONE,*locKX,-ONE,*MXT,*locR);
3627  }
3628  else
3629  {
3630  // Holds MX*T
3631  RCP<MV> MXT = MVT::CloneCopy(*MX_);
3632 
3633  // Holds the relevant part of theta
3634  std::vector<ScalarType> locTheta(blockSize_);
3635  for(int i=0; i<blockSize_; i++)
3636  locTheta[i] = theta_[i];
3637 
3638  // Compute MX*T
3639  MVT::MvScale(*MXT,locTheta);
3640 
3641  // form R <- KX - MX*T
3642  MVT::MvAddMv(ONE,*KX_,-ONE,*MXT,*R_);
3643  }
3644 
3645  // R has been updated; mark the norms as out-of-date
3646  Rnorms_current_ = false;
3647  R2norms_current_ = false;
3648  }
3649 
3650 
3652  // Check accuracy, orthogonality, and other debugging stuff
3653  //
3654  // bools specify which tests we want to run (instead of running more than we actually care about)
3655  //
3656  // we don't bother checking the following because they are computed explicitly:
3657  // H == Prec*R
3658  // KH == K*H
3659  //
3660  //
3661  // checkV : V orthonormal
3662  // orthogonal to auxvecs
3663  // checkX : X orthonormal
3664  // orthogonal to auxvecs
3665  // checkMX: check MX == M*X
3666  // checkKX: check KX == K*X
3667  // checkH : H orthonormal
3668  // orthogonal to V and H and auxvecs
3669  // checkMH: check MH == M*H
3670  // checkR : check R orthogonal to X
3671  // checkQ : check that auxiliary vectors are actually orthonormal
3672  // checkKK: check that KK is symmetric in memory
3673  //
3674  // TODO:
3675  // add checkTheta
3676  //
3677  template <class ScalarType, class MV, class OP>
3678  std::string TraceMinBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
3679  {
3680  using std::endl;
3681 
3682  std::stringstream os;
3683  os.precision(2);
3684  os.setf(std::ios::scientific, std::ios::floatfield);
3685 
3686  os << " Debugging checks: iteration " << iter_ << where << endl;
3687 
3688  // V and friends
3689  std::vector<int> lclind(curDim_);
3690  for (int i=0; i<curDim_; ++i) lclind[i] = i;
3691  RCP<const MV> lclV;
3692  if (initialized_) {
3693  lclV = MVT::CloneView(*V_,lclind);
3694  }
3695  if (chk.checkV && initialized_) {
3696  MagnitudeType err = orthman_->orthonormError(*lclV);
3697  os << " >> Error in V^H M V == I : " << err << endl;
3698  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3699  err = orthman_->orthogError(*lclV,*auxVecs_[i]);
3700  os << " >> Error in V^H M Q[" << i << "] == 0 : " << err << endl;
3701  }
3702  }
3703 
3704  // X and friends
3705  RCP<const MV> lclX;
3706  if(initialized_)
3707  {
3708  if(computeAllRes_)
3709  lclX = MVT::CloneView(*X_,lclind);
3710  else
3711  lclX = X_;
3712  }
3713 
3714  if (chk.checkX && initialized_) {
3715  MagnitudeType err = orthman_->orthonormError(*lclX);
3716  os << " >> Error in X^H M X == I : " << err << endl;
3717  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3718  err = orthman_->orthogError(*lclX,*auxVecs_[i]);
3719  os << " >> Error in X^H M Q[" << i << "] == 0 : " << err << endl;
3720  }
3721  }
3722  if (chk.checkMX && hasM_ && initialized_) {
3723  RCP<const MV> lclMX;
3724  if(computeAllRes_)
3725  lclMX = MVT::CloneView(*MX_,lclind);
3726  else
3727  lclMX = MX_;
3728 
3729  MagnitudeType err = Utils::errorEquality(*lclX, *lclMX, MOp_);
3730  os << " >> Error in MX == M*X : " << err << endl;
3731  }
3732  if (Op_ != Teuchos::null && chk.checkKX && initialized_) {
3733  RCP<const MV> lclKX;
3734  if(computeAllRes_)
3735  lclKX = MVT::CloneView(*KX_,lclind);
3736  else
3737  lclKX = KX_;
3738 
3739  MagnitudeType err = Utils::errorEquality(*lclX, *lclKX, Op_);
3740  os << " >> Error in KX == K*X : " << err << endl;
3741  }
3742 
3743  // KK
3744  if (chk.checkKK && initialized_) {
3745  Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
3746  if(Op_ != Teuchos::null) {
3747  RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
3748  OPT::Apply(*Op_,*lclV,*lclKV);
3749  MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
3750  }
3751  else {
3752  MVT::MvTransMv(ONE,*lclV,*lclV,curKK);
3753  }
3754  Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
3755  curKK -= subKK;
3756  os << " >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
3757 
3758  Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_);
3759  for (int j=0; j<curDim_; ++j) {
3760  for (int i=0; i<curDim_; ++i) {
3761  SDMerr(i,j) = subKK(i,j) - SCT::conjugate(subKK(j,i));
3762  }
3763  }
3764  os << " >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
3765  }
3766 
3767  // Q
3768  if (chk.checkQ) {
3769  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3770  MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
3771  os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << err << endl;
3772  for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
3773  err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
3774  os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << err << endl;
3775  }
3776  }
3777  }
3778 
3779  os << endl;
3780 
3781  return os.str();
3782  }
3783 
3784 }} // End of namespace Anasazi
3785 
3786 #endif
3787 
3788 // End of file AnasaziTraceMinBase.hpp
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
RCP< const MV > getRitzVectors()
Get access to the current Ritz vectors.
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 .
int NEV
Number of unconverged eigenvalues.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
RCP< const MV > X
The current eigenvectors.
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...
void setBlockSize(int blockSize)
Set the blocksize.
RCP< const MV > V
The current basis.
Declaration of basic traits for the multivector type.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms, computing the norms if they are not up-to-date with the current resid...
int getBlockSize() const
Get the blocksize used by the iterative solver.
An operator of the form [A Y; Y&#39; 0] where A is a sparse matrix and Y a multivector.
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
TraceMinBaseState< ScalarType, MV > getState() const
Get access to the current state of the eigensolver.
TraceMinBaseOrthoFailure is thrown when the orthogonalization manager is unable to orthogonalize the ...
Virtual base class which defines basic traits for the operator type.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms, computing the norms if they are not up-to-date with the current res...
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For the trace minimization methods...
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
static Teuchos::ScalarTraits< ScalarType >::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP< const OP > M=Teuchos::null)
Return the maximum coefficient of the matrix scaled by the maximum coefficient of MX...
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Basic implementation of the Anasazi::SortManager class.
TraceMinBaseInitFailure is thrown when the TraceMinBase solver is unable to generate an initial itera...
static void Assign(const MV &A, MV &mv)
mv := A
An exception class parent to all Anasazi exceptions.
std::vector< int > getRitzIndex()
Get the index used for extracting individual Ritz vectors from getRitzVectors().
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to the given output stream...
bool isOrtho
Whether V has been projected and orthonormalized already.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi&#39;s templated, static class providing utilities for the solvers.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
virtual ~TraceMinBase()
Anasazi::TraceMinBase destructor.
Output managers remove the need for the eigensolver to know any information about the required output...
static int directSolver(int size, const Teuchos::SerialDenseMatrix< int, ScalarType > &KK, Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > MM, Teuchos::SerialDenseMatrix< int, ScalarType > &EV, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &theta, int &nev, int esType=0)
Routine for computing the first NEV generalized eigenpairs of the Hermitian pencil (KK...
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
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).
void iterate()
This method performs trace minimization iterations until the status test indicates the need to stop o...
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Virtual base class which defines basic traits for the operator type.
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...
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).
RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
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 int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values for the previous iteration.
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
void sort(std::vector< MagnitudeType > &evals, Teuchos::RCP< std::vector< int > > perm=Teuchos::null, int n=-1) const
Sort real eigenvalues, optionally returning the permutation vector.
static void permuteVectors(const int n, const std::vector< int > &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *resids=0)
Permute the vectors in a multivector according to the permutation vector perm, and optionally the res...
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.
Types and exceptions used within Anasazi solvers and interfaces.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
This is an abstract base class for the trace minimization eigensolvers.
int getNumIters() const
Get the current iteration count.
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Stores a set of vectors of the form [V; L] where V is a distributed multivector and L is a serialdens...
Common interface of stopping criteria for Anasazi&#39;s solvers.
void setAuxVecs(const Teuchos::Array< RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
void resetNumIters()
Reset the iteration count.
TraceMinBase(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const RCP< OutputManager< ScalarType > > &printer, const RCP< StatusTest< ScalarType, MV, OP > > &tester, const RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
TraceMinBase constructor with eigenproblem, solver utilities, and parameter list of solver options...
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.