42#ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
43#define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
49#include "Teuchos_ParameterList.hpp"
50#include "Teuchos_RCPDecl.hpp"
91template <
class ScalarType,
class MV,
class OP>
128 Teuchos::ParameterList &pl );
151 typedef Teuchos::ScalarTraits<ScalarType> ST;
152 typedef typename ST::magnitudeType MagnitudeType;
153 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
155 RCP< Eigenproblem<ScalarType,MV,OP> > d_problem;
156 RCP< GeneralizedDavidson<ScalarType,MV,OP> > d_solver;
157 RCP< OutputManager<ScalarType> > d_outputMan;
158 RCP< OrthoManager<ScalarType,MV> > d_orthoMan;
159 RCP< SortManager<MagnitudeType> > d_sortMan;
160 RCP< StatusTest<ScalarType,MV,OP> > d_tester;
169template <
class MagnitudeType,
class MV,
class OP>
170class GeneralizedDavidsonSolMgr<std::complex<MagnitudeType>,MV,OP>
174 typedef std::complex<MagnitudeType> ScalarType;
175 GeneralizedDavidsonSolMgr(
176 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
177 Teuchos::ParameterList &pl )
180 MagnitudeType::this_class_is_missing_a_specialization();
191template <
class ScalarType,
class MV,
class OP>
194 Teuchos::ParameterList &pl )
197 TEUCHOS_TEST_FOR_EXCEPTION( d_problem == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager." );
198 TEUCHOS_TEST_FOR_EXCEPTION( !d_problem->isProblemSet(), std::invalid_argument,
"Problem not set." );
199 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getA() == Teuchos::null &&
200 d_problem->getOperator() == Teuchos::null, std::invalid_argument,
"A operator not supplied on Eigenproblem." );
201 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getInitVec() == Teuchos::null, std::invalid_argument,
"No vector to clone from on Eigenproblem." );
202 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV() <= 0, std::invalid_argument,
"Number of requested eigenvalues must be positive.");
204 if( !pl.isType<
int>(
"Block Size") )
206 pl.set<
int>(
"Block Size",1);
209 if( !pl.isType<
int>(
"Maximum Subspace Dimension") )
211 pl.set<
int>(
"Maximum Subspace Dimension",3*problem->getNEV()*pl.get<
int>(
"Block Size"));
214 if( !pl.isType<
int>(
"Print Number of Ritz Values") )
216 int numToPrint = std::max( pl.get<
int>(
"Block Size"), d_problem->getNEV() );
217 pl.set<
int>(
"Print Number of Ritz Values",numToPrint);
221 MagnitudeType tol = pl.get<MagnitudeType>(
"Convergence Tolerance", MT::eps() );
222 TEUCHOS_TEST_FOR_EXCEPTION( pl.get<MagnitudeType>(
"Convergence Tolerance") <= MT::zero(),
223 std::invalid_argument,
"Convergence Tolerance must be greater than zero." );
226 if( pl.isType<
int>(
"Maximum Restarts") )
228 d_maxRestarts = pl.get<
int>(
"Maximum Restarts");
229 TEUCHOS_TEST_FOR_EXCEPTION( d_maxRestarts < 0, std::invalid_argument,
"Maximum Restarts must be non-negative" );
237 d_restartDim = pl.get<
int>(
"Restart Dimension",d_problem->getNEV());
238 TEUCHOS_TEST_FOR_EXCEPTION( d_restartDim < d_problem->getNEV(),
239 std::invalid_argument,
"Restart Dimension must be at least NEV" );
242 std::string initType;
243 if( pl.isType<std::string>(
"Initial Guess") )
245 initType = pl.get<std::string>(
"Initial Guess");
246 TEUCHOS_TEST_FOR_EXCEPTION( initType!=
"User" && initType!=
"Random", std::invalid_argument,
247 "Initial Guess type must be 'User' or 'Random'." );
256 if( pl.isType<std::string>(
"Which") )
258 which = pl.get<std::string>(
"Which");
259 TEUCHOS_TEST_FOR_EXCEPTION( which!=
"LM" && which!=
"SM" && which!=
"LR" && which!=
"SR" && which!=
"LI" && which!=
"SI",
260 std::invalid_argument,
261 "Which must be one of LM,SM,LR,SR,LI,SI." );
272 std::string ortho = pl.get<std::string>(
"Orthogonalization",
"SVQB");
273 TEUCHOS_TEST_FOR_EXCEPTION( ortho!=
"DGKS" && ortho!=
"SVQB" && ortho!=
"ICGS", std::invalid_argument,
274 "Anasazi::GeneralizedDavidsonSolMgr::constructor: Invalid orthogonalization type" );
280 else if( ortho==
"SVQB" )
284 else if( ortho==
"ICGS" )
290 bool scaleRes =
false;
291 bool failOnNaN =
false;
292 RCP<StatusTest<ScalarType,MV,OP> > resNormTest = Teuchos::rcp(
294 RES_2NORM,scaleRes,failOnNaN) );
301 int osProc = pl.get(
"Output Processor", 0);
304 Teuchos::RCP<Teuchos::FancyOStream> osp;
306 if (pl.isParameter(
"Output Stream")) {
307 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
315 if (pl.isParameter(
"Verbosity")) {
316 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
317 verbosity = pl.get(
"Verbosity", verbosity);
319 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
326 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Building solver" << std::endl;
329 TEUCHOS_TEST_FOR_EXCEPTION(d_solver->getMaxSubspaceDim() < d_restartDim, std::invalid_argument,
330 "The maximum size of the subspace dimension (" << d_solver->getMaxSubspaceDim() <<
") must "
331 "not be smaller than the size of the restart space (" << d_restartDim <<
"). "
332 "Please adjust \"Restart Dimension\" and/or \"Maximum Subspace Dimension\" parameters.");
339template <
class ScalarType,
class MV,
class OP>
344 d_problem->setSolution(sol);
346 d_solver->initialize();
354 if( d_tester->getStatus() ==
Passed )
358 if( restarts == d_maxRestarts )
362 d_solver->sortProblem( d_restartDim );
364 getRestartState( state );
365 d_solver->initialize( state );
371 d_solver->currentStatus(d_outputMan->stream(
FinalSummary));
374 sol.
numVecs = d_tester->howMany();
377 std::vector<int> whichVecs = d_tester->whichVecs();
378 std::vector<int> origIndex = d_solver->getRitzIndex();
382 for(
int i=0; i<sol.
numVecs; ++i )
384 if( origIndex[ whichVecs[i] ] == 1 )
386 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]+1 ) == whichVecs.end() )
388 whichVecs.push_back( whichVecs[i]+1 );
392 else if( origIndex[ whichVecs[i] ] == -1 )
394 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]-1 ) == whichVecs.end() )
396 whichVecs.push_back( whichVecs[i]-1 );
402 if( d_outputMan->isVerbosity(
Debug) )
404 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: "
405 << sol.
numVecs <<
" eigenpairs converged" << std::endl;
409 std::vector< Value<ScalarType> > origVals = d_solver->getRitzValues();
410 std::vector<MagnitudeType> realParts;
411 std::vector<MagnitudeType> imagParts;
412 for(
int i=0; i<sol.
numVecs; ++i )
414 realParts.push_back( origVals[whichVecs[i]].realpart );
415 imagParts.push_back( origVals[whichVecs[i]].imagpart );
418 std::vector<int> permVec(sol.
numVecs);
419 d_sortMan->sort( realParts, imagParts, Teuchos::rcpFromRef(permVec), sol.
numVecs );
422 std::vector<int> newWhich;
423 for(
int i=0; i<sol.
numVecs; ++i )
424 newWhich.push_back( whichVecs[permVec[i]] );
428 for(
int i=0; i<sol.
numVecs; ++i )
440 sol.
index = origIndex;
442 sol.
Evals = d_solver->getRitzValues();
452 for(
int i=0; i<sol.
numVecs; ++i )
454 sol.
index[i] = origIndex[ newWhich[i] ];
455 sol.
Evals[i] = origVals[ newWhich[i] ];
458 sol.
Evecs = MVT::CloneCopy( *(d_solver->getRitzVectors()), newWhich );
460 d_problem->setSolution(sol);
463 if( sol.
numVecs < d_problem->getNEV() )
472template <
class ScalarType,
class MV,
class OP>
476 TEUCHOS_TEST_FOR_EXCEPTION( state.
curDim <= d_restartDim, std::runtime_error,
477 "Anasazi::GeneralizedDavidsonSolMgr: State dimension at restart is smaller than Restart Dimension" );
479 std::vector<int> ritzIndex = d_solver->getRitzIndex();
482 int restartDim = d_restartDim;
483 if( ritzIndex[d_restartDim-1]==1 )
486 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Restarting with "
487 << restartDim <<
" vectors" << std::endl;
496 const Teuchos::SerialDenseMatrix<int,ScalarType> Z_wanted =
497 Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*state.
Z,state.
curDim,restartDim);
500 std::vector<int> allIndices(state.
curDim);
501 for(
int i=0; i<state.
curDim; ++i )
504 RCP<const MV> V_orig = MVT::CloneView( *state.
V, allIndices );
507 std::vector<int> restartIndices(restartDim);
508 for(
int i=0; i<restartDim; ++i )
509 restartIndices[i] = i;
512 RCP<MV> V_restart = MVT::CloneViewNonConst( *state.
V, restartIndices );
515 RCP<MV> restartVecs = MVT::Clone(*state.
V,restartDim);
518 MVT::MvTimesMatAddMv(ST::one(),*V_orig,Z_wanted,ST::zero(),*restartVecs);
519 MVT::SetBlock(*restartVecs,restartIndices,*V_restart);
522 if( d_outputMan->isVerbosity(
Debug) )
524 MagnitudeType orthErr = d_orthoMan->orthonormError(*V_restart);
525 std::stringstream os;
526 os <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Error in V^T V == I after restart : " << orthErr << std::endl;
527 d_outputMan->print(
Debug,os.str());
531 RCP<MV> AV_restart = MVT::CloneViewNonConst( *state.
AV, restartIndices );
532 RCP<const MV> AV_orig = MVT::CloneView( *state.
AV, allIndices );
534 MVT::MvTimesMatAddMv(ST::one(),*AV_orig,Z_wanted,ST::zero(),*restartVecs);
535 MVT::SetBlock(*restartVecs,restartIndices,*AV_restart);
540 const Teuchos::SerialDenseMatrix<int,ScalarType> VAV_orig( Teuchos::View, *state.
VAV, state.
curDim, state.
curDim );
541 Teuchos::SerialDenseMatrix<int,ScalarType> tmpMat(state.
curDim, restartDim);
542 err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VAV_orig, Z_wanted, ST::zero() );
543 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
545 Teuchos::SerialDenseMatrix<int,ScalarType> VAV_restart( Teuchos::View, *state.
VAV, restartDim, restartDim );
546 err = VAV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
547 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
549 if( d_problem->getM() != Teuchos::null )
552 RCP<const MV> BV_orig = MVT::CloneView( *state.
BV, allIndices );
553 RCP<MV> BV_restart = MVT::CloneViewNonConst( *state.
BV, restartIndices );
555 MVT::MvTimesMatAddMv(ST::one(),*BV_orig,Z_wanted,ST::zero(),*restartVecs);
556 MVT::SetBlock(*restartVecs,restartIndices,*BV_restart);
560 const Teuchos::SerialDenseMatrix<int,ScalarType> VBV_orig( Teuchos::View, *state.
VBV, state.
curDim, state.
curDim );
561 err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VBV_orig, Z_wanted, ST::zero() );
562 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
564 Teuchos::SerialDenseMatrix<int,ScalarType> VBV_restart( Teuchos::View, *state.
VBV, restartDim, restartDim );
565 VBV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
566 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
570 state.
Q->putScalar( ST::zero() );
571 state.
Z->putScalar( ST::zero() );
572 for(
int ii=0; ii<restartDim; ii++ )
574 (*state.
Q)(ii,ii)= ST::one();
575 (*state.
Z)(ii,ii)= ST::one();
579 state.
curDim = restartDim;
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::SortManager class.
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...
Implementation of a block Generalized Davidson eigensolver.
Basic implementation of the Anasazi::OrthoManager class.
Abstract class definition for Anasazi Output Managers.
Abstract class definition for Anasazi output stream.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Pure virtual base class which describes the basic interface for a solver manager.
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Types and exceptions used within Anasazi solvers and interfaces.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
This class defines the interface required by an eigensolver and status test class to compute solution...
Solver Manager for GeneralizedDavidson.
GeneralizedDavidsonSolMgr(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for GeneralizedDavidsonSolMgr.
int getNumIters() const
Get the iteration count for the most recent call to solve()
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Solves eigenvalue problem using generalized Davidson method.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
Traits class which defines basic operations on multivectors.
Output managers remove the need for the eigensolver to know any information about the required output...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ReturnType
Enumerated type used to pass back information from a solver manager.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Structure to contain pointers to GeneralizedDavidson state variables.
RCP< MV > AV
Image of V under A.
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.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VBV
Projection of B onto V.
RCP< MV > BV
Image of V under B.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Z
Right generalized Schur vectors from QZ decomposition of (VAV,VBV)
RCP< MV > V
Orthonormal basis for search subspace.
Output managers remove the need for the eigensolver to know any information about the required output...