Anasazi  Version of the Day
AnasaziBlockDavidson.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_BLOCK_DAVIDSON_HPP
34 #define ANASAZI_BLOCK_DAVIDSON_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 
67 namespace Anasazi {
68 
70 
71 
76  template <class ScalarType, class MV>
82  int curDim;
87  Teuchos::RCP<const MV> V;
89  Teuchos::RCP<const MV> X;
91  Teuchos::RCP<const MV> KX;
93  Teuchos::RCP<const MV> MX;
95  Teuchos::RCP<const MV> R;
100  Teuchos::RCP<const MV> H;
102  Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
108  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > KK;
109  BlockDavidsonState() : curDim(0), V(Teuchos::null),
110  X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null),
111  R(Teuchos::null), H(Teuchos::null),
112  T(Teuchos::null), KK(Teuchos::null) {}
113  };
114 
116 
118 
119 
132  class BlockDavidsonInitFailure : public AnasaziError {public:
133  BlockDavidsonInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
134  {}};
135 
143  class BlockDavidsonOrthoFailure : public AnasaziError {public:
144  BlockDavidsonOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
145  {}};
146 
148 
149 
150  template <class ScalarType, class MV, class OP>
151  class BlockDavidson : public Eigensolver<ScalarType,MV,OP> {
152  public:
154 
155 
163  BlockDavidson( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
164  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
165  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
166  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
167  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
168  Teuchos::ParameterList &params
169  );
170 
172  virtual ~BlockDavidson();
174 
175 
177 
178 
202  void iterate();
203 
237  void initialize(BlockDavidsonState<ScalarType,MV>& newstate);
238 
242  void initialize();
243 
259  bool isInitialized() const;
260 
273  BlockDavidsonState<ScalarType,MV> getState() const;
274 
276 
277 
279 
280 
282  int getNumIters() const;
283 
285  void resetNumIters();
286 
294  Teuchos::RCP<const MV> getRitzVectors();
295 
301  std::vector<Value<ScalarType> > getRitzValues();
302 
303 
311  std::vector<int> getRitzIndex();
312 
313 
319  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
320 
321 
327  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
328 
329 
337  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
338 
344  int getCurSubspaceDim() const;
345 
347  int getMaxSubspaceDim() const;
348 
350 
351 
353 
354 
356  void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
357 
359  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
360 
362  const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
363 
373  void setBlockSize(int blockSize);
374 
376  int getBlockSize() const;
377 
390  void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
391 
393  Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
394 
396 
398 
399 
409  void setSize(int blockSize, int numBlocks);
410 
412 
414 
415 
417  void currentStatus(std::ostream &os);
418 
420 
421  private:
422  //
423  // Convenience typedefs
424  //
428  typedef Teuchos::ScalarTraits<ScalarType> SCT;
429  typedef typename SCT::magnitudeType MagnitudeType;
430  const MagnitudeType ONE;
431  const MagnitudeType ZERO;
432  const MagnitudeType NANVAL;
433  //
434  // Internal structs
435  //
436  struct CheckList {
437  bool checkV;
438  bool checkX, checkMX, checkKX;
439  bool checkH, checkMH, checkKH;
440  bool checkR, checkQ;
441  bool checkKK;
442  CheckList() : checkV(false),
443  checkX(false),checkMX(false),checkKX(false),
444  checkH(false),checkMH(false),checkKH(false),
445  checkR(false),checkQ(false),checkKK(false) {};
446  };
447  //
448  // Internal methods
449  //
450  std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
451  //
452  // Classes inputed through constructor that define the eigenproblem to be solved.
453  //
454  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
455  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
456  const Teuchos::RCP<OutputManager<ScalarType> > om_;
457  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
458  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
459  //
460  // Information obtained from the eigenproblem
461  //
462  Teuchos::RCP<const OP> Op_;
463  Teuchos::RCP<const OP> MOp_;
464  Teuchos::RCP<const OP> Prec_;
465  bool hasM_;
466  //
467  // Internal timers
468  //
469  Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
470  timerSortEval_, timerDS_,
471  timerLocal_, timerCompRes_,
472  timerOrtho_, timerInit_;
473  //
474  // Counters
475  //
476  int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
477 
478  //
479  // Algorithmic parameters.
480  //
481  // blockSize_ is the solver block size; it controls the number of eigenvectors that
482  // we compute, the number of residual vectors that we compute, and therefore the number
483  // of vectors added to the basis on each iteration.
484  int blockSize_;
485  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
486  int numBlocks_;
487 
488  //
489  // Current solver state
490  //
491  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
492  // is capable of running; _initialize is controlled by the initialize() member method
493  // For the implications of the state of initialized_, please see documentation for initialize()
494  bool initialized_;
495  //
496  // curDim_ reflects how much of the current basis is valid
497  // NOTE: 0 <= curDim_ <= blockSize_*numBlocks_
498  // this also tells us how many of the values in theta_ are valid Ritz values
499  int curDim_;
500  //
501  // State Multivecs
502  // H_,KH_,MH_ will not own any storage
503  // H_ will occasionally point at the current block of vectors in the basis V_
504  // MH_,KH_ will occasionally point at MX_,KX_ when they are used as temporary storage
505  Teuchos::RCP<MV> X_, KX_, MX_, R_,
506  H_, KH_, MH_,
507  V_;
508  //
509  // Projected matrices
510  //
511  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_;
512  //
513  // auxiliary vectors
514  Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
515  int numAuxVecs_;
516  //
517  // Number of iterations that have been performed.
518  int iter_;
519  //
520  // Current eigenvalues, residual norms
521  std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
522  //
523  // are the residual norms current with the residual?
524  bool Rnorms_current_, R2norms_current_;
525 
526  };
527 
530  //
531  // Implementations
532  //
535 
536 
538  // Constructor
539  template <class ScalarType, class MV, class OP>
541  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
542  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
543  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
544  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
545  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
546  Teuchos::ParameterList &params
547  ) :
548  ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
549  ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
550  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
551  // problem, tools
552  problem_(problem),
553  sm_(sorter),
554  om_(printer),
555  tester_(tester),
556  orthman_(ortho),
557  // timers, counters
558 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
559  timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation Op*x")),
560  timerMOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation M*x")),
561  timerPrec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation Prec*x")),
562  timerSortEval_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Sorting eigenvalues")),
563  timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Direct solve")),
564  timerLocal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Local update")),
565  timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Computing residuals")),
566  timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Orthogonalization")),
567  timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Initialization")),
568 #endif
569  count_ApplyOp_(0),
570  count_ApplyM_(0),
571  count_ApplyPrec_(0),
572  // internal data
573  blockSize_(0),
574  numBlocks_(0),
575  initialized_(false),
576  curDim_(0),
577  auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
578  numAuxVecs_(0),
579  iter_(0),
580  Rnorms_current_(false),
581  R2norms_current_(false)
582  {
583  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
584  "Anasazi::BlockDavidson::constructor: user passed null problem pointer.");
585  TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
586  "Anasazi::BlockDavidson::constructor: user passed null sort manager pointer.");
587  TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
588  "Anasazi::BlockDavidson::constructor: user passed null output manager pointer.");
589  TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
590  "Anasazi::BlockDavidson::constructor: user passed null status test pointer.");
591  TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
592  "Anasazi::BlockDavidson::constructor: user passed null orthogonalization manager pointer.");
593  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
594  "Anasazi::BlockDavidson::constructor: problem is not set.");
595  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
596  "Anasazi::BlockDavidson::constructor: problem is not hermitian.");
597 
598  // get the problem operators
599  Op_ = problem_->getOperator();
600  TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
601  "Anasazi::BlockDavidson::constructor: problem provides no operator.");
602  MOp_ = problem_->getM();
603  Prec_ = problem_->getPrec();
604  hasM_ = (MOp_ != Teuchos::null);
605 
606  // set the block size and allocate data
607  int bs = params.get("Block Size", problem_->getNEV());
608  int nb = params.get("Num Blocks", 2);
609  setSize(bs,nb);
610  }
611 
612 
614  // Destructor
615  template <class ScalarType, class MV, class OP>
617 
618 
620  // Set the block size
621  // This simply calls setSize(), modifying the block size while retaining the number of blocks.
622  template <class ScalarType, class MV, class OP>
624  {
625  setSize(blockSize,numBlocks_);
626  }
627 
628 
630  // Return the current auxiliary vectors
631  template <class ScalarType, class MV, class OP>
632  Teuchos::Array<Teuchos::RCP<const MV> > BlockDavidson<ScalarType,MV,OP>::getAuxVecs() const {
633  return auxVecs_;
634  }
635 
636 
637 
639  // return the current block size
640  template <class ScalarType, class MV, class OP>
642  return(blockSize_);
643  }
644 
645 
647  // return eigenproblem
648  template <class ScalarType, class MV, class OP>
650  return(*problem_);
651  }
652 
653 
655  // return max subspace dim
656  template <class ScalarType, class MV, class OP>
658  return blockSize_*numBlocks_;
659  }
660 
662  // return current subspace dim
663  template <class ScalarType, class MV, class OP>
665  if (!initialized_) return 0;
666  return curDim_;
667  }
668 
669 
671  // return ritz residual 2-norms
672  template <class ScalarType, class MV, class OP>
673  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
675  return this->getRes2Norms();
676  }
677 
678 
680  // return ritz index
681  template <class ScalarType, class MV, class OP>
683  std::vector<int> ret(curDim_,0);
684  return ret;
685  }
686 
687 
689  // return ritz values
690  template <class ScalarType, class MV, class OP>
691  std::vector<Value<ScalarType> > BlockDavidson<ScalarType,MV,OP>::getRitzValues() {
692  std::vector<Value<ScalarType> > ret(curDim_);
693  for (int i=0; i<curDim_; ++i) {
694  ret[i].realpart = theta_[i];
695  ret[i].imagpart = ZERO;
696  }
697  return ret;
698  }
699 
700 
702  // return pointer to ritz vectors
703  template <class ScalarType, class MV, class OP>
705  return X_;
706  }
707 
708 
710  // reset number of iterations
711  template <class ScalarType, class MV, class OP>
713  iter_=0;
714  }
715 
716 
718  // return number of iterations
719  template <class ScalarType, class MV, class OP>
721  return(iter_);
722  }
723 
724 
726  // return state pointers
727  template <class ScalarType, class MV, class OP>
730  state.curDim = curDim_;
731  state.V = V_;
732  state.X = X_;
733  state.KX = KX_;
734  if (hasM_) {
735  state.MX = MX_;
736  }
737  else {
738  state.MX = Teuchos::null;
739  }
740  state.R = R_;
741  state.H = H_;
742  state.KK = KK_;
743  if (curDim_ > 0) {
744  state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
745  } else {
746  state.T = Teuchos::rcp(new std::vector<MagnitudeType>(0));
747  }
748  return state;
749  }
750 
751 
753  // Return initialized state
754  template <class ScalarType, class MV, class OP>
755  bool BlockDavidson<ScalarType,MV,OP>::isInitialized() const { return initialized_; }
756 
757 
759  // Set the block size and make necessary adjustments.
760  template <class ScalarType, class MV, class OP>
761  void BlockDavidson<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
762  {
763  // time spent here counts towards timerInit_
764 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
765  Teuchos::TimeMonitor initimer( *timerInit_ );
766 #endif
767 
768  // This routine only allocates space; it doesn't not perform any computation
769  // any change in size will invalidate the state of the solver.
770 
771  TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): blocksize must be strictly positive.");
772  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 2, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): numblocks must be greater than one.");
773  if (blockSize == blockSize_ && numBlocks == numBlocks_) {
774  // do nothing
775  return;
776  }
777 
778  blockSize_ = blockSize;
779  numBlocks_ = numBlocks;
780 
781  Teuchos::RCP<const MV> tmp;
782  // grab some Multivector to Clone
783  // in practice, getInitVec() should always provide this, but it is possible to use a
784  // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
785  // in case of that strange scenario, we will try to Clone from X_ first, then resort to getInitVec()
786  if (X_ != Teuchos::null) { // this is equivalent to blockSize_ > 0
787  tmp = X_;
788  }
789  else {
790  tmp = problem_->getInitVec();
791  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
792  "Anasazi::BlockDavidson::setSize(): eigenproblem did not specify initial vectors to clone from.");
793  }
794 
795  TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*static_cast<ptrdiff_t>(numBlocks) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
796  "Anasazi::BlockDavidson::setSize(): max subspace dimension and auxilliary subspace too large.");
797 
798 
800  // blockSize dependent
801  //
802  // grow/allocate vectors
803  Rnorms_.resize(blockSize_,NANVAL);
804  R2norms_.resize(blockSize_,NANVAL);
805  //
806  // clone multivectors off of tmp
807  //
808  // free current allocation first, to make room for new allocation
809  X_ = Teuchos::null;
810  KX_ = Teuchos::null;
811  MX_ = Teuchos::null;
812  R_ = Teuchos::null;
813  V_ = Teuchos::null;
814 
815  om_->print(Debug," >> Allocating X_\n");
816  X_ = MVT::Clone(*tmp,blockSize_);
817  om_->print(Debug," >> Allocating KX_\n");
818  KX_ = MVT::Clone(*tmp,blockSize_);
819  if (hasM_) {
820  om_->print(Debug," >> Allocating MX_\n");
821  MX_ = MVT::Clone(*tmp,blockSize_);
822  }
823  else {
824  MX_ = X_;
825  }
826  om_->print(Debug," >> Allocating R_\n");
827  R_ = MVT::Clone(*tmp,blockSize_);
828 
830  // blockSize*numBlocks dependent
831  //
832  int newsd = blockSize_*numBlocks_;
833  theta_.resize(blockSize_*numBlocks_,NANVAL);
834  om_->print(Debug," >> Allocating V_\n");
835  V_ = MVT::Clone(*tmp,newsd);
836  KK_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
837 
838  om_->print(Debug," >> done allocating.\n");
839 
840  initialized_ = false;
841  curDim_ = 0;
842  }
843 
844 
846  // Set the auxiliary vectors
847  template <class ScalarType, class MV, class OP>
848  void BlockDavidson<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
849  typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
850 
851  // set new auxiliary vectors
852  auxVecs_ = auxvecs;
853  numAuxVecs_ = 0;
854  for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
855  numAuxVecs_ += MVT::GetNumberVecs(**i);
856  }
857 
858  // If the solver has been initialized, V is not necessarily orthogonal to new auxiliary vectors
859  if (numAuxVecs_ > 0 && initialized_) {
860  initialized_ = false;
861  }
862 
863  if (om_->isVerbosity( Debug ) ) {
864  CheckList chk;
865  chk.checkQ = true;
866  om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
867  }
868  }
869 
870 
872  /* Initialize the state of the solver
873  *
874  * POST-CONDITIONS:
875  *
876  * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
877  * theta_ contains Ritz w.r.t. V_(1:curDim_)
878  * X is Ritz vectors w.r.t. V_(1:curDim_)
879  * KX = Op*X
880  * MX = M*X if hasM_
881  * R = KX - MX*diag(theta_)
882  *
883  */
884  template <class ScalarType, class MV, class OP>
886  {
887  // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
888  // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
889 
890 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
891  Teuchos::TimeMonitor inittimer( *timerInit_ );
892 #endif
893 
894  std::vector<int> bsind(blockSize_);
895  for (int i=0; i<blockSize_; ++i) bsind[i] = i;
896 
897  Teuchos::BLAS<int,ScalarType> blas;
898 
899  // in BlockDavidson, V is primary
900  // the order of dependence follows like so.
901  // --init-> V,KK
902  // --ritz analysis-> theta,X
903  // --op apply-> KX,MX
904  // --compute-> R
905  //
906  // if the user specifies all data for a level, we will accept it.
907  // otherwise, we will generate the whole level, and all subsequent levels.
908  //
909  // the data members are ordered based on dependence, and the levels are
910  // partitioned according to the amount of work required to produce the
911  // items in a level.
912  //
913  // inconsistent multivectors widths and lengths will not be tolerated, and
914  // will be treated with exceptions.
915  //
916  // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to
917  // multivectors in the solver, the copy will not be affected.
918 
919  // set up V and KK: get them from newstate if user specified them
920  // otherwise, set them manually
921  Teuchos::RCP<MV> lclV;
922  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK;
923 
924  if (newstate.V != Teuchos::null && newstate.KK != Teuchos::null) {
925  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
926  "Anasazi::BlockDavidson::initialize(newstate): Vector length of V not correct." );
927  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
928  "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be at least blockSize().");
929  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
930  "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
931  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
932  "Anasazi::BlockDavidson::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
933 
934  curDim_ = newstate.curDim;
935  // pick an integral amount
936  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
937 
938  TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
939  "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
940 
941  // check size of KK
942  TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
943  "Anasazi::BlockDavidson::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
944 
945  // put data in V
946  std::vector<int> nevind(curDim_);
947  for (int i=0; i<curDim_; ++i) nevind[i] = i;
948  if (newstate.V != V_) {
949  if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
950  newstate.V = MVT::CloneView(*newstate.V,nevind);
951  }
952  MVT::SetBlock(*newstate.V,nevind,*V_);
953  }
954  lclV = MVT::CloneViewNonConst(*V_,nevind);
955 
956  // put data into KK_
957  lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
958  if (newstate.KK != KK_) {
959  if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
960  newstate.KK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
961  }
962  lclKK->assign(*newstate.KK);
963  }
964  //
965  // make lclKK Hermitian in memory (copy the upper half to the lower half)
966  for (int j=0; j<curDim_-1; ++j) {
967  for (int i=j+1; i<curDim_; ++i) {
968  (*lclKK)(i,j) = SCT::conjugate((*lclKK)(j,i));
969  }
970  }
971  }
972  else {
973  // user did not specify one of V or KK
974  // get vectors from problem or generate something, projectAndNormalize
975  Teuchos::RCP<const MV> ivec = problem_->getInitVec();
976  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
977  "Anasazi::BlockDavdison::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
978  // clear newstate so we won't use any data from it below
979  newstate.X = Teuchos::null;
980  newstate.MX = Teuchos::null;
981  newstate.KX = Teuchos::null;
982  newstate.R = Teuchos::null;
983  newstate.H = Teuchos::null;
984  newstate.T = Teuchos::null;
985  newstate.KK = Teuchos::null;
986  newstate.V = Teuchos::null;
987  newstate.curDim = 0;
988 
989  curDim_ = MVT::GetNumberVecs(*ivec);
990  // pick the largest multiple of blockSize_
991  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
992  if (curDim_ > blockSize_*numBlocks_) {
993  // user specified too many vectors... truncate
994  // this produces a full subspace, but that is okay
995  curDim_ = blockSize_*numBlocks_;
996  }
997  bool userand = false;
998  if (curDim_ == 0) {
999  // we need at least blockSize_ vectors
1000  // use a random multivec: ignore everything from InitVec
1001  userand = true;
1002  curDim_ = blockSize_;
1003  }
1004 
1005  // get pointers into V,KV,MV
1006  // tmpVecs will be used below for M*V and K*V (not simultaneously)
1007  // lclV has curDim vectors
1008  // if there is space for lclV and tmpVecs in V_, point tmpVecs into V_
1009  // otherwise, we must allocate space for these products
1010  //
1011  // get pointer to first curDim vector in V_
1012  std::vector<int> dimind(curDim_);
1013  for (int i=0; i<curDim_; ++i) dimind[i] = i;
1014  lclV = MVT::CloneViewNonConst(*V_,dimind);
1015  if (userand) {
1016  // generate random vector data
1017  MVT::MvRandom(*lclV);
1018  }
1019  else {
1020  if (MVT::GetNumberVecs(*ivec) > curDim_) {
1021  ivec = MVT::CloneView(*ivec,dimind);
1022  }
1023  // assign ivec to first part of lclV
1024  MVT::SetBlock(*ivec,dimind,*lclV);
1025  }
1026  Teuchos::RCP<MV> tmpVecs;
1027  if (curDim_*2 <= blockSize_*numBlocks_) {
1028  // partition V_ = [lclV tmpVecs _leftover_]
1029  std::vector<int> block2(curDim_);
1030  for (int i=0; i<curDim_; ++i) block2[i] = curDim_+i;
1031  tmpVecs = MVT::CloneViewNonConst(*V_,block2);
1032  }
1033  else {
1034  // allocate space for tmpVecs
1035  tmpVecs = MVT::Clone(*V_,curDim_);
1036  }
1037 
1038  // compute M*lclV if hasM_
1039  if (hasM_) {
1040 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1041  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1042 #endif
1043  OPT::Apply(*MOp_,*lclV,*tmpVecs);
1044  count_ApplyM_ += curDim_;
1045  }
1046 
1047  // remove auxVecs from lclV and normalize it
1048  if (auxVecs_.size() > 0) {
1049 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1050  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1051 #endif
1052 
1053  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1054  int rank = orthman_->projectAndNormalizeMat(*lclV,auxVecs_,dummyC,Teuchos::null,tmpVecs);
1055  TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
1056  "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1057  }
1058  else {
1059 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1060  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1061 #endif
1062 
1063  int rank = orthman_->normalizeMat(*lclV,Teuchos::null,tmpVecs);
1064  TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
1065  "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1066  }
1067 
1068  // compute K*lclV: we are re-using tmpVecs to store the result
1069  {
1070 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1071  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1072 #endif
1073  OPT::Apply(*Op_,*lclV,*tmpVecs);
1074  count_ApplyOp_ += curDim_;
1075  }
1076 
1077  // generate KK
1078  lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1079  MVT::MvTransMv(ONE,*lclV,*tmpVecs,*lclKK);
1080 
1081  // clear tmpVecs
1082  tmpVecs = Teuchos::null;
1083  }
1084 
1085  // X,theta require Ritz analysis; if we have to generate one of these, we might as well generate both
1086  if (newstate.X != Teuchos::null && newstate.T != Teuchos::null) {
1087  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1088  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of X must be consistent with block size and length of V.");
1089  TEUCHOS_TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
1090  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of T must be consistent with dimension of V.");
1091 
1092  if (newstate.X != X_) {
1093  MVT::SetBlock(*newstate.X,bsind,*X_);
1094  }
1095 
1096  std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
1097  }
1098  else {
1099  // compute ritz vecs/vals
1100  Teuchos::SerialDenseMatrix<int,ScalarType> S(curDim_,curDim_);
1101  {
1102 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1103  Teuchos::TimeMonitor lcltimer( *timerDS_ );
1104 #endif
1105  int rank = curDim_;
1106  Utils::directSolver(curDim_, *lclKK, Teuchos::null, S, theta_, rank, 10);
1107  // we want all ritz values back
1108  TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
1109  "Anasazi::BlockDavidson::initialize(newstate): Not enough Ritz vectors to initialize algorithm.");
1110  }
1111  // sort ritz pairs
1112  {
1113 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1114  Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1115 #endif
1116 
1117  std::vector<int> order(curDim_);
1118  //
1119  // sort the first curDim_ values in theta_
1120  sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_); // don't catch exception
1121  //
1122  // apply the same ordering to the primitive ritz vectors
1123  Utils::permuteVectors(order,S);
1124  }
1125 
1126  // compute eigenvectors
1127  Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,curDim_,blockSize_);
1128  {
1129 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1130  Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1131 #endif
1132 
1133  // X <- lclV*S
1134  MVT::MvTimesMatAddMv( ONE, *lclV, S1, ZERO, *X_ );
1135  }
1136  // we generated theta,X so we don't want to use the user's KX,MX
1137  newstate.KX = Teuchos::null;
1138  newstate.MX = Teuchos::null;
1139  }
1140 
1141  // done with local pointers
1142  lclV = Teuchos::null;
1143  lclKK = Teuchos::null;
1144 
1145  // set up KX
1146  if ( newstate.KX != Teuchos::null ) {
1147  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_,
1148  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
1149  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*X_),
1150  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
1151  if (newstate.KX != KX_) {
1152  MVT::SetBlock(*newstate.KX,bsind,*KX_);
1153  }
1154  }
1155  else {
1156  // generate KX
1157  {
1158 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1159  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1160 #endif
1161  OPT::Apply(*Op_,*X_,*KX_);
1162  count_ApplyOp_ += blockSize_;
1163  }
1164  // we generated KX; we will generate R as well
1165  newstate.R = Teuchos::null;
1166  }
1167 
1168  // set up MX
1169  if (hasM_) {
1170  if ( newstate.MX != Teuchos::null ) {
1171  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_,
1172  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
1173  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*X_),
1174  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
1175  if (newstate.MX != MX_) {
1176  MVT::SetBlock(*newstate.MX,bsind,*MX_);
1177  }
1178  }
1179  else {
1180  // generate MX
1181  {
1182 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1183  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1184 #endif
1185  OPT::Apply(*MOp_,*X_,*MX_);
1186  count_ApplyOp_ += blockSize_;
1187  }
1188  // we generated MX; we will generate R as well
1189  newstate.R = Teuchos::null;
1190  }
1191  }
1192  else {
1193  // the assignment MX_==X_ would be redundant; take advantage of this opportunity to debug a little
1194  TEUCHOS_TEST_FOR_EXCEPTION(MX_ != X_, std::logic_error, "Anasazi::BlockDavidson::initialize(): solver invariant not satisfied (MX==X).");
1195  }
1196 
1197  // set up R
1198  if (newstate.R != Teuchos::null) {
1199  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
1200  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
1201  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1202  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
1203  if (newstate.R != R_) {
1204  MVT::SetBlock(*newstate.R,bsind,*R_);
1205  }
1206  }
1207  else {
1208 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1209  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1210 #endif
1211 
1212  // form R <- KX - MX*T
1213  MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1214  Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1215  T.putScalar(ZERO);
1216  for (int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
1217  MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
1218 
1219  }
1220 
1221  // R has been updated; mark the norms as out-of-date
1222  Rnorms_current_ = false;
1223  R2norms_current_ = false;
1224 
1225  // finally, we are initialized
1226  initialized_ = true;
1227 
1228  if (om_->isVerbosity( Debug ) ) {
1229  // Check almost everything here
1230  CheckList chk;
1231  chk.checkV = true;
1232  chk.checkX = true;
1233  chk.checkKX = true;
1234  chk.checkMX = true;
1235  chk.checkR = true;
1236  chk.checkQ = true;
1237  chk.checkKK = true;
1238  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
1239  }
1240 
1241  // Print information on current status
1242  if (om_->isVerbosity(Debug)) {
1243  currentStatus( om_->stream(Debug) );
1244  }
1245  else if (om_->isVerbosity(IterationDetails)) {
1246  currentStatus( om_->stream(IterationDetails) );
1247  }
1248  }
1249 
1250 
1252  // initialize the solver with default state
1253  template <class ScalarType, class MV, class OP>
1255  {
1257  initialize(empty);
1258  }
1259 
1260 
1262  // Perform BlockDavidson iterations until the StatusTest tells us to stop.
1263  template <class ScalarType, class MV, class OP>
1265  {
1266  //
1267  // Initialize solver state
1268  if (initialized_ == false) {
1269  initialize();
1270  }
1271 
1272  // as a data member, this would be redundant and require synchronization with
1273  // blockSize_ and numBlocks_; we'll use a constant here.
1274  const int searchDim = blockSize_*numBlocks_;
1275 
1276  Teuchos::BLAS<int,ScalarType> blas;
1277 
1278  //
1279  // The projected matrices are part of the state, but the eigenvectors are defined locally.
1280  // S = Local eigenvectors (size: searchDim * searchDim
1281  Teuchos::SerialDenseMatrix<int,ScalarType> S( searchDim, searchDim );
1282 
1283 
1285  // iterate until the status test tells us to stop.
1286  // also break if our basis is full
1287  while (tester_->checkStatus(this) != Passed && curDim_ < searchDim) {
1288 
1289  // Print information on current iteration
1290  if (om_->isVerbosity(Debug)) {
1291  currentStatus( om_->stream(Debug) );
1292  }
1293  else if (om_->isVerbosity(IterationDetails)) {
1294  currentStatus( om_->stream(IterationDetails) );
1295  }
1296 
1297  ++iter_;
1298 
1299  // get the current part of the basis
1300  std::vector<int> curind(blockSize_);
1301  for (int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
1302  H_ = MVT::CloneViewNonConst(*V_,curind);
1303 
1304  // Apply the preconditioner on the residuals: H <- Prec*R
1305  // H = Prec*R
1306  if (Prec_ != Teuchos::null) {
1307 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1308  Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1309 #endif
1310  OPT::Apply( *Prec_, *R_, *H_ ); // don't catch the exception
1311  count_ApplyPrec_ += blockSize_;
1312  }
1313  else {
1314  std::vector<int> bsind(blockSize_);
1315  for (int i=0; i<blockSize_; ++i) { bsind[i] = i; }
1316  MVT::SetBlock(*R_,bsind,*H_);
1317  }
1318 
1319  // Apply the mass matrix on H
1320  if (hasM_) {
1321  // use memory at MX_ for temporary storage
1322  MH_ = MX_;
1323 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1324  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1325 #endif
1326  OPT::Apply( *MOp_, *H_, *MH_); // don't catch the exception
1327  count_ApplyM_ += blockSize_;
1328  }
1329  else {
1330  MH_ = H_;
1331  }
1332 
1333  // Get a view of the previous vectors
1334  // this is used for orthogonalization and for computing V^H K H
1335  std::vector<int> prevind(curDim_);
1336  for (int i=0; i<curDim_; ++i) prevind[i] = i;
1337  Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,prevind);
1338 
1339  // Orthogonalize H against the previous vectors and the auxiliary vectors, and normalize
1340  {
1341 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1342  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1343 #endif
1344 
1345  Teuchos::Array<Teuchos::RCP<const MV> > against = auxVecs_;
1346  against.push_back(Vprev);
1347  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1348  int rank = orthman_->projectAndNormalizeMat(*H_,against,
1349  dummyC,
1350  Teuchos::null,MH_);
1351  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockDavidsonOrthoFailure,
1352  "Anasazi::BlockDavidson::iterate(): unable to compute orthonormal basis for H.");
1353  }
1354 
1355  // Apply the stiffness matrix to H
1356  {
1357  // use memory at KX_ for temporary storage
1358  KH_ = KX_;
1359 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1360  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1361 #endif
1362  OPT::Apply( *Op_, *H_, *KH_); // don't catch the exception
1363  count_ApplyOp_ += blockSize_;
1364  }
1365 
1366  if (om_->isVerbosity( Debug ) ) {
1367  CheckList chk;
1368  chk.checkH = true;
1369  chk.checkMH = true;
1370  chk.checkKH = true;
1371  om_->print( Debug, accuracyCheck(chk, ": after ortho H") );
1372  }
1373  else if (om_->isVerbosity( OrthoDetails ) ) {
1374  CheckList chk;
1375  chk.checkH = true;
1376  chk.checkMH = true;
1377  chk.checkKH = true;
1378  om_->print( OrthoDetails, accuracyCheck(chk,": after ortho H") );
1379  }
1380 
1381  // compute next part of the projected matrices
1382  // this this in two parts
1383  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > nextKK;
1384  // Vprev*K*H
1385  nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_) );
1386  MVT::MvTransMv(ONE,*Vprev,*KH_,*nextKK);
1387  // H*K*H
1388  nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
1389  MVT::MvTransMv(ONE,*H_,*KH_,*nextKK);
1390  //
1391  // make sure that KK_ is Hermitian in memory
1392  nextKK = Teuchos::null;
1393  for (int i=curDim_; i<curDim_+blockSize_; ++i) {
1394  for (int j=0; j<i; ++j) {
1395  (*KK_)(i,j) = SCT::conjugate((*KK_)(j,i));
1396  }
1397  }
1398 
1399  // V has been extended, and KK has been extended. Update basis dim and release all pointers.
1400  curDim_ += blockSize_;
1401  H_ = KH_ = MH_ = Teuchos::null;
1402  Vprev = Teuchos::null;
1403 
1404  if (om_->isVerbosity( Debug ) ) {
1405  CheckList chk;
1406  chk.checkKK = true;
1407  om_->print( Debug, accuracyCheck(chk, ": after expanding KK") );
1408  }
1409 
1410  // Get pointer to complete basis
1411  curind.resize(curDim_);
1412  for (int i=0; i<curDim_; ++i) curind[i] = i;
1413  Teuchos::RCP<const MV> curV = MVT::CloneView(*V_,curind);
1414 
1415  // Perform spectral decomposition
1416  {
1417 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1418  Teuchos::TimeMonitor lcltimer(*timerDS_);
1419 #endif
1420  int nevlocal = curDim_;
1421  int info = Utils::directSolver(curDim_,*KK_,Teuchos::null,S,theta_,nevlocal,10);
1422  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve returned error code.");
1423  // we did not ask directSolver to perform deflation, so nevLocal better be curDim_
1424  TEUCHOS_TEST_FOR_EXCEPTION(nevlocal != curDim_,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve did not compute all eigenvectors."); // this should never happen
1425  }
1426 
1427  // Sort ritz pairs
1428  {
1429 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1430  Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1431 #endif
1432 
1433  std::vector<int> order(curDim_);
1434  //
1435  // sort the first curDim_ values in theta_
1436  sm_->sort(theta_, Teuchos::rcp(&order,false), curDim_); // don't catch exception
1437  //
1438  // apply the same ordering to the primitive ritz vectors
1439  Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,curDim_,curDim_);
1440  Utils::permuteVectors(order,curS);
1441  }
1442 
1443  // Create a view matrix of the first blockSize_ vectors
1444  Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
1445 
1446  // Compute the new Ritz vectors
1447  {
1448 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1449  Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1450 #endif
1451  MVT::MvTimesMatAddMv(ONE,*curV,S1,ZERO,*X_);
1452  }
1453 
1454  // Apply the stiffness matrix for the Ritz vectors
1455  {
1456 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1457  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1458 #endif
1459  OPT::Apply( *Op_, *X_, *KX_); // don't catch the exception
1460  count_ApplyOp_ += blockSize_;
1461  }
1462  // Apply the mass matrix for the Ritz vectors
1463  if (hasM_) {
1464 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1465  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1466 #endif
1467  OPT::Apply(*MOp_,*X_,*MX_);
1468  count_ApplyM_ += blockSize_;
1469  }
1470  else {
1471  MX_ = X_;
1472  }
1473 
1474  // Compute the residual
1475  // R = KX - MX*diag(theta)
1476  {
1477 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1478  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1479 #endif
1480 
1481  MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1482  Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1483  for (int i = 0; i < blockSize_; ++i) {
1484  T(i,i) = theta_[i];
1485  }
1486  MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
1487  }
1488 
1489  // R has been updated; mark the norms as out-of-date
1490  Rnorms_current_ = false;
1491  R2norms_current_ = false;
1492 
1493 
1494  // When required, monitor some orthogonalities
1495  if (om_->isVerbosity( Debug ) ) {
1496  // Check almost everything here
1497  CheckList chk;
1498  chk.checkV = true;
1499  chk.checkX = true;
1500  chk.checkKX = true;
1501  chk.checkMX = true;
1502  chk.checkR = true;
1503  om_->print( Debug, accuracyCheck(chk, ": after local update") );
1504  }
1505  else if (om_->isVerbosity( OrthoDetails )) {
1506  CheckList chk;
1507  chk.checkX = true;
1508  chk.checkKX = true;
1509  chk.checkMX = true;
1510  chk.checkR = true;
1511  om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
1512  }
1513  } // end while (statusTest == false)
1514 
1515  } // end of iterate()
1516 
1517 
1518 
1520  // compute/return residual M-norms
1521  template <class ScalarType, class MV, class OP>
1522  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1524  if (Rnorms_current_ == false) {
1525  // Update the residual norms
1526  orthman_->norm(*R_,Rnorms_);
1527  Rnorms_current_ = true;
1528  }
1529  return Rnorms_;
1530  }
1531 
1532 
1534  // compute/return residual 2-norms
1535  template <class ScalarType, class MV, class OP>
1536  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1538  if (R2norms_current_ == false) {
1539  // Update the residual 2-norms
1540  MVT::MvNorm(*R_,R2norms_);
1541  R2norms_current_ = true;
1542  }
1543  return R2norms_;
1544  }
1545 
1547  // Set a new StatusTest for the solver.
1548  template <class ScalarType, class MV, class OP>
1550  TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1551  "Anasazi::BlockDavidson::setStatusTest() was passed a null StatusTest.");
1552  tester_ = test;
1553  }
1554 
1556  // Get the current StatusTest used by the solver.
1557  template <class ScalarType, class MV, class OP>
1558  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockDavidson<ScalarType,MV,OP>::getStatusTest() const {
1559  return tester_;
1560  }
1561 
1562 
1564  // Check accuracy, orthogonality, and other debugging stuff
1565  //
1566  // bools specify which tests we want to run (instead of running more than we actually care about)
1567  //
1568  // we don't bother checking the following because they are computed explicitly:
1569  // H == Prec*R
1570  // KH == K*H
1571  //
1572  //
1573  // checkV : V orthonormal
1574  // orthogonal to auxvecs
1575  // checkX : X orthonormal
1576  // orthogonal to auxvecs
1577  // checkMX: check MX == M*X
1578  // checkKX: check KX == K*X
1579  // checkH : H orthonormal
1580  // orthogonal to V and H and auxvecs
1581  // checkMH: check MH == M*H
1582  // checkR : check R orthogonal to X
1583  // checkQ : check that auxiliary vectors are actually orthonormal
1584  // checkKK: check that KK is symmetric in memory
1585  //
1586  // TODO:
1587  // add checkTheta
1588  //
1589  template <class ScalarType, class MV, class OP>
1590  std::string BlockDavidson<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1591  {
1592  using std::endl;
1593 
1594  std::stringstream os;
1595  os.precision(2);
1596  os.setf(std::ios::scientific, std::ios::floatfield);
1597 
1598  os << " Debugging checks: iteration " << iter_ << where << endl;
1599 
1600  // V and friends
1601  std::vector<int> lclind(curDim_);
1602  for (int i=0; i<curDim_; ++i) lclind[i] = i;
1603  Teuchos::RCP<const MV> lclV;
1604  if (initialized_) {
1605  lclV = MVT::CloneView(*V_,lclind);
1606  }
1607  if (chk.checkV && initialized_) {
1608  MagnitudeType err = orthman_->orthonormError(*lclV);
1609  os << " >> Error in V^H M V == I : " << err << endl;
1610  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1611  err = orthman_->orthogError(*lclV,*auxVecs_[i]);
1612  os << " >> Error in V^H M Q[" << i << "] == 0 : " << err << endl;
1613  }
1614  Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
1615  Teuchos::RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
1616  OPT::Apply(*Op_,*lclV,*lclKV);
1617  MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
1618  Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
1619  curKK -= subKK;
1620  // dup the lower tri part
1621  for (int j=0; j<curDim_; ++j) {
1622  for (int i=j+1; i<curDim_; ++i) {
1623  curKK(i,j) = curKK(j,i);
1624  }
1625  }
1626  os << " >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
1627  }
1628 
1629  // X and friends
1630  if (chk.checkX && initialized_) {
1631  MagnitudeType err = orthman_->orthonormError(*X_);
1632  os << " >> Error in X^H M X == I : " << err << endl;
1633  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1634  err = orthman_->orthogError(*X_,*auxVecs_[i]);
1635  os << " >> Error in X^H M Q[" << i << "] == 0 : " << err << endl;
1636  }
1637  }
1638  if (chk.checkMX && hasM_ && initialized_) {
1639  MagnitudeType err = Utils::errorEquality(*X_, *MX_, MOp_);
1640  os << " >> Error in MX == M*X : " << err << endl;
1641  }
1642  if (chk.checkKX && initialized_) {
1643  MagnitudeType err = Utils::errorEquality(*X_, *KX_, Op_);
1644  os << " >> Error in KX == K*X : " << err << endl;
1645  }
1646 
1647  // H and friends
1648  if (chk.checkH && initialized_) {
1649  MagnitudeType err = orthman_->orthonormError(*H_);
1650  os << " >> Error in H^H M H == I : " << err << endl;
1651  err = orthman_->orthogError(*H_,*lclV);
1652  os << " >> Error in H^H M V == 0 : " << err << endl;
1653  err = orthman_->orthogError(*H_,*X_);
1654  os << " >> Error in H^H M X == 0 : " << err << endl;
1655  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1656  err = orthman_->orthogError(*H_,*auxVecs_[i]);
1657  os << " >> Error in H^H M Q[" << i << "] == 0 : " << err << endl;
1658  }
1659  }
1660  if (chk.checkKH && initialized_) {
1661  MagnitudeType err = Utils::errorEquality(*H_, *KH_, Op_);
1662  os << " >> Error in KH == K*H : " << err << endl;
1663  }
1664  if (chk.checkMH && hasM_ && initialized_) {
1665  MagnitudeType err = Utils::errorEquality(*H_, *MH_, MOp_);
1666  os << " >> Error in MH == M*H : " << err << endl;
1667  }
1668 
1669  // R: this is not M-orthogonality, but standard Euclidean orthogonality
1670  if (chk.checkR && initialized_) {
1671  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1672  MVT::MvTransMv(ONE,*X_,*R_,xTx);
1673  MagnitudeType err = xTx.normFrobenius();
1674  os << " >> Error in X^H R == 0 : " << err << endl;
1675  }
1676 
1677  // KK
1678  if (chk.checkKK && initialized_) {
1679  Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_), lclKK(Teuchos::View,*KK_,curDim_,curDim_);
1680  for (int j=0; j<curDim_; ++j) {
1681  for (int i=0; i<curDim_; ++i) {
1682  SDMerr(i,j) = lclKK(i,j) - SCT::conjugate(lclKK(j,i));
1683  }
1684  }
1685  os << " >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
1686  }
1687 
1688  // Q
1689  if (chk.checkQ) {
1690  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1691  MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
1692  os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << err << endl;
1693  for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
1694  err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1695  os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << err << endl;
1696  }
1697  }
1698  }
1699 
1700  os << endl;
1701 
1702  return os.str();
1703  }
1704 
1705 
1707  // Print the current status of the solver
1708  template <class ScalarType, class MV, class OP>
1709  void
1711  {
1712  using std::endl;
1713 
1714  os.setf(std::ios::scientific, std::ios::floatfield);
1715  os.precision(6);
1716  os <<endl;
1717  os <<"================================================================================" << endl;
1718  os << endl;
1719  os <<" BlockDavidson Solver Status" << endl;
1720  os << endl;
1721  os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1722  os <<"The number of iterations performed is " <<iter_<<endl;
1723  os <<"The block size is " << blockSize_<<endl;
1724  os <<"The number of blocks is " << numBlocks_<<endl;
1725  os <<"The current basis size is " << curDim_<<endl;
1726  os <<"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
1727  os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1728  os <<"The number of operations M*x is "<<count_ApplyM_<<endl;
1729  os <<"The number of operations Prec*x is "<<count_ApplyPrec_<<endl;
1730 
1731  os.setf(std::ios_base::right, std::ios_base::adjustfield);
1732 
1733  if (initialized_) {
1734  os << endl;
1735  os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
1736  os << std::setw(20) << "Eigenvalue"
1737  << std::setw(20) << "Residual(M)"
1738  << std::setw(20) << "Residual(2)"
1739  << endl;
1740  os <<"--------------------------------------------------------------------------------"<<endl;
1741  for (int i=0; i<blockSize_; ++i) {
1742  os << std::setw(20) << theta_[i];
1743  if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
1744  else os << std::setw(20) << "not current";
1745  if (R2norms_current_) os << std::setw(20) << R2norms_[i];
1746  else os << std::setw(20) << "not current";
1747  os << endl;
1748  }
1749  }
1750  os <<"================================================================================" << endl;
1751  os << endl;
1752  }
1753 
1754 } // End of namespace Anasazi
1755 
1756 #endif
1757 
1758 // End of file AnasaziBlockDavidson.hpp
virtual ~BlockDavidson()
Anasazi::BlockDavidson destructor.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
Teuchos::RCP< const MV > H
The current preconditioned residual 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 getNumIters() const
Get the current iteration count.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
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 class defines the interface required by an eigensolver and status test class to compute solution...
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
Declaration of basic traits for the multivector type.
BlockDavidsonOrthoFailure is thrown when the orthogonalization manager is unable to orthogonalize the...
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to the given output stream...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< const MV > KX
The image of the current eigenvectors under K.
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...
BlockDavidsonInitFailure is thrown when the BlockDavidson solver is unable to generate an initial ite...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
An exception class parent to all Anasazi exceptions.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values for the previous iteration.
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
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.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For BlockDavidson, this always returns n...
Output managers remove the need for the eigensolver to know any information about the required output...
Teuchos::RCP< const MV > V
The basis for the Krylov space.
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...
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Virtual base class which defines basic traits for the operator type.
int getBlockSize() const
Get the blocksize used by the iterative solver.
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).
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 ...
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
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).
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
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...
int curDim
The current dimension of the solver.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Structure to contain pointers to BlockDavidson state variables.
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...
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Types and exceptions used within Anasazi solvers and interfaces.
std::vector< int > getRitzIndex()
Get the index used for extracting individual Ritz vectors from getRitzVectors().
void iterate()
This method performs BlockDavidson iterations until the status test indicates the need to stop or an ...
BlockDavidsonState< ScalarType, MV > getState() const
Get access to the current state of the eigensolver.
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.
BlockDavidson(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< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
BlockDavidson constructor with eigenproblem, solver utilities, and parameter list of solver options...
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...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
void setBlockSize(int blockSize)
Set the blocksize.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
void resetNumIters()
Reset the iteration count.
Class which provides internal utilities for the Anasazi solvers.
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...
Teuchos::RCP< const MV > getRitzVectors()
Get access to the current Ritz vectors.