Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
46 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_HPP
47 #define ANASAZI_BLOCK_KRYLOV_SCHUR_HPP
48 
49 #include "AnasaziTypes.hpp"
50 
51 #include "AnasaziEigensolver.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 #include "AnasaziHelperTraits.hpp"
56 
57 #include "AnasaziOrthoManager.hpp"
58 
59 #include "Teuchos_LAPACK.hpp"
60 #include "Teuchos_BLAS.hpp"
63 #include "Teuchos_TimeMonitor.hpp"
64 
79 namespace Anasazi {
80 
82 
83 
88  template <class ScalarType, class MulVec>
94  int curDim;
107  BlockKrylovSchurState() : curDim(0), V(Teuchos::null),
108  H(Teuchos::null), S(Teuchos::null),
109  Q(Teuchos::null) {}
110  };
111 
113 
115 
116 
130  BlockKrylovSchurInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
131  {}};
132 
140  BlockKrylovSchurOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
141  {}};
142 
144 
145 
146  template <class ScalarType, class MV, class OP>
147  class BlockKrylovSchur : public Eigensolver<ScalarType,MV,OP> {
148  public:
150 
151 
165  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
168  Teuchos::ParameterList &params
169  );
170 
172  virtual ~BlockKrylovSchur() {};
174 
175 
177 
178 
200  void iterate();
201 
224 
228  void initialize();
229 
238  bool isInitialized() const { return initialized_; }
239 
249  state.curDim = curDim_;
250  state.V = V_;
251  state.H = H_;
252  state.Q = Q_;
253  state.S = schurH_;
254  return state;
255  }
256 
258 
259 
261 
262 
264  int getNumIters() const { return(iter_); }
265 
267  void resetNumIters() { iter_=0; }
268 
276  Teuchos::RCP<const MV> getRitzVectors() { return ritzVectors_; }
277 
285  std::vector<Value<ScalarType> > getRitzValues() {
286  std::vector<Value<ScalarType> > ret = ritzValues_;
287  ret.resize(ritzIndex_.size());
288  return ret;
289  }
290 
296  std::vector<int> getRitzIndex() { return ritzIndex_; }
297 
302  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms() {
303  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
304  return ret;
305  }
306 
311  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms() {
312  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
313  return ret;
314  }
315 
320  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms() {
321  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
322  ret.resize(ritzIndex_.size());
323  return ret;
324  }
325 
327 
329 
330 
333 
336 
338  const Eigenproblem<ScalarType,MV,OP>& getProblem() const { return(*problem_); };
339 
346  void setSize(int blockSize, int numBlocks);
347 
349  void setBlockSize(int blockSize);
350 
352  void setStepSize(int stepSize);
353 
355  void setNumRitzVectors(int numRitzVecs);
356 
358  int getStepSize() const { return(stepSize_); }
359 
361  int getBlockSize() const { return(blockSize_); }
362 
364  int getNumRitzVectors() const { return(numRitzVecs_); }
365 
371  int getCurSubspaceDim() const {
372  if (!initialized_) return 0;
373  return curDim_;
374  }
375 
377  int getMaxSubspaceDim() const { return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
378 
379 
392  void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
393 
396 
398 
400 
401 
403  void currentStatus(std::ostream &os);
404 
406 
408 
409 
411  bool isRitzVecsCurrent() const { return ritzVecsCurrent_; }
412 
414  bool isRitzValsCurrent() const { return ritzValsCurrent_; }
415 
417  bool isSchurCurrent() const { return schurCurrent_; }
418 
420 
422 
423 
425  void computeRitzVectors();
426 
428  void computeRitzValues();
429 
431  void computeSchurForm( const bool sort = true );
432 
434 
435  private:
436  //
437  // Convenience typedefs
438  //
439  typedef MultiVecTraits<ScalarType,MV> MVT;
442  typedef typename SCT::magnitudeType MagnitudeType;
443  typedef typename std::vector<ScalarType>::iterator STiter;
444  typedef typename std::vector<MagnitudeType>::iterator MTiter;
445  const MagnitudeType MT_ONE;
446  const MagnitudeType MT_ZERO;
447  const MagnitudeType NANVAL;
448  const ScalarType ST_ONE;
449  const ScalarType ST_ZERO;
450  //
451  // Internal structs
452  //
453  struct CheckList {
454  bool checkV;
455  bool checkArn;
456  bool checkAux;
457  CheckList() : checkV(false), checkArn(false), checkAux(false) {};
458  };
459  //
460  // Internal methods
461  //
462  std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
463  void sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
465  std::vector<int>& order );
466  //
467  // Classes inputed through constructor that define the eigenproblem to be solved.
468  //
474  //
475  // Information obtained from the eigenproblem
476  //
478  //
479  // Internal timers
480  //
481  Teuchos::RCP<Teuchos::Time> timerOp_, timerSortRitzVal_,
482  timerCompSF_, timerSortSF_,
483  timerCompRitzVec_, timerOrtho_;
484  //
485  // Counters
486  //
487  int count_ApplyOp_;
488 
489  //
490  // Algorithmic parameters.
491  //
492  // blockSize_ is the solver block size; it controls the number of eigenvectors that
493  // we compute, the number of residual vectors that we compute, and therefore the number
494  // of vectors added to the basis on each iteration.
495  int blockSize_;
496  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
497  int numBlocks_;
498  // stepSize_ dictates how many iterations are performed before eigenvectors and eigenvalues
499  // are computed again
500  int stepSize_;
501 
502  //
503  // Current solver state
504  //
505  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
506  // is capable of running; _initialize is controlled by the initialize() member method
507  // For the implications of the state of initialized_, please see documentation for initialize()
508  bool initialized_;
509  //
510  // curDim_ reflects how much of the current basis is valid
511  // NOTE: for Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_
512  // for non-Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_ + 1
513  // this also tells us how many of the values in _theta are valid Ritz values
514  int curDim_;
515  //
516  // State Multivecs
517  Teuchos::RCP<MV> ritzVectors_, V_;
518  int numRitzVecs_;
519  //
520  // Projected matrices
521  // H_ : Projected matrix from the Krylov-Schur factorization AV = VH + FB^T
522  //
524  //
525  // Schur form of Projected matrices (these are only updated when the Ritz values/vectors are updated).
526  // schurH_: Schur form reduction of H
527  // Q_: Schur vectors of H
530  //
531  // Auxiliary vectors
533  int numAuxVecs_;
534  //
535  // Number of iterations that have been performed.
536  int iter_;
537  //
538  // State flags
539  bool ritzVecsCurrent_, ritzValsCurrent_, schurCurrent_;
540  //
541  // Current eigenvalues, residual norms
542  std::vector<Value<ScalarType> > ritzValues_;
543  std::vector<MagnitudeType> ritzResiduals_;
544  //
545  // Current index vector for Ritz values and vectors
546  std::vector<int> ritzIndex_; // computed by BKS
547  std::vector<int> ritzOrder_; // returned from sort manager
548  //
549  // Number of Ritz pairs to be printed upon output, if possible
550  int numRitzPrint_;
551  };
552 
553 
555  // Constructor
556  template <class ScalarType, class MV, class OP>
560  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
563  Teuchos::ParameterList &params
564  ) :
565  MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
566  MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
567  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
568  ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
569  ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
570  // problem, tools
571  problem_(problem),
572  sm_(sorter),
573  om_(printer),
574  tester_(tester),
575  orthman_(ortho),
576  // timers, counters
577 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
578  timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Operation Op*x")),
579  timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Ritz values")),
580  timerCompSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Schur form")),
581  timerSortSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Schur form")),
582  timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Ritz vectors")),
583  timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Orthogonalization")),
584 #endif
585  count_ApplyOp_(0),
586  // internal data
587  blockSize_(0),
588  numBlocks_(0),
589  stepSize_(0),
590  initialized_(false),
591  curDim_(0),
592  numRitzVecs_(0),
593  auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
594  numAuxVecs_(0),
595  iter_(0),
596  ritzVecsCurrent_(false),
597  ritzValsCurrent_(false),
598  schurCurrent_(false),
599  numRitzPrint_(0)
600  {
601  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
602  "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
603  TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
604  "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
605  TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
606  "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
607  TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
608  "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
609  TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
610  "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
611  TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
612  "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
613  TEUCHOS_TEST_FOR_EXCEPTION(sorter == Teuchos::null,std::invalid_argument,
614  "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
615  TEUCHOS_TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
616  "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
617  TEUCHOS_TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
618  "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
619  TEUCHOS_TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
620  "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
621 
622  // Get problem operator
623  Op_ = problem_->getOperator();
624 
625  // get the step size
626  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Step Size"), std::invalid_argument,
627  "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
628  int ss = params.get("Step Size",numBlocks_);
629  setStepSize(ss);
630 
631  // set the block size and allocate data
632  int bs = params.get("Block Size", 1);
633  int nb = params.get("Num Blocks", 3*problem_->getNEV());
634  setSize(bs,nb);
635 
636  // get the number of Ritz vectors to compute and allocate data.
637  // --> if this parameter is not specified in the parameter list, then it's assumed that no Ritz vectors will be computed.
638  int numRitzVecs = params.get("Number of Ritz Vectors", 0);
639  setNumRitzVectors( numRitzVecs );
640 
641  // get the number of Ritz values to print out when currentStatus is called.
642  numRitzPrint_ = params.get("Print Number of Ritz Values", bs);
643  }
644 
645 
647  // Set the block size
648  // This simply calls setSize(), modifying the block size while retaining the number of blocks.
649  template <class ScalarType, class MV, class OP>
651  {
652  setSize(blockSize,numBlocks_);
653  }
654 
655 
657  // Set the step size.
658  template <class ScalarType, class MV, class OP>
660  {
661  TEUCHOS_TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
662  stepSize_ = stepSize;
663  }
664 
666  // Set the number of Ritz vectors to compute.
667  template <class ScalarType, class MV, class OP>
669  {
670  // This routine only allocates space; it doesn't not perform any computation
671  // any change in size will invalidate the state of the solver.
672 
673  TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
674 
675  // Check to see if the number of requested Ritz vectors has changed.
676  if (numRitzVecs != numRitzVecs_) {
677  if (numRitzVecs) {
678  ritzVectors_ = Teuchos::null;
679  ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
680  } else {
681  ritzVectors_ = Teuchos::null;
682  }
683  numRitzVecs_ = numRitzVecs;
684  ritzVecsCurrent_ = false;
685  }
686  }
687 
689  // Set the block size and make necessary adjustments.
690  template <class ScalarType, class MV, class OP>
691  void BlockKrylovSchur<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
692  {
693  // This routine only allocates space; it doesn't not perform any computation
694  // any change in size will invalidate the state of the solver.
695 
696  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
697  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
698  if (blockSize == blockSize_ && numBlocks == numBlocks_) {
699  // do nothing
700  return;
701  }
702 
703  blockSize_ = blockSize;
704  numBlocks_ = numBlocks;
705 
707  // grab some Multivector to Clone
708  // in practice, getInitVec() should always provide this, but it is possible to use a
709  // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
710  // in case of that strange scenario, we will try to Clone from V_; first resort to getInitVec(),
711  // because we would like to clear the storage associated with V_ so we have room for the new V_
712  if (problem_->getInitVec() != Teuchos::null) {
713  tmp = problem_->getInitVec();
714  }
715  else {
716  tmp = V_;
717  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
718  "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
719  }
720 
721 
723  // blockSize*numBlocks dependent
724  //
725  int newsd;
726  if (problem_->isHermitian()) {
727  newsd = blockSize_*numBlocks_;
728  } else {
729  newsd = blockSize_*numBlocks_+1;
730  }
731  // check that new size is valid
732  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(newsd) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
733  "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
734 
735  ritzValues_.resize(newsd);
736  ritzResiduals_.resize(newsd,MT_ONE);
737  ritzOrder_.resize(newsd);
738  V_ = Teuchos::null;
739  V_ = MVT::Clone(*tmp,newsd+blockSize_);
740  H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd+blockSize_,newsd) );
742 
743  initialized_ = false;
744  curDim_ = 0;
745  }
746 
747 
749  // Set the auxiliary vectors
750  template <class ScalarType, class MV, class OP>
752  typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
753 
754  // set new auxiliary vectors
755  auxVecs_ = auxvecs;
756 
757  if (om_->isVerbosity( Debug ) ) {
758  // Check almost everything here
759  CheckList chk;
760  chk.checkAux = true;
761  om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
762  }
763 
764  numAuxVecs_ = 0;
765  for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
766  numAuxVecs_ += MVT::GetNumberVecs(**i);
767  }
768 
769  // If the solver has been initialized, X and P are not necessarily orthogonal to new auxiliary vectors
770  if (numAuxVecs_ > 0 && initialized_) {
771  initialized_ = false;
772  }
773  }
774 
776  /* Initialize the state of the solver
777  *
778  * POST-CONDITIONS:
779  *
780  * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
781  *
782  */
783 
784  template <class ScalarType, class MV, class OP>
786  {
787  // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
788 
789  std::vector<int> bsind(blockSize_);
790  for (int i=0; i<blockSize_; i++) bsind[i] = i;
791 
792  // in BlockKrylovSchur, V and H are required.
793  // if either doesn't exist, then we will start with the initial vector.
794  //
795  // inconsistent multivectors widths and lengths will not be tolerated, and
796  // will be treated with exceptions.
797  //
798  std::string errstr("Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
799 
800  // set up V,H: if the user doesn't specify both of these these,
801  // we will start over with the initial vector.
802  if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
803 
804  // initialize V_,H_, and curDim_
805 
806  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
807  std::invalid_argument, errstr );
808  if (newstate.V != V_) {
809  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
810  std::invalid_argument, errstr );
811  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) > getMaxSubspaceDim(),
812  std::invalid_argument, errstr );
813  }
814  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > getMaxSubspaceDim(),
815  std::invalid_argument, errstr );
816 
817  curDim_ = newstate.curDim;
818  int lclDim = MVT::GetNumberVecs(*newstate.V);
819 
820  // check size of H
821  TEUCHOS_TEST_FOR_EXCEPTION(newstate.H->numRows() < curDim_ || newstate.H->numCols() < curDim_, std::invalid_argument, errstr);
822 
823  if (curDim_ == 0 && lclDim > blockSize_) {
824  om_->stream(Warnings) << "Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
825  << "The block size however is only " << blockSize_ << std::endl
826  << "The last " << lclDim - blockSize_ << " vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
827  }
828 
829 
830  // copy basis vectors from newstate into V
831  if (newstate.V != V_) {
832  std::vector<int> nevind(lclDim);
833  for (int i=0; i<lclDim; i++) nevind[i] = i;
834  MVT::SetBlock(*newstate.V,nevind,*V_);
835  }
836 
837  // put data into H_, make sure old information is not still hanging around.
838  if (newstate.H != H_) {
839  H_->putScalar( ST_ZERO );
840  Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H,curDim_+blockSize_,curDim_);
842  lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_,curDim_+blockSize_,curDim_) );
843  lclH->assign(newH);
844 
845  // done with local pointers
846  lclH = Teuchos::null;
847  }
848 
849  }
850  else {
851  // user did not specify a basis V
852  // get vectors from problem or generate something, projectAndNormalize, call initialize() recursively
853  Teuchos::RCP<const MV> ivec = problem_->getInitVec();
854  TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
855  "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
856 
857  int lclDim = MVT::GetNumberVecs(*ivec);
858  bool userand = false;
859  if (lclDim < blockSize_) {
860  // we need at least blockSize_ vectors
861  // use a random multivec
862  userand = true;
863  }
864 
865  if (userand) {
866  // make an index
867  std::vector<int> dimind2(lclDim);
868  for (int i=0; i<lclDim; i++) { dimind2[i] = i; }
869 
870  // alloc newV as a view of the first block of V
871  Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,dimind2);
872 
873  // copy the initial vectors into the first lclDim vectors of V
874  MVT::SetBlock(*ivec,dimind2,*newV1);
875 
876  // resize / reinitialize the index vector
877  dimind2.resize(blockSize_-lclDim);
878  for (int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
879 
880  // initialize the rest of the vectors with random vectors
881  Teuchos::RCP<MV> newV2 = MVT::CloneViewNonConst(*V_,dimind2);
882  MVT::MvRandom(*newV2);
883  }
884  else {
885  // alloc newV as a view of the first block of V
886  Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,bsind);
887 
888  // get a view of the first block of initial vectors
889  Teuchos::RCP<const MV> ivecV = MVT::CloneView(*ivec,bsind);
890 
891  // assign ivec to first part of newV
892  MVT::SetBlock(*ivecV,bsind,*newV1);
893  }
894 
895  // get pointer into first block of V
896  Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*V_,bsind);
897 
898  // remove auxVecs from newV and normalize newV
899  if (auxVecs_.size() > 0) {
900 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
901  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
902 #endif
903 
905  int rank = orthman_->projectAndNormalize(*newV,auxVecs_);
907  "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
908  }
909  else {
910 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
911  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
912 #endif
913 
914  int rank = orthman_->normalize(*newV);
916  "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
917  }
918 
919  // set curDim
920  curDim_ = 0;
921 
922  // clear pointer
923  newV = Teuchos::null;
924  }
925 
926  // The Ritz vectors/values and Schur form are no longer current.
927  ritzVecsCurrent_ = false;
928  ritzValsCurrent_ = false;
929  schurCurrent_ = false;
930 
931  // the solver is initialized
932  initialized_ = true;
933 
934  if (om_->isVerbosity( Debug ) ) {
935  // Check almost everything here
936  CheckList chk;
937  chk.checkV = true;
938  chk.checkArn = true;
939  chk.checkAux = true;
940  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
941  }
942 
943  // Print information on current status
944  if (om_->isVerbosity(Debug)) {
945  currentStatus( om_->stream(Debug) );
946  }
947  else if (om_->isVerbosity(IterationDetails)) {
948  currentStatus( om_->stream(IterationDetails) );
949  }
950  }
951 
952 
954  // initialize the solver with default state
955  template <class ScalarType, class MV, class OP>
957  {
959  initialize(empty);
960  }
961 
962 
964  // Perform BlockKrylovSchur iterations until the StatusTest tells us to stop.
965  template <class ScalarType, class MV, class OP>
967  {
968  //
969  // Allocate/initialize data structures
970  //
971  if (initialized_ == false) {
972  initialize();
973  }
974 
975  // Compute the current search dimension.
976  // If the problem is non-Hermitian and the blocksize is one, let the solver use the extra vector.
977  int searchDim = blockSize_*numBlocks_;
978  if (problem_->isHermitian() == false) {
979  searchDim++;
980  }
981 
983  // iterate until the status test tells us to stop.
984  //
985  // also break if our basis is full
986  //
987  while (tester_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
988 
989  iter_++;
990 
991  // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
992  int lclDim = curDim_ + blockSize_;
993 
994  // Get the current part of the basis.
995  std::vector<int> curind(blockSize_);
996  for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
997  Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
998 
999  // Get a view of the previous vectors
1000  // this is used for orthogonalization and for computing V^H K H
1001  for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
1002  Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
1003  // om_->stream(Debug) << "Vprev: " << std::endl;
1004  // MVT::MvPrint(*Vprev,om_->stream(Debug));
1005 
1006  // Compute the next vector in the Krylov basis: Vnext = Op*Vprev
1007  {
1008 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1009  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1010 #endif
1011  OPT::Apply(*Op_,*Vprev,*Vnext);
1012  count_ApplyOp_ += blockSize_;
1013  }
1014  // om_->stream(Debug) << "Vnext: " << std::endl;
1015  // MVT::MvPrint(*Vnext,om_->stream(Debug));
1016  Vprev = Teuchos::null;
1017 
1018  // Remove all previous Krylov-Schur basis vectors and auxVecs from Vnext
1019  {
1020 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1021  Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1022 #endif
1023 
1024  // Get a view of all the previous vectors
1025  std::vector<int> prevind(lclDim);
1026  for (int i=0; i<lclDim; i++) { prevind[i] = i; }
1027  Vprev = MVT::CloneView(*V_,prevind);
1028  Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
1029 
1030  // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
1033  ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
1035  AsubH.append( subH );
1036 
1037  // Add the auxiliary vectors to the current basis vectors if any exist
1038  if (auxVecs_.size() > 0) {
1039  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1040  AVprev.append( auxVecs_[i] );
1041  AsubH.append( Teuchos::null );
1042  }
1043  }
1044 
1045  // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
1046  // om_->stream(Debug) << "V before ortho: " << std::endl;
1047  // MVT::MvPrint(*Vprev,om_->stream(Debug));
1050  ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
1051  int rank = orthman_->projectAndNormalize(*Vnext,AVprev,AsubH,subR);
1052  // om_->stream(Debug) << "Vnext after ortho: " << std::endl;
1053  // MVT::MvPrint(*Vnext,om_->stream(Debug));
1054  // om_->stream(Debug) << "subH: " << std::endl << *subH << std::endl;
1055  // om_->stream(Debug) << "subR: " << std::endl << *subR << std::endl;
1056  // om_->stream(Debug) << "H: " << std::endl << *H_ << std::endl;
1058  "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
1059  }
1060  //
1061  // V has been extended, and H has been extended.
1062  //
1063  // Update basis dim and release all pointers.
1064  Vnext = Teuchos::null;
1065  curDim_ += blockSize_;
1066  // The Ritz vectors/values and Schur form are no longer current.
1067  ritzVecsCurrent_ = false;
1068  ritzValsCurrent_ = false;
1069  schurCurrent_ = false;
1070  //
1071  // Update Ritz values and residuals if needed
1072  if (!(iter_%stepSize_)) {
1073  computeRitzValues();
1074  }
1075 
1076  // When required, monitor some orthogonalities
1077  if (om_->isVerbosity( Debug ) ) {
1078  // Check almost everything here
1079  CheckList chk;
1080  chk.checkV = true;
1081  chk.checkArn = true;
1082  om_->print( Debug, accuracyCheck(chk, ": after local update") );
1083  }
1084  else if (om_->isVerbosity( OrthoDetails ) ) {
1085  CheckList chk;
1086  chk.checkV = true;
1087  om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
1088  }
1089 
1090  // Print information on current iteration
1091  if (om_->isVerbosity(Debug)) {
1092  currentStatus( om_->stream(Debug) );
1093  }
1094  else if (om_->isVerbosity(IterationDetails)) {
1095  currentStatus( om_->stream(IterationDetails) );
1096  }
1097 
1098  } // end while (statusTest == false)
1099 
1100  } // end of iterate()
1101 
1102 
1104  // Check accuracy, orthogonality, and other debugging stuff
1105  //
1106  // bools specify which tests we want to run (instead of running more than we actually care about)
1107  //
1108  // checkV : V orthonormal
1109  // orthogonal to auxvecs
1110  // checkAux: check that auxiliary vectors are actually orthonormal
1111  //
1112  // checkArn: check the Arnoldi factorization
1113  //
1114  // NOTE: This method needs to check the current dimension of the subspace, since it is possible to
1115  // call this method when curDim_ = 0 (after initialization).
1116  template <class ScalarType, class MV, class OP>
1117  std::string BlockKrylovSchur<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1118  {
1119  std::stringstream os;
1120  os.precision(2);
1121  os.setf(std::ios::scientific, std::ios::floatfield);
1122  MagnitudeType tmp;
1123 
1124  os << " Debugging checks: iteration " << iter_ << where << std::endl;
1125 
1126  // index vectors for V and F
1127  std::vector<int> lclind(curDim_);
1128  for (int i=0; i<curDim_; i++) lclind[i] = i;
1129  std::vector<int> bsind(blockSize_);
1130  for (int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
1131 
1132  Teuchos::RCP<const MV> lclV,lclF;
1133  Teuchos::RCP<MV> lclAV;
1134  if (curDim_)
1135  lclV = MVT::CloneView(*V_,lclind);
1136  lclF = MVT::CloneView(*V_,bsind);
1137 
1138  if (chk.checkV) {
1139  if (curDim_) {
1140  tmp = orthman_->orthonormError(*lclV);
1141  os << " >> Error in V^H M V == I : " << tmp << std::endl;
1142  }
1143  tmp = orthman_->orthonormError(*lclF);
1144  os << " >> Error in F^H M F == I : " << tmp << std::endl;
1145  if (curDim_) {
1146  tmp = orthman_->orthogError(*lclV,*lclF);
1147  os << " >> Error in V^H M F == 0 : " << tmp << std::endl;
1148  }
1149  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1150  if (curDim_) {
1151  tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
1152  os << " >> Error in V^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
1153  }
1154  tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
1155  os << " >> Error in F^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
1156  }
1157  }
1158 
1159  if (chk.checkArn) {
1160 
1161  if (curDim_) {
1162  // Compute AV
1163  lclAV = MVT::Clone(*V_,curDim_);
1164  {
1165 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1166  Teuchos::TimeMonitor lcltimer( *timerOp_ );
1167 #endif
1168  OPT::Apply(*Op_,*lclV,*lclAV);
1169  }
1170 
1171  // Compute AV - VH
1173  MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
1174 
1175  // Compute FB_k^T - (AV-VH)
1177  blockSize_,curDim_, curDim_ );
1178  MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
1179 
1180  // Compute || FE_k^T - (AV-VH) ||
1181  std::vector<MagnitudeType> arnNorms( curDim_ );
1182  orthman_->norm( *lclAV, arnNorms );
1183 
1184  for (int i=0; i<curDim_; i++) {
1185  os << " >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
1186  }
1187  }
1188  }
1189 
1190  if (chk.checkAux) {
1191  for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1192  tmp = orthman_->orthonormError(*auxVecs_[i]);
1193  os << " >> Error in Aux[" << i << "]^H M Aux[" << i << "] == I : " << tmp << std::endl;
1194  for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1195  tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1196  os << " >> Error in Aux[" << i << "]^H M Aux[" << j << "] == 0 : " << tmp << std::endl;
1197  }
1198  }
1199  }
1200 
1201  os << std::endl;
1202 
1203  return os.str();
1204  }
1205 
1207  /* Get the current approximate eigenvalues, i.e. Ritz values.
1208  *
1209  * POST-CONDITIONS:
1210  *
1211  * ritzValues_ contains Ritz w.r.t. V, H
1212  * Q_ contains the Schur vectors w.r.t. H
1213  * schurH_ contains the Schur matrix w.r.t. H
1214  * ritzOrder_ contains the current ordering from sort manager
1215  */
1216 
1217  template <class ScalarType, class MV, class OP>
1219  {
1220  // Can only call this if the solver is initialized
1221  if (initialized_) {
1222 
1223  // This just updates the Ritz values and residuals.
1224  // --> ritzValsCurrent_ will be set to 'true' by this method.
1225  if (!ritzValsCurrent_) {
1226  // Compute the current Ritz values, through computing the Schur form
1227  // without updating the current projection matrix or sorting the Schur form.
1228  computeSchurForm( false );
1229  }
1230  }
1231  }
1232 
1234  /* Get the current approximate eigenvectors, i.e. Ritz vectors.
1235  *
1236  * POST-CONDITIONS:
1237  *
1238  * ritzValues_ contains Ritz w.r.t. V, H
1239  * ritzVectors_ is first blockSize_ Ritz vectors w.r.t. V, H
1240  * Q_ contains the Schur vectors w.r.t. H
1241  * schurH_ contains the Schur matrix w.r.t. H
1242  * ritzOrder_ contains the current ordering from sort manager
1243  */
1244 
1245  template <class ScalarType, class MV, class OP>
1247  {
1248 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1249  Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
1250 #endif
1251 
1252  TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
1253  "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
1254 
1255  TEUCHOS_TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
1256  "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
1257 
1258 
1259  // Check to see if the current subspace dimension is non-trivial and the solver is initialized
1260  if (curDim_ && initialized_) {
1261 
1262  // Check to see if the Ritz vectors are current.
1263  if (!ritzVecsCurrent_) {
1264 
1265  // Check to see if the Schur factorization of H (schurH_, Q) is current and sorted.
1266  if (!schurCurrent_) {
1267  // Compute the Schur factorization of the current H, which will not directly change H,
1268  // the factorization will be sorted and placed in (schurH_, Q)
1269  computeSchurForm( true );
1270  }
1271 
1272  // After the Schur form is computed, then the Ritz values are current.
1273  // Thus, I can check the Ritz index vector to see if I have enough space for the Ritz vectors requested.
1274  TEUCHOS_TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
1275  "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
1276 
1279 
1280  // Compute the Ritz vectors.
1281  // --> For a Hermitian problem this is simply the current basis times the first numRitzVecs_ Schur vectors
1282  //
1283  // --> For a non-Hermitian problem, this involves solving the projected eigenproblem, then
1284  // placing the product of the current basis times the first numRitzVecs_ Schur vectors times the
1285  // eigenvectors of interest into the Ritz vectors.
1286 
1287  // Get a view of the current Krylov-Schur basis vectors and Schur vectors
1288  std::vector<int> curind( curDim_ );
1289  for (int i=0; i<curDim_; i++) { curind[i] = i; }
1290  Teuchos::RCP<const MV> Vtemp = MVT::CloneView( *V_, curind );
1291  if (problem_->isHermitian()) {
1292  // Get a view into the current Schur vectors
1293  Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
1294 
1295  // Compute the current Ritz vectors
1296  MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
1297 
1298  // Double check that no complex Ritz values have snuck into the set of converged nev.
1299  bool complexRitzVals = false;
1300  for (int i=0; i<numRitzVecs_; i++) {
1301  if (ritzIndex_[i] != 0) {
1302  complexRitzVals = true;
1303  break;
1304  }
1305  }
1306  if (complexRitzVals)
1307  om_->stream(Warnings) << " Eigenproblem is Hermitian and complex eigenvalues have converged, corresponding eigenvectors will be incorrect!!!"
1308  << std::endl;
1309  } else {
1310 
1311  // Get a view into the current Schur vectors.
1312  Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1313 
1314  // Get a set of work vectors to hold the current Ritz vectors.
1315  Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
1316 
1317  // Compute the current Krylov-Schur vectors.
1318  MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );
1319 
1320  // Now compute the eigenvectors of the Schur form
1321  // Reset the dense matrix and compute the eigenvalues of the Schur form.
1322  //
1323  // Allocate the work space. This space will be used below for calls to:
1324  // * TREVC (requires 3*N for real, 2*N for complex)
1325 
1326  int lwork = 3*curDim_;
1327  std::vector<ScalarType> work( lwork );
1328  std::vector<MagnitudeType> rwork( curDim_ );
1329  char side = 'R';
1330  int mm, info = 0;
1331  const int ldvl = 1;
1332  ScalarType vl[ ldvl ];
1333  Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
1334  lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1335  copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1336  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1337  "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
1338 
1339  // Get a view into the eigenvectors of the Schur form
1340  Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
1341 
1342  // Convert back to Ritz vectors of the operator.
1343  curind.resize( numRitzVecs_ ); // This is already initialized above
1344  Teuchos::RCP<MV> view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, curind );
1345  MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
1346 
1347  // Compute the norm of the new Ritz vectors
1348  std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
1349  MVT::MvNorm( *view_ritzVectors, ritzNrm );
1350 
1351  // Release memory used to compute Ritz vectors before scaling the current vectors.
1352  tmpritzVectors_ = Teuchos::null;
1353  view_ritzVectors = Teuchos::null;
1354 
1355  // Scale the Ritz vectors to have Euclidean norm.
1356  ScalarType ritzScale = ST_ONE;
1357  for (int i=0; i<numRitzVecs_; i++) {
1358 
1359  // If this is a conjugate pair then normalize by the real and imaginary parts.
1360  if (ritzIndex_[i] == 1 ) {
1361  ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
1362  std::vector<int> newind(2);
1363  newind[0] = i; newind[1] = i+1;
1364  tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1365  view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1366  MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1367 
1368  // Increment counter for imaginary part
1369  i++;
1370  } else {
1371 
1372  // This is a real Ritz value, normalize the vector
1373  std::vector<int> newind(1);
1374  newind[0] = i;
1375  tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1376  view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1377  MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1378  }
1379  }
1380 
1381  } // if (problem_->isHermitian())
1382 
1383  // The current Ritz vectors have been computed.
1384  ritzVecsCurrent_ = true;
1385 
1386  } // if (!ritzVecsCurrent_)
1387  } // if (curDim_)
1388  } // computeRitzVectors()
1389 
1390 
1392  // Set a new StatusTest for the solver.
1393  template <class ScalarType, class MV, class OP>
1395  TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1396  "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
1397  tester_ = test;
1398  }
1399 
1401  // Get the current StatusTest used by the solver.
1402  template <class ScalarType, class MV, class OP>
1404  return tester_;
1405  }
1406 
1407 
1409  /* Get the current approximate eigenvalues, i.e. Ritz values.
1410  *
1411  * POST-CONDITIONS:
1412  *
1413  * ritzValues_ contains Ritz w.r.t. V, H
1414  * Q_ contains the Schur vectors w.r.t. H
1415  * schurH_ contains the Schur matrix w.r.t. H
1416  * ritzOrder_ contains the current ordering from sort manager
1417  * schurCurrent_ = true if sort = true; i.e. the Schur form is sorted according to the index
1418  * vector returned by the sort manager.
1419  */
1420  template <class ScalarType, class MV, class OP>
1422  {
1423  // local timer
1424 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1425  Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
1426 #endif
1427 
1428  // Check to see if the dimension of the factorization is greater than zero.
1429  if (curDim_) {
1430 
1431  // Check to see if the Schur factorization is current.
1432  if (!schurCurrent_) {
1433 
1434  // Check to see if the Ritz values are current
1435  // --> If they are then the Schur factorization is current but not sorted.
1436  if (!ritzValsCurrent_) {
1441 
1442  // Get a view into Q, the storage for H's Schur vectors.
1443  Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1444 
1445  // Get a copy of H to compute/sort the Schur form.
1446  schurH_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
1447  //
1448  //---------------------------------------------------
1449  // Compute the Schur factorization of subH
1450  // ---> Use driver GEES to first reduce to upper Hessenberg
1451  // form and then compute Schur form, outputting Ritz values
1452  //---------------------------------------------------
1453  //
1454  // Allocate the work space. This space will be used below for calls to:
1455  // * GEES (requires 3*N for real, 2*N for complex)
1456  // * TREVC (requires 3*N for real, 2*N for complex)
1457  // * TREXC (requires N for real, none for complex)
1458  // Furthermore, GEES requires a real array of length curDim_ (for complex datatypes)
1459  //
1460  int lwork = 3*curDim_;
1461  std::vector<ScalarType> work( lwork );
1462  std::vector<MagnitudeType> rwork( curDim_ );
1463  std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
1464  std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
1465  std::vector<int> bwork( curDim_ );
1466  int info = 0, sdim = 0;
1467  char jobvs = 'V';
1468  lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
1469  &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork,
1470  &rwork[0], &bwork[0], &info );
1471 
1472  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1473  "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ << ") returned info " << info << " != 0.");
1474 
1475  // Check if imaginary part is detected by the Schur factorization of subH for Hermitian eigenproblems
1476  // NOTE: Because of full orthogonalization, there will be small entries above the block tridiagonal in the block Hessenberg matrix.
1477  // The spectrum of this matrix may include imaginary eigenvalues with small imaginary part, which will mess up the Schur
1478  // form sorting below.
1479  bool hermImagDetected = false;
1480  if (problem_->isHermitian()) {
1481  for (int i=0; i<curDim_; i++)
1482  {
1483  if (tmp_iRitzValues[i] != MT_ZERO)
1484  {
1485  hermImagDetected = true;
1486  break;
1487  }
1488  }
1489  if (hermImagDetected) {
1490  // Warn the user that complex eigenvalues have been detected.
1491  om_->stream(Warnings) << " Eigenproblem is Hermitian and complex eigenvalues have been detected!!!" << std::endl;
1492  // Compute || H - H' || to indicate how bad the symmetry is in the projected eigenproblem
1493  Teuchos::SerialDenseMatrix<int,ScalarType> localH( Teuchos::View, *H_, curDim_, curDim_ );
1497  else
1499  (*tLocalH) -= localH;
1500  MagnitudeType normF = tLocalH->normFrobenius();
1501  MagnitudeType norm1 = tLocalH->normOne();
1502  om_->stream(Warnings) << " Symmetry error in projected eigenproblem: || S - S' ||_F = " << normF
1503  << ", || S - S' ||_1 = " << norm1 << std::endl;
1504  }
1505  }
1506  //
1507  //---------------------------------------------------
1508  // Use the Krylov-Schur factorization to compute the current Ritz residuals
1509  // for ALL the eigenvalues estimates (Ritz values)
1510  // || Ax - x\theta || = || U_m+1*B_m+1^H*Q*s ||
1511  // = || B_m+1^H*Q*s ||
1512  //
1513  // where U_m+1 is the current Krylov-Schur basis, Q are the Schur vectors, and x = U_m+1*Q*s
1514  // NOTE: This means that s = e_i if the problem is hermitian, else the eigenvectors
1515  // of the Schur form need to be computed.
1516  //
1517  // First compute H_{m+1,m}*B_m^T, then determine what 's' is.
1518  //---------------------------------------------------
1519  //
1520  // Get current B_m+1
1522  blockSize_, curDim_, curDim_ );
1523  //
1524  // Compute B_m+1^H*Q
1525  Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
1526  blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1527  curB.values(), curB.stride(), subQ.values(), subQ.stride(),
1528  ST_ZERO, subB.values(), subB.stride() );
1529  //
1530  // Determine what 's' is and compute Ritz residuals.
1531  //
1532  ScalarType* b_ptr = subB.values();
1533  if (problem_->isHermitian() && !hermImagDetected) {
1534  //
1535  // 's' is the i-th canonical basis vector.
1536  //
1537  for (int i=0; i<curDim_ ; i++) {
1538  ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
1539  }
1540  } else {
1541  //
1542  // Compute S: the eigenvectors of the block upper triangular, Schur matrix.
1543  //
1544  char side = 'R';
1545  int mm;
1546  const int ldvl = 1;
1547  ScalarType vl[ ldvl ];
1548  Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
1549  lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1550  S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1551 
1552  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1553  "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
1554  //
1555  // Scale the eigenvectors so that their Euclidean norms are all one.
1556  //
1557  HelperTraits<ScalarType>::scaleRitzVectors( tmp_iRitzValues, &S );
1558  //
1559  // Compute ritzRes = *B_m+1^H*Q*S where the i-th column of S is 's' for the i-th Ritz-value
1560  //
1561  Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
1562  blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1563  subB.values(), subB.stride(), S.values(), S.stride(),
1564  ST_ZERO, ritzRes.values(), ritzRes.stride() );
1565 
1566  /* TO DO: There's be an incorrect assumption made in the computation of the Ritz residuals.
1567  This assumption is that the next vector in the Krylov subspace is Euclidean orthonormal.
1568  It may not be normalized using Euclidean norm.
1569  Teuchos::RCP<MV> ritzResVecs = MVT::Clone( *V_, curDim_ );
1570  std::vector<int> curind(blockSize_);
1571  for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
1572  Teuchos::RCP<MV> Vtemp = MVT::CloneView(*V_,curind);
1573 
1574  MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, ritzRes, ST_ZERO, *ritzResVecs );
1575  std::vector<MagnitudeType> ritzResNrms(curDim_);
1576  MVT::MvNorm( *ritzResVecs, ritzResNrms );
1577  i = 0;
1578  while( i < curDim_ ) {
1579  if ( tmp_ritzValues[curDim_+i] != MT_ZERO ) {
1580  ritzResiduals_[i] = lapack_mag.LAPY2( ritzResNrms[i], ritzResNrms[i+1] );
1581  ritzResiduals_[i+1] = ritzResiduals_[i];
1582  i = i+2;
1583  } else {
1584  ritzResiduals_[i] = ritzResNrms[i];
1585  i++;
1586  }
1587  }
1588  */
1589  //
1590  // Compute the Ritz residuals for each Ritz value.
1591  //
1592  HelperTraits<ScalarType>::computeRitzResiduals( tmp_iRitzValues, ritzRes, &ritzResiduals_ );
1593  }
1594  //
1595  // Sort the Ritz values.
1596  //
1597  {
1598 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1599  Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
1600 #endif
1601  int i=0;
1602  if (problem_->isHermitian() && !hermImagDetected) {
1603  //
1604  // Sort using just the real part of the Ritz values.
1605  sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_); // don't catch exception
1606  ritzIndex_.clear();
1607  while ( i < curDim_ ) {
1608  // The Ritz value is not complex.
1609  ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
1610  ritzIndex_.push_back(0);
1611  i++;
1612  }
1613  }
1614  else {
1615  //
1616  // Sort using both the real and imaginary parts of the Ritz values.
1617  sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
1618  HelperTraits<ScalarType>::sortRitzValues( tmp_rRitzValues, tmp_iRitzValues, &ritzValues_, &ritzOrder_, &ritzIndex_ );
1619  }
1620  //
1621  // Sort the ritzResiduals_ based on the ordering from the Sort Manager.
1622  std::vector<MagnitudeType> ritz2( curDim_ );
1623  for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
1624  blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
1625 
1626  // The Ritz values have now been updated.
1627  ritzValsCurrent_ = true;
1628  }
1629 
1630  } // if (!ritzValsCurrent_) ...
1631  //
1632  //---------------------------------------------------
1633  // * The Ritz values and residuals have been updated at this point.
1634  //
1635  // * The Schur factorization of the projected matrix has been computed,
1636  // and is stored in (schurH_, Q_).
1637  //
1638  // Now the Schur factorization needs to be sorted.
1639  //---------------------------------------------------
1640  //
1641  // Sort the Schur form using the ordering from the Sort Manager.
1642  if (sort) {
1643  sortSchurForm( *schurH_, *Q_, ritzOrder_ );
1644  //
1645  // Indicate the Schur form in (schurH_, Q_) is current and sorted
1646  schurCurrent_ = true;
1647  }
1648  } // if (!schurCurrent_) ...
1649 
1650  } // if (curDim_) ...
1651 
1652  } // computeSchurForm( ... )
1653 
1654 
1656  // Sort the Schur form of H stored in (H,Q) using the ordering vector.
1657  template <class ScalarType, class MV, class OP>
1660  std::vector<int>& order )
1661  {
1662  // local timer
1663 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1664  Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
1665 #endif
1666  //
1667  //---------------------------------------------------
1668  // Reorder real Schur factorization, remember to add one to the indices for the
1669  // fortran call and determine offset. The offset is necessary since the TREXC
1670  // method reorders in a nonsymmetric fashion, thus we use the reordering in
1671  // a stack-like fashion. Also take into account conjugate pairs, which may mess
1672  // up the reordering, since the pair is moved if one of the pair is moved.
1673  //---------------------------------------------------
1674  //
1675  int i = 0, nevtemp = 0;
1676  char compq = 'V';
1677  std::vector<int> offset2( curDim_ );
1678  std::vector<int> order2( curDim_ );
1679 
1680  // LAPACK objects.
1682  int lwork = 3*curDim_;
1683  std::vector<ScalarType> work( lwork );
1684 
1685  while (i < curDim_) {
1686  if ( ritzIndex_[i] != 0 ) { // This is the first value of a complex conjugate pair
1687  offset2[nevtemp] = 0;
1688  for (int j=i; j<curDim_; j++) {
1689  if (order[j] > order[i]) { offset2[nevtemp]++; }
1690  }
1691  order2[nevtemp] = order[i];
1692  i = i+2;
1693  } else {
1694  offset2[nevtemp] = 0;
1695  for (int j=i; j<curDim_; j++) {
1696  if (order[j] > order[i]) { offset2[nevtemp]++; }
1697  }
1698  order2[nevtemp] = order[i];
1699  i++;
1700  }
1701  nevtemp++;
1702  }
1703  ScalarType *ptr_h = H.values();
1704  ScalarType *ptr_q = Q.values();
1705  int ldh = H.stride(), ldq = Q.stride();
1706  int info = 0;
1707  for (i=nevtemp-1; i>=0; i--) {
1708  int ifst = order2[i]+1+offset2[i];
1709  int ilst = 1;
1710  lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, &ifst,
1711  &ilst, &work[0], &info );
1712  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1713  "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ << ") returned info " << info << " != 0.");
1714  }
1715  }
1716 
1718  // Print the current status of the solver
1719  template <class ScalarType, class MV, class OP>
1721  {
1722  using std::endl;
1723 
1724  os.setf(std::ios::scientific, std::ios::floatfield);
1725  os.precision(6);
1726  os <<"================================================================================" << endl;
1727  os << endl;
1728  os <<" BlockKrylovSchur Solver Status" << endl;
1729  os << endl;
1730  os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1731  os <<"The number of iterations performed is " <<iter_<<endl;
1732  os <<"The block size is " << blockSize_<<endl;
1733  os <<"The number of blocks is " << numBlocks_<<endl;
1734  os <<"The current basis size is " << curDim_<<endl;
1735  os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1736  os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1737 
1738  os.setf(std::ios_base::right, std::ios_base::adjustfield);
1739 
1740  os << endl;
1741  if (initialized_) {
1742  os <<"CURRENT RITZ VALUES "<<endl;
1743  if (ritzIndex_.size() != 0) {
1744  int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
1745  if (problem_->isHermitian()) {
1746  os << std::setw(20) << "Ritz Value"
1747  << std::setw(20) << "Ritz Residual"
1748  << endl;
1749  os <<"--------------------------------------------------------------------------------"<<endl;
1750  for (int i=0; i<numPrint; i++) {
1751  os << std::setw(20) << ritzValues_[i].realpart
1752  << std::setw(20) << ritzResiduals_[i]
1753  << endl;
1754  }
1755  } else {
1756  os << std::setw(24) << "Ritz Value"
1757  << std::setw(30) << "Ritz Residual"
1758  << endl;
1759  os <<"--------------------------------------------------------------------------------"<<endl;
1760  for (int i=0; i<numPrint; i++) {
1761  // Print out the real eigenvalue.
1762  os << std::setw(15) << ritzValues_[i].realpart;
1763  if (ritzValues_[i].imagpart < MT_ZERO) {
1764  os << " - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
1765  } else {
1766  os << " + i" << std::setw(15) << ritzValues_[i].imagpart;
1767  }
1768  os << std::setw(20) << ritzResiduals_[i] << endl;
1769  }
1770  }
1771  } else {
1772  os << std::setw(20) << "[ NONE COMPUTED ]" << endl;
1773  }
1774  }
1775  os << endl;
1776  os <<"================================================================================" << endl;
1777  os << endl;
1778  }
1779 
1780 } // End of namespace Anasazi
1781 
1782 #endif
1783 
1784 // End of file AnasaziBlockKrylovSchur.hpp
ScalarType * values() const
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.
BlockKrylovSchurState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
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.
Array< T > & append(const T &x)
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.
T & get(const std::string &name, T def_value)
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
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.
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
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...
void TREVC(const char &SIDE, const char &HOWMNY, OrdinalType *select, const OrdinalType &n, const ScalarType *T, const OrdinalType &ldt, ScalarType *VL, const OrdinalType &ldvl, ScalarType *VR, const OrdinalType &ldvr, const OrdinalType &mm, OrdinalType *m, ScalarType *WORK, OrdinalType *info) const
int getNumIters() const
Get the current iteration count.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
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.
TEUCHOS_DEPRECATED void TREXC(const char &COMPQ, const OrdinalType &n, ScalarType *T, const OrdinalType &ldt, ScalarType *Q, const OrdinalType &ldq, const OrdinalType &ifst, const OrdinalType &ilst, ScalarType *WORK, OrdinalType *info) const
void computeRitzVectors()
Compute the Ritz vectors using the current Krylov factorization.
int getStepSize() const
Get the step size.
bool isParameter(const std::string &name) const
void setNumRitzVectors(int numRitzVecs)
Set the number of Ritz vectors to compute.
void COPY(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isSchurCurrent() const
Get the status of the Schur form currently stored in the eigensolver.
void GEES(const char &JOBVS, const char &SORT, OrdinalType &(*ptr2func)(ScalarType *, ScalarType *), const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *sdim, ScalarType *WR, ScalarType *WI, ScalarType *VS, const OrdinalType &ldvs, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *BWORK, OrdinalType *info) const
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.
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.
ScalarType LAPY2(const ScalarType &x, const ScalarType &y) const
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values.
std::vector< int > getRitzIndex()
Get the Ritz index vector.
static magnitudeType magnitude(T a)
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.
OrdinalType stride() const