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