Anasazi Version of the Day
Loading...
Searching...
No Matches
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"
52#include "Teuchos_ParameterList.hpp"
53#include "Teuchos_SerialDenseMatrix.hpp"
54#include "Teuchos_SerialDenseVector.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"
67#include "AnasaziStatusTest.hpp"
68
69using Teuchos::RCP;
70
71namespace Anasazi {
72
76template <class ScalarType, class MV>
79 int curDim;
80
82 RCP<MV> V;
83
85 RCP<MV> AV;
86
88 RCP<MV> BV;
89
91 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > VAV;
92
94 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > VBV;
95
97 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > S;
98
100 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > T;
101
103 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > Q;
104
106 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > Z;
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
140template <class ScalarType, class MV, class OP>
141class GeneralizedDavidson : public Eigensolver<ScalarType,MV,OP>
142{
143 private:
144 // Convenience Typedefs
147 typedef Teuchos::ScalarTraits<ScalarType> ST;
148 typedef typename ST::magnitudeType MagnitudeType;
149 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
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,
178 Teuchos::ParameterList &pl);
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
214 RCP<const MV> getRitzVectors()
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,
362 Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
363 Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta );
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
381 RCP<Eigenproblem<ScalarType,MV,OP> > d_problem;
382 Teuchos::ParameterList d_pl;
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;
402 RCP<StatusTest<ScalarType,MV,OP> > d_tester;
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
420 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VAV;
421 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_VBV;
422 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_S;
423 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_T;
424 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Q;
425 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > d_Z;
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//---------------------------------------------------------------------------//
462template <class MagnitudeType, class MV, class OP>
463class GeneralizedDavidson<std::complex<MagnitudeType>,MV,OP>
464{
465 public:
466
467 typedef std::complex<MagnitudeType> ScalarType;
468 GeneralizedDavidson(
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,
474 Teuchos::ParameterList &pl)
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//---------------------------------------------------------------------------//
488template <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,
495 Teuchos::ParameterList &pl )
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//---------------------------------------------------------------------------//
566template <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//---------------------------------------------------------------------------//
623template <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//---------------------------------------------------------------------------//
644template <class ScalarType, class MV, class OP>
646{
647 setSize(blockSize,d_maxSubspaceDim);
648}
649
650//---------------------------------------------------------------------------//
651// Set block size and maximum subspace dimension.
652//---------------------------------------------------------------------------//
653template <class ScalarType, class MV, class OP>
654void 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
682 d_VAV = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
683 d_S = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
684 d_Q = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
685 d_Z = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
686
687 // If this is generalized eigenproblem, allocate B components
688 if( d_haveB )
689 {
690 d_BV = MVT::Clone(*initVec, msd);
691 d_VBV = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,msd) );
692 d_T = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(msd,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//---------------------------------------------------------------------------//
719template <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//---------------------------------------------------------------------------//
886template <class ScalarType, class MV, class OP>
888{
890 initialize( empty );
891}
892
893//---------------------------------------------------------------------------//
894// Get current residual norms
895//---------------------------------------------------------------------------//
896template <class ScalarType, class MV, class OP>
897std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
899{
900 return getResNorms(d_residualSize);
901}
902
903//---------------------------------------------------------------------------//
904// Get current residual norms
905//---------------------------------------------------------------------------//
906template <class ScalarType, class MV, class OP>
907std::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//---------------------------------------------------------------------------//
925template <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//---------------------------------------------------------------------------//
953template <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//---------------------------------------------------------------------------//
973template <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
1004 Teuchos::LAPACK<int,ScalarType> lapack;
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
1034 Teuchos::LAPACK<int,ScalarType> lapack;
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
1045 TEUCHOS_TEST_FOR_EXCEPTION(
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.
1055 TEUCHOS_TEST_FOR_EXCEPTION(
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
1061 TEUCHOS_TEST_FOR_EXCEPTION(
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//---------------------------------------------------------------------------//
1077template <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//---------------------------------------------------------------------------//
1134template <class ScalarType, class MV, class OP>
1135void 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//---------------------------------------------------------------------------//
1161template <class ScalarType, class MV, class OP>
1162void 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//---------------------------------------------------------------------------//
1232template <class ScalarType, class MV, class OP>
1233void 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;
1249 Teuchos::LAPACK<int,ScalarType> lapack;
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
1284 Teuchos::LAPACK<int,ScalarType> lapack;
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//---------------------------------------------------------------------------//
1307template <class ScalarType, class MV, class OP>
1308void 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//---------------------------------------------------------------------------//
1340template <class ScalarType, class MV, class OP>
1341void 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
1355 Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
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 );
1376 Teuchos::LAPACK<int,ScalarType> lapack;
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//---------------------------------------------------------------------------//
1399template <class ScalarType, class MV, class OP>
1400void 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//---------------------------------------------------------------------------//
1431template <class ScalarType, class MV, class OP>
1432void GeneralizedDavidson<ScalarType,MV,OP>::computeProjectedEigenvectors(
1433 Teuchos::SerialDenseMatrix<int,ScalarType> &X )
1434{
1435 int N = X.numRows();
1436 if( d_haveB )
1437 {
1438 Teuchos::SerialDenseMatrix<int,ScalarType> S(Teuchos::Copy, *d_S, N, N);
1439 Teuchos::SerialDenseMatrix<int,ScalarType> T(Teuchos::Copy, *d_T, N, N);
1440 Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Q, N, N);
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
1449 Teuchos::LAPACK<int,ScalarType> lapack;
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 {
1459 Teuchos::SerialDenseMatrix<int,ScalarType> S(Teuchos::Copy, *d_S, N, N);
1460 Teuchos::SerialDenseMatrix<int,ScalarType> VL(Teuchos::Copy, *d_Z, N, N);
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
1469 Teuchos::LAPACK<int,ScalarType> lapack;
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//---------------------------------------------------------------------------//
1483template <class ScalarType, class MV, class OP>
1484void GeneralizedDavidson<ScalarType,MV,OP>::scaleEigenvectors(
1485 const Teuchos::SerialDenseMatrix<int,ScalarType> &X,
1486 Teuchos::SerialDenseMatrix<int,ScalarType> &X_alpha,
1487 Teuchos::SerialDenseMatrix<int,ScalarType> &X_beta )
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
1507 Teuchos::SerialDenseMatrix<int,ScalarType> Alpha(Nc,Nc,true);
1508 Teuchos::SerialDenseMatrix<int,ScalarType> Beta(Nc,Nc,true);
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//---------------------------------------------------------------------------//
1551template <class ScalarType, class MV, class OP>
1552void 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
1573 Teuchos::SerialDenseMatrix<int,ScalarType> X(Teuchos::Copy,*d_Z,d_curDim,d_curDim);
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 */
1623 Teuchos::LAPACK<int,ScalarType> lapack;
1624 Teuchos::BLAS<int,ScalarType> blas;
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//---------------------------------------------------------------------------//
1721template <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
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Abstract class definition for Anasazi Output Managers.
Virtual base class which defines the interface between an eigensolver and a class whose job is the so...
Declaration and definition of Anasazi::StatusTest.
Types and exceptions used within Anasazi solvers and interfaces.
This class defines the interface required by an eigensolver and status test class to compute solution...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Solves eigenvalue problem using generalized Davidson method.
int getCurSubspaceDim() const
Get current subspace dimension.
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get status test.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get eigenproblem.
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > tester)
Set status test.
void currentStatus(std::ostream &myout)
Print current status of solver.
int getMaxSubspaceDim() const
Get maximum subspace dimension.
std::vector< MagnitudeType > getResNorms()
Get the current residual norms (w.r.t. norm defined by OrthoManager)
std::vector< int > getBlockIndex() const
Get indices of current block.
std::vector< int > getRitzIndex()
Get the current Ritz index vector.
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 initialize()
Initialize the eigenvalue problem.
GeneralizedDavidsonState< ScalarType, MV > getState()
Get the current state of the eigensolver.
bool isInitialized() const
Query if the solver is in an initialized state.
void setAuxVecs(const Teuchos::Array< RCP< const MV > > &auxVecs)
Set auxilliary vectors.
std::vector< Value< ScalarType > > getRitzValues()
Get the current Ritz values.
int getNumIters() const
Get number of iterations.
std::vector< MagnitudeType > getRes2Norms()
Get the current residual norms (2-norm)
void resetNumIters()
Reset the number of iterations.
RCP< const MV > getRitzVectors()
Get the current Ritz vectors.
int getBlockSize() const
Get block size.
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxilliary vectors.
void iterate()
Solves the eigenvalue problem.
std::vector< MagnitudeType > getRitzRes2Norms()
Get the current Ritz residual norms (2-norm)
void setSize(int blockSize, int maxSubDim)
Set problem size.
void setBlockSize(int blockSize)
Set block size.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Output managers remove the need for the eigensolver to know any information about the required output...
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
@ IterationDetails
Structure to contain pointers to GeneralizedDavidson state variables.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > T
Right quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VAV
Projection of A onto V.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Q
Left generalized Schur vectors from QZ decomposition of (VAV,VBV)
int curDim
The current subspace dimension.
std::vector< Value< ScalarType > > eVals
Vector of generalized eigenvalues.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VBV
Projection of B onto V.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Z
Right generalized Schur vectors from QZ decomposition of (VAV,VBV)
RCP< MV > V
Orthonormal basis for search subspace.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > S
Left quasi upper triangular matrix from QZ decomposition of (VAV,VBV)
This struct is used for storing eigenvalues and Ritz values, as a pair of real values.
Teuchos::ScalarTraits< ScalarType >::magnitudeType imagpart
The imaginary component of the eigenvalue.
Teuchos::ScalarTraits< ScalarType >::magnitudeType realpart
The real component of the eigenvalue.