Anasazi  Version of the Day
AnasaziBlockKrylovSchur.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
33 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_HPP
34 #define ANASAZI_BLOCK_KRYLOV_SCHUR_HPP
35 
36 #include "AnasaziTypes.hpp"
37 
38 #include "AnasaziEigensolver.hpp"
41 #include "Teuchos_ScalarTraits.hpp"
42 #include "AnasaziHelperTraits.hpp"
43 
44 #include "AnasaziOrthoManager.hpp"
45 
46 #include "Teuchos_LAPACK.hpp"
47 #include "Teuchos_BLAS.hpp"
48 #include "Teuchos_SerialDenseMatrix.hpp"
49 #include "Teuchos_ParameterList.hpp"
50 #include "Teuchos_TimeMonitor.hpp"
51 
66 namespace Anasazi {
67 
69 
70 
75  template <class ScalarType, class MulVec>
81  int curDim;
83  Teuchos::RCP<const MulVec> V;
89  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > H;
91  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > S;
93  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > Q;
94  BlockKrylovSchurState() : curDim(0), V(Teuchos::null),
95  H(Teuchos::null), S(Teuchos::null),
96  Q(Teuchos::null) {}
97  };
98 
100 
102 
103 
117  BlockKrylovSchurInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
118  {}};
119 
127  BlockKrylovSchurOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
128  {}};
129 
131 
132 
133  template <class ScalarType, class MV, class OP>
134  class BlockKrylovSchur : public Eigensolver<ScalarType,MV,OP> {
135  public:
137 
138 
150  BlockKrylovSchur( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
151  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
152  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
153  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
154  const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
155  Teuchos::ParameterList &params
156  );
157 
159  virtual ~BlockKrylovSchur() {};
161 
162 
164 
165 
187  void iterate();
188 
210  void initialize(BlockKrylovSchurState<ScalarType,MV>& state);
211 
215  void initialize();
216 
225  bool isInitialized() const { return initialized_; }
226 
236  state.curDim = curDim_;
237  state.V = V_;
238  state.H = H_;
239  state.Q = Q_;
240  state.S = schurH_;
241  return state;
242  }
243 
245 
246 
248 
249 
251  int getNumIters() const { return(iter_); }
252 
254  void resetNumIters() { iter_=0; }
255 
263  Teuchos::RCP<const MV> getRitzVectors() { return ritzVectors_; }
264 
272  std::vector<Value<ScalarType> > getRitzValues() {
273  std::vector<Value<ScalarType> > ret = ritzValues_;
274  ret.resize(ritzIndex_.size());
275  return ret;
276  }
277 
283  std::vector<int> getRitzIndex() { return ritzIndex_; }
284 
289  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms() {
290  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
291  return ret;
292  }
293 
298  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms() {
299  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
300  return ret;
301  }
302 
307  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms() {
308  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
309  ret.resize(ritzIndex_.size());
310  return ret;
311  }
312 
314 
316 
317 
319  void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
320 
322  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
323 
325  const Eigenproblem<ScalarType,MV,OP>& getProblem() const { return(*problem_); };
326 
333  void setSize(int blockSize, int numBlocks);
334 
336  void setBlockSize(int blockSize);
337 
339  void setStepSize(int stepSize);
340 
342  void setNumRitzVectors(int numRitzVecs);
343 
345  int getStepSize() const { return(stepSize_); }
346 
348  int getBlockSize() const { return(blockSize_); }
349 
351  int getNumRitzVectors() const { return(numRitzVecs_); }
352 
358  int getCurSubspaceDim() const {
359  if (!initialized_) return 0;
360  return curDim_;
361  }
362 
364  int getMaxSubspaceDim() const { return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
365 
366 
379  void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
380 
382  Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const {return auxVecs_;}
383 
385 
387 
388 
390  void currentStatus(std::ostream &os);
391 
393 
395 
396 
398  bool isRitzVecsCurrent() const { return ritzVecsCurrent_; }
399 
401  bool isRitzValsCurrent() const { return ritzValsCurrent_; }
402 
404  bool isSchurCurrent() const { return schurCurrent_; }
405 
407 
409 
410 
412  void computeRitzVectors();
413 
415  void computeRitzValues();
416 
418  void computeSchurForm( const bool sort = true );
419 
421 
422  private:
423  //
424  // Convenience typedefs
425  //
428  typedef Teuchos::ScalarTraits<ScalarType> SCT;
429  typedef typename SCT::magnitudeType MagnitudeType;
430  typedef typename std::vector<ScalarType>::iterator STiter;
431  typedef typename std::vector<MagnitudeType>::iterator MTiter;
432  const MagnitudeType MT_ONE;
433  const MagnitudeType MT_ZERO;
434  const MagnitudeType NANVAL;
435  const ScalarType ST_ONE;
436  const ScalarType ST_ZERO;
437  //
438  // Internal structs
439  //
440  struct CheckList {
441  bool checkV;
442  bool checkArn;
443  bool checkAux;
444  CheckList() : checkV(false), checkArn(false), checkAux(false) {};
445  };
446  //
447  // Internal methods
448  //
449  std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
450  void sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
451  Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
452  std::vector<int>& order );
453  //
454  // Classes inputed through constructor that define the eigenproblem to be solved.
455  //
456  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
457  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
458  const Teuchos::RCP<OutputManager<ScalarType> > om_;
459  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
460  const Teuchos::RCP<OrthoManager<ScalarType,MV> > orthman_;
461  //
462  // Information obtained from the eigenproblem
463  //
464  Teuchos::RCP<const OP> Op_;
465  //
466  // Internal timers
467  //
468  Teuchos::RCP<Teuchos::Time> timerOp_, timerSortRitzVal_,
469  timerCompSF_, timerSortSF_,
470  timerCompRitzVec_, timerOrtho_;
471  //
472  // Counters
473  //
474  int count_ApplyOp_;
475 
476  //
477  // Algorithmic parameters.
478  //
479  // blockSize_ is the solver block size; it controls the number of eigenvectors that
480  // we compute, the number of residual vectors that we compute, and therefore the number
481  // of vectors added to the basis on each iteration.
482  int blockSize_;
483  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
484  int numBlocks_;
485  // stepSize_ dictates how many iterations are performed before eigenvectors and eigenvalues
486  // are computed again
487  int stepSize_;
488 
489  //
490  // Current solver state
491  //
492  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
493  // is capable of running; _initialize is controlled by the initialize() member method
494  // For the implications of the state of initialized_, please see documentation for initialize()
495  bool initialized_;
496  //
497  // curDim_ reflects how much of the current basis is valid
498  // NOTE: for Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_
499  // for non-Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_ + 1
500  // this also tells us how many of the values in _theta are valid Ritz values
501  int curDim_;
502  //
503  // State Multivecs
504  Teuchos::RCP<MV> ritzVectors_, V_;
505  int numRitzVecs_;
506  //
507  // Projected matrices
508  // H_ : Projected matrix from the Krylov-Schur factorization AV = VH + FB^T
509  //
510  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
511  //
512  // Schur form of Projected matrices (these are only updated when the Ritz values/vectors are updated).
513  // schurH_: Schur form reduction of H
514  // Q_: Schur vectors of H
515  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > schurH_;
516  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Q_;
517  //
518  // Auxiliary vectors
519  Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
520  int numAuxVecs_;
521  //
522  // Number of iterations that have been performed.
523  int iter_;
524  //
525  // State flags
526  bool ritzVecsCurrent_, ritzValsCurrent_, schurCurrent_;
527  //
528  // Current eigenvalues, residual norms
529  std::vector<Value<ScalarType> > ritzValues_;
530  std::vector<MagnitudeType> ritzResiduals_;
531  //
532  // Current index vector for Ritz values and vectors
533  std::vector<int> ritzIndex_; // computed by BKS
534  std::vector<int> ritzOrder_; // returned from sort manager
535  //
536  // Number of Ritz pairs to be printed upon output, if possible
537  int numRitzPrint_;
538  };
539 
540 
542  // Constructor
543  template <class ScalarType, class MV, class OP>
545  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
546  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
547  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
548  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
549  const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
550  Teuchos::ParameterList &params
551  ) :
552  MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
553  MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
554  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
555  ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
556  ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
557  // problem, tools
558  problem_(problem),
559  sm_(sorter),
560  om_(printer),
561  tester_(tester),
562  orthman_(ortho),
563  // timers, counters
564 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
565  timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Operation Op*x")),
566  timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Ritz values")),
567  timerCompSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Schur form")),
568  timerSortSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Schur form")),
569  timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Ritz vectors")),
570  timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Orthogonalization")),
571 #endif
572  count_ApplyOp_(0),
573  // internal data
574  blockSize_(0),
575  numBlocks_(0),
576  stepSize_(0),
577  initialized_(false),
578  curDim_(0),
579  numRitzVecs_(0),
580  auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
581  numAuxVecs_(0),
582  iter_(0),
583  ritzVecsCurrent_(false),
584  ritzValsCurrent_(false),
585  schurCurrent_(false),
586  numRitzPrint_(0)
587  {
588  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
589  "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
590  TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
591  "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
592  TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
593  "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
594  TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
595  "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
596  TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
597  "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
598  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
599  "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
600  TEUCHOS_TEST_FOR_EXCEPTION(sorter == Teuchos::null,std::invalid_argument,
601  "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
602  TEUCHOS_TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
603  "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
604  TEUCHOS_TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
605  "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
606  TEUCHOS_TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
607  "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
608 
609  // Get problem operator
610  Op_ = problem_->getOperator();
611 
612  // get the step size
613  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Step Size"), std::invalid_argument,
614  "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
615  int ss = params.get("Step Size",numBlocks_);
616  setStepSize(ss);
617 
618  // set the block size and allocate data
619  int bs = params.get("Block Size", 1);
620  int nb = params.get("Num Blocks", 3*problem_->getNEV());
621  setSize(bs,nb);
622 
623  // get the number of Ritz vectors to compute and allocate data.
624  // --> if this parameter is not specified in the parameter list, then it's assumed that no Ritz vectors will be computed.
625  int numRitzVecs = params.get("Number of Ritz Vectors", 0);
626  setNumRitzVectors( numRitzVecs );
627 
628  // get the number of Ritz values to print out when currentStatus is called.
629  numRitzPrint_ = params.get("Print Number of Ritz Values", bs);
630  }
631 
632 
634  // Set the block size
635  // This simply calls setSize(), modifying the block size while retaining the number of blocks.
636  template <class ScalarType, class MV, class OP>
638  {
639  setSize(blockSize,numBlocks_);
640  }
641 
642 
644  // Set the step size.
645  template <class ScalarType, class MV, class OP>
647  {
648  TEUCHOS_TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
649  stepSize_ = stepSize;
650  }
651 
653  // Set the number of Ritz vectors to compute.
654  template <class ScalarType, class MV, class OP>
656  {
657  // This routine only allocates space; it doesn't not perform any computation
658  // any change in size will invalidate the state of the solver.
659 
660  TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
661 
662  // Check to see if the number of requested Ritz vectors has changed.
663  if (numRitzVecs != numRitzVecs_) {
664  if (numRitzVecs) {
665  ritzVectors_ = Teuchos::null;
666  ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
667  } else {
668  ritzVectors_ = Teuchos::null;
669  }
670  numRitzVecs_ = numRitzVecs;
671  ritzVecsCurrent_ = false;
672  }
673  }
674 
676  // Set the block size and make necessary adjustments.
677  template <class ScalarType, class MV, class OP>
678  void BlockKrylovSchur<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
679  {
680  // This routine only allocates space; it doesn't not perform any computation
681  // any change in size will invalidate the state of the solver.
682 
683  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
684  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
685  if (blockSize == blockSize_ && numBlocks == numBlocks_) {
686  // do nothing
687  return;
688  }
689 
690  blockSize_ = blockSize;
691  numBlocks_ = numBlocks;
692 
693  Teuchos::RCP<const MV> tmp;
694  // grab some Multivector to Clone
695  // in practice, getInitVec() should always provide this, but it is possible to use a
696  // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
697  // in case of that strange scenario, we will try to Clone from V_; first resort to getInitVec(),
698  // because we would like to clear the storage associated with V_ so we have room for the new V_
699  if (problem_->getInitVec() != Teuchos::null) {
700  tmp = problem_->getInitVec();
701  }
702  else {
703  tmp = V_;
704  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
705  "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
706  }
707 
708 
710  // blockSize*numBlocks dependent
711  //
712  int newsd;
713  if (problem_->isHermitian()) {
714  newsd = blockSize_*numBlocks_;
715  } else {
716  newsd = blockSize_*numBlocks_+1;
717  }
718  // check that new size is valid
719  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(newsd) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
720  "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
721 
722  ritzValues_.resize(newsd);
723  ritzResiduals_.resize(newsd,MT_ONE);
724  ritzOrder_.resize(newsd);
725  V_ = Teuchos::null;
726  V_ = MVT::Clone(*tmp,newsd+blockSize_);
727  H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd+blockSize_,newsd) );
728  Q_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
729 
730  initialized_ = false;
731  curDim_ = 0;
732  }
733 
734 
736  // Set the auxiliary vectors
737  template <class ScalarType, class MV, class OP>
738  void BlockKrylovSchur<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
739  typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
740 
741  // set new auxiliary vectors
742  auxVecs_ = auxvecs;
743 
744  if (om_->isVerbosity( Debug ) ) {
745  // Check almost everything here
746  CheckList chk;
747  chk.checkAux = true;
748  om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
749  }
750 
751  numAuxVecs_ = 0;
752  for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
753  numAuxVecs_ += MVT::GetNumberVecs(**i);
754  }
755 
756  // If the solver has been initialized, X and P are not necessarily orthogonal to new auxiliary vectors
757  if (numAuxVecs_ > 0 && initialized_) {
758  initialized_ = false;
759  }
760  }
761 
763  /* Initialize the state of the solver
764  *
765  * POST-CONDITIONS:
766  *
767  * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
768  *
769  */
770 
771  template <class ScalarType, class MV, class OP>
773  {
774  // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
775 
776  std::vector<int> bsind(blockSize_);
777  for (int i=0; i<blockSize_; i++) bsind[i] = i;
778 
779  // in BlockKrylovSchur, V and H are required.
780  // if either doesn't exist, then we will start with the initial vector.
781  //
782  // inconsistent multivectors widths and lengths will not be tolerated, and
783  // will be treated with exceptions.
784  //
785  std::string errstr("Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
786 
787  // set up V,H: if the user doesn't specify both of these these,
788  // we will start over with the initial vector.
789  if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
790 
791  // initialize V_,H_, and curDim_
792 
793  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
794  std::invalid_argument, errstr );
795  if (newstate.V != V_) {
796  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
797  std::invalid_argument, errstr );
798  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) > getMaxSubspaceDim(),
799  std::invalid_argument, errstr );
800  }
801  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > getMaxSubspaceDim(),
802  std::invalid_argument, errstr );
803 
804  curDim_ = newstate.curDim;
805  int lclDim = MVT::GetNumberVecs(*newstate.V);
806 
807  // check size of H
808  TEUCHOS_TEST_FOR_EXCEPTION(newstate.H->numRows() < curDim_ || newstate.H->numCols() < curDim_, std::invalid_argument, errstr);
809 
810  if (curDim_ == 0 && lclDim > blockSize_) {
811  om_->stream(Warnings) << "Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
812  << "The block size however is only " << blockSize_ << std::endl
813  << "The last " << lclDim - blockSize_ << " vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
814  }
815 
816 
817  // copy basis vectors from newstate into V
818  if (newstate.V != V_) {
819  std::vector<int> nevind(lclDim);
820  for (int i=0; i<lclDim; i++) nevind[i] = i;
821  MVT::SetBlock(*newstate.V,nevind,*V_);
822  }
823 
824  // put data into H_, make sure old information is not still hanging around.
825  if (newstate.H != H_) {
826  H_->putScalar( ST_ZERO );
827  Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H,curDim_+blockSize_,curDim_);
828  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
829  lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_,curDim_+blockSize_,curDim_) );
830  lclH->assign(newH);
831 
832  // done with local pointers
833  lclH = Teuchos::null;
834  }
835 
836  }
837  else {
838  // user did not specify a basis V
839  // get vectors from problem or generate something, projectAndNormalize, call initialize() recursively
840  Teuchos::RCP<const MV> ivec = problem_->getInitVec();
841  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
842  "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
843 
844  int lclDim = MVT::GetNumberVecs(*ivec);
845  bool userand = false;
846  if (lclDim < blockSize_) {
847  // we need at least blockSize_ vectors
848  // use a random multivec
849  userand = true;
850  }
851 
852  if (userand) {
853  // make an index
854  std::vector<int> dimind2(lclDim);
855  for (int i=0; i<lclDim; i++) { dimind2[i] = i; }
856 
857  // alloc newV as a view of the first block of V
858  Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,dimind2);
859 
860  // copy the initial vectors into the first lclDim vectors of V
861  MVT::SetBlock(*ivec,dimind2,*newV1);
862 
863  // resize / reinitialize the index vector
864  dimind2.resize(blockSize_-lclDim);
865  for (int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
866 
867  // initialize the rest of the vectors with random vectors
868  Teuchos::RCP<MV> newV2 = MVT::CloneViewNonConst(*V_,dimind2);
869  MVT::MvRandom(*newV2);
870  }
871  else {
872  // alloc newV as a view of the first block of V
873  Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,bsind);
874 
875  // get a view of the first block of initial vectors
876  Teuchos::RCP<const MV> ivecV = MVT::CloneView(*ivec,bsind);
877 
878  // assign ivec to first part of newV
879  MVT::SetBlock(*ivecV,bsind,*newV1);
880  }
881 
882  // get pointer into first block of V
883  Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*V_,bsind);
884 
885  // remove auxVecs from newV and normalize newV
886  if (auxVecs_.size() > 0) {
887 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
888  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
889 #endif
890 
891  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
892  int rank = orthman_->projectAndNormalize(*newV,auxVecs_);
893  TEUCHOS_TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
894  "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
895  }
896  else {
897 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
898  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
899 #endif
900 
901  int rank = orthman_->normalize(*newV);
902  TEUCHOS_TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
903  "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
904  }
905 
906  // set curDim
907  curDim_ = 0;
908 
909  // clear pointer
910  newV = Teuchos::null;
911  }
912 
913  // The Ritz vectors/values and Schur form are no longer current.
914  ritzVecsCurrent_ = false;
915  ritzValsCurrent_ = false;
916  schurCurrent_ = false;
917 
918  // the solver is initialized
919  initialized_ = true;
920 
921  if (om_->isVerbosity( Debug ) ) {
922  // Check almost everything here
923  CheckList chk;
924  chk.checkV = true;
925  chk.checkArn = true;
926  chk.checkAux = true;
927  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
928  }
929 
930  // Print information on current status
931  if (om_->isVerbosity(Debug)) {
932  currentStatus( om_->stream(Debug) );
933  }
934  else if (om_->isVerbosity(IterationDetails)) {
935  currentStatus( om_->stream(IterationDetails) );
936  }
937  }
938 
939 
941  // initialize the solver with default state
942  template <class ScalarType, class MV, class OP>
944  {
946  initialize(empty);
947  }
948 
949 
951  // Perform BlockKrylovSchur iterations until the StatusTest tells us to stop.
952  template <class ScalarType, class MV, class OP>
954  {
955  //
956  // Allocate/initialize data structures
957  //
958  if (initialized_ == false) {
959  initialize();
960  }
961 
962  // Compute the current search dimension.
963  // If the problem is non-Hermitian and the blocksize is one, let the solver use the extra vector.
964  int searchDim = blockSize_*numBlocks_;
965  if (problem_->isHermitian() == false) {
966  searchDim++;
967  }
968 
970  // iterate until the status test tells us to stop.
971  //
972  // also break if our basis is full
973  //
974  while (tester_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
975 
976  iter_++;
977 
978  // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
979  int lclDim = curDim_ + blockSize_;
980 
981  // Get the current part of the basis.
982  std::vector<int> curind(blockSize_);
983  for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
984  Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
985 
986  // Get a view of the previous vectors
987  // this is used for orthogonalization and for computing V^H K H
988  for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
989  Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
990  // om_->stream(Debug) << "Vprev: " << std::endl;
991  // MVT::MvPrint(*Vprev,om_->stream(Debug));
992 
993  // Compute the next vector in the Krylov basis: Vnext = Op*Vprev
994  {
995 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
996  Teuchos::TimeMonitor lcltimer( *timerOp_ );
997 #endif
998  OPT::Apply(*Op_,*Vprev,*Vnext);
999  count_ApplyOp_ += blockSize_;
1000  }
1001  // om_->stream(Debug) << "Vnext: " << std::endl;
1002  // MVT::MvPrint(*Vnext,om_->stream(Debug));
1003  Vprev = Teuchos::null;
1004 
1005  // Remove all previous Krylov-Schur basis vectors and auxVecs from Vnext
1006  {
1007 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1008  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1009 #endif
1010 
1011  // Get a view of all the previous vectors
1012  std::vector<int> prevind(lclDim);
1013  for (int i=0; i<lclDim; i++) { prevind[i] = i; }
1014  Vprev = MVT::CloneView(*V_,prevind);
1015  Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
1016 
1017  // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
1018  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
1019  subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
1020  ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
1021  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
1022  AsubH.append( subH );
1023 
1024  // Add the auxiliary vectors to the current basis vectors if any exist
1025  if (auxVecs_.size() > 0) {
1026  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1027  AVprev.append( auxVecs_[i] );
1028  AsubH.append( Teuchos::null );
1029  }
1030  }
1031 
1032  // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
1033  // om_->stream(Debug) << "V before ortho: " << std::endl;
1034  // MVT::MvPrint(*Vprev,om_->stream(Debug));
1035  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
1036  subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
1037  ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
1038  int rank = orthman_->projectAndNormalize(*Vnext,AVprev,AsubH,subR);
1039  // om_->stream(Debug) << "Vnext after ortho: " << std::endl;
1040  // MVT::MvPrint(*Vnext,om_->stream(Debug));
1041  // om_->stream(Debug) << "subH: " << std::endl << *subH << std::endl;
1042  // om_->stream(Debug) << "subR: " << std::endl << *subR << std::endl;
1043  // om_->stream(Debug) << "H: " << std::endl << *H_ << std::endl;
1044  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockKrylovSchurOrthoFailure,
1045  "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
1046  }
1047  //
1048  // V has been extended, and H has been extended.
1049  //
1050  // Update basis dim and release all pointers.
1051  Vnext = Teuchos::null;
1052  curDim_ += blockSize_;
1053  // The Ritz vectors/values and Schur form are no longer current.
1054  ritzVecsCurrent_ = false;
1055  ritzValsCurrent_ = false;
1056  schurCurrent_ = false;
1057  //
1058  // Update Ritz values and residuals if needed
1059  if (!(iter_%stepSize_)) {
1061  }
1062 
1063  // When required, monitor some orthogonalities
1064  if (om_->isVerbosity( Debug ) ) {
1065  // Check almost everything here
1066  CheckList chk;
1067  chk.checkV = true;
1068  chk.checkArn = true;
1069  om_->print( Debug, accuracyCheck(chk, ": after local update") );
1070  }
1071  else if (om_->isVerbosity( OrthoDetails ) ) {
1072  CheckList chk;
1073  chk.checkV = true;
1074  om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
1075  }
1076 
1077  // Print information on current iteration
1078  if (om_->isVerbosity(Debug)) {
1079  currentStatus( om_->stream(Debug) );
1080  }
1081  else if (om_->isVerbosity(IterationDetails)) {
1082  currentStatus( om_->stream(IterationDetails) );
1083  }
1084 
1085  } // end while (statusTest == false)
1086 
1087  } // end of iterate()
1088 
1089 
1091  // Check accuracy, orthogonality, and other debugging stuff
1092  //
1093  // bools specify which tests we want to run (instead of running more than we actually care about)
1094  //
1095  // checkV : V orthonormal
1096  // orthogonal to auxvecs
1097  // checkAux: check that auxiliary vectors are actually orthonormal
1098  //
1099  // checkArn: check the Arnoldi factorization
1100  //
1101  // NOTE: This method needs to check the current dimension of the subspace, since it is possible to
1102  // call this method when curDim_ = 0 (after initialization).
1103  template <class ScalarType, class MV, class OP>
1104  std::string BlockKrylovSchur<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1105  {
1106  std::stringstream os;
1107  os.precision(2);
1108  os.setf(std::ios::scientific, std::ios::floatfield);
1109  MagnitudeType tmp;
1110 
1111  os << " Debugging checks: iteration " << iter_ << where << std::endl;
1112 
1113  // index vectors for V and F
1114  std::vector<int> lclind(curDim_);
1115  for (int i=0; i<curDim_; i++) lclind[i] = i;
1116  std::vector<int> bsind(blockSize_);
1117  for (int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
1118 
1119  Teuchos::RCP<const MV> lclV,lclF;
1120  Teuchos::RCP<MV> lclAV;
1121  if (curDim_)
1122  lclV = MVT::CloneView(*V_,lclind);
1123  lclF = MVT::CloneView(*V_,bsind);
1124 
1125  if (chk.checkV) {
1126  if (curDim_) {
1127  tmp = orthman_->orthonormError(*lclV);
1128  os << " >> Error in V^H M V == I : " << tmp << std::endl;
1129  }
1130  tmp = orthman_->orthonormError(*lclF);
1131  os << " >> Error in F^H M F == I : " << tmp << std::endl;
1132  if (curDim_) {
1133  tmp = orthman_->orthogError(*lclV,*lclF);
1134  os << " >> Error in V^H M F == 0 : " << tmp << std::endl;
1135  }
1136  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1137  if (curDim_) {
1138  tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
1139  os << " >> Error in V^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
1140  }
1141  tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
1142  os << " >> Error in F^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
1143  }
1144  }
1145 
1146  if (chk.checkArn) {
1147 
1148  if (curDim_) {
1149  // Compute AV
1150  lclAV = MVT::Clone(*V_,curDim_);
1151  {
1152 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1153  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1154 #endif
1155  OPT::Apply(*Op_,*lclV,*lclAV);
1156  }
1157 
1158  // Compute AV - VH
1159  Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
1160  MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
1161 
1162  // Compute FB_k^T - (AV-VH)
1163  Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
1164  blockSize_,curDim_, curDim_ );
1165  MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
1166 
1167  // Compute || FE_k^T - (AV-VH) ||
1168  std::vector<MagnitudeType> arnNorms( curDim_ );
1169  orthman_->norm( *lclAV, arnNorms );
1170 
1171  for (int i=0; i<curDim_; i++) {
1172  os << " >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
1173  }
1174  }
1175  }
1176 
1177  if (chk.checkAux) {
1178  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1179  tmp = orthman_->orthonormError(*auxVecs_[i]);
1180  os << " >> Error in Aux[" << i << "]^H M Aux[" << i << "] == I : " << tmp << std::endl;
1181  for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1182  tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1183  os << " >> Error in Aux[" << i << "]^H M Aux[" << j << "] == 0 : " << tmp << std::endl;
1184  }
1185  }
1186  }
1187 
1188  os << std::endl;
1189 
1190  return os.str();
1191  }
1192 
1194  /* Get the current approximate eigenvalues, i.e. Ritz values.
1195  *
1196  * POST-CONDITIONS:
1197  *
1198  * ritzValues_ contains Ritz w.r.t. V, H
1199  * Q_ contains the Schur vectors w.r.t. H
1200  * schurH_ contains the Schur matrix w.r.t. H
1201  * ritzOrder_ contains the current ordering from sort manager
1202  */
1203 
1204  template <class ScalarType, class MV, class OP>
1206  {
1207  // Can only call this if the solver is initialized
1208  if (initialized_) {
1209 
1210  // This just updates the Ritz values and residuals.
1211  // --> ritzValsCurrent_ will be set to 'true' by this method.
1212  if (!ritzValsCurrent_) {
1213  // Compute the current Ritz values, through computing the Schur form
1214  // without updating the current projection matrix or sorting the Schur form.
1215  computeSchurForm( false );
1216  }
1217  }
1218  }
1219 
1221  /* Get the current approximate eigenvectors, i.e. Ritz vectors.
1222  *
1223  * POST-CONDITIONS:
1224  *
1225  * ritzValues_ contains Ritz w.r.t. V, H
1226  * ritzVectors_ is first blockSize_ Ritz vectors w.r.t. V, H
1227  * Q_ contains the Schur vectors w.r.t. H
1228  * schurH_ contains the Schur matrix w.r.t. H
1229  * ritzOrder_ contains the current ordering from sort manager
1230  */
1231 
1232  template <class ScalarType, class MV, class OP>
1234  {
1235 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1236  Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
1237 #endif
1238 
1239  TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
1240  "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
1241 
1242  TEUCHOS_TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
1243  "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
1244 
1245 
1246  // Check to see if the current subspace dimension is non-trivial and the solver is initialized
1247  if (curDim_ && initialized_) {
1248 
1249  // Check to see if the Ritz vectors are current.
1250  if (!ritzVecsCurrent_) {
1251 
1252  // Check to see if the Schur factorization of H (schurH_, Q) is current and sorted.
1253  if (!schurCurrent_) {
1254  // Compute the Schur factorization of the current H, which will not directly change H,
1255  // the factorization will be sorted and placed in (schurH_, Q)
1256  computeSchurForm( true );
1257  }
1258 
1259  // After the Schur form is computed, then the Ritz values are current.
1260  // Thus, I can check the Ritz index vector to see if I have enough space for the Ritz vectors requested.
1261  TEUCHOS_TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
1262  "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
1263 
1264  Teuchos::LAPACK<int,ScalarType> lapack;
1265  Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
1266 
1267  // Compute the Ritz vectors.
1268  // --> For a Hermitian problem this is simply the current basis times the first numRitzVecs_ Schur vectors
1269  //
1270  // --> For a non-Hermitian problem, this involves solving the projected eigenproblem, then
1271  // placing the product of the current basis times the first numRitzVecs_ Schur vectors times the
1272  // eigenvectors of interest into the Ritz vectors.
1273 
1274  // Get a view of the current Krylov-Schur basis vectors and Schur vectors
1275  std::vector<int> curind( curDim_ );
1276  for (int i=0; i<curDim_; i++) { curind[i] = i; }
1277  Teuchos::RCP<const MV> Vtemp = MVT::CloneView( *V_, curind );
1278  if (problem_->isHermitian()) {
1279  // Get a view into the current Schur vectors
1280  Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
1281 
1282  // Compute the current Ritz vectors
1283  MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
1284 
1285  // Double check that no complex Ritz values have snuck into the set of converged nev.
1286  bool complexRitzVals = false;
1287  for (int i=0; i<numRitzVecs_; i++) {
1288  if (ritzIndex_[i] != 0) {
1289  complexRitzVals = true;
1290  break;
1291  }
1292  }
1293  if (complexRitzVals)
1294  om_->stream(Warnings) << " Eigenproblem is Hermitian and complex eigenvalues have converged, corresponding eigenvectors will be incorrect!!!"
1295  << std::endl;
1296  } else {
1297 
1298  // Get a view into the current Schur vectors.
1299  Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1300 
1301  // Get a set of work vectors to hold the current Ritz vectors.
1302  Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
1303 
1304  // Compute the current Krylov-Schur vectors.
1305  MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );
1306 
1307  // Now compute the eigenvectors of the Schur form
1308  // Reset the dense matrix and compute the eigenvalues of the Schur form.
1309  //
1310  // Allocate the work space. This space will be used below for calls to:
1311  // * TREVC (requires 3*N for real, 2*N for complex)
1312 
1313  int lwork = 3*curDim_;
1314  std::vector<ScalarType> work( lwork );
1315  std::vector<MagnitudeType> rwork( curDim_ );
1316  char side = 'R';
1317  int mm, info = 0;
1318  const int ldvl = 1;
1319  ScalarType vl[ ldvl ];
1320  Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
1321  lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1322  copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1323  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1324  "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
1325 
1326  // Get a view into the eigenvectors of the Schur form
1327  Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
1328 
1329  // Convert back to Ritz vectors of the operator.
1330  curind.resize( numRitzVecs_ ); // This is already initialized above
1331  Teuchos::RCP<MV> view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, curind );
1332  MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
1333 
1334  // Compute the norm of the new Ritz vectors
1335  std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
1336  MVT::MvNorm( *view_ritzVectors, ritzNrm );
1337 
1338  // Release memory used to compute Ritz vectors before scaling the current vectors.
1339  tmpritzVectors_ = Teuchos::null;
1340  view_ritzVectors = Teuchos::null;
1341 
1342  // Scale the Ritz vectors to have Euclidean norm.
1343  ScalarType ritzScale = ST_ONE;
1344  for (int i=0; i<numRitzVecs_; i++) {
1345 
1346  // If this is a conjugate pair then normalize by the real and imaginary parts.
1347  if (ritzIndex_[i] == 1 ) {
1348  ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
1349  std::vector<int> newind(2);
1350  newind[0] = i; newind[1] = i+1;
1351  tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1352  view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1353  MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1354 
1355  // Increment counter for imaginary part
1356  i++;
1357  } else {
1358 
1359  // This is a real Ritz value, normalize the vector
1360  std::vector<int> newind(1);
1361  newind[0] = i;
1362  tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1363  view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1364  MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1365  }
1366  }
1367 
1368  } // if (problem_->isHermitian())
1369 
1370  // The current Ritz vectors have been computed.
1371  ritzVecsCurrent_ = true;
1372 
1373  } // if (!ritzVecsCurrent_)
1374  } // if (curDim_)
1375  } // computeRitzVectors()
1376 
1377 
1379  // Set a new StatusTest for the solver.
1380  template <class ScalarType, class MV, class OP>
1382  TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1383  "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
1384  tester_ = test;
1385  }
1386 
1388  // Get the current StatusTest used by the solver.
1389  template <class ScalarType, class MV, class OP>
1390  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockKrylovSchur<ScalarType,MV,OP>::getStatusTest() const {
1391  return tester_;
1392  }
1393 
1394 
1396  /* Get the current approximate eigenvalues, i.e. Ritz values.
1397  *
1398  * POST-CONDITIONS:
1399  *
1400  * ritzValues_ contains Ritz w.r.t. V, H
1401  * Q_ contains the Schur vectors w.r.t. H
1402  * schurH_ contains the Schur matrix w.r.t. H
1403  * ritzOrder_ contains the current ordering from sort manager
1404  * schurCurrent_ = true if sort = true; i.e. the Schur form is sorted according to the index
1405  * vector returned by the sort manager.
1406  */
1407  template <class ScalarType, class MV, class OP>
1409  {
1410  // local timer
1411 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1412  Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
1413 #endif
1414 
1415  // Check to see if the dimension of the factorization is greater than zero.
1416  if (curDim_) {
1417 
1418  // Check to see if the Schur factorization is current.
1419  if (!schurCurrent_) {
1420 
1421  // Check to see if the Ritz values are current
1422  // --> If they are then the Schur factorization is current but not sorted.
1423  if (!ritzValsCurrent_) {
1424  Teuchos::LAPACK<int,ScalarType> lapack;
1425  Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
1426  Teuchos::BLAS<int,ScalarType> blas;
1427  Teuchos::BLAS<int,MagnitudeType> blas_mag;
1428 
1429  // Get a view into Q, the storage for H's Schur vectors.
1430  Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1431 
1432  // Get a copy of H to compute/sort the Schur form.
1433  schurH_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
1434  //
1435  //---------------------------------------------------
1436  // Compute the Schur factorization of subH
1437  // ---> Use driver GEES to first reduce to upper Hessenberg
1438  // form and then compute Schur form, outputting Ritz values
1439  //---------------------------------------------------
1440  //
1441  // Allocate the work space. This space will be used below for calls to:
1442  // * GEES (requires 3*N for real, 2*N for complex)
1443  // * TREVC (requires 3*N for real, 2*N for complex)
1444  // * TREXC (requires N for real, none for complex)
1445  // Furthermore, GEES requires a real array of length curDim_ (for complex datatypes)
1446  //
1447  int lwork = 3*curDim_;
1448  std::vector<ScalarType> work( lwork );
1449  std::vector<MagnitudeType> rwork( curDim_ );
1450  std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
1451  std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
1452  std::vector<int> bwork( curDim_ );
1453  int info = 0, sdim = 0;
1454  char jobvs = 'V';
1455  lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
1456  &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork,
1457  &rwork[0], &bwork[0], &info );
1458 
1459  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1460  "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ << ") returned info " << info << " != 0.");
1461 
1462  // Check if imaginary part is detected by the Schur factorization of subH for Hermitian eigenproblems
1463  // NOTE: Because of full orthogonalization, there will be small entries above the block tridiagonal in the block Hessenberg matrix.
1464  // The spectrum of this matrix may include imaginary eigenvalues with small imaginary part, which will mess up the Schur
1465  // form sorting below.
1466  bool hermImagDetected = false;
1467  if (problem_->isHermitian()) {
1468  for (int i=0; i<curDim_; i++)
1469  {
1470  if (tmp_iRitzValues[i] != MT_ZERO)
1471  {
1472  hermImagDetected = true;
1473  break;
1474  }
1475  }
1476  if (hermImagDetected) {
1477  // Warn the user that complex eigenvalues have been detected.
1478  om_->stream(Warnings) << " Eigenproblem is Hermitian and complex eigenvalues have been detected!!!" << std::endl;
1479  // Compute || H - H' || to indicate how bad the symmetry is in the projected eigenproblem
1480  Teuchos::SerialDenseMatrix<int,ScalarType> localH( Teuchos::View, *H_, curDim_, curDim_ );
1481  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > tLocalH;
1482  if (Teuchos::ScalarTraits<ScalarType>::isComplex)
1483  tLocalH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( localH, Teuchos::CONJ_TRANS ) );
1484  else
1485  tLocalH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( localH, Teuchos::TRANS ) );
1486  (*tLocalH) -= localH;
1487  MagnitudeType normF = tLocalH->normFrobenius();
1488  MagnitudeType norm1 = tLocalH->normOne();
1489  om_->stream(Warnings) << " Symmetry error in projected eigenproblem: || S - S' ||_F = " << normF
1490  << ", || S - S' ||_1 = " << norm1 << std::endl;
1491  }
1492  }
1493  //
1494  //---------------------------------------------------
1495  // Use the Krylov-Schur factorization to compute the current Ritz residuals
1496  // for ALL the eigenvalues estimates (Ritz values)
1497  // || Ax - x\theta || = || U_m+1*B_m+1^H*Q*s ||
1498  // = || B_m+1^H*Q*s ||
1499  //
1500  // where U_m+1 is the current Krylov-Schur basis, Q are the Schur vectors, and x = U_m+1*Q*s
1501  // NOTE: This means that s = e_i if the problem is hermitian, else the eigenvectors
1502  // of the Schur form need to be computed.
1503  //
1504  // First compute H_{m+1,m}*B_m^T, then determine what 's' is.
1505  //---------------------------------------------------
1506  //
1507  // Get current B_m+1
1508  Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View, *H_,
1509  blockSize_, curDim_, curDim_ );
1510  //
1511  // Compute B_m+1^H*Q
1512  Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
1513  blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1514  curB.values(), curB.stride(), subQ.values(), subQ.stride(),
1515  ST_ZERO, subB.values(), subB.stride() );
1516  //
1517  // Determine what 's' is and compute Ritz residuals.
1518  //
1519  ScalarType* b_ptr = subB.values();
1520  if (problem_->isHermitian() && !hermImagDetected) {
1521  //
1522  // 's' is the i-th canonical basis vector.
1523  //
1524  for (int i=0; i<curDim_ ; i++) {
1525  ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
1526  }
1527  } else {
1528  //
1529  // Compute S: the eigenvectors of the block upper triangular, Schur matrix.
1530  //
1531  char side = 'R';
1532  int mm;
1533  const int ldvl = 1;
1534  ScalarType vl[ ldvl ];
1535  Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
1536  lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1537  S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1538 
1539  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1540  "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
1541  //
1542  // Scale the eigenvectors so that their Euclidean norms are all one.
1543  //
1544  HelperTraits<ScalarType>::scaleRitzVectors( tmp_iRitzValues, &S );
1545  //
1546  // Compute ritzRes = *B_m+1^H*Q*S where the i-th column of S is 's' for the i-th Ritz-value
1547  //
1548  Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
1549  blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1550  subB.values(), subB.stride(), S.values(), S.stride(),
1551  ST_ZERO, ritzRes.values(), ritzRes.stride() );
1552 
1553  /* TO DO: There's be an incorrect assumption made in the computation of the Ritz residuals.
1554  This assumption is that the next vector in the Krylov subspace is Euclidean orthonormal.
1555  It may not be normalized using Euclidean norm.
1556  Teuchos::RCP<MV> ritzResVecs = MVT::Clone( *V_, curDim_ );
1557  std::vector<int> curind(blockSize_);
1558  for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
1559  Teuchos::RCP<MV> Vtemp = MVT::CloneView(*V_,curind);
1560 
1561  MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, ritzRes, ST_ZERO, *ritzResVecs );
1562  std::vector<MagnitudeType> ritzResNrms(curDim_);
1563  MVT::MvNorm( *ritzResVecs, ritzResNrms );
1564  i = 0;
1565  while( i < curDim_ ) {
1566  if ( tmp_ritzValues[curDim_+i] != MT_ZERO ) {
1567  ritzResiduals_[i] = lapack_mag.LAPY2( ritzResNrms[i], ritzResNrms[i+1] );
1568  ritzResiduals_[i+1] = ritzResiduals_[i];
1569  i = i+2;
1570  } else {
1571  ritzResiduals_[i] = ritzResNrms[i];
1572  i++;
1573  }
1574  }
1575  */
1576  //
1577  // Compute the Ritz residuals for each Ritz value.
1578  //
1579  HelperTraits<ScalarType>::computeRitzResiduals( tmp_iRitzValues, ritzRes, &ritzResiduals_ );
1580  }
1581  //
1582  // Sort the Ritz values.
1583  //
1584  {
1585 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1586  Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
1587 #endif
1588  int i=0;
1589  if (problem_->isHermitian() && !hermImagDetected) {
1590  //
1591  // Sort using just the real part of the Ritz values.
1592  sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_); // don't catch exception
1593  ritzIndex_.clear();
1594  while ( i < curDim_ ) {
1595  // The Ritz value is not complex.
1596  ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
1597  ritzIndex_.push_back(0);
1598  i++;
1599  }
1600  }
1601  else {
1602  //
1603  // Sort using both the real and imaginary parts of the Ritz values.
1604  sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
1605  HelperTraits<ScalarType>::sortRitzValues( tmp_rRitzValues, tmp_iRitzValues, &ritzValues_, &ritzOrder_, &ritzIndex_ );
1606  }
1607  //
1608  // Sort the ritzResiduals_ based on the ordering from the Sort Manager.
1609  std::vector<MagnitudeType> ritz2( curDim_ );
1610  for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
1611  blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
1612 
1613  // The Ritz values have now been updated.
1614  ritzValsCurrent_ = true;
1615  }
1616 
1617  } // if (!ritzValsCurrent_) ...
1618  //
1619  //---------------------------------------------------
1620  // * The Ritz values and residuals have been updated at this point.
1621  //
1622  // * The Schur factorization of the projected matrix has been computed,
1623  // and is stored in (schurH_, Q_).
1624  //
1625  // Now the Schur factorization needs to be sorted.
1626  //---------------------------------------------------
1627  //
1628  // Sort the Schur form using the ordering from the Sort Manager.
1629  if (sort) {
1630  sortSchurForm( *schurH_, *Q_, ritzOrder_ );
1631  //
1632  // Indicate the Schur form in (schurH_, Q_) is current and sorted
1633  schurCurrent_ = true;
1634  }
1635  } // if (!schurCurrent_) ...
1636 
1637  } // if (curDim_) ...
1638 
1639  } // computeSchurForm( ... )
1640 
1641 
1643  // Sort the Schur form of H stored in (H,Q) using the ordering vector.
1644  template <class ScalarType, class MV, class OP>
1645  void BlockKrylovSchur<ScalarType,MV,OP>::sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
1646  Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
1647  std::vector<int>& order )
1648  {
1649  // local timer
1650 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1651  Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
1652 #endif
1653  //
1654  //---------------------------------------------------
1655  // Reorder real Schur factorization, remember to add one to the indices for the
1656  // fortran call and determine offset. The offset is necessary since the TREXC
1657  // method reorders in a nonsymmetric fashion, thus we use the reordering in
1658  // a stack-like fashion. Also take into account conjugate pairs, which may mess
1659  // up the reordering, since the pair is moved if one of the pair is moved.
1660  //---------------------------------------------------
1661  //
1662  int i = 0, nevtemp = 0;
1663  char compq = 'V';
1664  std::vector<int> offset2( curDim_ );
1665  std::vector<int> order2( curDim_ );
1666 
1667  // LAPACK objects.
1668  Teuchos::LAPACK<int,ScalarType> lapack;
1669  int lwork = 3*curDim_;
1670  std::vector<ScalarType> work( lwork );
1671 
1672  while (i < curDim_) {
1673  if ( ritzIndex_[i] != 0 ) { // This is the first value of a complex conjugate pair
1674  offset2[nevtemp] = 0;
1675  for (int j=i; j<curDim_; j++) {
1676  if (order[j] > order[i]) { offset2[nevtemp]++; }
1677  }
1678  order2[nevtemp] = order[i];
1679  i = i+2;
1680  } else {
1681  offset2[nevtemp] = 0;
1682  for (int j=i; j<curDim_; j++) {
1683  if (order[j] > order[i]) { offset2[nevtemp]++; }
1684  }
1685  order2[nevtemp] = order[i];
1686  i++;
1687  }
1688  nevtemp++;
1689  }
1690  ScalarType *ptr_h = H.values();
1691  ScalarType *ptr_q = Q.values();
1692  int ldh = H.stride(), ldq = Q.stride();
1693  int info = 0;
1694  for (i=nevtemp-1; i>=0; i--) {
1695  lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, order2[i]+1+offset2[i],
1696  1, &work[0], &info );
1697  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1698  "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ << ") returned info " << info << " != 0.");
1699  }
1700  }
1701 
1703  // Print the current status of the solver
1704  template <class ScalarType, class MV, class OP>
1706  {
1707  using std::endl;
1708 
1709  os.setf(std::ios::scientific, std::ios::floatfield);
1710  os.precision(6);
1711  os <<"================================================================================" << endl;
1712  os << endl;
1713  os <<" BlockKrylovSchur Solver Status" << endl;
1714  os << endl;
1715  os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1716  os <<"The number of iterations performed is " <<iter_<<endl;
1717  os <<"The block size is " << blockSize_<<endl;
1718  os <<"The number of blocks is " << numBlocks_<<endl;
1719  os <<"The current basis size is " << curDim_<<endl;
1720  os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1721  os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1722 
1723  os.setf(std::ios_base::right, std::ios_base::adjustfield);
1724 
1725  os << endl;
1726  if (initialized_) {
1727  os <<"CURRENT RITZ VALUES "<<endl;
1728  if (ritzIndex_.size() != 0) {
1729  int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
1730  if (problem_->isHermitian()) {
1731  os << std::setw(20) << "Ritz Value"
1732  << std::setw(20) << "Ritz Residual"
1733  << endl;
1734  os <<"--------------------------------------------------------------------------------"<<endl;
1735  for (int i=0; i<numPrint; i++) {
1736  os << std::setw(20) << ritzValues_[i].realpart
1737  << std::setw(20) << ritzResiduals_[i]
1738  << endl;
1739  }
1740  } else {
1741  os << std::setw(24) << "Ritz Value"
1742  << std::setw(30) << "Ritz Residual"
1743  << endl;
1744  os <<"--------------------------------------------------------------------------------"<<endl;
1745  for (int i=0; i<numPrint; i++) {
1746  // Print out the real eigenvalue.
1747  os << std::setw(15) << ritzValues_[i].realpart;
1748  if (ritzValues_[i].imagpart < MT_ZERO) {
1749  os << " - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
1750  } else {
1751  os << " + i" << std::setw(15) << ritzValues_[i].imagpart;
1752  }
1753  os << std::setw(20) << ritzResiduals_[i] << endl;
1754  }
1755  }
1756  } else {
1757  os << std::setw(20) << "[ NONE COMPUTED ]" << endl;
1758  }
1759  }
1760  os << endl;
1761  os <<"================================================================================" << endl;
1762  os << endl;
1763  }
1764 
1765 } // End of namespace Anasazi
1766 
1767 #endif
1768 
1769 // End of file AnasaziBlockKrylovSchur.hpp
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
void computeRitzValues()
Compute the Ritz values using the current Krylov factorization.
BlockKrylovSchurOrthoFailure is thrown when the orthogonalization manager is unable to generate ortho...
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
BlockKrylovSchurState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
void setStepSize(int stepSize)
Set the step size.
bool isRitzValsCurrent() const
Get the status of the Ritz values currently stored in the eigensolver.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the current Ritz residual 2-norms.
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors.
This class defines the interface required by an eigensolver and status test class to compute solution...
virtual ~BlockKrylovSchur()
BlockKrylovSchur destructor.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
Declaration of basic traits for the multivector type.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
void iterate()
This method performs Block Krylov-Schur iterations until the status test indicates the need to stop o...
Virtual base class which defines basic traits for the operator type.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
static void scaleRitzVectors(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, Teuchos::SerialDenseMatrix< int, ScalarType > *S)
Helper function for correctly scaling the eigenvectors of the projected eigenproblem.
Structure to contain pointers to BlockKrylovSchur state variables.
An exception class parent to all Anasazi exceptions.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
int getNumIters() const
Get the current iteration count.
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...
bool isRitzVecsCurrent() const
Get the status of the Ritz vectors currently stored in the eigensolver.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
void computeRitzVectors()
Compute the Ritz vectors using the current Krylov factorization.
int getStepSize() const
Get the step size.
void setNumRitzVectors(int numRitzVecs)
Set the number of Ritz vectors to compute.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
void computeSchurForm(const bool sort=true)
Compute the Schur form of the projected eigenproblem from the current Krylov factorization.
static void sortRitzValues(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &rRV, const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, std::vector< Value< ScalarType > > *RV, std::vector< int > *RO, std::vector< int > *RI)
Helper function for correctly storing the Ritz values when the eigenproblem is non-Hermitian.
bool isSchurCurrent() const
Get the status of the Schur form currently stored in the eigensolver.
BlockKrylovSchur(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< OrthoManager< ScalarType, MV > > &ortho, Teuchos::ParameterList &params)
BlockKrylovSchur constructor with eigenproblem, solver utilities, and parameter list of solver option...
Output managers remove the need for the eigensolver to know any information about the required output...
void resetNumIters()
Reset the iteration count.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Virtual base class which defines basic traits for the operator type.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
int getNumRitzVectors() const
Get the number of Ritz vectors to compute.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values.
std::vector< int > getRitzIndex()
Get the Ritz index vector.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms.
Teuchos::RCP< const MulVec > V
The current Krylov basis.
BlockKrylovSchurInitFailure is thrown when the BlockKrylovSchur solver is unable to generate an initi...
Types and exceptions used within Anasazi solvers and interfaces.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > Q
The current Schur vectors of the valid part of H.
static void computeRitzResiduals(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, const Teuchos::SerialDenseMatrix< int, ScalarType > &S, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *RR)
Helper function for correctly computing the Ritz residuals of the projected eigenproblem.
void setBlockSize(int blockSize)
Set the blocksize.
This class implements the block Krylov-Schur iteration, for solving linear eigenvalue problems...
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
int curDim
The current dimension of the reduction.
Common interface of stopping criteria for Anasazi&#39;s solvers.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.