Anasazi  Version of the Day
AnasaziRTRBase.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 
33 #ifndef ANASAZI_RTRBASE_HPP
34 #define ANASAZI_RTRBASE_HPP
35 
36 #include "AnasaziTypes.hpp"
37 
38 #include "AnasaziEigensolver.hpp"
41 #include "Teuchos_ScalarTraits.hpp"
42 
44 #include "AnasaziSolverUtils.hpp"
45 
46 #include "Teuchos_LAPACK.hpp"
47 #include "Teuchos_BLAS.hpp"
48 #include "Teuchos_SerialDenseMatrix.hpp"
49 #include "Teuchos_ParameterList.hpp"
50 #include "Teuchos_TimeMonitor.hpp"
51 
101 namespace Anasazi {
102 
104 
105 
110  template <class ScalarType, class MV>
111  struct RTRState {
113  Teuchos::RCP<const MV> X;
115  Teuchos::RCP<const MV> AX;
117  Teuchos::RCP<const MV> BX;
119  Teuchos::RCP<const MV> R;
121  Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
125  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType rho;
126  RTRState() : X(Teuchos::null),AX(Teuchos::null),BX(Teuchos::null),
127  R(Teuchos::null),T(Teuchos::null),rho(0) {};
128  };
129 
131 
133 
134 
148  class RTRRitzFailure : public AnasaziError {public:
149  RTRRitzFailure(const std::string& what_arg) : AnasaziError(what_arg)
150  {}};
151 
160  class RTRInitFailure : public AnasaziError {public:
161  RTRInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
162  {}};
163 
180  class RTROrthoFailure : public AnasaziError {public:
181  RTROrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
182  {}};
183 
184 
186 
187 
188  template <class ScalarType, class MV, class OP>
189  class RTRBase : public Eigensolver<ScalarType,MV,OP> {
190  public:
191 
193 
194 
200  RTRBase(const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
201  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
202  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
203  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
204  const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
205  Teuchos::ParameterList &params,
206  const std::string &solverLabel, bool skinnySolver
207  );
208 
210  virtual ~RTRBase() {};
211 
213 
215 
216 
238  virtual void iterate() = 0;
239 
264  void initialize(RTRState<ScalarType,MV>& newstate);
265 
269  void initialize();
270 
283  bool isInitialized() const;
284 
293 
295 
297 
298 
300  int getNumIters() const;
301 
303  void resetNumIters();
304 
312  Teuchos::RCP<const MV> getRitzVectors();
313 
319  std::vector<Value<ScalarType> > getRitzValues();
320 
328  std::vector<int> getRitzIndex();
329 
335  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
336 
337 
343  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
344 
345 
350  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
351 
352 
361  int getCurSubspaceDim() const;
362 
365  int getMaxSubspaceDim() const;
366 
368 
370 
371 
373  void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
374 
376  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
377 
380 
381 
390  void setBlockSize(int blockSize);
391 
392 
394  int getBlockSize() const;
395 
396 
417  void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
418 
420  Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
421 
423 
425 
426 
428  virtual void currentStatus(std::ostream &os);
429 
431 
432  protected:
433  //
434  // Convenience typedefs
435  //
436  typedef SolverUtils<ScalarType,MV,OP> Utils;
437  typedef MultiVecTraits<ScalarType,MV> MVT;
439  typedef Teuchos::ScalarTraits<ScalarType> SCT;
440  typedef typename SCT::magnitudeType MagnitudeType;
441  typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
442  const MagnitudeType ONE;
443  const MagnitudeType ZERO;
444  const MagnitudeType NANVAL;
445  typedef typename std::vector<MagnitudeType>::iterator vecMTiter;
446  typedef typename std::vector<ScalarType>::iterator vecSTiter;
447  //
448  // Internal structs
449  //
450  struct CheckList {
451  bool checkX, checkAX, checkBX;
452  bool checkEta, checkAEta, checkBEta;
453  bool checkR, checkQ, checkBR;
454  bool checkZ, checkPBX;
455  CheckList() : checkX(false),checkAX(false),checkBX(false),
456  checkEta(false),checkAEta(false),checkBEta(false),
457  checkR(false),checkQ(false),checkBR(false),
458  checkZ(false), checkPBX(false) {};
459  };
460  //
461  // Internal methods
462  //
463  std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
464  // solves the model minimization
465  virtual void solveTRSubproblem() = 0;
466  // Riemannian metric
467  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi) const;
468  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi, const MV &zeta) const;
469  void ginnersep(const MV &xi, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
470  void ginnersep(const MV &xi, const MV &zeta, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
471  //
472  // Classes input through constructor that define the eigenproblem to be solved.
473  //
474  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
475  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
476  const Teuchos::RCP<OutputManager<ScalarType> > om_;
477  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
478  const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > orthman_;
479  //
480  // Information obtained from the eigenproblem
481  //
482  Teuchos::RCP<const OP> AOp_;
483  Teuchos::RCP<const OP> BOp_;
484  Teuchos::RCP<const OP> Prec_;
485  bool hasBOp_, hasPrec_, olsenPrec_;
486  //
487  // Internal timers
488  //
489  Teuchos::RCP<Teuchos::Time> timerAOp_, timerBOp_, timerPrec_,
490  timerSort_,
491  timerLocalProj_, timerDS_,
492  timerLocalUpdate_, timerCompRes_,
493  timerOrtho_, timerInit_;
494  //
495  // Counters
496  //
497  // Number of operator applications
498  int counterAOp_, counterBOp_, counterPrec_;
499 
500  //
501  // Algorithmic parameters.
502  //
503  // blockSize_ is the solver block size
504  int blockSize_;
505  //
506  // Current solver state
507  //
508  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
509  // is capable of running; _initialize is controlled by the initialize() member method
510  // For the implications of the state of initialized_, please see documentation for initialize()
511  bool initialized_;
512  //
513  // nevLocal_ reflects how much of the current basis is valid (0 <= nevLocal_ <= blockSize_)
514  // this tells us how many of the values in theta_ are valid Ritz values
515  int nevLocal_;
516  //
517  // are we implementing a skinny solver? (SkinnyIRTR)
518  bool isSkinny_;
519  //
520  // are we computing leftmost or rightmost eigenvalue?
521  bool leftMost_;
522  //
523  // State Multivecs
524  //
525  // if we are implementing a skinny solver (SkinnyIRTR),
526  // then some of these will never be allocated
527  //
528  // In order to handle auxiliary vectors, we need to handle the projector
529  // P_{[BQ BX],[BQ BX]}
530  // Using an orthomanager with B-inner product, this requires calling with multivectors
531  // [BQ,BX] and [Q,X].
532  // These multivectors must be combined because <[BQ,BX],[Q,X]>_B != I
533  // Hence, we will create two multivectors V and BV, which store
534  // V = [Q,X]
535  // BV = [BQ,BX]
536  //
537  // In the context of preconditioning, we may need to apply the projector
538  // P_{prec*[BQ,BX],[BQ,BX]}
539  // Because [BQ,BX] do not change during the supproblem solver, we will apply
540  // the preconditioner to [BQ,BX] only once. This result is stored in PBV.
541  //
542  // X,BX are views into V,BV
543  // We don't need views for Q,BQ
544  // Inside the subproblem solver, X,BX are constant, so we can allow these
545  // views to exist alongside the full view of V,BV without violating
546  // view semantics.
547  //
548  // Skinny solver allocates 6/7/8 multivectors:
549  // V_, BV_ (if hasB)
550  // PBV_ (if hasPrec and olsenPrec)
551  // R_, Z_ (regardless of hasPrec)
552  // eta_, delta_, Hdelta_
553  //
554  // Hefty solver allocates 10/11/12/13 multivectors:
555  // V_, AX_, BV_ (if hasB)
556  // PBV_ (if hasPrec and olsenPrec)
557  // R_, Z_ (if hasPrec)
558  // eta_, Aeta_, Beta_
559  // delta_, Adelta_, Bdelta_, Hdelta_
560  //
561  Teuchos::RCP<MV> V_, BV_, PBV_, // V = [Q,X]; B*V; Prec*B*V
562  AX_, R_, // A*X_; residual, gradient, and residual of model minimization
563  eta_, Aeta_, Beta_, // update vector from model minimization
564  delta_, Adelta_, Bdelta_, Hdelta_, // search direction in model minimization
565  Z_; // preconditioned residual
566  Teuchos::RCP<const MV> X_, BX_;
567  //
568  // auxiliary vectors
569  Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
570  int numAuxVecs_;
571  //
572  // Number of iterations that have been performed.
573  int iter_;
574  //
575  // Current eigenvalues, residual norms
576  std::vector<MagnitudeType> theta_, Rnorms_, R2norms_, ritz2norms_;
577  //
578  // are the residual norms current with the residual?
579  bool Rnorms_current_, R2norms_current_;
580  //
581  // parameters solver and inner solver
582  MagnitudeType conv_kappa_, conv_theta_;
583  MagnitudeType rho_;
584  //
585  // current objective function value
586  MagnitudeType fx_;
587  };
588 
589 
590 
591 
593  // Constructor
594  template <class ScalarType, class MV, class OP>
596  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
597  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
598  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
599  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
600  const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
601  Teuchos::ParameterList &params,
602  const std::string &solverLabel, bool skinnySolver
603  ) :
604  ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
605  ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
606  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
607  // problem, tools
608  problem_(problem),
609  sm_(sorter),
610  om_(printer),
611  tester_(tester),
612  orthman_(ortho),
613 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
614  timerAOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Operation A*x")),
615  timerBOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Operation B*x")),
616  timerPrec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Operation Prec*x")),
617  timerSort_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Sorting eigenvalues")),
618  timerLocalProj_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Local projection")),
619  timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Direct solve")),
620  timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Local update")),
621  timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Computing residuals")),
622  timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Orthogonalization")),
623  timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Initialization")),
624 #endif
625  counterAOp_(0),
626  counterBOp_(0),
627  counterPrec_(0),
628  // internal data
629  blockSize_(0),
630  initialized_(false),
631  nevLocal_(0),
632  isSkinny_(skinnySolver),
633  leftMost_(true),
634  auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
635  numAuxVecs_(0),
636  iter_(0),
637  Rnorms_current_(false),
638  R2norms_current_(false),
639  conv_kappa_(.1),
640  conv_theta_(1),
641  rho_( MAT::nan() ),
642  fx_( MAT::nan() )
643  {
644  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
645  "Anasazi::RTRBase::constructor: user passed null problem pointer.");
646  TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
647  "Anasazi::RTRBase::constructor: user passed null sort manager pointer.");
648  TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
649  "Anasazi::RTRBase::constructor: user passed null output manager pointer.");
650  TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
651  "Anasazi::RTRBase::constructor: user passed null status test pointer.");
652  TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
653  "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer.");
654  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
655  "Anasazi::RTRBase::constructor: problem is not set.");
656  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
657  "Anasazi::RTRBase::constructor: problem is not Hermitian.");
658 
659  // get the problem operators
660  AOp_ = problem_->getOperator();
661  TEUCHOS_TEST_FOR_EXCEPTION(AOp_ == Teuchos::null, std::invalid_argument,
662  "Anasazi::RTRBase::constructor: problem provides no A matrix.");
663  BOp_ = problem_->getM();
664  Prec_ = problem_->getPrec();
665  hasBOp_ = (BOp_ != Teuchos::null);
666  hasPrec_ = (Prec_ != Teuchos::null);
667  olsenPrec_ = params.get<bool>("Olsen Prec", true);
668 
669  TEUCHOS_TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument,
670  "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product.");
671 
672  // set the block size and allocate data
673  int bs = params.get("Block Size", problem_->getNEV());
674  setBlockSize(bs);
675 
676  // leftmost or rightmost?
677  leftMost_ = params.get("Leftmost",leftMost_);
678 
679  conv_kappa_ = params.get("Kappa Convergence",conv_kappa_);
680  TEUCHOS_TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument,
681  "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
682  conv_theta_ = params.get("Theta Convergence",conv_theta_);
683  TEUCHOS_TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument,
684  "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
685  }
686 
687 
689  // Set the block size and make necessary adjustments.
690  template <class ScalarType, class MV, class OP>
692  {
693  // time spent here counts towards timerInit_
694 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
695  Teuchos::TimeMonitor lcltimer( *timerInit_ );
696 #endif
697 
698  // This routine only allocates space; it doesn't not perform any computation
699  // if solver is initialized and size is to be decreased, take the first blockSize vectors of all to preserve state
700  // otherwise, shrink/grow/allocate space and set solver to unitialized
701 
702  Teuchos::RCP<const MV> tmp;
703  // grab some Multivector to Clone
704  // in practice, getInitVec() should always provide this, but it is possible to use a
705  // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
706  // in case of that strange scenario, we will try to Clone from R_
707  // we like R_ for this, because it has minimal size (blockSize_), as opposed to V_ (blockSize_+numAuxVecs_)
708  if (blockSize_ > 0) {
709  tmp = R_;
710  }
711  else {
712  tmp = problem_->getInitVec();
713  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
714  "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
715  }
716 
717  TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetGlobalLength(*tmp), std::invalid_argument,
718  "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
719 
720  // last chance to quit before causing side-effects
721  if (blockSize == blockSize_) {
722  // do nothing
723  return;
724  }
725 
726  // clear views
727  X_ = Teuchos::null;
728  BX_ = Teuchos::null;
729 
730  // regardless of whether we preserve any data, any content in V, BV and PBV corresponding to the
731  // auxiliary vectors must be retained
732  // go ahead and do these first
733  //
734  // two cases here:
735  // * we are growing (possibly, from empty)
736  // any aux data must be copied over, nothing from X need copying
737  // * we are shrinking
738  // any aux data must be copied over, go ahead and copy X material (initialized or not)
739  //
740  if (blockSize > blockSize_)
741  {
742  // GROWING
743  // get a pointer for Q-related material, and an index for extracting/setting it
744  Teuchos::RCP<const MV> Q;
745  std::vector<int> indQ(numAuxVecs_);
746  for (int i=0; i<numAuxVecs_; i++) indQ[i] = i;
747  // if numAuxVecs_ > 0, then necessarily blockSize_ > 0 (we have already been allocated once)
748  TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error,
749  "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team.");
750  // V
751  if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ);
752  V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
753  if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_);
754  // BV
755  if (hasBOp_) {
756  if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ);
757  BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
758  if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_);
759  }
760  else {
761  BV_ = V_;
762  }
763  // PBV
764  if (hasPrec_ && olsenPrec_) {
765  if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ);
766  PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
767  if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_);
768  }
769  }
770  else
771  {
772  // SHRINKING
773  std::vector<int> indV(numAuxVecs_+blockSize);
774  for (int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
775  // V
776  V_ = MVT::CloneCopy(*V_,indV);
777  // BV
778  if (hasBOp_) {
779  BV_ = MVT::CloneCopy(*BV_,indV);
780  }
781  else {
782  BV_ = V_;
783  }
784  // PBV
785  if (hasPrec_ && olsenPrec_) {
786  PBV_ = MVT::CloneCopy(*PBV_,indV);
787  }
788  }
789 
790  if (blockSize < blockSize_) {
791  // shrink vectors
792  blockSize_ = blockSize;
793 
794  theta_.resize(blockSize_);
795  ritz2norms_.resize(blockSize_);
796  Rnorms_.resize(blockSize_);
797  R2norms_.resize(blockSize_);
798 
799  if (initialized_) {
800  // shrink multivectors with copy
801  std::vector<int> ind(blockSize_);
802  for (int i=0; i<blockSize_; i++) ind[i] = i;
803 
804  // Z can be deleted, no important info there
805  Z_ = Teuchos::null;
806 
807  // we will not use tmp for cloning, clear it and free some space
808  tmp = Teuchos::null;
809 
810  R_ = MVT::CloneCopy(*R_ ,ind);
811  eta_ = MVT::CloneCopy(*eta_ ,ind);
812  delta_ = MVT::CloneCopy(*delta_ ,ind);
813  Hdelta_ = MVT::CloneCopy(*Hdelta_,ind);
814  if (!isSkinny_) {
815  AX_ = MVT::CloneCopy(*AX_ ,ind);
816  Aeta_ = MVT::CloneCopy(*Aeta_ ,ind);
817  Adelta_ = MVT::CloneCopy(*Adelta_,ind);
818  }
819  else {
820  AX_ = Teuchos::null;
821  Aeta_ = Teuchos::null;
822  Adelta_ = Teuchos::null;
823  }
824 
825  if (hasBOp_) {
826  if (!isSkinny_) {
827  Beta_ = MVT::CloneCopy(*Beta_,ind);
828  Bdelta_ = MVT::CloneCopy(*Bdelta_,ind);
829  }
830  else {
831  Beta_ = Teuchos::null;
832  Bdelta_ = Teuchos::null;
833  }
834  }
835  else {
836  Beta_ = eta_;
837  Bdelta_ = delta_;
838  }
839 
840  // we need Z if we have a preconditioner
841  // we also use Z for temp storage in the skinny solvers
842  if (hasPrec_ || isSkinny_) {
843  Z_ = MVT::Clone(*V_,blockSize_);
844  }
845  else {
846  Z_ = R_;
847  }
848 
849  }
850  else {
851  // release previous multivectors
852  // shrink multivectors without copying
853  AX_ = Teuchos::null;
854  R_ = Teuchos::null;
855  eta_ = Teuchos::null;
856  Aeta_ = Teuchos::null;
857  delta_ = Teuchos::null;
858  Adelta_ = Teuchos::null;
859  Hdelta_ = Teuchos::null;
860  Beta_ = Teuchos::null;
861  Bdelta_ = Teuchos::null;
862  Z_ = Teuchos::null;
863 
864  R_ = MVT::Clone(*tmp,blockSize_);
865  eta_ = MVT::Clone(*tmp,blockSize_);
866  delta_ = MVT::Clone(*tmp,blockSize_);
867  Hdelta_ = MVT::Clone(*tmp,blockSize_);
868  if (!isSkinny_) {
869  AX_ = MVT::Clone(*tmp,blockSize_);
870  Aeta_ = MVT::Clone(*tmp,blockSize_);
871  Adelta_ = MVT::Clone(*tmp,blockSize_);
872  }
873 
874  if (hasBOp_) {
875  if (!isSkinny_) {
876  Beta_ = MVT::Clone(*tmp,blockSize_);
877  Bdelta_ = MVT::Clone(*tmp,blockSize_);
878  }
879  }
880  else {
881  Beta_ = eta_;
882  Bdelta_ = delta_;
883  }
884 
885  // we need Z if we have a preconditioner
886  // we also use Z for temp storage in the skinny solvers
887  if (hasPrec_ || isSkinny_) {
888  Z_ = MVT::Clone(*tmp,blockSize_);
889  }
890  else {
891  Z_ = R_;
892  }
893  } // if initialized_
894  } // if blockSize is shrinking
895  else { // if blockSize > blockSize_
896  // this is also the scenario for our initial call to setBlockSize(), in the constructor
897  initialized_ = false;
898 
899  // grow/allocate vectors
900  theta_.resize(blockSize,NANVAL);
901  ritz2norms_.resize(blockSize,NANVAL);
902  Rnorms_.resize(blockSize,NANVAL);
903  R2norms_.resize(blockSize,NANVAL);
904 
905  // deallocate old multivectors
906  AX_ = Teuchos::null;
907  R_ = Teuchos::null;
908  eta_ = Teuchos::null;
909  Aeta_ = Teuchos::null;
910  delta_ = Teuchos::null;
911  Adelta_ = Teuchos::null;
912  Hdelta_ = Teuchos::null;
913  Beta_ = Teuchos::null;
914  Bdelta_ = Teuchos::null;
915  Z_ = Teuchos::null;
916 
917  // clone multivectors off of tmp
918  R_ = MVT::Clone(*tmp,blockSize);
919  eta_ = MVT::Clone(*tmp,blockSize);
920  delta_ = MVT::Clone(*tmp,blockSize);
921  Hdelta_ = MVT::Clone(*tmp,blockSize);
922  if (!isSkinny_) {
923  AX_ = MVT::Clone(*tmp,blockSize);
924  Aeta_ = MVT::Clone(*tmp,blockSize);
925  Adelta_ = MVT::Clone(*tmp,blockSize);
926  }
927 
928  if (hasBOp_) {
929  if (!isSkinny_) {
930  Beta_ = MVT::Clone(*tmp,blockSize);
931  Bdelta_ = MVT::Clone(*tmp,blockSize);
932  }
933  }
934  else {
935  Beta_ = eta_;
936  Bdelta_ = delta_;
937  }
938  if (hasPrec_ || isSkinny_) {
939  Z_ = MVT::Clone(*tmp,blockSize);
940  }
941  else {
942  Z_ = R_;
943  }
944  blockSize_ = blockSize;
945  }
946 
947  // get view of X from V, BX from BV
948  // these are located after the first numAuxVecs columns
949  {
950  std::vector<int> indX(blockSize_);
951  for (int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
952  X_ = MVT::CloneView(*V_,indX);
953  if (hasBOp_) {
954  BX_ = MVT::CloneView(*BV_,indX);
955  }
956  else {
957  BX_ = X_;
958  }
959  }
960  }
961 
962 
964  // Set a new StatusTest for the solver.
965  template <class ScalarType, class MV, class OP>
967  TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
968  "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest.");
969  tester_ = test;
970  }
971 
972 
974  // Get the current StatusTest used by the solver.
975  template <class ScalarType, class MV, class OP>
976  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > RTRBase<ScalarType,MV,OP>::getStatusTest() const {
977  return tester_;
978  }
979 
980 
982  // Set the auxiliary vectors
983  template <class ScalarType, class MV, class OP>
984  void RTRBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
985  typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv;
986 
987  // set new auxiliary vectors
988  auxVecs_.resize(0);
989  auxVecs_.reserve(auxvecs.size());
990 
991  numAuxVecs_ = 0;
992  for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
993  numAuxVecs_ += MVT::GetNumberVecs(**v);
994  }
995 
996  // If the solver has been initialized, X is not necessarily orthogonal to new auxiliary vectors
997  if (numAuxVecs_ > 0 && initialized_) {
998  initialized_ = false;
999  }
1000 
1001  // clear X,BX views
1002  X_ = Teuchos::null;
1003  BX_ = Teuchos::null;
1004  // deallocate, we'll clone off R_ below
1005  V_ = Teuchos::null;
1006  BV_ = Teuchos::null;
1007  PBV_ = Teuchos::null;
1008 
1009  // put auxvecs into V, update BV and PBV if necessary
1010  if (numAuxVecs_ > 0) {
1011  V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1012  int numsofar = 0;
1013  for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
1014  std::vector<int> ind(MVT::GetNumberVecs(**v));
1015  for (size_t j=0; j<ind.size(); j++) ind[j] = numsofar++;
1016  MVT::SetBlock(**v,ind,*V_);
1017  auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind));
1018  }
1019  TEUCHOS_TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error,
1020  "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team.");
1021  // compute B*V, Prec*B*V
1022  if (hasBOp_) {
1023  BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1024  OPT::Apply(*BOp_,*V_,*BV_);
1025  }
1026  else {
1027  BV_ = V_;
1028  }
1029  if (hasPrec_ && olsenPrec_) {
1030  PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1031  OPT::Apply(*Prec_,*BV_,*V_);
1032  }
1033  }
1034  //
1035 
1036  if (om_->isVerbosity( Debug ) ) {
1037  // Check almost everything here
1038  CheckList chk;
1039  chk.checkQ = true;
1040  om_->print( Debug, accuracyCheck(chk, "in setAuxVecs()") );
1041  }
1042  }
1043 
1044 
1046  /* Initialize the state of the solver
1047  *
1048  * POST-CONDITIONS:
1049  *
1050  * initialized_ == true
1051  * X is orthonormal, orthogonal to auxVecs_
1052  * AX = A*X if not skinnny
1053  * BX = B*X if hasBOp_
1054  * theta_ contains Ritz values of X
1055  * R = AX - BX*diag(theta_)
1056  */
1057  template <class ScalarType, class MV, class OP>
1059  {
1060  // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
1061  // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
1062 
1063  // clear const views to X,BX in V,BV
1064  // work with temporary non-const views
1065  X_ = Teuchos::null;
1066  BX_ = Teuchos::null;
1067  Teuchos::RCP<MV> X, BX;
1068  {
1069  std::vector<int> ind(blockSize_);
1070  for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1071  X = MVT::CloneViewNonConst(*V_,ind);
1072  if (hasBOp_) {
1073  BX = MVT::CloneViewNonConst(*BV_,ind);
1074  }
1075  else {
1076  BX = X;
1077  }
1078  }
1079 
1080 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1081  Teuchos::TimeMonitor inittimer( *timerInit_ );
1082 #endif
1083 
1084  std::vector<int> bsind(blockSize_);
1085  for (int i=0; i<blockSize_; i++) bsind[i] = i;
1086 
1087  // in RTR, X (the subspace iterate) is primary
1088  // the order of dependence follows like so.
1089  // --init-> X
1090  // --op apply-> AX,BX
1091  // --ritz analysis-> theta
1092  //
1093  // if the user specifies all data for a level, we will accept it.
1094  // otherwise, we will generate the whole level, and all subsequent levels.
1095  //
1096  // the data members are ordered based on dependence, and the levels are
1097  // partitioned according to the amount of work required to produce the
1098  // items in a level.
1099  //
1100  // inconsitent multivectors widths and lengths will not be tolerated, and
1101  // will be treated with exceptions.
1102 
1103  // set up X, AX, BX: get them from "state" if user specified them
1104  if (newstate.X != Teuchos::null) {
1105  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X),
1106  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
1107  // newstate.X must have blockSize_ vectors; any more will be ignored
1108  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
1109  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
1110 
1111  // put data in X
1112  MVT::SetBlock(*newstate.X,bsind,*X);
1113 
1114  // put data in AX
1115  // if we are implementing a skinny solver, then we don't have memory allocated for AX
1116  // in this case, point AX at Z (skinny solvers allocate Z) and use it for temporary storage
1117  // we will clear the pointer later
1118  if (isSkinny_) {
1119  AX_ = Z_;
1120  }
1121  if (newstate.AX != Teuchos::null) {
1122  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.AX) != MVT::GetGlobalLength(*X),
1123  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
1124  // newstate.AX must have blockSize_ vectors; any more will be ignored
1125  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.AX) < blockSize_,
1126  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
1127  MVT::SetBlock(*newstate.AX,bsind,*AX_);
1128  }
1129  else {
1130  {
1131 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1132  Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1133 #endif
1134  OPT::Apply(*AOp_,*X,*AX_);
1135  counterAOp_ += blockSize_;
1136  }
1137  // we generated AX; we will generate R as well
1138  newstate.R = Teuchos::null;
1139  }
1140 
1141  // put data in BX
1142  // skinny solvers always allocate BX if hasB, so this is unconditionally appropriate
1143  if (hasBOp_) {
1144  if (newstate.BX != Teuchos::null) {
1145  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.BX) != MVT::GetGlobalLength(*X),
1146  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
1147  // newstate.BX must have blockSize_ vectors; any more will be ignored
1148  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.BX) < blockSize_,
1149  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
1150  MVT::SetBlock(*newstate.BX,bsind,*BX);
1151  }
1152  else {
1153  {
1154 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1155  Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1156 #endif
1157  OPT::Apply(*BOp_,*X,*BX);
1158  counterBOp_ += blockSize_;
1159  }
1160  // we generated BX; we will generate R as well
1161  newstate.R = Teuchos::null;
1162  }
1163  }
1164  else {
1165  // the assignment BX_==X_ would be redundant; take advantage of this opportunity to debug a little
1166  TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1167  }
1168 
1169  }
1170  else {
1171  // user did not specify X
1172 
1173  // clear state so we won't use any data from it below
1174  newstate.R = Teuchos::null;
1175  newstate.T = Teuchos::null;
1176 
1177  // generate something and projectAndNormalize
1178  Teuchos::RCP<const MV> ivec = problem_->getInitVec();
1179  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
1180  "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1181 
1182  int initSize = MVT::GetNumberVecs(*ivec);
1183  if (initSize > blockSize_) {
1184  // we need only the first blockSize_ vectors from ivec; get a view of them
1185  initSize = blockSize_;
1186  std::vector<int> ind(blockSize_);
1187  for (int i=0; i<blockSize_; i++) ind[i] = i;
1188  ivec = MVT::CloneView(*ivec,ind);
1189  }
1190 
1191  // assign ivec to first part of X
1192  if (initSize > 0) {
1193  std::vector<int> ind(initSize);
1194  for (int i=0; i<initSize; i++) ind[i] = i;
1195  MVT::SetBlock(*ivec,ind,*X);
1196  }
1197  // fill the rest of X with random
1198  if (blockSize_ > initSize) {
1199  std::vector<int> ind(blockSize_ - initSize);
1200  for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1201  Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X,ind);
1202  MVT::MvRandom(*rX);
1203  rX = Teuchos::null;
1204  }
1205 
1206  // put data in BX
1207  if (hasBOp_) {
1208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1209  Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1210 #endif
1211  OPT::Apply(*BOp_,*X,*BX);
1212  counterBOp_ += blockSize_;
1213  }
1214  else {
1215  // the assignment BX==X would be redundant; take advantage of this opportunity to debug a little
1216  TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1217  }
1218 
1219  // remove auxVecs from X and normalize it
1220  if (numAuxVecs_ > 0) {
1221 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1222  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1223 #endif
1224  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1225  int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX);
1226  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
1227  "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1228  }
1229  else {
1230 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1231  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1232 #endif
1233  int rank = orthman_->normalizeMat(*X,Teuchos::null,BX);
1234  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
1235  "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1236  }
1237 
1238  // put data in AX
1239  if (isSkinny_) {
1240  AX_ = Z_;
1241  }
1242  {
1243 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1244  Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1245 #endif
1246  OPT::Apply(*AOp_,*X,*AX_);
1247  counterAOp_ += blockSize_;
1248  }
1249 
1250  } // end if (newstate.X != Teuchos::null)
1251 
1252 
1253  // set up Ritz values
1254  if (newstate.T != Teuchos::null) {
1255  TEUCHOS_TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_,
1256  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values.");
1257  for (int i=0; i<blockSize_; i++) {
1258  theta_[i] = (*newstate.T)[i];
1259  }
1260  }
1261  else {
1262  // get ritz vecs/vals
1263  Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_),
1264  BB(blockSize_,blockSize_),
1265  S(blockSize_,blockSize_);
1266  // project A
1267  {
1268 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1269  Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1270 #endif
1271  MVT::MvTransMv(ONE,*X,*AX_,AA);
1272  if (hasBOp_) {
1273  MVT::MvTransMv(ONE,*X,*BX,BB);
1274  }
1275  }
1276  nevLocal_ = blockSize_;
1277 
1278  // solve the projected problem
1279  // project B
1280  int ret;
1281  if (hasBOp_) {
1282 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1283  Teuchos::TimeMonitor lcltimer( *timerDS_ );
1284 #endif
1285  ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
1286  }
1287  else {
1288 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1289  Teuchos::TimeMonitor lcltimer( *timerDS_ );
1290 #endif
1291  ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10);
1292  }
1293  TEUCHOS_TEST_FOR_EXCEPTION(ret != 0,RTRInitFailure,
1294  "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret);
1295  TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,RTRInitFailure,"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis.");
1296 
1297  // We only have blockSize_ ritz pairs, ergo we do not need to select.
1298  // However, we still require them to be ordered correctly
1299  {
1300 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1301  Teuchos::TimeMonitor lcltimer( *timerSort_ );
1302 #endif
1303 
1304  std::vector<int> order(blockSize_);
1305  //
1306  // sort the first blockSize_ values in theta_
1307  sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_); // don't catch exception
1308  //
1309  // apply the same ordering to the primitive ritz vectors
1310  Utils::permuteVectors(order,S);
1311  }
1312 
1313  // compute ritz residual norms
1314  {
1315  Teuchos::BLAS<int,ScalarType> blas;
1316  Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_);
1317  // RR = AA*S - BB*S*diag(theta)
1318  int info;
1319  if (hasBOp_) {
1320  info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO);
1321  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1322  }
1323  else {
1324  RR.assign(S);
1325  }
1326  for (int i=0; i<blockSize_; i++) {
1327  blas.SCAL(blockSize_,theta_[i],RR[i],1);
1328  }
1329  info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE);
1330  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1331  for (int i=0; i<blockSize_; i++) {
1332  ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1);
1333  }
1334  }
1335 
1336  // update the solution
1337  // use R as local work space
1338  // Z may already be in use as work space (holding AX if isSkinny)
1339  {
1340 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1341  Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1342 #endif
1343  // X <- X*S
1344  MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ );
1345  MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X );
1346  // AX <- AX*S
1347  MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ );
1348  MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ );
1349  if (hasBOp_) {
1350  // BX <- BX*S
1351  MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ );
1352  MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX );
1353  }
1354  }
1355  }
1356 
1357  // done modifying X,BX; clear temp views and setup const views
1358  X = Teuchos::null;
1359  BX = Teuchos::null;
1360  {
1361  std::vector<int> ind(blockSize_);
1362  for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1363  this->X_ = MVT::CloneView(static_cast<const MV&>(*V_),ind);
1364  if (hasBOp_) {
1365  this->BX_ = MVT::CloneView(static_cast<const MV&>(*BV_),ind);
1366  }
1367  else {
1368  this->BX_ = this->X_;
1369  }
1370  }
1371 
1372 
1373  // get objective function value
1374  fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
1375 
1376  // set up R
1377  if (newstate.R != Teuchos::null) {
1378  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
1379  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
1380  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
1381  std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
1382  MVT::SetBlock(*newstate.R,bsind,*R_);
1383  }
1384  else {
1385 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1386  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1387 #endif
1388  // form R <- AX - BX*T
1389  MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
1390  Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1391  T.putScalar(ZERO);
1392  for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
1393  MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_);
1394  }
1395 
1396  // R has been updated; mark the norms as out-of-date
1397  Rnorms_current_ = false;
1398  R2norms_current_ = false;
1399 
1400  // if isSkinny, then AX currently points to Z for temp storage
1401  // set AX back to null
1402  if (isSkinny_) {
1403  AX_ = Teuchos::null;
1404  }
1405 
1406  // finally, we are initialized
1407  initialized_ = true;
1408 
1409  if (om_->isVerbosity( Debug ) ) {
1410  // Check almost everything here
1411  CheckList chk;
1412  chk.checkX = true;
1413  chk.checkAX = true;
1414  chk.checkBX = true;
1415  chk.checkR = true;
1416  chk.checkQ = true;
1417  om_->print( Debug, accuracyCheck(chk, "after initialize()") );
1418  }
1419  }
1420 
1421  template <class ScalarType, class MV, class OP>
1423  {
1425  initialize(empty);
1426  }
1427 
1428 
1429 
1430 
1432  // compute/return residual B-norms
1433  template <class ScalarType, class MV, class OP>
1434  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1436  if (Rnorms_current_ == false) {
1437  // Update the residual norms
1438  orthman_->norm(*R_,Rnorms_);
1439  Rnorms_current_ = true;
1440  }
1441  return Rnorms_;
1442  }
1443 
1444 
1446  // compute/return residual 2-norms
1447  template <class ScalarType, class MV, class OP>
1448  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1450  if (R2norms_current_ == false) {
1451  // Update the residual 2-norms
1452  MVT::MvNorm(*R_,R2norms_);
1453  R2norms_current_ = true;
1454  }
1455  return R2norms_;
1456  }
1457 
1458 
1459 
1460 
1462  // Check accuracy, orthogonality, and other debugging stuff
1463  //
1464  // bools specify which tests we want to run (instead of running more than we actually care about)
1465  //
1466  // we don't bother checking the following because they are computed explicitly:
1467  // AH == A*H
1468  //
1469  //
1470  // checkX : X orthonormal
1471  // orthogonal to auxvecs
1472  // checkAX : check AX == A*X
1473  // checkBX : check BX == B*X
1474  // checkEta : check that Eta'*B*X == 0
1475  // orthogonal to auxvecs
1476  // checkAEta : check that AEta = A*Eta
1477  // checkBEta : check that BEta = B*Eta
1478  // checkR : check R orthogonal to X
1479  // checkBR : check R B-orthogonal to X, valid in and immediately after solveTRSubproblem
1480  // checkQ : check that auxiliary vectors are actually orthonormal
1481  //
1482  // TODO:
1483  // add checkTheta
1484  //
1485  template <class ScalarType, class MV, class OP>
1486  std::string RTRBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1487  {
1488  using std::setprecision;
1489  using std::scientific;
1490  using std::setw;
1491  using std::endl;
1492  std::stringstream os;
1493  MagnitudeType tmp;
1494 
1495  os << " Debugging checks: " << where << endl;
1496 
1497  // X and friends
1498  if (chk.checkX && initialized_) {
1499  tmp = orthman_->orthonormError(*X_);
1500  os << " >> Error in X^H B X == I : " << scientific << setprecision(10) << tmp << endl;
1501  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1502  tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
1503  os << " >> Error in X^H B Q[" << i << "] == 0 : " << scientific << setprecision(10) << tmp << endl;
1504  }
1505  }
1506  if (chk.checkAX && !isSkinny_ && initialized_) {
1507  tmp = Utils::errorEquality(*X_, *AX_, AOp_);
1508  os << " >> Error in AX == A*X : " << scientific << setprecision(10) << tmp << endl;
1509  }
1510  if (chk.checkBX && hasBOp_ && initialized_) {
1511  Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_);
1512  OPT::Apply(*BOp_,*this->X_,*BX);
1513  tmp = Utils::errorEquality(*BX, *BX_);
1514  os << " >> Error in BX == B*X : " << scientific << setprecision(10) << tmp << endl;
1515  }
1516 
1517  // Eta and friends
1518  if (chk.checkEta && initialized_) {
1519  tmp = orthman_->orthogError(*X_,*eta_);
1520  os << " >> Error in X^H B Eta == 0 : " << scientific << setprecision(10) << tmp << endl;
1521  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1522  tmp = orthman_->orthogError(*eta_,*auxVecs_[i]);
1523  os << " >> Error in Eta^H B Q[" << i << "]==0 : " << scientific << setprecision(10) << tmp << endl;
1524  }
1525  }
1526  if (chk.checkAEta && !isSkinny_ && initialized_) {
1527  tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_);
1528  os << " >> Error in AEta == A*Eta : " << scientific << setprecision(10) << tmp << endl;
1529  }
1530  if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
1531  tmp = Utils::errorEquality(*eta_, *Beta_, BOp_);
1532  os << " >> Error in BEta == B*Eta : " << scientific << setprecision(10) << tmp << endl;
1533  }
1534 
1535  // R: this is not B-orthogonality, but standard euclidean orthogonality
1536  if (chk.checkR && initialized_) {
1537  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1538  MVT::MvTransMv(ONE,*X_,*R_,xTx);
1539  tmp = xTx.normFrobenius();
1540  os << " >> Error in X^H R == 0 : " << scientific << setprecision(10) << tmp << endl;
1541  }
1542 
1543  // BR: this is B-orthogonality: this is only valid inside and immediately after solveTRSubproblem
1544  // only check if B != I, otherwise, it is equivalent to the above test
1545  if (chk.checkBR && hasBOp_ && initialized_) {
1546  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1547  MVT::MvTransMv(ONE,*BX_,*R_,xTx);
1548  tmp = xTx.normFrobenius();
1549  os << " >> Error in X^H B R == 0 : " << scientific << setprecision(10) << tmp << endl;
1550  }
1551 
1552  // Z: Z is preconditioned R, should be on tangent plane
1553  if (chk.checkZ && initialized_) {
1554  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1555  MVT::MvTransMv(ONE,*BX_,*Z_,xTx);
1556  tmp = xTx.normFrobenius();
1557  os << " >> Error in X^H B Z == 0 : " << scientific << setprecision(10) << tmp << endl;
1558  }
1559 
1560  // Q
1561  if (chk.checkQ) {
1562  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1563  tmp = orthman_->orthonormError(*auxVecs_[i]);
1564  os << " >> Error in Q[" << i << "]^H B Q[" << i << "]==I: " << scientific << setprecision(10) << tmp << endl;
1565  for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1566  tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1567  os << " >> Error in Q[" << i << "]^H B Q[" << j << "]==0: " << scientific << setprecision(10) << tmp << endl;
1568  }
1569  }
1570  }
1571  os << endl;
1572  return os.str();
1573  }
1574 
1575 
1577  // Print the current status of the solver
1578  template <class ScalarType, class MV, class OP>
1579  void
1581  {
1582  using std::setprecision;
1583  using std::scientific;
1584  using std::setw;
1585  using std::endl;
1586 
1587  os <<endl;
1588  os <<"================================================================================" << endl;
1589  os << endl;
1590  os <<" RTRBase Solver Status" << endl;
1591  os << endl;
1592  os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1593  os <<"The number of iterations performed is " << iter_ << endl;
1594  os <<"The current block size is " << blockSize_ << endl;
1595  os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1596  os <<"The number of operations A*x is " << counterAOp_ << endl;
1597  os <<"The number of operations B*x is " << counterBOp_ << endl;
1598  os <<"The number of operations Prec*x is " << counterPrec_ << endl;
1599  os <<"The most recent rho was " << scientific << setprecision(10) << rho_ << endl;
1600  os <<"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl;
1601 
1602  if (initialized_) {
1603  os << endl;
1604  os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
1605  os << setw(20) << "Eigenvalue"
1606  << setw(20) << "Residual(B)"
1607  << setw(20) << "Residual(2)"
1608  << endl;
1609  os <<"--------------------------------------------------------------------------------"<<endl;
1610  for (int i=0; i<blockSize_; i++) {
1611  os << scientific << setprecision(10) << setw(20) << theta_[i];
1612  if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i];
1613  else os << scientific << setprecision(10) << setw(20) << "not current";
1614  if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i];
1615  else os << scientific << setprecision(10) << setw(20) << "not current";
1616  os << endl;
1617  }
1618  }
1619  os <<"================================================================================" << endl;
1620  os << endl;
1621  }
1622 
1623 
1625  // Inner product 1
1626  template <class ScalarType, class MV, class OP>
1627  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1628  RTRBase<ScalarType,MV,OP>::ginner(const MV &xi) const
1629  {
1630  std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi));
1631  MVT::MvNorm(xi,d);
1632  MagnitudeType ret = 0;
1633  for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1634  ret += (*i)*(*i);
1635  }
1636  return ret;
1637  }
1638 
1639 
1641  // Inner product 2
1642  template <class ScalarType, class MV, class OP>
1643  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1644  RTRBase<ScalarType,MV,OP>::ginner(const MV &xi, const MV &zeta) const
1645  {
1646  std::vector<ScalarType> d(MVT::GetNumberVecs(xi));
1647  MVT::MvDot(xi,zeta,d);
1648  return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
1649  }
1650 
1651 
1653  // Inner product 1 without trace accumulation
1654  template <class ScalarType, class MV, class OP>
1655  void RTRBase<ScalarType,MV,OP>::ginnersep(
1656  const MV &xi,
1657  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const
1658  {
1659  MVT::MvNorm(xi,d);
1660  for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1661  (*i) = (*i)*(*i);
1662  }
1663  }
1664 
1665 
1667  // Inner product 2 without trace accumulation
1668  template <class ScalarType, class MV, class OP>
1669  void RTRBase<ScalarType,MV,OP>::ginnersep(
1670  const MV &xi, const MV &zeta,
1671  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const
1672  {
1673  std::vector<ScalarType> dC(MVT::GetNumberVecs(xi));
1674  MVT::MvDot(xi,zeta,dC);
1675  vecMTiter iR=d.begin();
1676  vecSTiter iS=dC.begin();
1677  for (; iR != d.end(); iR++, iS++) {
1678  (*iR) = SCT::real(*iS);
1679  }
1680  }
1681 
1682  template <class ScalarType, class MV, class OP>
1683  Teuchos::Array<Teuchos::RCP<const MV> > RTRBase<ScalarType,MV,OP>::getAuxVecs() const {
1684  return auxVecs_;
1685  }
1686 
1687  template <class ScalarType, class MV, class OP>
1689  return(blockSize_);
1690  }
1691 
1692  template <class ScalarType, class MV, class OP>
1694  return(*problem_);
1695  }
1696 
1697  template <class ScalarType, class MV, class OP>
1699  return blockSize_;
1700  }
1701 
1702  template <class ScalarType, class MV, class OP>
1704  {
1705  if (!initialized_) return 0;
1706  return nevLocal_;
1707  }
1708 
1709  template <class ScalarType, class MV, class OP>
1710  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1712  {
1713  std::vector<MagnitudeType> ret = ritz2norms_;
1714  ret.resize(nevLocal_);
1715  return ret;
1716  }
1717 
1718  template <class ScalarType, class MV, class OP>
1719  std::vector<Value<ScalarType> >
1721  {
1722  std::vector<Value<ScalarType> > ret(nevLocal_);
1723  for (int i=0; i<nevLocal_; i++) {
1724  ret[i].realpart = theta_[i];
1725  ret[i].imagpart = ZERO;
1726  }
1727  return ret;
1728  }
1729 
1730  template <class ScalarType, class MV, class OP>
1731  Teuchos::RCP<const MV>
1733  {
1734  return X_;
1735  }
1736 
1737  template <class ScalarType, class MV, class OP>
1739  {
1740  iter_=0;
1741  }
1742 
1743  template <class ScalarType, class MV, class OP>
1745  {
1746  return iter_;
1747  }
1748 
1749  template <class ScalarType, class MV, class OP>
1751  {
1753  state.X = X_;
1754  state.AX = AX_;
1755  if (hasBOp_) {
1756  state.BX = BX_;
1757  }
1758  else {
1759  state.BX = Teuchos::null;
1760  }
1761  state.rho = rho_;
1762  state.R = R_;
1763  state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_));
1764  return state;
1765  }
1766 
1767  template <class ScalarType, class MV, class OP>
1769  {
1770  return initialized_;
1771  }
1772 
1773  template <class ScalarType, class MV, class OP>
1775  {
1776  std::vector<int> ret(nevLocal_,0);
1777  return ret;
1778  }
1779 
1780 
1781 } // end Anasazi namespace
1782 
1783 #endif // ANASAZI_RTRBASE_HPP
RTROrthoFailure is thrown when an orthogonalization attempt fails.
void resetNumIters()
Reset the iteration count.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values from the previous iteration.
Structure to contain pointers to RTR state variables.
virtual void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
Teuchos::RCP< const MV > AX
The image of the current eigenvectors under A, or Teuchos::null is we implement a skinny solver...
This class defines the interface required by an eigensolver and status test class to compute solution...
Teuchos::RCP< const MV > BX
The image of the current eigenvectors under B, or Teuchos::null if B was not specified.
Declaration of basic traits for the multivector type.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors from the previous iteration.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the Ritz residuals.
RTRState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
This class is an abstract base class for Implicit Riemannian Trust-Region based eigensolvers. The class provides the interfaces shared by the IRTR solvers (e.g., getState() and initialize()) as well as the shared implementations (e.g., inner products).
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
virtual void iterate()=0
This method performs RTR iterations until the status test indicates the need to stop or an error occu...
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For RTR, this always returns getBlockSiz...
Teuchos::RCP< const MV > X
The current eigenvectors.
An exception class parent to all Anasazi exceptions.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
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...
virtual ~RTRBase()
RTRBase destructor
Anasazi&#39;s templated, static class providing utilities for the solvers.
Output managers remove the need for the eigensolver to know any information about the required output...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this eigenproblem.
std::vector< int > getRitzIndex()
Get the index used for extracting Ritz vectors from getRitzVectors().
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
RTRInitFailure is thrown when the RTR solver is unable to generate an initial iterate in the RTRBase:...
Teuchos::ScalarTraits< ScalarType >::magnitudeType rho
The current rho value. This is only valid if the debugging level of verbosity is enabled.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Teuchos::RCP< const MV > R
The current residual vectors.
void setAuxVecs(const Teuchos::Array< Teuchos::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...
Types and exceptions used within Anasazi solvers and interfaces.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi&#39;s solvers.
RTRRitzFailure is thrown when the RTR solver is unable to continue a call to RTRBase::iterate() due t...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
RTRBase(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< GenOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params, const std::string &solverLabel, bool skinnySolver)
RTRBase constructor with eigenproblem, solver utilities, and parameter list of solver options...
Class which provides internal utilities for the Anasazi solvers.