Anasazi  Version of the Day
AnasaziLOBPCG.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 
34 /*
35  LOBPCG contains local storage of up to 10*blockSize_ vectors, representing 10 entities
36  X,H,P,R
37  KX,KH,KP (product of K and the above)
38  MX,MH,MP (product of M and the above, not allocated if we don't have an M matrix)
39  If full orthogonalization is enabled, one extra multivector of blockSize_ vectors is required to
40  compute the local update of X and P.
41 
42  A solver is bound to an eigenproblem at declaration.
43  Other solver parameters (e.g., block size, auxiliary vectors) can be changed dynamically.
44 
45  The orthogonalization manager is used to project away from the auxiliary vectors.
46  If full orthogonalization is enabled, the orthogonalization manager is also used to construct an M orthonormal basis.
47  The orthogonalization manager is subclass of MatOrthoManager, which LOBPCG assumes to be defined by the M inner product.
48  LOBPCG will not work correctly if the orthomanager uses a different inner product.
49  */
50 
51 
52 #ifndef ANASAZI_LOBPCG_HPP
53 #define ANASAZI_LOBPCG_HPP
54 
55 #include "AnasaziTypes.hpp"
56 
57 #include "AnasaziEigensolver.hpp"
60 #include "Teuchos_ScalarTraits.hpp"
61 
63 #include "AnasaziSolverUtils.hpp"
64 
65 #include "Teuchos_LAPACK.hpp"
66 #include "Teuchos_BLAS.hpp"
69 #include "Teuchos_TimeMonitor.hpp"
70 
91 namespace Anasazi {
92 
94 
95 
100  template <class ScalarType, class MultiVector>
101  struct LOBPCGState {
108 
115 
122 
132 
135 
138 
139  LOBPCGState() :
140  V(Teuchos::null),KV(Teuchos::null),MV(Teuchos::null),
141  X(Teuchos::null),KX(Teuchos::null),MX(Teuchos::null),
142  P(Teuchos::null),KP(Teuchos::null),MP(Teuchos::null),
143  H(Teuchos::null),KH(Teuchos::null),MH(Teuchos::null),
144  R(Teuchos::null),T(Teuchos::null) {};
145  };
146 
148 
150 
151 
165  class LOBPCGRitzFailure : public AnasaziError {public:
166  LOBPCGRitzFailure(const std::string& what_arg) : AnasaziError(what_arg)
167  {}};
168 
181  class LOBPCGInitFailure : public AnasaziError {public:
182  LOBPCGInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
183  {}};
184 
201  class LOBPCGOrthoFailure : public AnasaziError {public:
202  LOBPCGOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
203  {}};
204 
206 
207 
208  template <class ScalarType, class MV, class OP>
209  class LOBPCG : public Eigensolver<ScalarType,MV,OP> {
210  public:
211 
213 
214 
224  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
227  Teuchos::ParameterList &params
228  );
229 
231  virtual ~LOBPCG() {};
232 
234 
236 
237 
264  void iterate();
265 
291  void initialize(LOBPCGState<ScalarType,MV>& newstate);
292 
296  void initialize();
297 
317  bool isInitialized() const;
318 
330 
332 
334 
335 
337  int getNumIters() const;
338 
340  void resetNumIters();
341 
350 
356  std::vector<Value<ScalarType> > getRitzValues();
357 
365  std::vector<int> getRitzIndex();
366 
367 
373  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
374 
375 
381  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
382 
383 
391  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
392 
393 
402  int getCurSubspaceDim() const;
403 
407  int getMaxSubspaceDim() const;
408 
410 
412 
413 
416 
419 
422 
423 
432  void setBlockSize(int blockSize);
433 
434 
436  int getBlockSize() const;
437 
438 
450  void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
451 
454 
456 
458 
459 
465  void setFullOrtho(bool fullOrtho);
466 
468  bool getFullOrtho() const;
469 
471  bool hasP();
472 
474 
476 
477 
479  void currentStatus(std::ostream &os);
480 
482 
483  private:
484  //
485  //
486  //
487  void setupViews();
488  //
489  // Convenience typedefs
490  //
491  typedef SolverUtils<ScalarType,MV,OP> Utils;
492  typedef MultiVecTraits<ScalarType,MV> MVT;
495  typedef typename SCT::magnitudeType MagnitudeType;
496  const MagnitudeType ONE;
497  const MagnitudeType ZERO;
498  const MagnitudeType NANVAL;
499  //
500  // Internal structs
501  //
502  struct CheckList {
503  bool checkX, checkMX, checkKX;
504  bool checkH, checkMH;
505  bool checkP, checkMP, checkKP;
506  bool checkR, checkQ;
507  CheckList() : checkX(false),checkMX(false),checkKX(false),
508  checkH(false),checkMH(false),
509  checkP(false),checkMP(false),checkKP(false),
510  checkR(false),checkQ(false) {};
511  };
512  //
513  // Internal methods
514  //
515  std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
516  //
517  // Classes inputed through constructor that define the eigenproblem to be solved.
518  //
524  //
525  // Information obtained from the eigenproblem
526  //
530  bool hasM_;
531  //
532  // Internal timers
533  //
534  Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
535  timerSort_,
536  timerLocalProj_, timerDS_,
537  timerLocalUpdate_, timerCompRes_,
538  timerOrtho_, timerInit_;
539  //
540  // Counters
541  //
542  // Number of operator applications
543  int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
544 
545  //
546  // Algorithmic parameters.
547  //
548  // blockSize_ is the solver block size
549  int blockSize_;
550  //
551  // fullOrtho_ dictates whether the orthogonalization procedures specified by Hetmaniuk and Lehoucq should
552  // be activated (see citations at the top of this file)
553  bool fullOrtho_;
554 
555  //
556  // Current solver state
557  //
558  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
559  // is capable of running; _initialize is controlled by the initialize() member method
560  // For the implications of the state of initialized_, please see documentation for initialize()
561  bool initialized_;
562  //
563  // nevLocal_ reflects how much of the current basis is valid (0 <= nevLocal_ <= 3*blockSize_)
564  // this tells us how many of the values in theta_ are valid Ritz values
565  int nevLocal_;
566  //
567  // hasP_ tells us whether there is valid data in P (and KP,MP)
568  bool hasP_;
569  //
570  // State Multivecs
571  // V_, KV_ MV_ and R_ are primary pointers to allocated multivectors
572  Teuchos::RCP<MV> V_, KV_, MV_, R_;
573  // the rest are multivector views into V_, KV_ and MV_
574  Teuchos::RCP<MV> X_, KX_, MX_,
575  H_, KH_, MH_,
576  P_, KP_, MP_;
577 
578  //
579  // if fullOrtho_ == true, then we must produce the following on every iteration:
580  // [newX newP] = [X H P] [CX;CP]
581  // the structure of [CX;CP] when using full orthogonalization does not allow us to
582  // do this in situ, and R_ does not have enough storage for newX and newP. therefore,
583  // we must allocate additional storage for this.
584  // otherwise, when not using full orthogonalization, the structure
585  // [newX newP] = [X H P] [CX1 0 ]
586  // [CX2 CP2] allows us to work using only R as work space
587  // [CX3 CP3]
588  Teuchos::RCP<MV> tmpmvec_;
589  //
590  // auxiliary vectors
592  int numAuxVecs_;
593  //
594  // Number of iterations that have been performed.
595  int iter_;
596  //
597  // Current eigenvalues, residual norms
598  std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
599  //
600  // are the residual norms current with the residual?
601  bool Rnorms_current_, R2norms_current_;
602 
603  };
604 
605 
606 
607 
609  // Constructor
610  template <class ScalarType, class MV, class OP>
614  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
617  Teuchos::ParameterList &params
618  ) :
619  ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
620  ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
621  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
622  // problem, tools
623  problem_(problem),
624  sm_(sorter),
625  om_(printer),
626  tester_(tester),
627  orthman_(ortho),
628  // timers, counters
629 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
630  timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Operation Op*x")),
631  timerMOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Operation M*x")),
632  timerPrec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Operation Prec*x")),
633  timerSort_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Sorting eigenvalues")),
634  timerLocalProj_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Local projection")),
635  timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Direct solve")),
636  timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Local update")),
637  timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Computing residuals")),
638  timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Orthogonalization")),
639  timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Initialization")),
640 #endif
641  count_ApplyOp_(0),
642  count_ApplyM_(0),
643  count_ApplyPrec_(0),
644  // internal data
645  blockSize_(0),
646  fullOrtho_(params.get("Full Ortho", true)),
647  initialized_(false),
648  nevLocal_(0),
649  hasP_(false),
650  auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
651  numAuxVecs_(0),
652  iter_(0),
653  Rnorms_current_(false),
654  R2norms_current_(false)
655  {
656  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
657  "Anasazi::LOBPCG::constructor: user passed null problem pointer.");
658  TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
659  "Anasazi::LOBPCG::constructor: user passed null sort manager pointer.");
660  TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
661  "Anasazi::LOBPCG::constructor: user passed null output manager pointer.");
662  TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
663  "Anasazi::LOBPCG::constructor: user passed null status test pointer.");
664  TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
665  "Anasazi::LOBPCG::constructor: user passed null orthogonalization manager pointer.");
666  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
667  "Anasazi::LOBPCG::constructor: problem is not set.");
668  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
669  "Anasazi::LOBPCG::constructor: problem is not Hermitian; LOBPCG requires Hermitian problem.");
670 
671  // get the problem operators
672  Op_ = problem_->getOperator();
673  TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
674  "Anasazi::LOBPCG::constructor: problem provides no operator.");
675  MOp_ = problem_->getM();
676  Prec_ = problem_->getPrec();
677  hasM_ = (MOp_ != Teuchos::null);
678 
679  // set the block size and allocate data
680  int bs = params.get("Block Size", problem_->getNEV());
681  setBlockSize(bs);
682  }
683 
684 
686  // Set the block size and make necessary adjustments.
687  template <class ScalarType, class MV, class OP>
689  {
690  // time spent here counts towards timerInit_
691 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
692  Teuchos::TimeMonitor lcltimer( *timerInit_ );
693 #endif
694 
695  // This routine only allocates space; it doesn't not perform any computation
696  // if size is decreased, take the first newBS vectors of all and leave state as is
697  // otherwise, grow/allocate space and set solver to unitialized
698 
700  // grab some Multivector to Clone
701  // in practice, getInitVec() should always provide this, but it is possible to use a
702  // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
703  // in case of that strange scenario, we will try to Clone from R_ because it is smaller
704  // than V_, and we don't want to keep V_ around longer than necessary
705  if (blockSize_ > 0) {
706  tmp = R_;
707  }
708  else {
709  tmp = problem_->getInitVec();
710  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
711  "Anasazi::LOBPCG::setBlockSize(): eigenproblem did not specify initial vectors to clone from.");
712  }
713 
714  TEUCHOS_TEST_FOR_EXCEPTION(newBS <= 0 || static_cast<ptrdiff_t>(newBS) > MVT::GetGlobalLength(*tmp),
715  std::invalid_argument, "Anasazi::LOBPCG::setBlockSize(): block size must be strictly positive.");
716  if (newBS == blockSize_) {
717  // do nothing
718  return;
719  }
720  else if (newBS < blockSize_ && initialized_) {
721  //
722  // shrink vectors
723 
724  // release views so we can modify the bases
725  X_ = Teuchos::null;
726  KX_ = Teuchos::null;
727  MX_ = Teuchos::null;
728  H_ = Teuchos::null;
729  KH_ = Teuchos::null;
730  MH_ = Teuchos::null;
731  P_ = Teuchos::null;
732  KP_ = Teuchos::null;
733  MP_ = Teuchos::null;
734 
735  // make new indices vectors
736  std::vector<int> newind(newBS), oldind(newBS);
737  for (int i=0; i<newBS; i++) {
738  newind[i] = i;
739  oldind[i] = i;
740  }
741 
742  Teuchos::RCP<MV> newV, newMV, newKV, newR;
744  // allocate R and newV
745  newR = MVT::Clone(*tmp,newBS);
746  newV = MVT::Clone(*tmp,newBS*3);
747  newKV = MVT::Clone(*tmp,newBS*3);
748  if (hasM_) {
749  newMV = MVT::Clone(*tmp,newBS*3);
750  }
751 
752  //
753  // if we are initialized, we want to pull the data from V_ into newV:
754  // bs | bs | bs
755  // newV = [newX | **** |newP ]
756  // newKV = [newKX| **** |newKP]
757  // newMV = [newMX| **** |newMP]
758  // where
759  // oldbs | oldbs | oldbs
760  // V_ = [newX *** | ******* | newP ***]
761  // KV_ = [newKX *** | ******* | newKP ***]
762  // MV_ = [newMX *** | ******* | newMP ***]
763  //
764  // we don't care to copy the data corresponding to H
765  // we will not copy the M data if !hasM_, because it doesn't exist
766  //
767 
768  // these are shrink operations which preserve their data
769  theta_.resize(3*newBS);
770  Rnorms_.resize(newBS);
771  R2norms_.resize(newBS);
772 
773  // copy residual vectors: oldind,newind currently contains [0,...,newBS-1]
774  src = MVT::CloneView(*R_,newind);
775  MVT::SetBlock(*src,newind,*newR);
776  // free old memory and point to new memory
777  R_ = newR;
778 
779  // copy in order: newX newKX newMX, then newP newKP newMP
780  // for X: [0,bs-1] <-- [0,bs-1]
781  src = MVT::CloneView(*V_,oldind);
782  MVT::SetBlock(*src,newind,*newV);
783  src = MVT::CloneView(*KV_,oldind);
784  MVT::SetBlock(*src,newind,*newKV);
785  if (hasM_) {
786  src = MVT::CloneView(*MV_,oldind);
787  MVT::SetBlock(*src,newind,*newMV);
788  }
789  // for P: [2*bs, 3*bs-1] <-- [2*oldbs, 2*oldbs+bs-1]
790  for (int i=0; i<newBS; i++) {
791  newind[i] += 2*newBS;
792  oldind[i] += 2*blockSize_;
793  }
794  src = MVT::CloneView(*V_,oldind);
795  MVT::SetBlock(*src,newind,*newV);
796  src = MVT::CloneView(*KV_,oldind);
797  MVT::SetBlock(*src,newind,*newKV);
798  if (hasM_) {
799  src = MVT::CloneView(*MV_,oldind);
800  MVT::SetBlock(*src,newind,*newMV);
801  }
802 
803  // release temp view
804  src = Teuchos::null;
805 
806  // release old allocations and point at new memory
807  V_ = newV;
808  KV_ = newKV;
809  if (hasM_) {
810  MV_ = newMV;
811  }
812  else {
813  MV_ = V_;
814  }
815  }
816  else {
817  // newBS > blockSize_ or not initialized
818  // this is also the scenario for our initial call to setBlockSize(), in the constructor
819  initialized_ = false;
820  hasP_ = false;
821 
822  // release views
823  X_ = Teuchos::null;
824  KX_ = Teuchos::null;
825  MX_ = Teuchos::null;
826  H_ = Teuchos::null;
827  KH_ = Teuchos::null;
828  MH_ = Teuchos::null;
829  P_ = Teuchos::null;
830  KP_ = Teuchos::null;
831  MP_ = Teuchos::null;
832 
833  // free allocated storage
834  R_ = Teuchos::null;
835  V_ = Teuchos::null;
836 
837  // allocate scalar vectors
838  theta_.resize(3*newBS,NANVAL);
839  Rnorms_.resize(newBS,NANVAL);
840  R2norms_.resize(newBS,NANVAL);
841 
842  // clone multivectors off of tmp
843  R_ = MVT::Clone(*tmp,newBS);
844  V_ = MVT::Clone(*tmp,newBS*3);
845  KV_ = MVT::Clone(*tmp,newBS*3);
846  if (hasM_) {
847  MV_ = MVT::Clone(*tmp,newBS*3);
848  }
849  else {
850  MV_ = V_;
851  }
852  }
853 
854  // allocate tmp space
855  tmpmvec_ = Teuchos::null;
856  if (fullOrtho_) {
857  tmpmvec_ = MVT::Clone(*tmp,newBS);
858  }
859 
860  // set new block size
861  blockSize_ = newBS;
862 
863  // setup new views
864  setupViews();
865  }
866 
867 
869  // Setup views into V,KV,MV
870  template <class ScalarType, class MV, class OP>
872  {
873  std::vector<int> ind(blockSize_);
874 
875  for (int i=0; i<blockSize_; i++) {
876  ind[i] = i;
877  }
878  X_ = MVT::CloneViewNonConst(*V_,ind);
879  KX_ = MVT::CloneViewNonConst(*KV_,ind);
880  if (hasM_) {
881  MX_ = MVT::CloneViewNonConst(*MV_,ind);
882  }
883  else {
884  MX_ = X_;
885  }
886 
887  for (int i=0; i<blockSize_; i++) {
888  ind[i] += blockSize_;
889  }
890  H_ = MVT::CloneViewNonConst(*V_,ind);
891  KH_ = MVT::CloneViewNonConst(*KV_,ind);
892  if (hasM_) {
893  MH_ = MVT::CloneViewNonConst(*MV_,ind);
894  }
895  else {
896  MH_ = H_;
897  }
898 
899  for (int i=0; i<blockSize_; i++) {
900  ind[i] += blockSize_;
901  }
902  P_ = MVT::CloneViewNonConst(*V_,ind);
903  KP_ = MVT::CloneViewNonConst(*KV_,ind);
904  if (hasM_) {
905  MP_ = MVT::CloneViewNonConst(*MV_,ind);
906  }
907  else {
908  MP_ = P_;
909  }
910  }
911 
912 
914  // Set the auxiliary vectors
915  template <class ScalarType, class MV, class OP>
917  typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
918 
919  // set new auxiliary vectors
920  auxVecs_ = auxvecs;
921 
922  numAuxVecs_ = 0;
923  for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
924  numAuxVecs_ += MVT::GetNumberVecs(**i);
925  }
926 
927  // If the solver has been initialized, X and P are not necessarily orthogonal to new auxiliary vectors
928  if (numAuxVecs_ > 0 && initialized_) {
929  initialized_ = false;
930  hasP_ = false;
931  }
932 
933  if (om_->isVerbosity( Debug ) ) {
934  // Check almost everything here
935  CheckList chk;
936  chk.checkQ = true;
937  om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
938  }
939  }
940 
941 
943  /* Initialize the state of the solver
944  *
945  * POST-CONDITIONS:
946  *
947  * initialized_ == true
948  * X is orthonormal, orthogonal to auxVecs_
949  * KX = Op*X
950  * MX = M*X if hasM_
951  * theta_ contains Ritz values of X
952  * R = KX - MX*diag(theta_)
953  * if hasP() == true,
954  * P orthogonal to auxVecs_
955  * if fullOrtho_ == true,
956  * P orthonormal and orthogonal to X
957  * KP = Op*P
958  * MP = M*P
959  */
960  template <class ScalarType, class MV, class OP>
962  {
963  // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
964  // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
965 
966 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
967  Teuchos::TimeMonitor inittimer( *timerInit_ );
968 #endif
969 
970  std::vector<int> bsind(blockSize_);
971  for (int i=0; i<blockSize_; i++) bsind[i] = i;
972 
973  // in LOBPCG, X (the subspace iterate) is primary
974  // the order of dependence follows like so.
975  // --init-> X
976  // --op apply-> MX,KX
977  // --ritz analysis-> theta
978  // --optional-> P,MP,KP
979  //
980  // if the user specifies all data for a level, we will accept it.
981  // otherwise, we will generate the whole level, and all subsequent levels.
982  //
983  // the data members are ordered based on dependence, and the levels are
984  // partitioned according to the amount of work required to produce the
985  // items in a level.
986  //
987  // inconsitent multivectors widths and lengths will not be tolerated, and
988  // will be treated with exceptions.
989 
990  // set up X, KX, MX: get them from "state" if user specified them
991 
992  //----------------------------------------
993  // set up X, MX, KX
994  //----------------------------------------
995  if (newstate.X != Teuchos::null) {
996  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
997  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.X not correct." );
998  // newstate.X must have blockSize_ vectors; any more will be ignored
999  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
1000  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.X must have at least block size vectors.");
1001 
1002  // put X data in X_
1003  MVT::SetBlock(*newstate.X,bsind,*X_);
1004 
1005  // put MX data in MX_
1006  if (hasM_) {
1007  if (newstate.MX != Teuchos::null) {
1008  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
1009  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MX not correct." );
1010  // newstate.MX must have blockSize_ vectors; any more will be ignored
1011  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MX) < blockSize_,
1012  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.MX must have at least block size vectors.");
1013  MVT::SetBlock(*newstate.MX,bsind,*MX_);
1014  }
1015  else {
1016  // user didn't specify MX, compute it
1017  {
1018 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1019  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1020 #endif
1021  OPT::Apply(*MOp_,*X_,*MX_);
1022  count_ApplyM_ += blockSize_;
1023  }
1024  // we generated MX; we will generate R as well
1025  newstate.R = Teuchos::null;
1026  }
1027  }
1028 
1029  // put data in KX
1030  if (newstate.KX != Teuchos::null) {
1031  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
1032  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KX not correct." );
1033  // newstate.KX must have blockSize_ vectors; any more will be ignored
1034  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KX) < blockSize_,
1035  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.KX must have at least block size vectors.");
1036  MVT::SetBlock(*newstate.KX,bsind,*KX_);
1037  }
1038  else {
1039  // user didn't specify KX, compute it
1040  {
1041 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1042  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1043 #endif
1044  OPT::Apply(*Op_,*X_,*KX_);
1045  count_ApplyOp_ += blockSize_;
1046  }
1047  // we generated KX; we will generate R as well
1048  newstate.R = Teuchos::null;
1049  }
1050  }
1051  else {
1052  // user did not specify X
1053  // we will initialize X, compute KX and MX, and compute R
1054  //
1055  // clear state so we won't use any data from it below
1056  newstate.P = Teuchos::null;
1057  newstate.KP = Teuchos::null;
1058  newstate.MP = Teuchos::null;
1059  newstate.R = Teuchos::null;
1060  newstate.T = Teuchos::null;
1061 
1062  // generate a basis and projectAndNormalize
1063  Teuchos::RCP<const MV> ivec = problem_->getInitVec();
1064  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
1065  "Anasazi::LOBPCG::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1066 
1067  int initSize = MVT::GetNumberVecs(*ivec);
1068  if (initSize > blockSize_) {
1069  // we need only the first blockSize_ vectors from ivec; get a view of them
1070  initSize = blockSize_;
1071  std::vector<int> ind(blockSize_);
1072  for (int i=0; i<blockSize_; i++) ind[i] = i;
1073  ivec = MVT::CloneView(*ivec,ind);
1074  }
1075 
1076  // assign ivec to first part of X_
1077  if (initSize > 0) {
1078  std::vector<int> ind(initSize);
1079  for (int i=0; i<initSize; i++) ind[i] = i;
1080  MVT::SetBlock(*ivec,ind,*X_);
1081  }
1082  // fill the rest of X_ with random
1083  if (blockSize_ > initSize) {
1084  std::vector<int> ind(blockSize_ - initSize);
1085  for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1086  Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X_,ind);
1087  MVT::MvRandom(*rX);
1088  rX = Teuchos::null;
1089  }
1090 
1091  // put data in MX
1092  if (hasM_) {
1093 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1094  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1095 #endif
1096  OPT::Apply(*MOp_,*X_,*MX_);
1097  count_ApplyM_ += blockSize_;
1098  }
1099 
1100  // remove auxVecs from X_ and normalize it
1101  if (numAuxVecs_ > 0) {
1102 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1103  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1104 #endif
1106  int rank = orthman_->projectAndNormalizeMat(*X_,auxVecs_,dummy,Teuchos::null,MX_);
1107  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, LOBPCGInitFailure,
1108  "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
1109  }
1110  else {
1111 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1112  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1113 #endif
1114  int rank = orthman_->normalizeMat(*X_,Teuchos::null,MX_);
1115  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, LOBPCGInitFailure,
1116  "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
1117  }
1118 
1119  // put data in KX
1120  {
1121 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1122  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1123 #endif
1124  OPT::Apply(*Op_,*X_,*KX_);
1125  count_ApplyOp_ += blockSize_;
1126  }
1127  } // end if (newstate.X != Teuchos::null)
1128 
1129 
1130  //----------------------------------------
1131  // set up Ritz values
1132  //----------------------------------------
1133  theta_.resize(3*blockSize_,NANVAL);
1134  if (newstate.T != Teuchos::null) {
1135  TEUCHOS_TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_,
1136  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.T must contain at least block size Ritz values.");
1137  for (int i=0; i<blockSize_; i++) {
1138  theta_[i] = (*newstate.T)[i];
1139  }
1140  nevLocal_ = blockSize_;
1141  }
1142  else {
1143  // get ritz vecs/vals
1144  Teuchos::SerialDenseMatrix<int,ScalarType> KK(blockSize_,blockSize_),
1145  MM(blockSize_,blockSize_),
1146  S(blockSize_,blockSize_);
1147  {
1148 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1149  Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1150 #endif
1151  // project K
1152  MVT::MvTransMv(ONE,*X_,*KX_,KK);
1153  // project M
1154  MVT::MvTransMv(ONE,*X_,*MX_,MM);
1155  nevLocal_ = blockSize_;
1156  }
1157 
1158  // solve the projected problem
1159  {
1160 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1161  Teuchos::TimeMonitor lcltimer( *timerDS_ );
1162 #endif
1163  Utils::directSolver(blockSize_, KK, Teuchos::rcpFromRef(MM), S, theta_, nevLocal_, 1);
1164  TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,LOBPCGInitFailure,
1165  "Anasazi::LOBPCG::initialize(): Initial Ritz analysis did not produce enough Ritz pairs to initialize algorithm.");
1166  }
1167 
1168  // We only have blockSize_ ritz pairs, ergo we do not need to select.
1169  // However, we still require them to be ordered correctly
1170  {
1171 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1172  Teuchos::TimeMonitor lcltimer( *timerSort_ );
1173 #endif
1174 
1175  std::vector<int> order(blockSize_);
1176  //
1177  // sort the first blockSize_ values in theta_
1178  sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_); // don't catch exception
1179  //
1180  // apply the same ordering to the primitive ritz vectors
1181  Utils::permuteVectors(order,S);
1182  }
1183 
1184  // update the solution, use R for storage
1185  {
1186 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1187  Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1188 #endif
1189  // X <- X*S
1190  MVT::MvAddMv( ONE, *X_, ZERO, *X_, *R_ );
1191  MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X_ );
1192  // KX <- KX*S
1193  MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1194  MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *KX_ );
1195  if (hasM_) {
1196  // MX <- MX*S
1197  MVT::MvAddMv( ONE, *MX_, ZERO, *MX_, *R_ );
1198  MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *MX_ );
1199  }
1200  }
1201  }
1202 
1203  //----------------------------------------
1204  // compute R
1205  //----------------------------------------
1206  if (newstate.R != Teuchos::null) {
1207  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
1208  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.R not correct." );
1209  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
1210  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.R must have blockSize number of vectors." );
1211  MVT::SetBlock(*newstate.R,bsind,*R_);
1212  }
1213  else {
1214 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1215  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1216 #endif
1217  // form R <- KX - MX*T
1218  MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1219  Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1220  for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
1221  MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
1222  }
1223 
1224  // R has been updated; mark the norms as out-of-date
1225  Rnorms_current_ = false;
1226  R2norms_current_ = false;
1227 
1228  // put data in P,KP,MP: P is not used to set theta
1229  if (newstate.P != Teuchos::null) {
1230  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.P) < blockSize_ ,
1231  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.P must have blockSize number of vectors." );
1232  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.P) != MVT::GetGlobalLength(*P_),
1233  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.P not correct." );
1234  hasP_ = true;
1235 
1236  // set P_
1237  MVT::SetBlock(*newstate.P,bsind,*P_);
1238 
1239  // set/compute KP_
1240  if (newstate.KP != Teuchos::null) {
1241  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KP) < blockSize_,
1242  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.KP must have blockSize number of vectors." );
1243  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KP) != MVT::GetGlobalLength(*KP_),
1244  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KP not correct." );
1245  MVT::SetBlock(*newstate.KP,bsind,*KP_);
1246  }
1247  else {
1248 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1249  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1250 #endif
1251  OPT::Apply(*Op_,*P_,*KP_);
1252  count_ApplyOp_ += blockSize_;
1253  }
1254 
1255  // set/compute MP_
1256  if (hasM_) {
1257  if (newstate.MP != Teuchos::null) {
1258  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MP) < blockSize_,
1259  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.MP must have blockSize number of vectors." );
1260  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MP) != MVT::GetGlobalLength(*MP_),
1261  std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MP not correct." );
1262  MVT::SetBlock(*newstate.MP,bsind,*MP_);
1263  }
1264  else {
1265 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1266  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1267 #endif
1268  OPT::Apply(*MOp_,*P_,*MP_);
1269  count_ApplyM_ += blockSize_;
1270  }
1271  }
1272  }
1273  else {
1274  hasP_ = false;
1275  }
1276 
1277  // finally, we are initialized
1278  initialized_ = true;
1279 
1280  if (om_->isVerbosity( Debug ) ) {
1281  // Check almost everything here
1282  CheckList chk;
1283  chk.checkX = true;
1284  chk.checkKX = true;
1285  chk.checkMX = true;
1286  chk.checkP = true;
1287  chk.checkKP = true;
1288  chk.checkMP = true;
1289  chk.checkR = true;
1290  chk.checkQ = true;
1291  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
1292  }
1293 
1294  }
1295 
1296  template <class ScalarType, class MV, class OP>
1298  {
1300  initialize(empty);
1301  }
1302 
1303 
1305  // Instruct the solver to use full orthogonalization
1306  template <class ScalarType, class MV, class OP>
1308  {
1309  if ( fullOrtho_ == true || initialized_ == false || fullOrtho == fullOrtho_ ) {
1310  // state is already orthogonalized or solver is not initialized
1311  fullOrtho_ = fullOrtho;
1312  }
1313  else {
1314  // solver is initialized, state is not fully orthogonalized, and user has requested full orthogonalization
1315  // ergo, we must throw away data in P
1316  fullOrtho_ = true;
1317  hasP_ = false;
1318  }
1319 
1320  // the user has called setFullOrtho, so the class has been instantiated
1321  // ergo, the data has already been allocated, i.e., setBlockSize() has been called
1322  // if it is already allocated, it should be the proper size
1323  if (fullOrtho_ && tmpmvec_ == Teuchos::null) {
1324  // allocated the workspace
1325  tmpmvec_ = MVT::Clone(*X_,blockSize_);
1326  }
1327  else if (fullOrtho_==false) {
1328  // free the workspace
1329  tmpmvec_ = Teuchos::null;
1330  }
1331  }
1332 
1333 
1335  // Perform LOBPCG iterations until the StatusTest tells us to stop.
1336  template <class ScalarType, class MV, class OP>
1338  {
1339  //
1340  // Allocate/initialize data structures
1341  //
1342  if (initialized_ == false) {
1343  initialize();
1344  }
1345 
1346  //
1347  // Miscellaneous definitions
1348  const int oneBlock = blockSize_;
1349  const int twoBlocks = 2*blockSize_;
1350  const int threeBlocks = 3*blockSize_;
1351 
1352  std::vector<int> indblock1(blockSize_), indblock2(blockSize_), indblock3(blockSize_);
1353  for (int i=0; i<blockSize_; i++) {
1354  indblock1[i] = i;
1355  indblock2[i] = i + blockSize_;
1356  indblock3[i] = i + 2*blockSize_;
1357  }
1358 
1359  //
1360  // Define dense projected/local matrices
1361  // KK = Local stiffness matrix (size: 3*blockSize_ x 3*blockSize_)
1362  // MM = Local mass matrix (size: 3*blockSize_ x 3*blockSize_)
1363  // S = Local eigenvectors (size: 3*blockSize_ x 3*blockSize_)
1364  Teuchos::SerialDenseMatrix<int,ScalarType> KK( threeBlocks, threeBlocks ),
1365  MM( threeBlocks, threeBlocks ),
1366  S( threeBlocks, threeBlocks );
1367 
1368  while (tester_->checkStatus(this) != Passed) {
1369 
1370  // Print information on current status
1371  if (om_->isVerbosity(Debug)) {
1372  currentStatus( om_->stream(Debug) );
1373  }
1374  else if (om_->isVerbosity(IterationDetails)) {
1375  currentStatus( om_->stream(IterationDetails) );
1376  }
1377 
1378  // increment iteration counter
1379  iter_++;
1380 
1381  // Apply the preconditioner on the residuals: H <- Prec*R
1382  if (Prec_ != Teuchos::null) {
1383 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1384  Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1385 #endif
1386  OPT::Apply( *Prec_, *R_, *H_ ); // don't catch the exception
1387  count_ApplyPrec_ += blockSize_;
1388  }
1389  else {
1390  MVT::MvAddMv(ONE,*R_,ZERO,*R_,*H_);
1391  }
1392 
1393  // Apply the mass matrix on H
1394  if (hasM_) {
1395 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1396  Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1397 #endif
1398  OPT::Apply( *MOp_, *H_, *MH_); // don't catch the exception
1399  count_ApplyM_ += blockSize_;
1400  }
1401 
1402  // orthogonalize H against the auxiliary vectors
1403  // optionally: orthogonalize H against X and P ([X P] is already orthonormal)
1405  Q = auxVecs_;
1406  if (fullOrtho_) {
1407  // X and P are not contiguous, so there is not much point in putting them under
1408  // a single multivector view
1409  Q.push_back(X_);
1410  if (hasP_) {
1411  Q.push_back(P_);
1412  }
1413  }
1414  {
1415 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1416  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1417 #endif
1419  Teuchos::tuple<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > >(Teuchos::null);
1420  int rank = orthman_->projectAndNormalizeMat(*H_,Q,dummyC,Teuchos::null,MH_);
1421  // our views are currently in place; it is safe to throw an exception
1423  "Anasazi::LOBPCG::iterate(): unable to compute orthonormal basis for H.");
1424  }
1425 
1426  if (om_->isVerbosity( Debug ) ) {
1427  CheckList chk;
1428  chk.checkH = true;
1429  chk.checkMH = true;
1430  om_->print( Debug, accuracyCheck(chk, ": after ortho H") );
1431  }
1432  else if (om_->isVerbosity( OrthoDetails ) ) {
1433  CheckList chk;
1434  chk.checkH = true;
1435  chk.checkMH = true;
1436  om_->print( OrthoDetails, accuracyCheck(chk,": after ortho H") );
1437  }
1438 
1439  // Apply the stiffness matrix to H
1440  {
1441 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1442  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1443 #endif
1444  OPT::Apply( *Op_, *H_, *KH_); // don't catch the exception
1445  count_ApplyOp_ += blockSize_;
1446  }
1447 
1448  if (hasP_) {
1449  nevLocal_ = threeBlocks;
1450  }
1451  else {
1452  nevLocal_ = twoBlocks;
1453  }
1454 
1455  //
1456  // we need bases: [X H P] and [H P] (only need the latter if fullOrtho == false)
1457  // we need to perform the following operations:
1458  // X' [KX KH KP]
1459  // X' [MX MH MP]
1460  // H' [KH KP]
1461  // H' [MH MP]
1462  // P' [KP]
1463  // P' [MP]
1464  // [X H P] CX
1465  // [X H P] CP if fullOrtho
1466  // [H P] CP if !fullOrtho
1467  //
1468  // since M[X H P] is potentially the same memory as [X H P], and
1469  // because we are not allowed to have overlapping non-const views of
1470  // a multivector, we will now abandon our non-const views in favor of
1471  // const views
1472  //
1473  X_ = Teuchos::null;
1474  KX_ = Teuchos::null;
1475  MX_ = Teuchos::null;
1476  H_ = Teuchos::null;
1477  KH_ = Teuchos::null;
1478  MH_ = Teuchos::null;
1479  P_ = Teuchos::null;
1480  KP_ = Teuchos::null;
1481  MP_ = Teuchos::null;
1482  Teuchos::RCP<const MV> cX, cH, cXHP, cHP, cK_XHP, cK_HP, cM_XHP, cM_HP, cP, cK_P, cM_P;
1483  {
1484  cX = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock1);
1485  cH = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock2);
1486 
1487  std::vector<int> indXHP(nevLocal_);
1488  for (int i=0; i<nevLocal_; i++) {
1489  indXHP[i] = i;
1490  }
1491  cXHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indXHP);
1492  cK_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indXHP);
1493  if (hasM_) {
1494  cM_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indXHP);
1495  }
1496  else {
1497  cM_XHP = cXHP;
1498  }
1499 
1500  std::vector<int> indHP(nevLocal_-blockSize_);
1501  for (int i=blockSize_; i<nevLocal_; i++) {
1502  indHP[i-blockSize_] = i;
1503  }
1504  cHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indHP);
1505  cK_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indHP);
1506  if (hasM_) {
1507  cM_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indHP);
1508  }
1509  else {
1510  cM_HP = cHP;
1511  }
1512 
1513  if (nevLocal_ == threeBlocks) {
1514  cP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock3);
1515  cK_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indblock3);
1516  if (hasM_) {
1517  cM_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indblock3);
1518  }
1519  else {
1520  cM_P = cP;
1521  }
1522  }
1523  }
1524 
1525  //
1526  //----------------------------------------
1527  // Form "local" mass and stiffness matrices
1528  //----------------------------------------
1529  {
1530  // We will form only the block upper triangular part of
1531  // [X H P]' K [X H P] and [X H P]' M [X H P]
1532  // Get the necessary views into KK and MM:
1533  // [--K1--] [--M1--]
1534  // KK = [ -K2-] MM = [ -M2-]
1535  // [ K3] [ M3]
1536  //
1537  // It is okay to declare a zero-area view of a Teuchos::SerialDenseMatrix
1538  //
1540  K1(Teuchos::View,KK,blockSize_,nevLocal_ ,0*blockSize_,0*blockSize_),
1541  K2(Teuchos::View,KK,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
1542  K3(Teuchos::View,KK,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_),
1543  M1(Teuchos::View,MM,blockSize_,nevLocal_ ,0*blockSize_,0*blockSize_),
1544  M2(Teuchos::View,MM,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
1545  M3(Teuchos::View,MM,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_);
1546  {
1547 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1548  Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1549 #endif
1550  MVT::MvTransMv( ONE, *cX, *cK_XHP, K1 );
1551  MVT::MvTransMv( ONE, *cX, *cM_XHP, M1 );
1552  MVT::MvTransMv( ONE, *cH, *cK_HP , K2 );
1553  MVT::MvTransMv( ONE, *cH, *cM_HP , M2 );
1554  if (nevLocal_ == threeBlocks) {
1555  MVT::MvTransMv( ONE, *cP, *cK_P, K3 );
1556  MVT::MvTransMv( ONE, *cP, *cM_P, M3 );
1557  }
1558  }
1559  }
1560  // below, we only need bases [X H P] and [H P] and friends
1561  // furthermore, we only need [H P] and friends if fullOrtho == false
1562  // clear the others now
1563  cX = Teuchos::null;
1564  cH = Teuchos::null;
1565  cP = Teuchos::null;
1566  cK_P = Teuchos::null;
1567  cM_P = Teuchos::null;
1568  if (fullOrtho_ == true) {
1569  cHP = Teuchos::null;
1570  cK_HP = Teuchos::null;
1571  cM_HP = Teuchos::null;
1572  }
1573 
1574  //
1575  //---------------------------------------------------
1576  // Perform a spectral decomposition of (KK,MM)
1577  //---------------------------------------------------
1578  //
1579  // Get pointers to relevant part of KK, MM and S for Rayleigh-Ritz analysis
1580  Teuchos::SerialDenseMatrix<int,ScalarType> lclKK(Teuchos::View,KK,nevLocal_,nevLocal_),
1581  lclMM(Teuchos::View,MM,nevLocal_,nevLocal_),
1582  lclS(Teuchos::View, S,nevLocal_,nevLocal_);
1583  {
1584 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1585  Teuchos::TimeMonitor lcltimer( *timerDS_ );
1586 #endif
1587  int localSize = nevLocal_;
1588  Utils::directSolver(localSize, lclKK, Teuchos::rcpFromRef(lclMM), lclS, theta_, nevLocal_, 0);
1589  // localSize tells directSolver() how big KK,MM are
1590  // however, directSolver() may choose to use only the principle submatrices of KK,MM
1591  // because of loss of MM-orthogonality in the projected eigenvectors
1592  // nevLocal_ tells us how much it used, telling us the effective localSize
1593  // (i.e., how much of KK,MM used by directSolver)
1594  // we will not tolerate any indefiniteness, and will throw an exception if it was
1595  // detected by directSolver
1596  //
1597  if (nevLocal_ != localSize) {
1598  // before throwing the exception, and thereby leaving iterate(), setup the views again
1599  // first, clear the const views
1600  cXHP = Teuchos::null;
1601  cK_XHP = Teuchos::null;
1602  cM_XHP = Teuchos::null;
1603  cHP = Teuchos::null;
1604  cK_HP = Teuchos::null;
1605  cM_HP = Teuchos::null;
1606  setupViews();
1607  }
1608  TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != localSize, LOBPCGRitzFailure,
1609  "Anasazi::LOBPCG::iterate(): indefiniteness detected in projected mass matrix." );
1610  }
1611 
1612  //
1613  //---------------------------------------------------
1614  // Sort the ritz values using the sort manager
1615  //---------------------------------------------------
1618  {
1619 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1620  Teuchos::TimeMonitor lcltimer( *timerSort_ );
1621 #endif
1622 
1623  std::vector<int> order(nevLocal_);
1624  //
1625  // Sort the first nevLocal_ values in theta_
1626  sm_->sort(theta_, Teuchos::rcpFromRef(order), nevLocal_); // don't catch exception
1627  //
1628  // Sort the primitive ritz vectors
1629  Utils::permuteVectors(order,lclS);
1630  }
1631 
1632  //
1633  //----------------------------------------
1634  // Compute coefficients for X and P under [X H P]
1635  //----------------------------------------
1636  // Before computing X,P, optionally perform orthogonalization per Hetmaniuk,Lehoucq paper
1637  // CX will be the coefficients of [X,H,P] for new X, CP for new P
1638  // The paper suggests orthogonalizing CP against CX and orthonormalizing CP, w.r.t. MM
1639  // Here, we will also orthonormalize CX.
1640  // This is accomplished using the Cholesky factorization of [CX CP]^H lclMM [CX CP]
1642  if (fullOrtho_) {
1643  // build orthonormal basis for ( 0 ) that is MM orthogonal to ( S11 )
1644  // ( S21 ) ( S21 )
1645  // ( S31 ) ( S31 )
1646  // Do this using Cholesky factorization of ( S11 0 )
1647  // ( S21 S21 )
1648  // ( S31 S31 )
1649  // ( S11 0 )
1650  // Build C = ( S21 S21 )
1651  // ( S31 S31 )
1652  Teuchos::SerialDenseMatrix<int,ScalarType> C(nevLocal_,twoBlocks),
1653  MMC(nevLocal_,twoBlocks),
1654  CMMC(twoBlocks ,twoBlocks);
1655 
1656  // first block of rows: ( S11 0 )
1657  for (int j=0; j<oneBlock; j++) {
1658  for (int i=0; i<oneBlock; i++) {
1659  // CX
1660  C(i,j) = lclS(i,j);
1661  // CP
1662  C(i,j+oneBlock) = ZERO;
1663  }
1664  }
1665  // second block of rows: (S21 S21)
1666  for (int j=0; j<oneBlock; j++) {
1667  for (int i=oneBlock; i<twoBlocks; i++) {
1668  // CX
1669  C(i,j) = lclS(i,j);
1670  // CP
1671  C(i,j+oneBlock) = lclS(i,j);
1672  }
1673  }
1674  // third block of rows
1675  if (nevLocal_ == threeBlocks) {
1676  for (int j=0; j<oneBlock; j++) {
1677  for (int i=twoBlocks; i<threeBlocks; i++) {
1678  // CX
1679  C(i,j) = lclS(i,j);
1680  // CP
1681  C(i,j+oneBlock) = lclS(i,j);
1682  }
1683  }
1684  }
1685 
1686  // compute MMC = lclMM*C
1687  {
1688  int teuchosret;
1689  teuchosret = MMC.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,C,ZERO);
1690  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1691  "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1692  // compute CMMC = C^H*MMC == C^H*lclMM*C
1693  teuchosret = CMMC.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,C,MMC,ZERO);
1694  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1695  "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1696  }
1697 
1698  // compute R (cholesky) of CMMC
1699  {
1700  int info;
1701  lapack.POTRF('U',twoBlocks,CMMC.values(),CMMC.stride(),&info);
1702  // our views ARE NOT currently in place; we must reestablish them before throwing an exception
1703  if (info != 0) {
1704  // Check symmetry of H'*K*H, this might be the first place a non-Hermitian operator will cause a problem.
1705  Teuchos::SerialDenseMatrix<int,ScalarType> K22(Teuchos::View,KK,blockSize_,blockSize_,1*blockSize_,1*blockSize_);
1707  K22_trans -= K22;
1708  MagnitudeType symNorm = K22_trans.normOne();
1709  MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
1710 
1711  cXHP = Teuchos::null;
1712  cHP = Teuchos::null;
1713  cK_XHP = Teuchos::null;
1714  cK_HP = Teuchos::null;
1715  cM_XHP = Teuchos::null;
1716  cM_HP = Teuchos::null;
1717  setupViews();
1718 
1719  std::string errMsg = "Anasazi::LOBPCG::iterate(): Cholesky factorization failed during full orthogonalization.";
1720  if ( symNorm < tol )
1721  {
1722  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, LOBPCGOrthoFailure, errMsg );
1723  }
1724  else
1725  {
1726  errMsg += " Check the operator A (or K), it may not be Hermitian!";
1727  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, LOBPCGOrthoFailure, errMsg );
1728  }
1729  }
1730  }
1731  // compute C = C inv(R)
1733  nevLocal_,twoBlocks,ONE,CMMC.values(),CMMC.stride(),C.values(),C.stride());
1734  // put C(:,0:oneBlock-1) into CX
1735  CX = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,0) );
1736  // put C(:,oneBlock:twoBlocks-1) into CP
1737  CP = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,oneBlock) );
1738 
1739  // check the results
1740  if (om_->isVerbosity( Debug ) ) {
1741  Teuchos::SerialDenseMatrix<int,ScalarType> tmp1(nevLocal_,oneBlock),
1742  tmp2(oneBlock,oneBlock);
1743  MagnitudeType tmp;
1744  int teuchosret;
1745  std::stringstream os;
1746  os.precision(2);
1747  os.setf(std::ios::scientific, std::ios::floatfield);
1748 
1749  os << " Checking Full Ortho: iteration " << iter_ << std::endl;
1750 
1751  // check CX^T MM CX == I
1752  // compute tmp1 = MM*CX
1753  teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CX,ZERO);
1754  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1755  "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1756  // compute tmp2 = CX^H*tmp1 == CX^H*MM*CX
1757  teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
1758  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1759  "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1760  // subtrace tmp2 - I == CX^H * MM * CX - I
1761  for (int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1762  tmp = tmp2.normFrobenius();
1763  os << " >> Error in CX^H MM CX == I : " << tmp << std::endl;
1764 
1765  // check CP^T MM CP == I
1766  // compute tmp1 = MM*CP
1767  teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
1768  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1769  "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1770  // compute tmp2 = CP^H*tmp1 == CP^H*MM*CP
1771  teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CP,tmp1,ZERO);
1772  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
1773  "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1774  // subtrace tmp2 - I == CP^H * MM * CP - I
1775  for (int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
1776  tmp = tmp2.normFrobenius();
1777  os << " >> Error in CP^H MM CP == I : " << tmp << std::endl;
1778 
1779  // check CX^T MM CP == 0
1780  // compute tmp1 = MM*CP
1781  teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
1782  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1783  // compute tmp2 = CX^H*tmp1 == CX^H*MM*CP
1784  teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
1785  TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
1786  // subtrace tmp2 == CX^H * MM * CP
1787  tmp = tmp2.normFrobenius();
1788  os << " >> Error in CX^H MM CP == 0 : " << tmp << std::endl;
1789 
1790  os << std::endl;
1791  om_->print(Debug,os.str());
1792  }
1793  }
1794  else {
1795  // [S11 ... ...]
1796  // S = [S21 ... ...]
1797  // [S31 ... ...]
1798  //
1799  // CX = [S11]
1800  // [S21]
1801  // [S31] -> X = [X H P] CX
1802  //
1803  // CP = [S21] -> P = [H P] CP
1804  // [S31]
1805  //
1806  CX = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_ ,oneBlock,0 ,0) );
1807  CP = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_-oneBlock,oneBlock,oneBlock,0) );
1808  }
1809 
1810  //
1811  //----------------------------------------
1812  // Compute new X and new P
1813  //----------------------------------------
1814  // Note: Use R as a temporary work space and (if full ortho) tmpMV as well
1815  {
1816 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1817  Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1818 #endif
1819 
1820  // if full ortho, then CX and CP are dense
1821  // we multiply [X H P]*CX into tmpMV
1822  // [X H P]*CP into R
1823  // then put V(:,firstblock) <- tmpMV
1824  // V(:,thirdblock) <- R
1825  //
1826  // if no full ortho, then [H P]*CP doesn't reference first block (X)
1827  // of V, so that we can modify it before computing P
1828  // so we multiply [X H P]*CX into R
1829  // V(:,firstblock) <- R
1830  // multiply [H P]*CP into R
1831  // V(:,thirdblock) <- R
1832  //
1833  // mutatis mutandis for K[XP] and M[XP]
1834  //
1835  // use SetBlock to do the assignments into V_
1836  //
1837  // in either case, views are only allowed to be overlapping
1838  // if they are const, and it should be assume that SetBlock
1839  // creates a view of the associated part
1840  //
1841  // we have from above const-pointers to [KM]XHP, [KM]HP and (if hasP) [KM]P
1842  //
1843  if (fullOrtho_) {
1844  // X,P
1845  MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*tmpmvec_);
1846  MVT::MvTimesMatAddMv(ONE,*cXHP,*CP,ZERO,*R_);
1847  cXHP = Teuchos::null;
1848  MVT::SetBlock(*tmpmvec_,indblock1,*V_);
1849  MVT::SetBlock(*R_ ,indblock3,*V_);
1850  // KX,KP
1851  MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*tmpmvec_);
1852  MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CP,ZERO,*R_);
1853  cK_XHP = Teuchos::null;
1854  MVT::SetBlock(*tmpmvec_,indblock1,*KV_);
1855  MVT::SetBlock(*R_ ,indblock3,*KV_);
1856  // MX,MP
1857  if (hasM_) {
1858  MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*tmpmvec_);
1859  MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CP,ZERO,*R_);
1860  cM_XHP = Teuchos::null;
1861  MVT::SetBlock(*tmpmvec_,indblock1,*MV_);
1862  MVT::SetBlock(*R_ ,indblock3,*MV_);
1863  }
1864  else {
1865  cM_XHP = Teuchos::null;
1866  }
1867  }
1868  else {
1869  // X,P
1870  MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*R_);
1871  cXHP = Teuchos::null;
1872  MVT::SetBlock(*R_,indblock1,*V_);
1873  MVT::MvTimesMatAddMv(ONE,*cHP,*CP,ZERO,*R_);
1874  cHP = Teuchos::null;
1875  MVT::SetBlock(*R_,indblock3,*V_);
1876  // KX,KP
1877  MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*R_);
1878  cK_XHP = Teuchos::null;
1879  MVT::SetBlock(*R_,indblock1,*KV_);
1880  MVT::MvTimesMatAddMv(ONE,*cK_HP,*CP,ZERO,*R_);
1881  cK_HP = Teuchos::null;
1882  MVT::SetBlock(*R_,indblock3,*KV_);
1883  // MX,MP
1884  if (hasM_) {
1885  MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*R_);
1886  cM_XHP = Teuchos::null;
1887  MVT::SetBlock(*R_,indblock1,*MV_);
1888  MVT::MvTimesMatAddMv(ONE,*cM_HP,*CP,ZERO,*R_);
1889  cM_HP = Teuchos::null;
1890  MVT::SetBlock(*R_,indblock3,*MV_);
1891  }
1892  else {
1893  cM_XHP = Teuchos::null;
1894  cM_HP = Teuchos::null;
1895  }
1896  }
1897  } // end timing block
1898  // done with coefficient matrices
1899  CX = Teuchos::null;
1900  CP = Teuchos::null;
1901 
1902  //
1903  // we now have a P direction
1904  hasP_ = true;
1905 
1906  // debugging check: all of our const views should have been cleared by now
1907  // if not, we have a logic error above
1908  TEUCHOS_TEST_FOR_EXCEPTION( cXHP != Teuchos::null || cK_XHP != Teuchos::null || cM_XHP != Teuchos::null
1909  || cHP != Teuchos::null || cK_HP != Teuchos::null || cM_HP != Teuchos::null
1910  || cP != Teuchos::null || cK_P != Teuchos::null || cM_P != Teuchos::null,
1911  std::logic_error,
1912  "Anasazi::BlockKrylovSchur::iterate(): const views were not all cleared! Something went wrong!" );
1913 
1914  //
1915  // recreate our const MV views of X,H,P and friends
1916  setupViews();
1917 
1918  //
1919  // Compute the new residuals, explicitly
1920  {
1921 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1922  Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1923 #endif
1924  MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1925  Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1926  for (int i = 0; i < blockSize_; i++) {
1927  T(i,i) = theta_[i];
1928  }
1929  MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
1930  }
1931 
1932  // R has been updated; mark the norms as out-of-date
1933  Rnorms_current_ = false;
1934  R2norms_current_ = false;
1935 
1936  // When required, monitor some orthogonalities
1937  if (om_->isVerbosity( Debug ) ) {
1938  // Check almost everything here
1939  CheckList chk;
1940  chk.checkX = true;
1941  chk.checkKX = true;
1942  chk.checkMX = true;
1943  chk.checkP = true;
1944  chk.checkKP = true;
1945  chk.checkMP = true;
1946  chk.checkR = true;
1947  om_->print( Debug, accuracyCheck(chk, ": after local update") );
1948  }
1949  else if (om_->isVerbosity( OrthoDetails )) {
1950  CheckList chk;
1951  chk.checkX = true;
1952  chk.checkP = true;
1953  chk.checkR = true;
1954  om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
1955  }
1956  } // end while (statusTest == false)
1957  }
1958 
1959 
1961  // compute/return residual M-norms
1962  template <class ScalarType, class MV, class OP>
1963  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1965  if (Rnorms_current_ == false) {
1966  // Update the residual norms
1967  orthman_->norm(*R_,Rnorms_);
1968  Rnorms_current_ = true;
1969  }
1970  return Rnorms_;
1971  }
1972 
1973 
1975  // compute/return residual 2-norms
1976  template <class ScalarType, class MV, class OP>
1977  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1979  if (R2norms_current_ == false) {
1980  // Update the residual 2-norms
1981  MVT::MvNorm(*R_,R2norms_);
1982  R2norms_current_ = true;
1983  }
1984  return R2norms_;
1985  }
1986 
1987 
1988 
1989 
1991  // Check accuracy, orthogonality, and other debugging stuff
1992  //
1993  // bools specify which tests we want to run (instead of running more than we actually care about)
1994  //
1995  // we don't bother checking the following because they are computed explicitly:
1996  // H == Prec*R
1997  // KH == K*H
1998  //
1999  //
2000  // checkX : X orthonormal
2001  // orthogonal to auxvecs
2002  // checkMX: check MX == M*X
2003  // checkKX: check KX == K*X
2004  // checkP : if fullortho P orthonormal and orthogonal to X
2005  // orthogonal to auxvecs
2006  // checkMP: check MP == M*P
2007  // checkKP: check KP == K*P
2008  // checkH : if fullortho H orthonormal and orthogonal to X and P
2009  // orthogonal to auxvecs
2010  // checkMH: check MH == M*H
2011  // checkR : check R orthogonal to X
2012  // checkQ : check that auxiliary vectors are actually orthonormal
2013  //
2014  // TODO:
2015  // add checkTheta
2016  //
2017  template <class ScalarType, class MV, class OP>
2018  std::string LOBPCG<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
2019  {
2020  using std::endl;
2021 
2022  std::stringstream os;
2023  os.precision(2);
2024  os.setf(std::ios::scientific, std::ios::floatfield);
2025  MagnitudeType tmp;
2026 
2027  os << " Debugging checks: iteration " << iter_ << where << endl;
2028 
2029  // X and friends
2030  if (chk.checkX && initialized_) {
2031  tmp = orthman_->orthonormError(*X_);
2032  os << " >> Error in X^H M X == I : " << tmp << endl;
2033  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2034  tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
2035  os << " >> Error in X^H M Q[" << i << "] == 0 : " << tmp << endl;
2036  }
2037  }
2038  if (chk.checkMX && hasM_ && initialized_) {
2039  tmp = Utils::errorEquality(*X_, *MX_, MOp_);
2040  os << " >> Error in MX == M*X : " << tmp << endl;
2041  }
2042  if (chk.checkKX && initialized_) {
2043  tmp = Utils::errorEquality(*X_, *KX_, Op_);
2044  os << " >> Error in KX == K*X : " << tmp << endl;
2045  }
2046 
2047  // P and friends
2048  if (chk.checkP && hasP_ && initialized_) {
2049  if (fullOrtho_) {
2050  tmp = orthman_->orthonormError(*P_);
2051  os << " >> Error in P^H M P == I : " << tmp << endl;
2052  tmp = orthman_->orthogError(*P_,*X_);
2053  os << " >> Error in P^H M X == 0 : " << tmp << endl;
2054  }
2055  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2056  tmp = orthman_->orthogError(*P_,*auxVecs_[i]);
2057  os << " >> Error in P^H M Q[" << i << "] == 0 : " << tmp << endl;
2058  }
2059  }
2060  if (chk.checkMP && hasM_ && hasP_ && initialized_) {
2061  tmp = Utils::errorEquality(*P_, *MP_, MOp_);
2062  os << " >> Error in MP == M*P : " << tmp << endl;
2063  }
2064  if (chk.checkKP && hasP_ && initialized_) {
2065  tmp = Utils::errorEquality(*P_, *KP_, Op_);
2066  os << " >> Error in KP == K*P : " << tmp << endl;
2067  }
2068 
2069  // H and friends
2070  if (chk.checkH && initialized_) {
2071  if (fullOrtho_) {
2072  tmp = orthman_->orthonormError(*H_);
2073  os << " >> Error in H^H M H == I : " << tmp << endl;
2074  tmp = orthman_->orthogError(*H_,*X_);
2075  os << " >> Error in H^H M X == 0 : " << tmp << endl;
2076  if (hasP_) {
2077  tmp = orthman_->orthogError(*H_,*P_);
2078  os << " >> Error in H^H M P == 0 : " << tmp << endl;
2079  }
2080  }
2081  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2082  tmp = orthman_->orthogError(*H_,*auxVecs_[i]);
2083  os << " >> Error in H^H M Q[" << i << "] == 0 : " << tmp << endl;
2084  }
2085  }
2086  if (chk.checkMH && hasM_ && initialized_) {
2087  tmp = Utils::errorEquality(*H_, *MH_, MOp_);
2088  os << " >> Error in MH == M*H : " << tmp << endl;
2089  }
2090 
2091  // R: this is not M-orthogonality, but standard euclidean orthogonality
2092  if (chk.checkR && initialized_) {
2093  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
2094  MVT::MvTransMv(ONE,*X_,*R_,xTx);
2095  tmp = xTx.normFrobenius();
2096  MVT::MvTransMv(ONE,*R_,*R_,xTx);
2097  double normR = xTx.normFrobenius();
2098  os << " >> RelError in X^H R == 0: " << tmp/normR << endl;
2099  }
2100 
2101  // Q
2102  if (chk.checkQ) {
2103  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
2104  tmp = orthman_->orthonormError(*auxVecs_[i]);
2105  os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << tmp << endl;
2106  for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
2107  tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
2108  os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << tmp << endl;
2109  }
2110  }
2111  }
2112 
2113  os << endl;
2114 
2115  return os.str();
2116  }
2117 
2118 
2120  // Print the current status of the solver
2121  template <class ScalarType, class MV, class OP>
2122  void
2124  {
2125  using std::endl;
2126 
2127  os.setf(std::ios::scientific, std::ios::floatfield);
2128  os.precision(6);
2129  os <<endl;
2130  os <<"================================================================================" << endl;
2131  os << endl;
2132  os <<" LOBPCG Solver Status" << endl;
2133  os << endl;
2134  os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
2135  os <<"The number of iterations performed is " << iter_ << endl;
2136  os <<"The current block size is " << blockSize_ << endl;
2137  os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
2138  os <<"The number of operations Op*x is " << count_ApplyOp_ << endl;
2139  os <<"The number of operations M*x is " << count_ApplyM_ << endl;
2140  os <<"The number of operations Prec*x is " << count_ApplyPrec_ << endl;
2141 
2142  os.setf(std::ios_base::right, std::ios_base::adjustfield);
2143 
2144  if (initialized_) {
2145  os << endl;
2146  os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
2147  os << std::setw(20) << "Eigenvalue"
2148  << std::setw(20) << "Residual(M)"
2149  << std::setw(20) << "Residual(2)"
2150  << endl;
2151  os <<"--------------------------------------------------------------------------------"<<endl;
2152  for (int i=0; i<blockSize_; i++) {
2153  os << std::setw(20) << theta_[i];
2154  if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
2155  else os << std::setw(20) << "not current";
2156  if (R2norms_current_) os << std::setw(20) << R2norms_[i];
2157  else os << std::setw(20) << "not current";
2158  os << endl;
2159  }
2160  }
2161  os <<"================================================================================" << endl;
2162  os << endl;
2163  }
2164 
2166  // are we initialized or not?
2167  template <class ScalarType, class MV, class OP>
2169  return initialized_;
2170  }
2171 
2172 
2174  // is P valid or not?
2175  template <class ScalarType, class MV, class OP>
2177  return hasP_;
2178  }
2179 
2181  // is full orthogonalization enabled or not?
2182  template <class ScalarType, class MV, class OP>
2184  return(fullOrtho_);
2185  }
2186 
2187 
2189  // return the current auxilliary vectors
2190  template <class ScalarType, class MV, class OP>
2192  return auxVecs_;
2193  }
2194 
2196  // return the current block size
2197  template <class ScalarType, class MV, class OP>
2199  return(blockSize_);
2200  }
2201 
2203  // return the current eigenproblem
2204  template <class ScalarType, class MV, class OP>
2206  return(*problem_);
2207  }
2208 
2209 
2211  // return the max subspace dimension
2212  template <class ScalarType, class MV, class OP>
2214  return 3*blockSize_;
2215  }
2216 
2218  // return the current subspace dimension
2219  template <class ScalarType, class MV, class OP>
2221  if (!initialized_) return 0;
2222  return nevLocal_;
2223  }
2224 
2225 
2227  // return the current ritz residual norms
2228  template <class ScalarType, class MV, class OP>
2229  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2231  {
2232  return this->getRes2Norms();
2233  }
2234 
2235 
2237  // return the current compression indices
2238  template <class ScalarType, class MV, class OP>
2240  std::vector<int> ret(nevLocal_,0);
2241  return ret;
2242  }
2243 
2244 
2246  // return the current ritz values
2247  template <class ScalarType, class MV, class OP>
2248  std::vector<Value<ScalarType> > LOBPCG<ScalarType,MV,OP>::getRitzValues() {
2249  std::vector<Value<ScalarType> > ret(nevLocal_);
2250  for (int i=0; i<nevLocal_; i++) {
2251  ret[i].realpart = theta_[i];
2252  ret[i].imagpart = ZERO;
2253  }
2254  return ret;
2255  }
2256 
2258  // Set a new StatusTest for the solver.
2259  template <class ScalarType, class MV, class OP>
2261  TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
2262  "Anasazi::LOBPCG::setStatusTest() was passed a null StatusTest.");
2263  tester_ = test;
2264  }
2265 
2267  // Get the current StatusTest used by the solver.
2268  template <class ScalarType, class MV, class OP>
2270  return tester_;
2271  }
2272 
2274  // return the current ritz vectors
2275  template <class ScalarType, class MV, class OP>
2277  return X_;
2278  }
2279 
2280 
2282  // reset the iteration counter
2283  template <class ScalarType, class MV, class OP>
2285  iter_=0;
2286  }
2287 
2288 
2290  // return the number of iterations
2291  template <class ScalarType, class MV, class OP>
2293  return(iter_);
2294  }
2295 
2296 
2298  // return the state
2299  template <class ScalarType, class MV, class OP>
2302  state.V = V_;
2303  state.KV = KV_;
2304  state.X = X_;
2305  state.KX = KX_;
2306  state.P = P_;
2307  state.KP = KP_;
2308  state.H = H_;
2309  state.KH = KH_;
2310  state.R = R_;
2311  state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_));
2312  if (hasM_) {
2313  state.MV = MV_;
2314  state.MX = MX_;
2315  state.MP = MP_;
2316  state.MH = MH_;
2317  }
2318  else {
2319  state.MX = Teuchos::null;
2320  state.MP = Teuchos::null;
2321  state.MH = Teuchos::null;
2322  }
2323  return state;
2324  }
2325 
2326 } // end Anasazi namespace
2327 
2328 #endif // ANASAZI_LOBPCG_HPP
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
ScalarType * values() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
LOBPCG(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)
LOBPCG constructor with eigenproblem, solver utilities, and parameter list of solver options...
This class defines the interface required by an eigensolver and status test class to compute solution...
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
Declaration of basic traits for the multivector type.
T & get(const std::string &name, T def_value)
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb) const
LOBPCGState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration...
Teuchos::RCP< const MultiVector > V
The current test basis.
virtual ~LOBPCG()
LOBPCG destructor
Virtual base class which defines basic traits for the operator type.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For LOBPCG, this always returns 3*getBlo...
Teuchos::RCP< const MultiVector > H
The current preconditioned residual vectors.
Teuchos::RCP< const MultiVector > P
The current search direction.
An exception class parent to all Anasazi exceptions.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
Teuchos::RCP< const MultiVector > MV
The image of the current test basis under M, or Teuchos::null if M was not specified.
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...
Teuchos::RCP< const MultiVector > R
The current residual vectors.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Anasazi&#39;s templated, static class providing utilities for the solvers.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this eigenproblem.
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Output managers remove the need for the eigensolver to know any information about the required output...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms.
Teuchos::RCP< const MultiVector > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Teuchos::RCP< const MultiVector > KX
The image of the current eigenvectors under K.
bool getFullOrtho() const
Determine if the LOBPCG iteration is using full orthogonality.
Teuchos::RCP< const MultiVector > KH
The image of the current preconditioned residual vectors under K.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
std::vector< int > getRitzIndex()
Get the index used for extracting Ritz vectors from getRitzVectors().
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Teuchos::RCP< const MultiVector > MH
The image of the current preconditioned residual vectors under M, or Teuchos::null if M was not speci...
Teuchos::RCP< const MultiVector > KP
The image of the current search direction under K.
ScalarTraits< ScalarType >::magnitudeType normOne() const
void push_back(const value_type &x)
int getNumIters() const
Get the current iteration count.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
void iterate()
This method performs LOBPCG iterations until the status test indicates the need to stop or an error o...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
LOBPCGRitzFailure is thrown when the LOBPCG solver is unable to continue a call to LOBPCG::iterate() ...
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values from the previous iteration.
Teuchos::RCP< const MultiVector > KV
The image of the current test basis under K.
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
Types and exceptions used within Anasazi solvers and interfaces.
void resetNumIters()
Reset the iteration count.
OrdinalType stride() const
Teuchos::RCP< const MultiVector > MP
The image of the current search direction under M, or Teuchos::null if M was not specified.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
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.
LOBPCGOrthoFailure is thrown when an orthogonalization attempt fails.
LOBPCGInitFailure is thrown when the LOBPCG solver is unable to generate an initial iterate in the LO...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
void setFullOrtho(bool fullOrtho)
Instruct the LOBPCG iteration to use full orthogonality.
Structure to contain pointers to Anasazi state variables.
Teuchos::RCP< const MultiVector > X
The current eigenvectors.
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors from the previous iteration.
bool hasP()
Indicates whether the search direction given by getState() is valid.
void POTRF(const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *info) const
Class which provides internal utilities for the Anasazi solvers.