Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
46 #ifndef ANASAZI_BLOCK_DAVIDSON_HPP
47 #define ANASAZI_BLOCK_DAVIDSON_HPP
48 
49 #include "AnasaziTypes.hpp"
50 
51 #include "AnasaziEigensolver.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 
57 #include "AnasaziSolverUtils.hpp"
58 
59 #include "Teuchos_LAPACK.hpp"
60 #include "Teuchos_BLAS.hpp"
63 #include "Teuchos_TimeMonitor.hpp"
64 
80 namespace Anasazi {
81 
83 
84 
89  template <class ScalarType, class MV>
95  int curDim;
122  BlockDavidsonState() : curDim(0), V(Teuchos::null),
123  X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null),
124  R(Teuchos::null), H(Teuchos::null),
125  T(Teuchos::null), KK(Teuchos::null) {}
126  };
127 
129 
131 
132 
145  class BlockDavidsonInitFailure : public AnasaziError {public:
146  BlockDavidsonInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
147  {}};
148 
156  class BlockDavidsonOrthoFailure : public AnasaziError {public:
157  BlockDavidsonOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
158  {}};
159 
161 
162 
163  template <class ScalarType, class MV, class OP>
164  class BlockDavidson : public Eigensolver<ScalarType,MV,OP> {
165  public:
167 
168 
178  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
181  Teuchos::ParameterList &params
182  );
183 
185  virtual ~BlockDavidson();
187 
188 
190 
191 
215  void iterate();
216 
251 
255  void initialize();
256 
272  bool isInitialized() const;
273 
287 
289 
290 
292 
293 
295  int getNumIters() const;
296 
298  void resetNumIters();
299 
308 
314  std::vector<Value<ScalarType> > getRitzValues();
315 
316 
324  std::vector<int> getRitzIndex();
325 
326 
332  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
333 
334 
340  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
341 
342 
350  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
351 
357  int getCurSubspaceDim() const;
358 
360  int getMaxSubspaceDim() const;
361 
363 
364 
366 
367 
370 
373 
376 
386  void setBlockSize(int blockSize);
387 
389  int getBlockSize() const;
390 
403  void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
404 
407 
409 
411 
412 
422  void setSize(int blockSize, int numBlocks);
423 
425 
427 
428 
430  void currentStatus(std::ostream &os);
431 
433 
434  private:
435  //
436  // Convenience typedefs
437  //
442  typedef typename SCT::magnitudeType MagnitudeType;
443  const MagnitudeType ONE;
444  const MagnitudeType ZERO;
445  const MagnitudeType NANVAL;
446  //
447  // Internal structs
448  //
449  struct CheckList {
450  bool checkV;
451  bool checkX, checkMX, checkKX;
452  bool checkH, checkMH, checkKH;
453  bool checkR, checkQ;
454  bool checkKK;
455  CheckList() : checkV(false),
456  checkX(false),checkMX(false),checkKX(false),
457  checkH(false),checkMH(false),checkKH(false),
458  checkR(false),checkQ(false),checkKK(false) {};
459  };
460  //
461  // Internal methods
462  //
463  std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
464  //
465  // Classes inputed through constructor that define the eigenproblem to be solved.
466  //
472  //
473  // Information obtained from the eigenproblem
474  //
478  bool hasM_;
479  //
480  // Internal timers
481  //
482  Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
483  timerSortEval_, timerDS_,
484  timerLocal_, timerCompRes_,
485  timerOrtho_, timerInit_;
486  //
487  // Counters
488  //
489  int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
490 
491  //
492  // Algorithmic parameters.
493  //
494  // blockSize_ is the solver block size; it controls the number of eigenvectors that
495  // we compute, the number of residual vectors that we compute, and therefore the number
496  // of vectors added to the basis on each iteration.
497  int blockSize_;
498  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
499  int numBlocks_;
500 
501  //
502  // Current solver state
503  //
504  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
505  // is capable of running; _initialize is controlled by the initialize() member method
506  // For the implications of the state of initialized_, please see documentation for initialize()
507  bool initialized_;
508  //
509  // curDim_ reflects how much of the current basis is valid
510  // NOTE: 0 <= curDim_ <= blockSize_*numBlocks_
511  // this also tells us how many of the values in theta_ are valid Ritz values
512  int curDim_;
513  //
514  // State Multivecs
515  // H_,KH_,MH_ will not own any storage
516  // H_ will occasionally point at the current block of vectors in the basis V_
517  // MH_,KH_ will occasionally point at MX_,KX_ when they are used as temporary storage
518  Teuchos::RCP<MV> X_, KX_, MX_, R_,
519  H_, KH_, MH_,
520  V_;
521  //
522  // Projected matrices
523  //
525  //
526  // auxiliary vectors
528  int numAuxVecs_;
529  //
530  // Number of iterations that have been performed.
531  int iter_;
532  //
533  // Current eigenvalues, residual norms
534  std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
535  //
536  // are the residual norms current with the residual?
537  bool Rnorms_current_, R2norms_current_;
538 
539  };
540 
543  //
544  // Implementations
545  //
548 
549 
551  // Constructor
552  template <class ScalarType, class MV, class OP>
556  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
559  Teuchos::ParameterList &params
560  ) :
561  ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
562  ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
563  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
564  // problem, tools
565  problem_(problem),
566  sm_(sorter),
567  om_(printer),
568  tester_(tester),
569  orthman_(ortho),
570  // timers, counters
571 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
572  timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation Op*x")),
573  timerMOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation M*x")),
574  timerPrec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation Prec*x")),
575  timerSortEval_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Sorting eigenvalues")),
576  timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Direct solve")),
577  timerLocal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Local update")),
578  timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Computing residuals")),
579  timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Orthogonalization")),
580  timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Initialization")),
581 #endif
582  count_ApplyOp_(0),
583  count_ApplyM_(0),
584  count_ApplyPrec_(0),
585  // internal data
586  blockSize_(0),
587  numBlocks_(0),
588  initialized_(false),
589  curDim_(0),
590  auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
591  numAuxVecs_(0),
592  iter_(0),
593  Rnorms_current_(false),
594  R2norms_current_(false)
595  {
596  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
597  "Anasazi::BlockDavidson::constructor: user passed null problem pointer.");
598  TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
599  "Anasazi::BlockDavidson::constructor: user passed null sort manager pointer.");
600  TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
601  "Anasazi::BlockDavidson::constructor: user passed null output manager pointer.");
602  TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
603  "Anasazi::BlockDavidson::constructor: user passed null status test pointer.");
604  TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
605  "Anasazi::BlockDavidson::constructor: user passed null orthogonalization manager pointer.");
606  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
607  "Anasazi::BlockDavidson::constructor: problem is not set.");
608  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
609  "Anasazi::BlockDavidson::constructor: problem is not hermitian.");
610 
611  // get the problem operators
612  Op_ = problem_->getOperator();
613  TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
614  "Anasazi::BlockDavidson::constructor: problem provides no operator.");
615  MOp_ = problem_->getM();
616  Prec_ = problem_->getPrec();
617  hasM_ = (MOp_ != Teuchos::null);
618 
619  // set the block size and allocate data
620  int bs = params.get("Block Size", problem_->getNEV());
621  int nb = params.get("Num Blocks", 2);
622  setSize(bs,nb);
623  }
624 
625 
627  // Destructor
628  template <class ScalarType, class MV, class OP>
630 
631 
633  // Set the block size
634  // This simply calls setSize(), modifying the block size while retaining the number of blocks.
635  template <class ScalarType, class MV, class OP>
637  {
638  setSize(blockSize,numBlocks_);
639  }
640 
641 
643  // Return the current auxiliary vectors
644  template <class ScalarType, class MV, class OP>
646  return auxVecs_;
647  }
648 
649 
650 
652  // return the current block size
653  template <class ScalarType, class MV, class OP>
655  return(blockSize_);
656  }
657 
658 
660  // return eigenproblem
661  template <class ScalarType, class MV, class OP>
663  return(*problem_);
664  }
665 
666 
668  // return max subspace dim
669  template <class ScalarType, class MV, class OP>
671  return blockSize_*numBlocks_;
672  }
673 
675  // return current subspace dim
676  template <class ScalarType, class MV, class OP>
678  if (!initialized_) return 0;
679  return curDim_;
680  }
681 
682 
684  // return ritz residual 2-norms
685  template <class ScalarType, class MV, class OP>
686  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
688  return this->getRes2Norms();
689  }
690 
691 
693  // return ritz index
694  template <class ScalarType, class MV, class OP>
696  std::vector<int> ret(curDim_,0);
697  return ret;
698  }
699 
700 
702  // return ritz values
703  template <class ScalarType, class MV, class OP>
704  std::vector<Value<ScalarType> > BlockDavidson<ScalarType,MV,OP>::getRitzValues() {
705  std::vector<Value<ScalarType> > ret(curDim_);
706  for (int i=0; i<curDim_; ++i) {
707  ret[i].realpart = theta_[i];
708  ret[i].imagpart = ZERO;
709  }
710  return ret;
711  }
712 
713 
715  // return pointer to ritz vectors
716  template <class ScalarType, class MV, class OP>
718  return X_;
719  }
720 
721 
723  // reset number of iterations
724  template <class ScalarType, class MV, class OP>
726  iter_=0;
727  }
728 
729 
731  // return number of iterations
732  template <class ScalarType, class MV, class OP>
734  return(iter_);
735  }
736 
737 
739  // return state pointers
740  template <class ScalarType, class MV, class OP>
743  state.curDim = curDim_;
744  state.V = V_;
745  state.X = X_;
746  state.KX = KX_;
747  if (hasM_) {
748  state.MX = MX_;
749  }
750  else {
751  state.MX = Teuchos::null;
752  }
753  state.R = R_;
754  state.H = H_;
755  state.KK = KK_;
756  if (curDim_ > 0) {
757  state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
758  } else {
759  state.T = Teuchos::rcp(new std::vector<MagnitudeType>(0));
760  }
761  return state;
762  }
763 
764 
766  // Return initialized state
767  template <class ScalarType, class MV, class OP>
768  bool BlockDavidson<ScalarType,MV,OP>::isInitialized() const { return initialized_; }
769 
770 
772  // Set the block size and make necessary adjustments.
773  template <class ScalarType, class MV, class OP>
774  void BlockDavidson<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
775  {
776  // time spent here counts towards timerInit_
777 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
778  Teuchos::TimeMonitor initimer( *timerInit_ );
779 #endif
780 
781  // This routine only allocates space; it doesn't not perform any computation
782  // any change in size will invalidate the state of the solver.
783 
784  TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): blocksize must be strictly positive.");
785  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 2, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): numblocks must be greater than one.");
786  if (blockSize == blockSize_ && numBlocks == numBlocks_) {
787  // do nothing
788  return;
789  }
790 
791  blockSize_ = blockSize;
792  numBlocks_ = numBlocks;
793 
795  // grab some Multivector to Clone
796  // in practice, getInitVec() should always provide this, but it is possible to use a
797  // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
798  // in case of that strange scenario, we will try to Clone from X_ first, then resort to getInitVec()
799  if (X_ != Teuchos::null) { // this is equivalent to blockSize_ > 0
800  tmp = X_;
801  }
802  else {
803  tmp = problem_->getInitVec();
804  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
805  "Anasazi::BlockDavidson::setSize(): eigenproblem did not specify initial vectors to clone from.");
806  }
807 
808  TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*static_cast<ptrdiff_t>(numBlocks) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
809  "Anasazi::BlockDavidson::setSize(): max subspace dimension and auxilliary subspace too large.");
810 
811 
813  // blockSize dependent
814  //
815  // grow/allocate vectors
816  Rnorms_.resize(blockSize_,NANVAL);
817  R2norms_.resize(blockSize_,NANVAL);
818  //
819  // clone multivectors off of tmp
820  //
821  // free current allocation first, to make room for new allocation
822  X_ = Teuchos::null;
823  KX_ = Teuchos::null;
824  MX_ = Teuchos::null;
825  R_ = Teuchos::null;
826  V_ = Teuchos::null;
827 
828  om_->print(Debug," >> Allocating X_\n");
829  X_ = MVT::Clone(*tmp,blockSize_);
830  om_->print(Debug," >> Allocating KX_\n");
831  KX_ = MVT::Clone(*tmp,blockSize_);
832  if (hasM_) {
833  om_->print(Debug," >> Allocating MX_\n");
834  MX_ = MVT::Clone(*tmp,blockSize_);
835  }
836  else {
837  MX_ = X_;
838  }
839  om_->print(Debug," >> Allocating R_\n");
840  R_ = MVT::Clone(*tmp,blockSize_);
841 
843  // blockSize*numBlocks dependent
844  //
845  int newsd = blockSize_*numBlocks_;
846  theta_.resize(blockSize_*numBlocks_,NANVAL);
847  om_->print(Debug," >> Allocating V_\n");
848  V_ = MVT::Clone(*tmp,newsd);
850 
851  om_->print(Debug," >> done allocating.\n");
852 
853  initialized_ = false;
854  curDim_ = 0;
855  }
856 
857 
859  // Set the auxiliary vectors
860  template <class ScalarType, class MV, class OP>
862  typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
863 
864  // set new auxiliary vectors
865  auxVecs_ = auxvecs;
866  numAuxVecs_ = 0;
867  for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
868  numAuxVecs_ += MVT::GetNumberVecs(**i);
869  }
870 
871  // If the solver has been initialized, V is not necessarily orthogonal to new auxiliary vectors
872  if (numAuxVecs_ > 0 && initialized_) {
873  initialized_ = false;
874  }
875 
876  if (om_->isVerbosity( Debug ) ) {
877  CheckList chk;
878  chk.checkQ = true;
879  om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
880  }
881  }
882 
883 
885  /* Initialize the state of the solver
886  *
887  * POST-CONDITIONS:
888  *
889  * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
890  * theta_ contains Ritz w.r.t. V_(1:curDim_)
891  * X is Ritz vectors w.r.t. V_(1:curDim_)
892  * KX = Op*X
893  * MX = M*X if hasM_
894  * R = KX - MX*diag(theta_)
895  *
896  */
897  template <class ScalarType, class MV, class OP>
899  {
900  // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
901  // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
902 
903 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
904  Teuchos::TimeMonitor inittimer( *timerInit_ );
905 #endif
906 
907  std::vector<int> bsind(blockSize_);
908  for (int i=0; i<blockSize_; ++i) bsind[i] = i;
909 
911 
912  // in BlockDavidson, V is primary
913  // the order of dependence follows like so.
914  // --init-> V,KK
915  // --ritz analysis-> theta,X
916  // --op apply-> KX,MX
917  // --compute-> R
918  //
919  // if the user specifies all data for a level, we will accept it.
920  // otherwise, we will generate the whole level, and all subsequent levels.
921  //
922  // the data members are ordered based on dependence, and the levels are
923  // partitioned according to the amount of work required to produce the
924  // items in a level.
925  //
926  // inconsistent multivectors widths and lengths will not be tolerated, and
927  // will be treated with exceptions.
928  //
929  // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to
930  // multivectors in the solver, the copy will not be affected.
931 
932  // set up V and KK: get them from newstate if user specified them
933  // otherwise, set them manually
934  Teuchos::RCP<MV> lclV;
936 
937  if (newstate.V != Teuchos::null && newstate.KK != Teuchos::null) {
938  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
939  "Anasazi::BlockDavidson::initialize(newstate): Vector length of V not correct." );
940  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
941  "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be at least blockSize().");
942  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
943  "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
944  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
945  "Anasazi::BlockDavidson::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
946 
947  curDim_ = newstate.curDim;
948  // pick an integral amount
949  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
950 
951  TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
952  "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
953 
954  // check size of KK
955  TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
956  "Anasazi::BlockDavidson::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
957 
958  // put data in V
959  std::vector<int> nevind(curDim_);
960  for (int i=0; i<curDim_; ++i) nevind[i] = i;
961  if (newstate.V != V_) {
962  if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
963  newstate.V = MVT::CloneView(*newstate.V,nevind);
964  }
965  MVT::SetBlock(*newstate.V,nevind,*V_);
966  }
967  lclV = MVT::CloneViewNonConst(*V_,nevind);
968 
969  // put data into KK_
970  lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
971  if (newstate.KK != KK_) {
972  if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
973  newstate.KK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
974  }
975  lclKK->assign(*newstate.KK);
976  }
977  //
978  // make lclKK Hermitian in memory (copy the upper half to the lower half)
979  for (int j=0; j<curDim_-1; ++j) {
980  for (int i=j+1; i<curDim_; ++i) {
981  (*lclKK)(i,j) = SCT::conjugate((*lclKK)(j,i));
982  }
983  }
984  }
985  else {
986  // user did not specify one of V or KK
987  // get vectors from problem or generate something, projectAndNormalize
988  Teuchos::RCP<const MV> ivec = problem_->getInitVec();
989  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
990  "Anasazi::BlockDavdison::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
991  // clear newstate so we won't use any data from it below
992  newstate.X = Teuchos::null;
993  newstate.MX = Teuchos::null;
994  newstate.KX = Teuchos::null;
995  newstate.R = Teuchos::null;
996  newstate.H = Teuchos::null;
997  newstate.T = Teuchos::null;
998  newstate.KK = Teuchos::null;
999  newstate.V = Teuchos::null;
1000  newstate.curDim = 0;
1001 
1002  curDim_ = MVT::GetNumberVecs(*ivec);
1003  // pick the largest multiple of blockSize_
1004  curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1005  if (curDim_ > blockSize_*numBlocks_) {
1006  // user specified too many vectors... truncate
1007  // this produces a full subspace, but that is okay
1008  curDim_ = blockSize_*numBlocks_;
1009  }
1010  bool userand = false;
1011  if (curDim_ == 0) {
1012  // we need at least blockSize_ vectors
1013  // use a random multivec: ignore everything from InitVec
1014  userand = true;
1015  curDim_ = blockSize_;
1016  }
1017 
1018  // get pointers into V,KV,MV
1019  // tmpVecs will be used below for M*V and K*V (not simultaneously)
1020  // lclV has curDim vectors
1021  // if there is space for lclV and tmpVecs in V_, point tmpVecs into V_
1022  // otherwise, we must allocate space for these products
1023  //
1024  // get pointer to first curDim vector in V_
1025  std::vector<int> dimind(curDim_);
1026  for (int i=0; i<curDim_; ++i) dimind[i] = i;
1027  lclV = MVT::CloneViewNonConst(*V_,dimind);
1028  if (userand) {
1029  // generate random vector data
1030  MVT::MvRandom(*lclV);
1031  }
1032  else {
1033  if (MVT::GetNumberVecs(*ivec) > curDim_) {
1034  ivec = MVT::CloneView(*ivec,dimind);
1035  }
1036  // assign ivec to first part of lclV
1037  MVT::SetBlock(*ivec,dimind,*lclV);
1038  }
1039  Teuchos::RCP<MV> tmpVecs;
1040  if (curDim_*2 <= blockSize_*numBlocks_) {
1041  // partition V_ = [lclV tmpVecs _leftover_]
1042  std::vector<int> block2(curDim_);
1043  for (int i=0; i<curDim_; ++i) block2[i] = curDim_+i;
1044  tmpVecs = MVT::CloneViewNonConst(*V_,block2);
1045  }
1046  else {
1047  // allocate space for tmpVecs
1048  tmpVecs = MVT::Clone(*V_,curDim_);
1049  }
1050 
1051  // compute M*lclV if hasM_
1052  if (hasM_) {
1053 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1054  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1055 #endif
1056  OPT::Apply(*MOp_,*lclV,*tmpVecs);
1057  count_ApplyM_ += curDim_;
1058  }
1059 
1060  // remove auxVecs from lclV and normalize it
1061  if (auxVecs_.size() > 0) {
1062 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1063  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1064 #endif
1065 
1067  int rank = orthman_->projectAndNormalizeMat(*lclV,auxVecs_,dummyC,Teuchos::null,tmpVecs);
1069  "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1070  }
1071  else {
1072 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1073  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1074 #endif
1075 
1076  int rank = orthman_->normalizeMat(*lclV,Teuchos::null,tmpVecs);
1078  "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1079  }
1080 
1081  // compute K*lclV: we are re-using tmpVecs to store the result
1082  {
1083 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1084  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1085 #endif
1086  OPT::Apply(*Op_,*lclV,*tmpVecs);
1087  count_ApplyOp_ += curDim_;
1088  }
1089 
1090  // generate KK
1091  lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1092  MVT::MvTransMv(ONE,*lclV,*tmpVecs,*lclKK);
1093 
1094  // clear tmpVecs
1095  tmpVecs = Teuchos::null;
1096  }
1097 
1098  // X,theta require Ritz analysis; if we have to generate one of these, we might as well generate both
1099  if (newstate.X != Teuchos::null && newstate.T != Teuchos::null) {
1100  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1101  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of X must be consistent with block size and length of V.");
1102  TEUCHOS_TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
1103  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of T must be consistent with dimension of V.");
1104 
1105  if (newstate.X != X_) {
1106  MVT::SetBlock(*newstate.X,bsind,*X_);
1107  }
1108 
1109  std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
1110  }
1111  else {
1112  // compute ritz vecs/vals
1113  Teuchos::SerialDenseMatrix<int,ScalarType> S(curDim_,curDim_);
1114  {
1115 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1116  Teuchos::TimeMonitor lcltimer( *timerDS_ );
1117 #endif
1118  int rank = curDim_;
1119  Utils::directSolver(curDim_, *lclKK, Teuchos::null, S, theta_, rank, 10);
1120  // we want all ritz values back
1122  "Anasazi::BlockDavidson::initialize(newstate): Not enough Ritz vectors to initialize algorithm.");
1123  }
1124  // sort ritz pairs
1125  {
1126 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1127  Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1128 #endif
1129 
1130  std::vector<int> order(curDim_);
1131  //
1132  // sort the first curDim_ values in theta_
1133  sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_); // don't catch exception
1134  //
1135  // apply the same ordering to the primitive ritz vectors
1136  Utils::permuteVectors(order,S);
1137  }
1138 
1139  // compute eigenvectors
1141  {
1142 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1143  Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1144 #endif
1145 
1146  // X <- lclV*S
1147  MVT::MvTimesMatAddMv( ONE, *lclV, S1, ZERO, *X_ );
1148  }
1149  // we generated theta,X so we don't want to use the user's KX,MX
1150  newstate.KX = Teuchos::null;
1151  newstate.MX = Teuchos::null;
1152  }
1153 
1154  // done with local pointers
1155  lclV = Teuchos::null;
1156  lclKK = Teuchos::null;
1157 
1158  // set up KX
1159  if ( newstate.KX != Teuchos::null ) {
1160  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_,
1161  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
1162  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*X_),
1163  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
1164  if (newstate.KX != KX_) {
1165  MVT::SetBlock(*newstate.KX,bsind,*KX_);
1166  }
1167  }
1168  else {
1169  // generate KX
1170  {
1171 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1172  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1173 #endif
1174  OPT::Apply(*Op_,*X_,*KX_);
1175  count_ApplyOp_ += blockSize_;
1176  }
1177  // we generated KX; we will generate R as well
1178  newstate.R = Teuchos::null;
1179  }
1180 
1181  // set up MX
1182  if (hasM_) {
1183  if ( newstate.MX != Teuchos::null ) {
1184  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_,
1185  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
1186  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*X_),
1187  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
1188  if (newstate.MX != MX_) {
1189  MVT::SetBlock(*newstate.MX,bsind,*MX_);
1190  }
1191  }
1192  else {
1193  // generate MX
1194  {
1195 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1196  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1197 #endif
1198  OPT::Apply(*MOp_,*X_,*MX_);
1199  count_ApplyOp_ += blockSize_;
1200  }
1201  // we generated MX; we will generate R as well
1202  newstate.R = Teuchos::null;
1203  }
1204  }
1205  else {
1206  // the assignment MX_==X_ would be redundant; take advantage of this opportunity to debug a little
1207  TEUCHOS_TEST_FOR_EXCEPTION(MX_ != X_, std::logic_error, "Anasazi::BlockDavidson::initialize(): solver invariant not satisfied (MX==X).");
1208  }
1209 
1210  // set up R
1211  if (newstate.R != Teuchos::null) {
1212  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
1213  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
1214  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1215  std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
1216  if (newstate.R != R_) {
1217  MVT::SetBlock(*newstate.R,bsind,*R_);
1218  }
1219  }
1220  else {
1221 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1222  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1223 #endif
1224 
1225  // form R <- KX - MX*T
1226  MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1227  Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1228  T.putScalar(ZERO);
1229  for (int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
1230  MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
1231 
1232  }
1233 
1234  // R has been updated; mark the norms as out-of-date
1235  Rnorms_current_ = false;
1236  R2norms_current_ = false;
1237 
1238  // finally, we are initialized
1239  initialized_ = true;
1240 
1241  if (om_->isVerbosity( Debug ) ) {
1242  // Check almost everything here
1243  CheckList chk;
1244  chk.checkV = true;
1245  chk.checkX = true;
1246  chk.checkKX = true;
1247  chk.checkMX = true;
1248  chk.checkR = true;
1249  chk.checkQ = true;
1250  chk.checkKK = true;
1251  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
1252  }
1253 
1254  // Print information on current status
1255  if (om_->isVerbosity(Debug)) {
1256  currentStatus( om_->stream(Debug) );
1257  }
1258  else if (om_->isVerbosity(IterationDetails)) {
1259  currentStatus( om_->stream(IterationDetails) );
1260  }
1261  }
1262 
1263 
1265  // initialize the solver with default state
1266  template <class ScalarType, class MV, class OP>
1268  {
1270  initialize(empty);
1271  }
1272 
1273 
1275  // Perform BlockDavidson iterations until the StatusTest tells us to stop.
1276  template <class ScalarType, class MV, class OP>
1278  {
1279  //
1280  // Initialize solver state
1281  if (initialized_ == false) {
1282  initialize();
1283  }
1284 
1285  // as a data member, this would be redundant and require synchronization with
1286  // blockSize_ and numBlocks_; we'll use a constant here.
1287  const int searchDim = blockSize_*numBlocks_;
1288 
1290 
1291  //
1292  // The projected matrices are part of the state, but the eigenvectors are defined locally.
1293  // S = Local eigenvectors (size: searchDim * searchDim
1294  Teuchos::SerialDenseMatrix<int,ScalarType> S( searchDim, searchDim );
1295 
1296 
1298  // iterate until the status test tells us to stop.
1299  // also break if our basis is full
1300  while (tester_->checkStatus(this) != Passed && curDim_ < searchDim) {
1301 
1302  // Print information on current iteration
1303  if (om_->isVerbosity(Debug)) {
1304  currentStatus( om_->stream(Debug) );
1305  }
1306  else if (om_->isVerbosity(IterationDetails)) {
1307  currentStatus( om_->stream(IterationDetails) );
1308  }
1309 
1310  ++iter_;
1311 
1312  // get the current part of the basis
1313  std::vector<int> curind(blockSize_);
1314  for (int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
1315  H_ = MVT::CloneViewNonConst(*V_,curind);
1316 
1317  // Apply the preconditioner on the residuals: H <- Prec*R
1318  // H = Prec*R
1319  if (Prec_ != Teuchos::null) {
1320 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1321  Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1322 #endif
1323  OPT::Apply( *Prec_, *R_, *H_ ); // don't catch the exception
1324  count_ApplyPrec_ += blockSize_;
1325  }
1326  else {
1327  std::vector<int> bsind(blockSize_);
1328  for (int i=0; i<blockSize_; ++i) { bsind[i] = i; }
1329  MVT::SetBlock(*R_,bsind,*H_);
1330  }
1331 
1332  // Apply the mass matrix on H
1333  if (hasM_) {
1334  // use memory at MX_ for temporary storage
1335  MH_ = MX_;
1336 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1337  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1338 #endif
1339  OPT::Apply( *MOp_, *H_, *MH_); // don't catch the exception
1340  count_ApplyM_ += blockSize_;
1341  }
1342  else {
1343  MH_ = H_;
1344  }
1345 
1346  // Get a view of the previous vectors
1347  // this is used for orthogonalization and for computing V^H K H
1348  std::vector<int> prevind(curDim_);
1349  for (int i=0; i<curDim_; ++i) prevind[i] = i;
1350  Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,prevind);
1351 
1352  // Orthogonalize H against the previous vectors and the auxiliary vectors, and normalize
1353  {
1354 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1355  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1356 #endif
1357 
1358  Teuchos::Array<Teuchos::RCP<const MV> > against = auxVecs_;
1359  against.push_back(Vprev);
1361  int rank = orthman_->projectAndNormalizeMat(*H_,against,
1362  dummyC,
1363  Teuchos::null,MH_);
1365  "Anasazi::BlockDavidson::iterate(): unable to compute orthonormal basis for H.");
1366  }
1367 
1368  // Apply the stiffness matrix to H
1369  {
1370  // use memory at KX_ for temporary storage
1371  KH_ = KX_;
1372 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1373  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1374 #endif
1375  OPT::Apply( *Op_, *H_, *KH_); // don't catch the exception
1376  count_ApplyOp_ += blockSize_;
1377  }
1378 
1379  if (om_->isVerbosity( Debug ) ) {
1380  CheckList chk;
1381  chk.checkH = true;
1382  chk.checkMH = true;
1383  chk.checkKH = true;
1384  om_->print( Debug, accuracyCheck(chk, ": after ortho H") );
1385  }
1386  else if (om_->isVerbosity( OrthoDetails ) ) {
1387  CheckList chk;
1388  chk.checkH = true;
1389  chk.checkMH = true;
1390  chk.checkKH = true;
1391  om_->print( OrthoDetails, accuracyCheck(chk,": after ortho H") );
1392  }
1393 
1394  // compute next part of the projected matrices
1395  // this this in two parts
1397  // Vprev*K*H
1398  nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_) );
1399  MVT::MvTransMv(ONE,*Vprev,*KH_,*nextKK);
1400  // H*K*H
1401  nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
1402  MVT::MvTransMv(ONE,*H_,*KH_,*nextKK);
1403  //
1404  // make sure that KK_ is Hermitian in memory
1405  nextKK = Teuchos::null;
1406  for (int i=curDim_; i<curDim_+blockSize_; ++i) {
1407  for (int j=0; j<i; ++j) {
1408  (*KK_)(i,j) = SCT::conjugate((*KK_)(j,i));
1409  }
1410  }
1411 
1412  // V has been extended, and KK has been extended. Update basis dim and release all pointers.
1413  curDim_ += blockSize_;
1414  H_ = KH_ = MH_ = Teuchos::null;
1415  Vprev = Teuchos::null;
1416 
1417  if (om_->isVerbosity( Debug ) ) {
1418  CheckList chk;
1419  chk.checkKK = true;
1420  om_->print( Debug, accuracyCheck(chk, ": after expanding KK") );
1421  }
1422 
1423  // Get pointer to complete basis
1424  curind.resize(curDim_);
1425  for (int i=0; i<curDim_; ++i) curind[i] = i;
1426  Teuchos::RCP<const MV> curV = MVT::CloneView(*V_,curind);
1427 
1428  // Perform spectral decomposition
1429  {
1430 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1431  Teuchos::TimeMonitor lcltimer(*timerDS_);
1432 #endif
1433  int nevlocal = curDim_;
1434  int info = Utils::directSolver(curDim_,*KK_,Teuchos::null,S,theta_,nevlocal,10);
1435  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve returned error code.");
1436  // we did not ask directSolver to perform deflation, so nevLocal better be curDim_
1437  TEUCHOS_TEST_FOR_EXCEPTION(nevlocal != curDim_,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve did not compute all eigenvectors."); // this should never happen
1438  }
1439 
1440  // Sort ritz pairs
1441  {
1442 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1443  Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1444 #endif
1445 
1446  std::vector<int> order(curDim_);
1447  //
1448  // sort the first curDim_ values in theta_
1449  sm_->sort(theta_, Teuchos::rcp(&order,false), curDim_); // don't catch exception
1450  //
1451  // apply the same ordering to the primitive ritz vectors
1453  Utils::permuteVectors(order,curS);
1454  }
1455 
1456  // Create a view matrix of the first blockSize_ vectors
1457  Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
1458 
1459  // Compute the new Ritz vectors
1460  {
1461 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1462  Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1463 #endif
1464  MVT::MvTimesMatAddMv(ONE,*curV,S1,ZERO,*X_);
1465  }
1466 
1467  // Apply the stiffness matrix for the Ritz vectors
1468  {
1469 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1470  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1471 #endif
1472  OPT::Apply( *Op_, *X_, *KX_); // don't catch the exception
1473  count_ApplyOp_ += blockSize_;
1474  }
1475  // Apply the mass matrix for the Ritz vectors
1476  if (hasM_) {
1477 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1478  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1479 #endif
1480  OPT::Apply(*MOp_,*X_,*MX_);
1481  count_ApplyM_ += blockSize_;
1482  }
1483  else {
1484  MX_ = X_;
1485  }
1486 
1487  // Compute the residual
1488  // R = KX - MX*diag(theta)
1489  {
1490 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1491  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1492 #endif
1493 
1494  MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1495  Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1496  for (int i = 0; i < blockSize_; ++i) {
1497  T(i,i) = theta_[i];
1498  }
1499  MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
1500  }
1501 
1502  // R has been updated; mark the norms as out-of-date
1503  Rnorms_current_ = false;
1504  R2norms_current_ = false;
1505 
1506 
1507  // When required, monitor some orthogonalities
1508  if (om_->isVerbosity( Debug ) ) {
1509  // Check almost everything here
1510  CheckList chk;
1511  chk.checkV = true;
1512  chk.checkX = true;
1513  chk.checkKX = true;
1514  chk.checkMX = true;
1515  chk.checkR = true;
1516  om_->print( Debug, accuracyCheck(chk, ": after local update") );
1517  }
1518  else if (om_->isVerbosity( OrthoDetails )) {
1519  CheckList chk;
1520  chk.checkX = true;
1521  chk.checkKX = true;
1522  chk.checkMX = true;
1523  chk.checkR = true;
1524  om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
1525  }
1526  } // end while (statusTest == false)
1527 
1528  } // end of iterate()
1529 
1530 
1531 
1533  // compute/return residual M-norms
1534  template <class ScalarType, class MV, class OP>
1535  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1537  if (Rnorms_current_ == false) {
1538  // Update the residual norms
1539  orthman_->norm(*R_,Rnorms_);
1540  Rnorms_current_ = true;
1541  }
1542  return Rnorms_;
1543  }
1544 
1545 
1547  // compute/return residual 2-norms
1548  template <class ScalarType, class MV, class OP>
1549  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1551  if (R2norms_current_ == false) {
1552  // Update the residual 2-norms
1553  MVT::MvNorm(*R_,R2norms_);
1554  R2norms_current_ = true;
1555  }
1556  return R2norms_;
1557  }
1558 
1560  // Set a new StatusTest for the solver.
1561  template <class ScalarType, class MV, class OP>
1563  TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1564  "Anasazi::BlockDavidson::setStatusTest() was passed a null StatusTest.");
1565  tester_ = test;
1566  }
1567 
1569  // Get the current StatusTest used by the solver.
1570  template <class ScalarType, class MV, class OP>
1572  return tester_;
1573  }
1574 
1575 
1577  // Check accuracy, orthogonality, and other debugging stuff
1578  //
1579  // bools specify which tests we want to run (instead of running more than we actually care about)
1580  //
1581  // we don't bother checking the following because they are computed explicitly:
1582  // H == Prec*R
1583  // KH == K*H
1584  //
1585  //
1586  // checkV : V orthonormal
1587  // orthogonal to auxvecs
1588  // checkX : X orthonormal
1589  // orthogonal to auxvecs
1590  // checkMX: check MX == M*X
1591  // checkKX: check KX == K*X
1592  // checkH : H orthonormal
1593  // orthogonal to V and H and auxvecs
1594  // checkMH: check MH == M*H
1595  // checkR : check R orthogonal to X
1596  // checkQ : check that auxiliary vectors are actually orthonormal
1597  // checkKK: check that KK is symmetric in memory
1598  //
1599  // TODO:
1600  // add checkTheta
1601  //
1602  template <class ScalarType, class MV, class OP>
1603  std::string BlockDavidson<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1604  {
1605  using std::endl;
1606 
1607  std::stringstream os;
1608  os.precision(2);
1609  os.setf(std::ios::scientific, std::ios::floatfield);
1610 
1611  os << " Debugging checks: iteration " << iter_ << where << endl;
1612 
1613  // V and friends
1614  std::vector<int> lclind(curDim_);
1615  for (int i=0; i<curDim_; ++i) lclind[i] = i;
1617  if (initialized_) {
1618  lclV = MVT::CloneView(*V_,lclind);
1619  }
1620  if (chk.checkV && initialized_) {
1621  MagnitudeType err = orthman_->orthonormError(*lclV);
1622  os << " >> Error in V^H M V == I : " << err << endl;
1623  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1624  err = orthman_->orthogError(*lclV,*auxVecs_[i]);
1625  os << " >> Error in V^H M Q[" << i << "] == 0 : " << err << endl;
1626  }
1627  Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
1628  Teuchos::RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
1629  OPT::Apply(*Op_,*lclV,*lclKV);
1630  MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
1631  Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
1632  curKK -= subKK;
1633  // dup the lower tri part
1634  for (int j=0; j<curDim_; ++j) {
1635  for (int i=j+1; i<curDim_; ++i) {
1636  curKK(i,j) = curKK(j,i);
1637  }
1638  }
1639  os << " >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
1640  }
1641 
1642  // X and friends
1643  if (chk.checkX && initialized_) {
1644  MagnitudeType err = orthman_->orthonormError(*X_);
1645  os << " >> Error in X^H M X == I : " << err << endl;
1646  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1647  err = orthman_->orthogError(*X_,*auxVecs_[i]);
1648  os << " >> Error in X^H M Q[" << i << "] == 0 : " << err << endl;
1649  }
1650  }
1651  if (chk.checkMX && hasM_ && initialized_) {
1652  MagnitudeType err = Utils::errorEquality(*X_, *MX_, MOp_);
1653  os << " >> Error in MX == M*X : " << err << endl;
1654  }
1655  if (chk.checkKX && initialized_) {
1656  MagnitudeType err = Utils::errorEquality(*X_, *KX_, Op_);
1657  os << " >> Error in KX == K*X : " << err << endl;
1658  }
1659 
1660  // H and friends
1661  if (chk.checkH && initialized_) {
1662  MagnitudeType err = orthman_->orthonormError(*H_);
1663  os << " >> Error in H^H M H == I : " << err << endl;
1664  err = orthman_->orthogError(*H_,*lclV);
1665  os << " >> Error in H^H M V == 0 : " << err << endl;
1666  err = orthman_->orthogError(*H_,*X_);
1667  os << " >> Error in H^H M X == 0 : " << err << endl;
1668  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1669  err = orthman_->orthogError(*H_,*auxVecs_[i]);
1670  os << " >> Error in H^H M Q[" << i << "] == 0 : " << err << endl;
1671  }
1672  }
1673  if (chk.checkKH && initialized_) {
1674  MagnitudeType err = Utils::errorEquality(*H_, *KH_, Op_);
1675  os << " >> Error in KH == K*H : " << err << endl;
1676  }
1677  if (chk.checkMH && hasM_ && initialized_) {
1678  MagnitudeType err = Utils::errorEquality(*H_, *MH_, MOp_);
1679  os << " >> Error in MH == M*H : " << err << endl;
1680  }
1681 
1682  // R: this is not M-orthogonality, but standard Euclidean orthogonality
1683  if (chk.checkR && initialized_) {
1684  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1685  MVT::MvTransMv(ONE,*X_,*R_,xTx);
1686  MagnitudeType err = xTx.normFrobenius();
1687  os << " >> Error in X^H R == 0 : " << err << endl;
1688  }
1689 
1690  // KK
1691  if (chk.checkKK && initialized_) {
1692  Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_), lclKK(Teuchos::View,*KK_,curDim_,curDim_);
1693  for (int j=0; j<curDim_; ++j) {
1694  for (int i=0; i<curDim_; ++i) {
1695  SDMerr(i,j) = lclKK(i,j) - SCT::conjugate(lclKK(j,i));
1696  }
1697  }
1698  os << " >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
1699  }
1700 
1701  // Q
1702  if (chk.checkQ) {
1703  for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1704  MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
1705  os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << err << endl;
1706  for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
1707  err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1708  os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << err << endl;
1709  }
1710  }
1711  }
1712 
1713  os << endl;
1714 
1715  return os.str();
1716  }
1717 
1718 
1720  // Print the current status of the solver
1721  template <class ScalarType, class MV, class OP>
1722  void
1724  {
1725  using std::endl;
1726 
1727  os.setf(std::ios::scientific, std::ios::floatfield);
1728  os.precision(6);
1729  os <<endl;
1730  os <<"================================================================================" << endl;
1731  os << endl;
1732  os <<" BlockDavidson Solver Status" << endl;
1733  os << endl;
1734  os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1735  os <<"The number of iterations performed is " <<iter_<<endl;
1736  os <<"The block size is " << blockSize_<<endl;
1737  os <<"The number of blocks is " << numBlocks_<<endl;
1738  os <<"The current basis size is " << curDim_<<endl;
1739  os <<"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
1740  os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1741  os <<"The number of operations M*x is "<<count_ApplyM_<<endl;
1742  os <<"The number of operations Prec*x is "<<count_ApplyPrec_<<endl;
1743 
1744  os.setf(std::ios_base::right, std::ios_base::adjustfield);
1745 
1746  if (initialized_) {
1747  os << endl;
1748  os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
1749  os << std::setw(20) << "Eigenvalue"
1750  << std::setw(20) << "Residual(M)"
1751  << std::setw(20) << "Residual(2)"
1752  << endl;
1753  os <<"--------------------------------------------------------------------------------"<<endl;
1754  for (int i=0; i<blockSize_; ++i) {
1755  os << std::setw(20) << theta_[i];
1756  if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
1757  else os << std::setw(20) << "not current";
1758  if (R2norms_current_) os << std::setw(20) << R2norms_[i];
1759  else os << std::setw(20) << "not current";
1760  os << endl;
1761  }
1762  }
1763  os <<"================================================================================" << endl;
1764  os << endl;
1765  }
1766 
1767 } // End of namespace Anasazi
1768 
1769 #endif
1770 
1771 // End of file AnasaziBlockDavidson.hpp
virtual ~BlockDavidson()
Anasazi::BlockDavidson destructor.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Teuchos::RCP< const MV > H
The current preconditioned residual vectors.
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.
T & get(const std::string &name, T def_value)
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...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
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.
Anasazi&#39;s templated, static class providing utilities for the solvers.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
int getBlockSize() const
Get the blocksize used by the iterative solver.
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
int curDim
The current dimension of the solver.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
void push_back(const value_type &x)
Structure to contain pointers to BlockDavidson state variables.
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.
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.