Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
ScalarType * values() const
void TGEVC(const char &SIDE, const char &HOWMNY, const OrdinalType *SELECT, const OrdinalType &n, const ScalarType *S, const OrdinalType &lds, const 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
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxilliary vectors.
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > tester)
Set status test.
RCP< MV > V
Orthonormal basis for search subspace.
std::vector< MagnitudeType > getRitzRes2Norms()
Get the current Ritz residual norms (2-norm)
Teuchos::ScalarTraits< ScalarType >::magnitudeType imagpart
The imaginary component of the eigenvalue.
std::vector< int > getBlockIndex() const
Get indices of current block.
This class defines the interface required by an eigensolver and status test class to compute solution...
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.
int getMaxSubspaceDim() const
Get maximum subspace dimension.
void setBlockSize(int blockSize)
Set block size.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
int curDim
The current subspace dimension.
std::vector< int > getRitzIndex()
Get the current Ritz index vector.
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
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get eigenproblem.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
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 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
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.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
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
Output managers remove the need for the eigensolver to know any information about the required output...
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< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get status test.
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.
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...
ScalarType LAPY2(const ScalarType &x, const ScalarType &y) const
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.
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
iterator end()
OrdinalType numCols() const
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)
bool isInitialized() const
Query if the solver is in an initialized state.
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.
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
int getBlockSize() const
Get block size.
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
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.
Common interface of stopping criteria for Anasazi&#39;s solvers.
iterator begin()
int getNumIters() const
Get number of iterations.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
int getCurSubspaceDim() const
Get current subspace dimension.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
OrdinalType stride() const
OrdinalType numRows() const
Declaration and definition of Anasazi::StatusTest.
void resetNumIters()
Reset the number of iterations.