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