Anasazi  Version of the Day
AnasaziGeneralizedDavidson.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_GENERALIZED_DAVIDSON_HPP
30 #define ANASAZI_GENERALIZED_DAVIDSON_HPP
31 
38 #include "Teuchos_RCPDecl.hpp"
42 #include "Teuchos_Array.hpp"
43 #include "Teuchos_BLAS.hpp"
44 #include "Teuchos_LAPACK.hpp"
45 
46 #include "AnasaziConfigDefs.hpp"
47 #include "AnasaziTypes.hpp"
48 #include "AnasaziEigenproblem.hpp"
49 #include "AnasaziEigensolver.hpp"
50 #include "AnasaziOrthoManager.hpp"
51 #include "AnasaziOutputManager.hpp"
52 #include "AnasaziSortManager.hpp"
53 #include "AnasaziStatusTest.hpp"
54 
55 using Teuchos::RCP;
56 
57 namespace Anasazi {
58 
62 template <class ScalarType, class MV>
65  int curDim;
66 
69 
72 
75 
78 
81 
84 
87 
90 
93 
95  std::vector< Value<ScalarType> > eVals;
96 
97  GeneralizedDavidsonState() : curDim(0), V(Teuchos::null), AV(Teuchos::null),
98  BV(Teuchos::null), VAV(Teuchos::null),
99  VBV(Teuchos::null), S(Teuchos::null),
100  T(Teuchos::null), Q(Teuchos::null),
101  Z(Teuchos::null), eVals(0) {}
102 
103 };
104 
105 
126 template <class ScalarType, class MV, class OP>
127 class GeneralizedDavidson : public Eigensolver<ScalarType,MV,OP>
128 {
129  private:
130  // Convenience Typedefs
134  typedef typename ST::magnitudeType MagnitudeType;
136 
137  public:
138 
160  const RCP<SortManager<MagnitudeType> > &sortman,
161  const RCP<OutputManager<ScalarType> > &outputman,
162  const RCP<StatusTest<ScalarType,MV,OP> > &tester,
163  const RCP<OrthoManager<ScalarType,MV> > &orthoman,
165 
169  void iterate();
170 
180  void initialize();
181 
186 
190  int getNumIters() const { return d_iteration; }
191 
195  void resetNumIters() { d_iteration=0; d_opApplies=0; }
196 
201  {
202  if( !d_ritzVectorsValid )
203  computeRitzVectors();
204  return d_ritzVecs;
205  }
206 
210  std::vector< Value<ScalarType> > getRitzValues();
211 
215  std::vector<int> getRitzIndex()
216  {
217  if( !d_ritzIndexValid )
218  computeRitzIndex();
219  return d_ritzIndex;
220  }
221 
227  std::vector<int> getBlockIndex() const
228  {
229  return d_expansionIndices;
230  }
231 
235  std::vector<MagnitudeType> getResNorms();
236 
240  std::vector<MagnitudeType> getResNorms(int numWanted);
241 
245  std::vector<MagnitudeType> getRes2Norms() { return d_resNorms; }
246 
253  std::vector<MagnitudeType> getRitzRes2Norms() { return d_resNorms; }
254 
258  int getCurSubspaceDim() const { return d_curDim; }
259 
263  int getMaxSubspaceDim() const { return d_maxSubspaceDim; }
264 
268  void setStatusTest( RCP<StatusTest<ScalarType,MV,OP> > tester) { d_tester = tester; }
269 
273  RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const { return d_tester; }
274 
278  const Eigenproblem<ScalarType,MV,OP> & getProblem() const { return *d_problem; }
279 
283  int getBlockSize() const { return d_expansionSize; }
284 
288  void setBlockSize(int blockSize);
289 
293  void setSize(int blockSize, int maxSubDim);
294 
298  Teuchos::Array< RCP<const MV> > getAuxVecs() const { return d_auxVecs; }
299 
307  void setAuxVecs( const Teuchos::Array< RCP<const MV> > &auxVecs );
308 
312  bool isInitialized() const { return d_initialized; }
313 
317  void currentStatus( std::ostream &myout );
318 
323 
327  void sortProblem( int numWanted );
328 
329  private:
330 
331  // Expand subspace
332  void expandSearchSpace();
333 
334  // Apply Operators
335  void applyOperators();
336 
337  // Update projections
338  void updateProjections();
339 
340  // Solve projected eigenproblem
341  void solveProjectedEigenproblem();
342 
343  // Compute eigenvectors of matrix pair
344  void computeProjectedEigenvectors( Teuchos::SerialDenseMatrix<int,ScalarType> &X );
345 
346  // Scale projected eigenvectors by alpha/beta
347  void scaleEigenvectors( const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
350 
351  // Sort vectors of pairs
352  void sortValues( std::vector<MagnitudeType> &realParts,
353  std::vector<MagnitudeType> &imagParts,
354  std::vector<int> &permVec,
355  int N);
356 
357  // Compute Residual
358  void computeResidual();
359 
360  // Update the current Ritz index vector
361  void computeRitzIndex();
362 
363  // Compute the current Ritz vectors
364  void computeRitzVectors();
365 
366  // Operators
369  RCP<const OP> d_A;
370  RCP<const OP> d_B;
371  RCP<const OP> d_P;
372  bool d_haveB;
373  bool d_haveP;
374 
375  // Parameters
376  int d_blockSize;
377  int d_maxSubspaceDim;
378  int d_NEV;
379  int d_numToPrint;
380  std::string d_initType;
381  int d_verbosity;
382  bool d_relativeConvergence;
383 
384  // Managers
385  RCP<OutputManager<ScalarType> > d_outputMan;
386  RCP<OrthoManager<ScalarType,MV> > d_orthoMan;
387  RCP<SortManager<MagnitudeType> > d_sortMan;
389 
390  // Eigenvalues
391  std::vector< Value<ScalarType> > d_eigenvalues;
392 
393  // Residual Vector
394  RCP<MV> d_R;
395  std::vector<MagnitudeType> d_resNorms;
396 
397  // Subspace Vectors
398  RCP<MV> d_V;
399  RCP<MV> d_AV;
400  RCP<MV> d_BV;
401  RCP<MV> d_ritzVecSpace;
402  RCP<MV> d_ritzVecs;
403  Teuchos::Array< RCP<const MV> > d_auxVecs;
404 
405  // Serial Matrices
412 
413  // Arrays for holding Ritz values
414  std::vector<MagnitudeType> d_alphar;
415  std::vector<MagnitudeType> d_alphai;
416  std::vector<MagnitudeType> d_betar;
417  std::vector<int> d_ritzIndex;
418  std::vector<int> d_convergedIndices;
419  std::vector<int> d_expansionIndices;
420 
421  // Current subspace dimension
422  int d_curDim;
423 
424  // How many vectors are to be added to the subspace
425  int d_expansionSize;
426 
427  // Should subspace expansion use leading vectors
428  // (if false, will use leading unconverged vectors)
429  bool d_useLeading;
430 
431  // What should be used for test subspace (V, AV, or BV)
432  std::string d_testSpace;
433 
434  // How many residual vectors are valid
435  int d_residualSize;
436 
437  int d_iteration;
438  int d_opApplies;
439  bool d_initialized;
440  bool d_ritzIndexValid;
441  bool d_ritzVectorsValid;
442 
443 };
444 
445 //---------------------------------------------------------------------------//
446 // Prevent instantiation on complex scalar type
447 //---------------------------------------------------------------------------//
448 template <class MagnitudeType, class MV, class OP>
449 class GeneralizedDavidson<std::complex<MagnitudeType>,MV,OP>
450 {
451  public:
452 
453  typedef std::complex<MagnitudeType> ScalarType;
455  const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
456  const RCP<SortManager<MagnitudeType> > &sortman,
457  const RCP<OutputManager<ScalarType> > &outputman,
458  const RCP<StatusTest<ScalarType,MV,OP> > &tester,
459  const RCP<OrthoManager<ScalarType,MV> > &orthoman,
461  {
462  // Provide a compile error when attempting to instantiate on complex type
463  MagnitudeType::this_class_is_missing_a_specialization();
464  }
465 };
466 
467 //---------------------------------------------------------------------------//
468 // PUBLIC METHODS
469 //---------------------------------------------------------------------------//
470 
471 //---------------------------------------------------------------------------//
472 // Constructor
473 //---------------------------------------------------------------------------//
474 template <class ScalarType, class MV, class OP>
476  const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
477  const RCP<SortManager<MagnitudeType> > &sortman,
478  const RCP<OutputManager<ScalarType> > &outputman,
479  const RCP<StatusTest<ScalarType,MV,OP> > &tester,
480  const RCP<OrthoManager<ScalarType,MV> > &orthoman,
482 {
483  TEUCHOS_TEST_FOR_EXCEPTION( problem == Teuchos::null, std::invalid_argument, "No Eigenproblem given to solver." );
484  TEUCHOS_TEST_FOR_EXCEPTION( outputman == Teuchos::null, std::invalid_argument, "No OutputManager given to solver." );
485  TEUCHOS_TEST_FOR_EXCEPTION( orthoman == Teuchos::null, std::invalid_argument, "No OrthoManager given to solver." );
486  TEUCHOS_TEST_FOR_EXCEPTION( sortman == Teuchos::null, std::invalid_argument, "No SortManager given to solver." );
487  TEUCHOS_TEST_FOR_EXCEPTION( tester == Teuchos::null, std::invalid_argument, "No StatusTest given to solver." );
488  TEUCHOS_TEST_FOR_EXCEPTION( !problem->isProblemSet(), std::invalid_argument, "Problem has not been set." );
489 
490  d_problem = problem;
491  d_pl = pl;
492  TEUCHOS_TEST_FOR_EXCEPTION( problem->getA()==Teuchos::null && problem->getOperator()==Teuchos::null,
493  std::invalid_argument, "Either A or Operator must be non-null on Eigenproblem");
494  d_A = problem->getA()!=Teuchos::null ? problem->getA() : problem->getOperator();
495  d_B = problem->getM();
496  d_P = problem->getPrec();
497  d_sortMan = sortman;
498  d_outputMan = outputman;
499  d_tester = tester;
500  d_orthoMan = orthoman;
501 
502  // Pull entries from the ParameterList and Eigenproblem
503  d_NEV = d_problem->getNEV();
504  d_initType = d_pl.get<std::string>("Initial Guess","Random");
505  d_numToPrint = d_pl.get<int>("Print Number of Ritz Values",-1);
506  d_useLeading = d_pl.get<bool>("Use Leading Vectors",false);
507 
508  if( d_B != Teuchos::null )
509  d_haveB = true;
510  else
511  d_haveB = false;
512 
513  if( d_P != Teuchos::null )
514  d_haveP = true;
515  else
516  d_haveP = false;
517 
518  d_testSpace = d_pl.get<std::string>("Test Space","V");
519  TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace!="V" && d_testSpace!="AV" && d_testSpace!="BV", std::invalid_argument,
520  "Anasazi::GeneralizedDavidson: Test Space must be V, AV, or BV" );
521  TEUCHOS_TEST_FOR_EXCEPTION( d_testSpace=="V" ? false : !d_haveB, std::invalid_argument,
522  "Anasazi::GeneralizedDavidson: Test Space must be V for standard eigenvalue problem" );
523 
524  // Allocate space for subspace vectors, projected matrices
525  int blockSize = d_pl.get<int>("Block Size",1);
526  int maxSubDim = d_pl.get<int>("Maximum Subspace Dimension",3*d_NEV*blockSize);
527  d_blockSize = -1;
528  d_maxSubspaceDim = -1;
529  setSize( blockSize, maxSubDim );
530  d_relativeConvergence = d_pl.get<bool>("Relative Convergence Tolerance",false);
531 
532  // Make sure subspace size is consistent with requested eigenvalues
533  TEUCHOS_TEST_FOR_EXCEPTION( d_blockSize <= 0, std::invalid_argument, "Block size must be positive");
534  TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim <= 0, std::invalid_argument, "Maximum Subspace Dimension must be positive" );
535  TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV()+2 > pl.get<int>("Maximum Subspace Dimension"),
536  std::invalid_argument, "Maximum Subspace Dimension must be strictly greater than NEV");
537  TEUCHOS_TEST_FOR_EXCEPTION( d_maxSubspaceDim > MVT::GetGlobalLength(*problem->getInitVec()),
538  std::invalid_argument, "Maximum Subspace Dimension cannot exceed problem size");
539 
540 
541  d_curDim = 0;
542  d_iteration = 0;
543  d_opApplies = 0;
544  d_ritzIndexValid = false;
545  d_ritzVectorsValid = false;
546 }
547 
548 
549 //---------------------------------------------------------------------------//
550 // Iterate
551 //---------------------------------------------------------------------------//
552 template <class ScalarType, class MV, class OP>
554 {
555  // Initialize Problem
556  if( !d_initialized )
557  {
558  d_outputMan->stream(Warnings) << "WARNING: GeneralizedDavidson::iterate called without first calling initialize" << std::endl;
559  d_outputMan->stream(Warnings) << " Default initialization will be performed" << std::endl;
560  initialize();
561  }
562 
563  // Print current status
564  if( d_outputMan->isVerbosity(Debug) )
565  {
566  currentStatus( d_outputMan->stream(Debug) );
567  }
568  else if( d_outputMan->isVerbosity(IterationDetails) )
569  {
570  currentStatus( d_outputMan->stream(IterationDetails) );
571  }
572 
573  while( d_tester->getStatus() != Passed && d_curDim+d_expansionSize <= d_maxSubspaceDim )
574  {
575  d_iteration++;
576 
577  expandSearchSpace();
578 
579  applyOperators();
580 
581  updateProjections();
582 
583  solveProjectedEigenproblem();
584 
585  // Make sure the most significant Ritz values are in front
586  // We want the greater of the block size and the number of
587  // requested values, but can't exceed the current dimension
588  int numToSort = std::max(d_blockSize,d_NEV);
589  numToSort = std::min(numToSort,d_curDim);
590  sortProblem( numToSort );
591 
592  computeResidual();
593 
594  // Print current status
595  if( d_outputMan->isVerbosity(Debug) )
596  {
597  currentStatus( d_outputMan->stream(Debug) );
598  }
599  else if( d_outputMan->isVerbosity(IterationDetails) )
600  {
601  currentStatus( d_outputMan->stream(IterationDetails) );
602  }
603  }
604 }
605 
606 //---------------------------------------------------------------------------//
607 // Return the current state struct
608 //---------------------------------------------------------------------------//
609 template <class ScalarType, class MV, class OP>
611 {
613  state.curDim = d_curDim;
614  state.V = d_V;
615  state.AV = d_AV;
616  state.BV = d_BV;
617  state.VAV = d_VAV;
618  state.VBV = d_VBV;
619  state.S = d_S;
620  state.T = d_T;
621  state.Q = d_Q;
622  state.Z = d_Z;
623  state.eVals = getRitzValues();
624  return state;
625 }
626 
627 //---------------------------------------------------------------------------//
628 // Set block size
629 //---------------------------------------------------------------------------//
630 template <class ScalarType, class MV, class OP>
632 {
633  setSize(blockSize,d_maxSubspaceDim);
634 }
635 
636 //---------------------------------------------------------------------------//
637 // Set block size and maximum subspace dimension.
638 //---------------------------------------------------------------------------//
639 template <class ScalarType, class MV, class OP>
640 void GeneralizedDavidson<ScalarType,MV,OP>::setSize(int blockSize, int maxSubDim )
641 {
642  if( blockSize != d_blockSize || maxSubDim != d_maxSubspaceDim )
643  {
644  d_blockSize = blockSize;
645  d_maxSubspaceDim = maxSubDim;
646  d_initialized = false;
647 
648  d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Allocating eigenproblem"
649  << " state with block size of " << d_blockSize
650  << " and maximum subspace dimension of " << d_maxSubspaceDim << std::endl;
651 
652  // Resize arrays for Ritz values
653  d_alphar.resize(d_maxSubspaceDim);
654  d_alphai.resize(d_maxSubspaceDim);
655  d_betar.resize(d_maxSubspaceDim);
656 
657  // Shorten for convenience here
658  int msd = d_maxSubspaceDim;
659 
660  // Temporarily save initialization vector to clone needed vectors
661  RCP<const MV> initVec = d_problem->getInitVec();
662 
663  // Allocate subspace vectors
664  d_V = MVT::Clone(*initVec, msd);
665  d_AV = MVT::Clone(*initVec, msd);
666 
667  // Allocate serial matrices
672 
673  // If this is generalized eigenproblem, allocate B components
674  if( d_haveB )
675  {
676  d_BV = MVT::Clone(*initVec, msd);
679  }
680 
681  /* Allocate space for residual and Ritz vectors
682  * The residual serves two purposes in the Davidson algorithm --
683  * subspace expansion (via the preconditioner) and convergence checking.
684  * We need "Block Size" vectors for subspace expantion and NEV vectors
685  * for convergence checking. Allocate space for max of these, one
686  * extra to avoid splitting conjugate pairs
687  * Allocate one more than "Block Size" to avoid splitting a conjugate pair
688  */
689  d_R = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
690  d_ritzVecSpace = MVT::Clone(*initVec,std::max(d_blockSize,d_NEV)+1);
691  }
692 }
693 
694 //---------------------------------------------------------------------------//
695 /*
696  * Initialize the eigenvalue problem
697  *
698  * Anything on the state that is not null is assumed to be valid.
699  * Anything not present on the state will be generated
700  * Very limited error checking can be performed to ensure the validity of
701  * state components (e.g. we cannot verify that state.AV actually corresponds
702  * to A*state.V), so this function should be used carefully.
703  */
704 //---------------------------------------------------------------------------//
705 template <class ScalarType, class MV, class OP>
707 {
708  // If state has nonzero dimension, we initialize from that, otherwise
709  // we'll pick d_blockSize vectors to start with
710  d_curDim = (state.curDim > 0 ? state.curDim : d_blockSize );
711 
712  d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Initializing state"
713  << " with subspace dimension " << d_curDim << std::endl;
714 
715  // Index for 1st block_size vectors
716  std::vector<int> initInds(d_curDim);
717  for( int i=0; i<d_curDim; ++i )
718  initInds[i] = i;
719 
720  // View of vectors that need to be initialized
721  RCP<MV> V1 = MVT::CloneViewNonConst(*d_V,initInds);
722 
723  // If state's dimension is large enough, use state.V to initialize
724  bool reset_V = false;;
725  if( state.curDim > 0 && state.V != Teuchos::null && MVT::GetNumberVecs(*state.V) >= d_curDim )
726  {
727  if( state.V != d_V )
728  MVT::SetBlock(*state.V,initInds,*V1);
729  }
730  // If there aren't enough vectors in problem->getInitVec() or the user specifically
731  // wants to use random data, set V to random
732  else if( MVT::GetNumberVecs(*d_problem->getInitVec()) < d_blockSize || d_initType == "Random" )
733  {
734  MVT::MvRandom(*V1);
735  reset_V = true;
736  }
737  // Use vectors in problem->getInitVec()
738  else
739  {
740  RCP<const MV> initVec = MVT::CloneView(*d_problem->getInitVec(),initInds);
741  MVT::SetBlock(*initVec,initInds,*V1);
742  reset_V = true;
743  }
744 
745  // If we reset V, it needs to be orthonormalized
746  if( reset_V )
747  {
748  int rank = d_orthoMan->projectAndNormalize( *V1, d_auxVecs );
749  TEUCHOS_TEST_FOR_EXCEPTION( rank < d_blockSize, std::logic_error,
750  "Anasazi::GeneralizedDavidson::initialize(): Error generating initial orthonormal basis" );
751  }
752 
753  if( d_outputMan->isVerbosity(Debug) )
754  {
755  d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Error in V^T V == I: "
756  << d_orthoMan->orthonormError( *V1 ) << std::endl;
757  }
758 
759  // Now process AV
760  RCP<MV> AV1 = MVT::CloneViewNonConst(*d_AV,initInds);
761 
762  // If AV in the state is valid and of appropriate size, use it
763  // We have no way to check that AV is actually A*V
764  if( !reset_V && state.AV != Teuchos::null && MVT::GetNumberVecs(*state.AV) >= d_curDim )
765  {
766  if( state.AV != d_AV )
767  MVT::SetBlock(*state.AV,initInds,*AV1);
768  }
769  // Otherwise apply A to V
770  else
771  {
772  OPT::Apply( *d_A, *V1, *AV1 );
773  d_opApplies += MVT::GetNumberVecs( *V1 );
774  }
775 
776  // Views of matrix to be updated
777  Teuchos::SerialDenseMatrix<int,ScalarType> VAV1( Teuchos::View, *d_VAV, d_curDim, d_curDim );
778 
779  // If the state has a valid VAV, use it
780  if( !reset_V && state.VAV != Teuchos::null && state.VAV->numRows() >= d_curDim && state.VAV->numCols() >= d_curDim )
781  {
782  if( state.VAV != d_VAV )
783  {
784  Teuchos::SerialDenseMatrix<int,ScalarType> state_VAV( Teuchos::View, *state.VAV, d_curDim, d_curDim );
785  VAV1.assign( state_VAV );
786  }
787  }
788  // Otherwise compute VAV from V,AV
789  else
790  {
791  if( d_testSpace == "V" )
792  {
793  MVT::MvTransMv( ST::one(), *V1, *AV1, VAV1 );
794  }
795  else if( d_testSpace == "AV" )
796  {
797  MVT::MvTransMv( ST::one(), *AV1, *AV1, VAV1 );
798  }
799  else if( d_testSpace == "BV" )
800  {
801  RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
802  MVT::MvTransMv( ST::one(), *BV1, *AV1, VAV1 );
803  }
804  }
805 
806  // Process BV if we have it
807  if( d_haveB )
808  {
809  RCP<MV> BV1 = MVT::CloneViewNonConst(*d_BV,initInds);
810 
811  // If BV in the state is valid and of appropriate size, use it
812  // We have no way to check that BV is actually B*V
813  if( !reset_V && state.BV != Teuchos::null && MVT::GetNumberVecs(*state.BV) >= d_curDim )
814  {
815  if( state.BV != d_BV )
816  MVT::SetBlock(*state.BV,initInds,*BV1);
817  }
818  // Otherwise apply B to V
819  else
820  {
821  OPT::Apply( *d_B, *V1, *BV1 );
822  }
823 
824  // Views of matrix to be updated
825  Teuchos::SerialDenseMatrix<int,ScalarType> VBV1( Teuchos::View, *d_VBV, d_curDim, d_curDim );
826 
827  // If the state has a valid VBV, use it
828  if( !reset_V && state.VBV != Teuchos::null && state.VBV->numRows() >= d_curDim && state.VBV->numCols() >= d_curDim )
829  {
830  if( state.VBV != d_VBV )
831  {
832  Teuchos::SerialDenseMatrix<int,ScalarType> state_VBV( Teuchos::View, *state.VBV, d_curDim, d_curDim );
833  VBV1.assign( state_VBV );
834  }
835  }
836  // Otherwise compute VBV from V,BV
837  else
838  {
839  if( d_testSpace == "V" )
840  {
841  MVT::MvTransMv( ST::one(), *V1, *BV1, VBV1 );
842  }
843  else if( d_testSpace == "AV" )
844  {
845  MVT::MvTransMv( ST::one(), *AV1, *BV1, VBV1 );
846  }
847  else if( d_testSpace == "BV" )
848  {
849  MVT::MvTransMv( ST::one(), *BV1, *BV1, VBV1 );
850  }
851  }
852  }
853 
854  // Update Ritz values
855  solveProjectedEigenproblem();
856 
857  // Sort
858  int numToSort = std::max(d_blockSize,d_NEV);
859  numToSort = std::min(numToSort,d_curDim);
860  sortProblem( numToSort );
861 
862  // Get valid residual
863  computeResidual();
864 
865  // Set solver to initialized
866  d_initialized = true;
867 }
868 
869 //---------------------------------------------------------------------------//
870 // Initialize the eigenvalue problem with empty state
871 //---------------------------------------------------------------------------//
872 template <class ScalarType, class MV, class OP>
874 {
876  initialize( empty );
877 }
878 
879 //---------------------------------------------------------------------------//
880 // Get current residual norms
881 //---------------------------------------------------------------------------//
882 template <class ScalarType, class MV, class OP>
883 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
885 {
886  return getResNorms(d_residualSize);
887 }
888 
889 //---------------------------------------------------------------------------//
890 // Get current residual norms
891 //---------------------------------------------------------------------------//
892 template <class ScalarType, class MV, class OP>
893 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
895 {
896  std::vector<int> resIndices(numWanted);
897  for( int i=0; i<numWanted; ++i )
898  resIndices[i]=i;
899 
900  RCP<const MV> R_checked = MVT::CloneView( *d_R, resIndices );
901 
902  std::vector<MagnitudeType> resNorms;
903  d_orthoMan->norm( *R_checked, resNorms );
904 
905  return resNorms;
906 }
907 
908 //---------------------------------------------------------------------------//
909 // Get current Ritz values
910 //---------------------------------------------------------------------------//
911 template <class ScalarType, class MV, class OP>
913 {
914  std::vector< Value<ScalarType> > ritzValues;
915  for( int ival=0; ival<d_curDim; ++ival )
916  {
917  Value<ScalarType> thisVal;
918  thisVal.realpart = d_alphar[ival] / d_betar[ival];
919  if( d_betar[ival] != MT::zero() )
920  thisVal.imagpart = d_alphai[ival] / d_betar[ival];
921  else
922  thisVal.imagpart = MT::zero();
923 
924  ritzValues.push_back( thisVal );
925  }
926 
927  return ritzValues;
928 }
929 
930 //---------------------------------------------------------------------------//
931 /*
932  * Set auxilliary vectors
933  *
934  * Manually setting the auxilliary vectors invalidates the current state
935  * of the solver. Reuse of any components of the solver requires extracting
936  * the state, orthogonalizing V against the aux vecs and reinitializing.
937  */
938 //---------------------------------------------------------------------------//
939 template <class ScalarType, class MV, class OP>
941  const Teuchos::Array< RCP<const MV> > &auxVecs )
942 {
943  d_auxVecs = auxVecs;
944 
945  // Set state to uninitialized if any vectors were set here
946  typename Teuchos::Array< RCP<const MV> >::const_iterator arrItr;
947  int numAuxVecs=0;
948  for( arrItr=auxVecs.begin(); arrItr!=auxVecs.end(); ++arrItr )
949  {
950  numAuxVecs += MVT::GetNumberVecs( *(*arrItr) );
951  }
952  if( numAuxVecs > 0 )
953  d_initialized = false;
954 }
955 
956 //---------------------------------------------------------------------------//
957 // Reorder Schur form, bringing wanted values to front
958 //---------------------------------------------------------------------------//
959 template <class ScalarType, class MV, class OP>
961 {
962  // Get permutation vector
963  std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
964  std::vector< Value<ScalarType> > ritzVals = getRitzValues();
965  for( int i=0; i<d_curDim; ++i )
966  {
967  realRitz[i] = ritzVals[i].realpart;
968  imagRitz[i] = ritzVals[i].imagpart;
969  }
970 
971  std::vector<int> permVec;
972  sortValues( realRitz, imagRitz, permVec, d_curDim );
973 
974  std::vector<int> sel(d_curDim,0);
975  for( int ii=0; ii<numWanted; ++ii )
976  sel[ permVec[ii] ]=1;
977 
978  if( d_haveB )
979  {
980  int ijob = 0; // reorder only, no condition number estimates
981  int wantq = 1; // keep left Schur vectors
982  int wantz = 1; // keep right Schur vectors
983  int work_size=10*d_maxSubspaceDim+16;
984  std::vector<ScalarType> work(work_size);
985  int sdim = 0;
986  int iwork_size = 1;
987  int iwork;
988  int info = 0;
989 
991  lapack.TGSEN( ijob, wantq, wantz, &sel[0], d_curDim, d_S->values(), d_S->stride(), d_T->values(), d_T->stride(),
992  &d_alphar[0], &d_alphai[0], &d_betar[0], d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(),
993  &sdim, 0, 0, 0, &work[0], work_size, &iwork, iwork_size, &info );
994 
995  d_ritzIndexValid = false;
996  d_ritzVectorsValid = false;
997 
998  std::stringstream ss;
999  ss << "Anasazi::GeneralizedDavidson: TGSEN returned error code " << info << std::endl;
1000  TEUCHOS_TEST_FOR_EXCEPTION( info<0, std::runtime_error, ss.str() );
1001  if( info > 0 )
1002  {
1003  // Only issue a warning for positive error code, this usually indicates
1004  // that the system has not been fully reordered, presumably due to ill-conditioning.
1005  // This is usually not detrimental to the calculation.
1006  d_outputMan->stream(Warnings) << "WARNING: " << ss.str() << std::endl;
1007  d_outputMan->stream(Warnings) << " Problem may not be correctly sorted" << std::endl;
1008  }
1009  }
1010  else {
1011  char getCondNum = 'N'; // no condition number estimates
1012  char getQ = 'V'; // keep Schur vectors
1013  int subDim = 0;
1014  int work_size = d_curDim;
1015  std::vector<ScalarType> work(work_size);
1016  int iwork_size = 1;
1017  int iwork = 0;
1018  int info = 0;
1019 
1021  lapack.TRSEN (getCondNum, getQ, &sel[0], d_curDim, d_S->values (),
1022  d_S->stride (), d_Z->values (), d_Z->stride (),
1023  &d_alphar[0], &d_alphai[0], &subDim, 0, 0, &work[0],
1024  work_size, &iwork, iwork_size, &info );
1025 
1026  std::fill( d_betar.begin(), d_betar.end(), ST::one() );
1027 
1028  d_ritzIndexValid = false;
1029  d_ritzVectorsValid = false;
1030 
1032  info < 0, std::runtime_error, "Anasazi::GeneralizedDavidson::"
1033  "sortProblem: LAPACK routine TRSEN returned error code INFO = "
1034  << info << " < 0. This indicates that argument " << -info
1035  << " (counting starts with one) of TRSEN had an illegal value.");
1036 
1037  // LAPACK's documentation suggests that this should only happen
1038  // in the real-arithmetic case, because I only see INFO == 1
1039  // possible for DTRSEN, not for ZTRSEN. Nevertheless, it's
1040  // harmless to check regardless.
1042  info == 1, std::runtime_error, "Anasazi::GeneralizedDavidson::"
1043  "sortProblem: LAPACK routine TRSEN returned error code INFO = 1. "
1044  "This indicates that the reordering failed because some eigenvalues "
1045  "are too close to separate (the problem is very ill-conditioned).");
1046 
1048  info > 1, std::logic_error, "Anasazi::GeneralizedDavidson::"
1049  "sortProblem: LAPACK routine TRSEN returned error code INFO = "
1050  << info << " > 1. This should not be possible. It may indicate an "
1051  "error either in LAPACK itself, or in Teuchos' LAPACK wrapper.");
1052  }
1053 }
1054 
1055 
1056 //---------------------------------------------------------------------------//
1057 // PRIVATE IMPLEMENTATION
1058 //---------------------------------------------------------------------------//
1059 
1060 //---------------------------------------------------------------------------//
1061 // Expand subspace using preconditioner and orthogonalize
1062 //---------------------------------------------------------------------------//
1063 template <class ScalarType, class MV, class OP>
1065 {
1066  // Get indices into relevant portion of residual and
1067  // location to be added to search space
1068  std::vector<int> newIndices(d_expansionSize);
1069  for( int i=0; i<d_expansionSize; ++i )
1070  {
1071  newIndices[i] = d_curDim+i;
1072  }
1073 
1074  // Get indices into pre-existing search space
1075  std::vector<int> curIndices(d_curDim);
1076  for( int i=0; i<d_curDim; ++i )
1077  curIndices[i] = i;
1078 
1079  // Get View of vectors
1080  RCP<MV> V_new = MVT::CloneViewNonConst( *d_V, newIndices);
1081  RCP<const MV> V_cur = MVT::CloneView( *d_V, curIndices);
1082  RCP<const MV> R_active = MVT::CloneView( *d_R, d_expansionIndices);
1083 
1084  if( d_haveP )
1085  {
1086  // Apply Preconditioner to Residual
1087  OPT::Apply( *d_P, *R_active, *V_new );
1088  }
1089  else
1090  {
1091  // Just copy the residual
1092  MVT::SetBlock( *R_active, newIndices, *d_V );
1093  }
1094 
1095  // Normalize new vector against existing vectors in V plus auxVecs
1096  Teuchos::Array< RCP<const MV> > against = d_auxVecs;
1097  against.push_back( V_cur );
1098  int rank = d_orthoMan->projectAndNormalize(*V_new,against);
1099 
1100  if( d_outputMan->isVerbosity(Debug) )
1101  {
1102  std::vector<int> allIndices(d_curDim+d_expansionSize);
1103  for( int i=0; i<d_curDim+d_expansionSize; ++i )
1104  allIndices[i]=i;
1105 
1106  RCP<const MV> V_all = MVT::CloneView( *d_V, allIndices );
1107 
1108  d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidson: Error in V^T V == I: "
1109  << d_orthoMan->orthonormError( *V_all ) << std::endl;
1110  }
1111 
1112  TEUCHOS_TEST_FOR_EXCEPTION( rank != d_expansionSize, std::runtime_error,
1113  "Anasazi::GeneralizedDavidson::ExpandSearchSpace(): Orthonormalization of new vectors failed" );
1114 
1115 }
1116 
1117 //---------------------------------------------------------------------------//
1118 // Apply operators
1119 //---------------------------------------------------------------------------//
1120 template <class ScalarType, class MV, class OP>
1121 void GeneralizedDavidson<ScalarType,MV,OP>::applyOperators()
1122 {
1123  // Get indices for different components
1124  std::vector<int> newIndices(d_expansionSize);
1125  for( int i=0; i<d_expansionSize; ++i )
1126  newIndices[i] = d_curDim+i;
1127 
1128  // Get Views
1129  RCP<const MV> V_new = MVT::CloneView( *d_V, newIndices );
1130  RCP<MV> AV_new = MVT::CloneViewNonConst( *d_AV, newIndices );
1131 
1132  // Multiply by A
1133  OPT::Apply( *d_A, *V_new, *AV_new );
1134  d_opApplies += MVT::GetNumberVecs( *V_new );
1135 
1136  // Multiply by B
1137  if( d_haveB )
1138  {
1139  RCP<MV> BV_new = MVT::CloneViewNonConst( *d_BV, newIndices );
1140  OPT::Apply( *d_B, *V_new, *BV_new );
1141  }
1142 }
1143 
1144 //---------------------------------------------------------------------------//
1145 // Update projected matrices.
1146 //---------------------------------------------------------------------------//
1147 template <class ScalarType, class MV, class OP>
1148 void GeneralizedDavidson<ScalarType,MV,OP>::updateProjections()
1149 {
1150  // Get indices for different components
1151  std::vector<int> newIndices(d_expansionSize);
1152  for( int i=0; i<d_expansionSize; ++i )
1153  newIndices[i] = d_curDim+i;
1154 
1155  std::vector<int> curIndices(d_curDim);
1156  for( int i=0; i<d_curDim; ++i )
1157  curIndices[i] = i;
1158 
1159  std::vector<int> allIndices(d_curDim+d_expansionSize);
1160  for( int i=0; i<d_curDim+d_expansionSize; ++i )
1161  allIndices[i] = i;
1162 
1163  // Test subspace can be V, AV, or BV
1164  RCP<const MV> W_new, W_all;
1165  if( d_testSpace == "V" )
1166  {
1167  W_new = MVT::CloneView(*d_V, newIndices );
1168  W_all = MVT::CloneView(*d_V, allIndices );
1169  }
1170  else if( d_testSpace == "AV" )
1171  {
1172  W_new = MVT::CloneView(*d_AV, newIndices );
1173  W_all = MVT::CloneView(*d_AV, allIndices );
1174  }
1175  else if( d_testSpace == "BV" )
1176  {
1177  W_new = MVT::CloneView(*d_BV, newIndices );
1178  W_all = MVT::CloneView(*d_BV, allIndices );
1179  }
1180 
1181  // Get views of AV
1182  RCP<const MV> AV_new = MVT::CloneView(*d_AV, newIndices);
1183  RCP<const MV> AV_current = MVT::CloneView(*d_AV, curIndices);
1184 
1185  // Last block_size rows of VAV (minus final entry)
1186  Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastrow( Teuchos::View, *d_VAV, d_expansionSize, d_curDim, d_curDim, 0 );
1187  MVT::MvTransMv( ST::one(), *W_new, *AV_current, VAV_lastrow );
1188 
1189  // Last block_size columns of VAV
1190  Teuchos::SerialDenseMatrix<int,ScalarType> VAV_lastcol( Teuchos::View, *d_VAV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
1191  MVT::MvTransMv( ST::one(), *W_all, *AV_new, VAV_lastcol );
1192 
1193  if( d_haveB )
1194  {
1195  // Get views of BV
1196  RCP<const MV> BV_new = MVT::CloneView(*d_BV, newIndices);
1197  RCP<const MV> BV_current = MVT::CloneView(*d_BV, curIndices);
1198 
1199  // Last block_size rows of VBV (minus final entry)
1200  Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastrow( Teuchos::View, *d_VBV, d_expansionSize, d_curDim, d_curDim, 0 );
1201  MVT::MvTransMv( ST::one(), *W_new, *BV_current, VBV_lastrow );
1202 
1203  // Last block_size columns of VBV
1204  Teuchos::SerialDenseMatrix<int,ScalarType> VBV_lastcol( Teuchos::View, *d_VBV, d_curDim+d_expansionSize, d_expansionSize, 0, d_curDim );
1205  MVT::MvTransMv( ST::one(), *W_all, *BV_new, VBV_lastcol );
1206  }
1207 
1208  // All bases are expanded, increase current subspace dimension
1209  d_curDim += d_expansionSize;
1210 
1211  d_ritzIndexValid = false;
1212  d_ritzVectorsValid = false;
1213 }
1214 
1215 //---------------------------------------------------------------------------//
1216 // Solve low dimensional eigenproblem using LAPACK
1217 //---------------------------------------------------------------------------//
1218 template <class ScalarType, class MV, class OP>
1219 void GeneralizedDavidson<ScalarType,MV,OP>::solveProjectedEigenproblem()
1220 {
1221  if( d_haveB )
1222  {
1223  // VAV and VBV need to stay unchanged, GGES will overwrite
1224  // S and T with the triangular matrices from the generalized
1225  // Schur form
1226  d_S->assign(*d_VAV);
1227  d_T->assign(*d_VBV);
1228 
1229  // Get QZ Decomposition of Projected Problem
1230  char leftVecs = 'V'; // compute left vectors
1231  char rightVecs = 'V'; // compute right vectors
1232  char sortVals = 'N'; // don't sort
1233  int sdim;
1234  // int work_size = 10*d_curDim+16;
1236  int info;
1237  // workspace query
1238  int work_size = -1;
1239  std::vector<ScalarType> work(1);
1240  lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
1241  d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
1242  d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
1243  // actual call
1244  work_size = work[0];
1245  work.resize(work_size);
1246  lapack.GGES( leftVecs, rightVecs, sortVals, NULL, d_curDim, d_S->values(), d_S->stride(),
1247  d_T->values(), d_T->stride(), &sdim, &d_alphar[0], &d_alphai[0], &d_betar[0],
1248  d_Q->values(), d_Q->stride(), d_Z->values(), d_Z->stride(), &work[0], work_size, 0, &info );
1249 
1250  d_ritzIndexValid = false;
1251  d_ritzVectorsValid = false;
1252 
1253  std::stringstream ss;
1254  ss << "Anasazi::GeneralizedDavidson: GGES returned error code " << info << std::endl;
1255  TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1256  }
1257  else
1258  {
1259  // VAV needs to stay unchanged, GGES will overwrite
1260  // S with the triangular matrix from the Schur form
1261  d_S->assign(*d_VAV);
1262 
1263  // Get QR Decomposition of Projected Problem
1264  char vecs = 'V'; // compute Schur vectors
1265  int sdim;
1266  int work_size = 3*d_curDim;
1267  std::vector<ScalarType> work(work_size);
1268  int info;
1269 
1271  lapack.GEES( vecs, d_curDim, d_S->values(), d_S->stride(), &sdim, &d_alphar[0], &d_alphai[0],
1272  d_Z->values(), d_Z->stride(), &work[0], work_size, 0, 0, &info);
1273 
1274  std::fill( d_betar.begin(), d_betar.end(), ST::one() );
1275 
1276  d_ritzIndexValid = false;
1277  d_ritzVectorsValid = false;
1278 
1279  std::stringstream ss;
1280  ss << "Anasazi::GeneralizedDavidson: GEES returned error code " << info << std::endl;
1281  TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1282  }
1283 }
1284 
1285 //---------------------------------------------------------------------------//
1286 /*
1287  * Get index vector into current Ritz values/vectors
1288  *
1289  * The current ordering of d_alphar, d_alphai, d_betar will be used.
1290  * Reordering those vectors will invalidate the vector returned here.
1291  */
1292 //---------------------------------------------------------------------------//
1293 template <class ScalarType, class MV, class OP>
1294 void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzIndex()
1295 {
1296  if( d_ritzIndexValid )
1297  return;
1298 
1299  d_ritzIndex.resize( d_curDim );
1300  int i=0;
1301  while( i < d_curDim )
1302  {
1303  if( d_alphai[i] == ST::zero() )
1304  {
1305  d_ritzIndex[i] = 0;
1306  i++;
1307  }
1308  else
1309  {
1310  d_ritzIndex[i] = 1;
1311  d_ritzIndex[i+1] = -1;
1312  i+=2;
1313  }
1314  }
1315  d_ritzIndexValid = true;
1316 }
1317 
1318 //---------------------------------------------------------------------------//
1319 /*
1320  * Compute current Ritz vectors
1321  *
1322  * The current ordering of d_alphar, d_alphai, d_betar will be used.
1323  * Reordering those vectors will invalidate the vector returned here.
1324  */
1325 //---------------------------------------------------------------------------//
1326 template <class ScalarType, class MV, class OP>
1327 void GeneralizedDavidson<ScalarType,MV,OP>::computeRitzVectors()
1328 {
1329  if( d_ritzVectorsValid )
1330  return;
1331 
1332  // Make Ritz indices current
1333  computeRitzIndex();
1334 
1335  // Get indices of converged vector
1336  std::vector<int> checkedIndices(d_residualSize);
1337  for( int ii=0; ii<d_residualSize; ++ii )
1338  checkedIndices[ii] = ii;
1339 
1340  // Get eigenvectors of projected system
1342  computeProjectedEigenvectors( X );
1343 
1344  // Get view of wanted vectors
1345  Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View,X,d_curDim,d_residualSize);
1346 
1347  // Get views of relevant portion of V, evecs
1348  d_ritzVecs = MVT::CloneViewNonConst( *d_ritzVecSpace, checkedIndices );
1349 
1350  std::vector<int> curIndices(d_curDim);
1351  for( int i=0; i<d_curDim; ++i )
1352  curIndices[i] = i;
1353 
1354  RCP<const MV> V_current = MVT::CloneView( *d_V, curIndices );
1355 
1356  // Now form Ritz vector
1357  MVT::MvTimesMatAddMv(ST::one(),*V_current,X_wanted,ST::zero(),*d_ritzVecs);
1358 
1359  // Normalize vectors, conjugate pairs get normalized together
1360  std::vector<MagnitudeType> scale(d_residualSize);
1361  MVT::MvNorm( *d_ritzVecs, scale );
1363  for( int i=0; i<d_residualSize; ++i )
1364  {
1365  if( d_ritzIndex[i] == 0 )
1366  {
1367  scale[i] = 1.0/scale[i];
1368  }
1369  else if( d_ritzIndex[i] == 1 )
1370  {
1371  MagnitudeType nrm = lapack.LAPY2(scale[i],scale[i+1]);
1372  scale[i] = 1.0/nrm;
1373  scale[i+1] = 1.0/nrm;
1374  }
1375  }
1376  MVT::MvScale( *d_ritzVecs, scale );
1377 
1378  d_ritzVectorsValid = true;
1379 
1380 }
1381 
1382 //---------------------------------------------------------------------------//
1383 // Use sort manager to sort generalized eigenvalues
1384 //---------------------------------------------------------------------------//
1385 template <class ScalarType, class MV, class OP>
1386 void GeneralizedDavidson<ScalarType,MV,OP>::sortValues( std::vector<MagnitudeType> &realParts,
1387  std::vector<MagnitudeType> &imagParts,
1388  std::vector<int> &permVec,
1389  int N)
1390 {
1391  permVec.resize(N);
1392 
1393  TEUCHOS_TEST_FOR_EXCEPTION( (int) realParts.size()<N, std::runtime_error,
1394  "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1395  TEUCHOS_TEST_FOR_EXCEPTION( (int) imagParts.size()<N, std::runtime_error,
1396  "Anasazi::GeneralizedDavidson::SortValues: Number of requested sorted values greater than vector length." );
1397 
1398  RCP< std::vector<int> > rcpPermVec = Teuchos::rcpFromRef(permVec);
1399 
1400  d_sortMan->sort( realParts, imagParts, rcpPermVec, N );
1401 
1402  d_ritzIndexValid = false;
1403  d_ritzVectorsValid = false;
1404 }
1405 
1406 //---------------------------------------------------------------------------//
1407 /*
1408  * Compute (right) scaled eigenvectors of a pair of dense matrices
1409  *
1410  * This routine computes the eigenvectors for the generalized eigenvalue
1411  * problem \f$ \beta A x = \alpha B x \f$. The input matrices are the upper
1412  * quasi-triangular matrices S and T from a real QZ decomposition,
1413  * the routine dtgevc will back-calculate the eigenvectors of the original
1414  * pencil (A,B) using the orthogonal matrices Q and Z.
1415  */
1416 //---------------------------------------------------------------------------//
1417 template <class ScalarType, class MV, class OP>
1418 void GeneralizedDavidson<ScalarType,MV,OP>::computeProjectedEigenvectors(
1420 {
1421  int N = X.numRows();
1422  if( d_haveB )
1423  {
1427 
1428  char whichVecs = 'R'; // only need right eigenvectors
1429  char howMany = 'B'; // back-compute eigenvectors of original A,B (we have Schur decomposition here)
1430  int work_size = 6*d_maxSubspaceDim;
1431  std::vector<ScalarType> work(work_size,ST::zero());
1432  int info;
1433  int M;
1434 
1436  lapack.TGEVC( whichVecs, howMany, 0, N, S.values(), S.stride(), T.values(), T.stride(),
1437  VL.values(), VL.stride(), X.values(), X.stride(), N, &M, &work[0], &info );
1438 
1439  std::stringstream ss;
1440  ss << "Anasazi::GeneralizedDavidson: TGEVC returned error code " << info << std::endl;
1441  TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1442  }
1443  else
1444  {
1447 
1448  char whichVecs = 'R'; // only need right eigenvectors
1449  char howMany = 'B'; // back-compute eigenvectors of original A (we have Schur decomposition here)
1450  int sel = 0;
1451  std::vector<ScalarType> work(3*N);
1452  int m;
1453  int info;
1454 
1456 
1457  lapack.TREVC( whichVecs, howMany, &sel, N, S.values(), S.stride(), VL.values(), VL.stride(),
1458  X.values(), X.stride(), N, &m, &work[0], &info );
1459 
1460  std::stringstream ss;
1461  ss << "Anasazi::GeneralizedDavidson: TREVC returned error code " << info << std::endl;
1462  TEUCHOS_TEST_FOR_EXCEPTION( info!=0, std::runtime_error, ss.str() );
1463  }
1464 }
1465 
1466 //---------------------------------------------------------------------------//
1467 // Scale eigenvectors by quasi-diagonal matrices alpha and beta
1468 //---------------------------------------------------------------------------//
1469 template <class ScalarType, class MV, class OP>
1470 void GeneralizedDavidson<ScalarType,MV,OP>::scaleEigenvectors(
1474 {
1475  int Nr = X.numRows();
1476  int Nc = X.numCols();
1477 
1478  TEUCHOS_TEST_FOR_EXCEPTION( Nr>d_curDim, std::logic_error,
1479  "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1480  TEUCHOS_TEST_FOR_EXCEPTION( Nc>d_curDim, std::logic_error,
1481  "Anasazi::GeneralizedDavidson::ScaleEigenvectors: Matrix size exceeds current dimension");
1482  TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numRows()!=Nr, std::logic_error,
1483  "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xalpha does not match X");
1484  TEUCHOS_TEST_FOR_EXCEPTION( X_alpha.numCols()!=Nc, std::logic_error,
1485  "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xalpha does not match X");
1486  TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numRows()!=Nr, std::logic_error,
1487  "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of rows in Xbeta does not match X");
1488  TEUCHOS_TEST_FOR_EXCEPTION( X_beta.numCols()!=Nc, std::logic_error,
1489  "Anasazi::GeneralizedDavidson::ScaleEigenvectors: number of cols in Xbeta does not match X");
1490 
1491  // Now form quasi-diagonal matrices
1492  // containing alpha and beta
1495 
1496  computeRitzIndex();
1497 
1498  for( int i=0; i<Nc; ++i )
1499  {
1500  if( d_ritzIndex[i] == 0 )
1501  {
1502  Alpha(i,i) = d_alphar[i];
1503  Beta(i,i) = d_betar[i];
1504  }
1505  else if( d_ritzIndex[i] == 1 )
1506  {
1507  Alpha(i,i) = d_alphar[i];
1508  Alpha(i,i+1) = d_alphai[i];
1509  Beta(i,i) = d_betar[i];
1510  }
1511  else
1512  {
1513  Alpha(i,i-1) = d_alphai[i];
1514  Alpha(i,i) = d_alphar[i];
1515  Beta(i,i) = d_betar[i];
1516  }
1517  }
1518 
1519  int err;
1520 
1521  // Multiply the eigenvectors by alpha
1522  err = X_alpha.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Alpha, ST::zero() );
1523  std::stringstream astream;
1524  astream << "GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1525  TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, astream.str() );
1526 
1527  // Multiply the eigenvectors by beta
1528  err = X_beta.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), X, Beta, ST::zero() );
1529  std::stringstream bstream;
1530  bstream << "GeneralizedDavidson::ScaleEigenvectors: multiply returned error code " << err;
1531  TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, bstream.str() );
1532 }
1533 
1534 //---------------------------------------------------------------------------//
1535 // Compute residual
1536 //---------------------------------------------------------------------------//
1537 template <class ScalarType, class MV, class OP>
1538 void GeneralizedDavidson<ScalarType,MV,OP>::computeResidual()
1539 {
1540  computeRitzIndex();
1541 
1542  // Determine how many residual vectors need to be computed
1543  d_residualSize = std::max( d_blockSize, d_NEV );
1544  if( d_curDim < d_residualSize )
1545  {
1546  d_residualSize = d_curDim;
1547  }
1548  else if( d_ritzIndex[d_residualSize-1] == 1 )
1549  {
1550  d_residualSize++;
1551  }
1552 
1553  // Get indices of all valid residual vectors
1554  std::vector<int> residualIndices(d_residualSize);
1555  for( int i=0; i<d_residualSize; ++i )
1556  residualIndices[i] = i;
1557 
1558  // X will store (right) eigenvectors of projected system
1560 
1561  // Get eigenvectors of projected problem -- computed from previous Schur decomposition
1562  computeProjectedEigenvectors( X );
1563 
1564  // X_alpha and X_beta will be eigenvectors right-multiplied by alpha, beta (which are quasi-diagonal portions of S,T)
1565  Teuchos::SerialDenseMatrix<int,ScalarType> X_alpha(d_curDim,d_residualSize);
1566  Teuchos::SerialDenseMatrix<int,ScalarType> X_beta(d_curDim,d_residualSize);
1567 
1568  // X_wanted is the wanted portion of X
1569  Teuchos::SerialDenseMatrix<int,ScalarType> X_wanted(Teuchos::View, X, d_curDim, d_residualSize);
1570 
1571  // Scale Eigenvectors by alpha or beta
1572  scaleEigenvectors( X_wanted, X_alpha, X_beta );
1573 
1574  // Get view of residual vector(s)
1575  RCP<MV> R_active = MVT::CloneViewNonConst( *d_R, residualIndices );
1576 
1577  // View of active portion of AV,BV
1578  std::vector<int> activeIndices(d_curDim);
1579  for( int i=0; i<d_curDim; ++i )
1580  activeIndices[i]=i;
1581 
1582  // Compute residual
1583  RCP<const MV> AV_active = MVT::CloneView( *d_AV, activeIndices );
1584  MVT::MvTimesMatAddMv(ST::one(),*AV_active, X_beta, ST::zero(),*R_active);
1585 
1586  if( d_haveB )
1587  {
1588  RCP<const MV> BV_active = MVT::CloneView( *d_BV, activeIndices );
1589  MVT::MvTimesMatAddMv(ST::one(),*BV_active, X_alpha,-ST::one(), *R_active);
1590  }
1591  else
1592  {
1593  RCP<const MV> V_active = MVT::CloneView( *d_V, activeIndices );
1594  MVT::MvTimesMatAddMv(ST::one(),*V_active, X_alpha,-ST::one(), *R_active);
1595  }
1596 
1597  /* Apply a scaling to the residual
1598  * For generalized eigenvalue problems, LAPACK scales eigenvectors
1599  * to have unit length in the infinity norm, we want them to have unit
1600  * length in the 2-norm. For conjugate pairs, the scaling is such that
1601  * |xr|^2 + |xi|^2 = 1
1602  * Additionally, the residual is currently computed as r=beta*A*x-alpha*B*x
1603  * but the "standard" residual is r=A*x-(alpha/beta)*B*x, or if we want
1604  * to scale the residual by the Ritz value then it is r=(beta/alpha)*A*x-B*x
1605  * Performing the scaling this way allows us to avoid the possibility of
1606  * diving by infinity or zero if the StatusTest were allowed to handle the
1607  * scaling.
1608  */
1611  std::vector<MagnitudeType> resScaling(d_residualSize);
1612  for( int icol=0; icol<d_residualSize; ++icol )
1613  {
1614  if( d_ritzIndex[icol] == 0 )
1615  {
1616  MagnitudeType Xnrm = blas.NRM2( d_curDim, X_wanted[icol], 1);
1617  MagnitudeType ABscaling = d_relativeConvergence ? d_alphar[icol] : d_betar[icol];
1618  resScaling[icol] = MT::one() / (Xnrm * ABscaling);
1619  }
1620  else if( d_ritzIndex[icol] == 1 )
1621  {
1622  MagnitudeType Xnrm1 = blas.NRM2( d_curDim, X_wanted[icol], 1 );
1623  MagnitudeType Xnrm2 = blas.NRM2( d_curDim, X_wanted[icol+1], 1 );
1624  MagnitudeType Xnrm = lapack.LAPY2(Xnrm1,Xnrm2);
1625  MagnitudeType ABscaling = d_relativeConvergence ? lapack.LAPY2(d_alphar[icol],d_alphai[icol])
1626  : d_betar[icol];
1627  resScaling[icol] = MT::one() / (Xnrm * ABscaling);
1628  resScaling[icol+1] = MT::one() / (Xnrm * ABscaling);
1629  }
1630  }
1631  MVT::MvScale( *R_active, resScaling );
1632 
1633  // Compute residual norms
1634  d_resNorms.resize(d_residualSize);
1635  MVT::MvNorm(*R_active,d_resNorms);
1636 
1637  // If Ritz value i is real, then the corresponding residual vector
1638  // is the true residual
1639  // If Ritz values i and i+1 form a conjugate pair, then the
1640  // corresponding residual vectors are the real and imaginary components
1641  // of the residual. Adjust the residual norms appropriately...
1642  for( int i=0; i<d_residualSize; ++i )
1643  {
1644  if( d_ritzIndex[i] == 1 )
1645  {
1646  MagnitudeType nrm = lapack.LAPY2(d_resNorms[i],d_resNorms[i+1]);
1647  d_resNorms[i] = nrm;
1648  d_resNorms[i+1] = nrm;
1649  }
1650  }
1651 
1652  // Evaluate with status test
1653  d_tester->checkStatus(this);
1654 
1655  // Determine which residual vectors should be used for subspace expansion
1656  if( d_useLeading || d_blockSize >= d_NEV )
1657  {
1658  d_expansionSize=d_blockSize;
1659  if( d_ritzIndex[d_blockSize-1]==1 )
1660  d_expansionSize++;
1661 
1662  d_expansionIndices.resize(d_expansionSize);
1663  for( int i=0; i<d_expansionSize; ++i )
1664  d_expansionIndices[i] = i;
1665  }
1666  else
1667  {
1668  std::vector<int> convergedVectors = d_tester->whichVecs();
1669 
1670  // Get index of first unconverged vector
1671  int startVec;
1672  for( startVec=0; startVec<d_residualSize; ++startVec )
1673  {
1674  if( std::find(convergedVectors.begin(),convergedVectors.end(),startVec)==convergedVectors.end() )
1675  break;
1676  }
1677 
1678  // Now get a contiguous block of indices starting at startVec
1679  // If this crosses the end of our residual vectors, take the final d_blockSize vectors
1680  int endVec = startVec + d_blockSize - 1;
1681  if( endVec > (d_residualSize-1) )
1682  {
1683  endVec = d_residualSize-1;
1684  startVec = d_residualSize-d_blockSize;
1685  }
1686 
1687  // Don't split conjugate pairs on either end of the range
1688  if( d_ritzIndex[startVec]==-1 )
1689  {
1690  startVec--;
1691  endVec--;
1692  }
1693 
1694  if( d_ritzIndex[endVec]==1 )
1695  endVec++;
1696 
1697  d_expansionSize = 1+endVec-startVec;
1698  d_expansionIndices.resize(d_expansionSize);
1699  for( int i=0; i<d_expansionSize; ++i )
1700  d_expansionIndices[i] = startVec+i;
1701  }
1702 }
1703 
1704 //---------------------------------------------------------------------------//
1705 // Print current status.
1706 //---------------------------------------------------------------------------//
1707 template <class ScalarType, class MV, class OP>
1709 {
1710  using std::endl;
1711 
1712  myout.setf(std::ios::scientific, std::ios::floatfield);
1713  myout.precision(6);
1714  myout <<endl;
1715  myout <<"================================================================================" << endl;
1716  myout << endl;
1717  myout <<" GeneralizedDavidson Solver Solver Status" << endl;
1718  myout << endl;
1719  myout <<"The solver is "<<(d_initialized ? "initialized." : "not initialized.") << endl;
1720  myout <<"The number of iterations performed is " << d_iteration << endl;
1721  myout <<"The number of operator applies performed is " << d_opApplies << endl;
1722  myout <<"The block size is " << d_expansionSize << endl;
1723  myout <<"The current basis size is " << d_curDim << endl;
1724  myout <<"The number of requested eigenvalues is " << d_NEV << endl;
1725  myout <<"The number of converged values is " << d_tester->howMany() << endl;
1726  myout << endl;
1727 
1728  myout.setf(std::ios_base::right, std::ios_base::adjustfield);
1729 
1730  if( d_initialized )
1731  {
1732  myout << "CURRENT RITZ VALUES" << endl;
1733 
1734  myout << std::setw(24) << "Ritz Value"
1735  << std::setw(30) << "Residual Norm" << endl;
1736  myout << "--------------------------------------------------------------------------------" << endl;
1737  if( d_residualSize > 0 )
1738  {
1739  std::vector<MagnitudeType> realRitz(d_curDim), imagRitz(d_curDim);
1740  std::vector< Value<ScalarType> > ritzVals = getRitzValues();
1741  for( int i=0; i<d_curDim; ++i )
1742  {
1743  realRitz[i] = ritzVals[i].realpart;
1744  imagRitz[i] = ritzVals[i].imagpart;
1745  }
1746  std::vector<int> permvec;
1747  sortValues( realRitz, imagRitz, permvec, d_curDim );
1748 
1749  int numToPrint = std::max( d_numToPrint, d_NEV );
1750  numToPrint = std::min( d_curDim, numToPrint );
1751 
1752  // Because the sort manager does not use a stable sort, occasionally
1753  // the portion of a conjugate pair with negative imaginary part will be placed
1754  // first...in that case the following will not give the usual expected behavior
1755  // and an extra value will be printed. This is only an issue with the output
1756  // format because the actually sorting of Schur forms is guaranteed to be stable.
1757  if( d_ritzIndex[permvec[numToPrint-1]] != 0 )
1758  numToPrint++;
1759 
1760  int i=0;
1761  while( i<numToPrint )
1762  {
1763  if( imagRitz[i] == ST::zero() )
1764  {
1765  myout << std::setw(15) << realRitz[i];
1766  myout << " + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1767  if( i < d_residualSize )
1768  myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1769  else
1770  myout << " Not Computed" << endl;
1771 
1772  i++;
1773  }
1774  else
1775  {
1776  // Positive imaginary part
1777  myout << std::setw(15) << realRitz[i];
1778  myout << " + i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1779  if( i < d_residualSize )
1780  myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1781  else
1782  myout << " Not Computed" << endl;
1783 
1784  // Negative imaginary part
1785  myout << std::setw(15) << realRitz[i];
1786  myout << " - i" << std::setw(15) << ST::magnitude( imagRitz[i] );
1787  if( i < d_residualSize )
1788  myout << std::setw(20) << d_resNorms[permvec[i]] << endl;
1789  else
1790  myout << " Not Computed" << endl;
1791 
1792  i+=2;
1793  }
1794  }
1795  }
1796  else
1797  {
1798  myout << std::setw(20) << "[ NONE COMPUTED ]" << endl;
1799  }
1800  }
1801  myout << endl;
1802  myout << "================================================================================" << endl;
1803  myout << endl;
1804 }
1805 
1806 } // namespace Anasazi
1807 
1808 #endif // ANASAZI_GENERALIZED_DAVIDSON_HPP
1809 
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > tester)
Set status test.
RCP< MV > V
Orthonormal basis for search subspace.
ScalarType LAPY2(const ScalarType x, const ScalarType y) const
std::vector< MagnitudeType > getRitzRes2Norms()
Get the current Ritz residual norms (2-norm)
ScalarType * values() const
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get status test.
Teuchos::ScalarTraits< ScalarType >::magnitudeType imagpart
The imaginary component of the eigenvalue.
This class defines the interface required by an eigensolver and status test class to compute solution...
int getBlockSize() const
Get block size.
void TGEVC(const char SIDE, const char HOWMNY, const OrdinalType *SELECT, const OrdinalType n, ScalarType *S, const OrdinalType lds, ScalarType *P, const OrdinalType ldp, ScalarType *VL, const OrdinalType ldvl, ScalarType *VR, const OrdinalType ldvr, const OrdinalType mm, OrdinalType *M, ScalarType *WORK, OrdinalType *info) const
Solves eigenvalue problem using generalized Davidson method.
std::vector< Value< ScalarType > > getRitzValues()
Get the current Ritz values.
T & get(const std::string &name, T def_value)
void TRSEN(const char JOB, const char COMPQ, const OrdinalType *SELECT, const OrdinalType n, ScalarType *T, const OrdinalType ldt, ScalarType *Q, const OrdinalType ldq, MagnitudeType *WR, MagnitudeType *WI, OrdinalType *M, ScalarType *S, MagnitudeType *SEP, ScalarType *WORK, const OrdinalType lwork, OrdinalType *IWORK, const OrdinalType liwork, OrdinalType *info) const
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > T
Right quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
GeneralizedDavidsonState< ScalarType, MV > getState()
Get the current state of the eigensolver.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VAV
Projection of A onto V.
void setBlockSize(int blockSize)
Set block size.
int curDim
The current subspace dimension.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType n, const ScalarType *x, const OrdinalType incx) const
std::vector< int > getRitzIndex()
Get the current Ritz index vector.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get eigenproblem.
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...
void iterate()
Solves the eigenvalue problem.
std::vector< MagnitudeType > getRes2Norms()
Get the current residual norms (2-norm)
Structure to contain pointers to GeneralizedDavidson state variables.
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
Abstract class definition for Anasazi Output Managers.
std::vector< Value< ScalarType > > eVals
Vector of generalized eigenvalues.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setSize(int blockSize, int maxSubDim)
Set problem size.
OrdinalType numRows() const
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Output managers remove the need for the eigensolver to know any information about the required output...
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxilliary vectors.
std::vector< MagnitudeType > getResNorms()
Get the current residual norms (w.r.t. norm defined by OrthoManager)
Templated virtual class for providing orthogonalization/orthonormalization methods.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Q
Left generalized Schur vectors from QZ decomposition of (VAV,VBV)
Traits class which defines basic operations on multivectors.
void GGES(const char JOBVL, const char JOBVR, const char SORT, OrdinalType(*ptr2func)(ScalarType *, ScalarType *, ScalarType *), const OrdinalType n, ScalarType *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb, OrdinalType *sdim, MagnitudeType *ALPHAR, MagnitudeType *ALPHAI, MagnitudeType *BETA, ScalarType *VL, const OrdinalType ldvl, ScalarType *VR, const OrdinalType ldvr, ScalarType *WORK, const OrdinalType lwork, OrdinalType *BWORK, OrdinalType *info) const
This struct is used for storing eigenvalues and Ritz values, as a pair of real values.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Virtual base class which defines the interface between an eigensolver and a class whose job is the so...
bool isInitialized() const
Query if the solver is in an initialized state.
GeneralizedDavidson(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const RCP< SortManager< MagnitudeType > > &sortman, const RCP< OutputManager< ScalarType > > &outputman, const RCP< StatusTest< ScalarType, MV, OP > > &tester, const RCP< OrthoManager< ScalarType, MV > > &orthoman, Teuchos::ParameterList &pl)
Constructor.
iterator end()
int getCurSubspaceDim() const
Get current subspace dimension.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > S
Left quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VBV
Projection of B onto V.
void initialize()
Initialize the eigenvalue problem.
void push_back(const value_type &x)
int getNumIters() const
Get number of iterations.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Z
Right generalized Schur vectors from QZ decomposition of (VAV,VBV)
void setAuxVecs(const Teuchos::Array< RCP< const MV > > &auxVecs)
Set auxilliary vectors.
void TGSEN(const OrdinalType ijob, const OrdinalType wantq, const OrdinalType wantz, const OrdinalType *SELECT, const OrdinalType n, ScalarType *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb, MagnitudeType *ALPHAR, MagnitudeType *ALPHAI, MagnitudeType *BETA, ScalarType *Q, const OrdinalType ldq, ScalarType *Z, const OrdinalType ldz, OrdinalType *M, MagnitudeType *PL, MagnitudeType *PR, MagnitudeType *DIF, ScalarType *WORK, const OrdinalType lwork, OrdinalType *IWORK, const OrdinalType liwork, OrdinalType *info) const
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
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
OrdinalType stride() const
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
RCP< const MV > getRitzVectors()
Get the current Ritz vectors.
void currentStatus(std::ostream &myout)
Print current status of solver.
Teuchos::ScalarTraits< ScalarType >::magnitudeType realpart
The real component of the eigenvalue.
OrdinalType numCols() const
Common interface of stopping criteria for Anasazi&#39;s solvers.
std::vector< int > getBlockIndex() const
Get indices of current block.
iterator begin()
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
int getMaxSubspaceDim() const
Get maximum subspace dimension.
Declaration and definition of Anasazi::StatusTest.
void resetNumIters()
Reset the number of iterations.