42 #ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 43 #define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 91 template <
class ScalarType,
class MV,
class OP>
169 template <
class MagnitudeType,
class MV,
class OP>
170 class GeneralizedDavidsonSolMgr<std::complex<MagnitudeType>,MV,OP>
174 typedef std::complex<MagnitudeType> ScalarType;
176 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
180 MagnitudeType::this_class_is_missing_a_specialization();
191 template <
class ScalarType,
class MV,
class OP>
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() );
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");
237 d_restartDim = pl.
get<
int>(
"Restart Dimension",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");
247 "Initial Guess type must be 'User' or 'Random'." );
256 if( pl.
isType<std::string>(
"Which") )
258 which = pl.
get<std::string>(
"Which");
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");
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;
294 RES_2NORM,scaleRes,failOnNaN) );
301 int osProc = pl.
get(
"Output Processor", 0);
307 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
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;
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.");
339 template <
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() )
472 template <
class ScalarType,
class MV,
class OP>
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;
500 std::vector<int> allIndices(state.
curDim);
501 for(
int i=0; i<state.
curDim; ++i )
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 );
534 MVT::MvTimesMatAddMv(ST::one(),*AV_orig,Z_wanted,ST::zero(),*restartVecs);
535 MVT::SetBlock(*restartVecs,restartIndices,*AV_restart);
543 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
547 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
549 if( d_problem->getM() != Teuchos::null )
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);
562 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
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;
584 #endif // ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Pure virtual base class which describes the basic interface for a solver manager. ...
RCP< MV > V
Orthonormal basis for search subspace.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
RCP< MV > AV
Image of V under A.
This class defines the interface required by an eigensolver and status test class to compute solution...
static magnitudeType eps()
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
Solves eigenvalue problem using generalized Davidson method.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VAV
Projection of A onto V.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
int curDim
The current subspace dimension.
Basic implementation of the Anasazi::SortManager class.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
Output managers remove the need for the eigensolver to know any information about the required output...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
RCP< MV > BV
Image of V under B.
Structure to contain pointers to GeneralizedDavidson state variables.
int numVecs
The number of computed eigenpairs.
Abstract class definition for Anasazi Output Managers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Abstract base class which defines the interface required by an eigensolver and status test class to c...
ReturnType
Enumerated type used to pass back information from a solver manager.
Output managers remove the need for the eigensolver to know any information about the required output...
A status test for testing the norm of the eigenvectors residuals.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Q
Left generalized Schur vectors from QZ decomposition of (VAV,VBV)
Traits class which defines basic operations on multivectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
bool isType(const std::string &name) const
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
GeneralizedDavidsonSolMgr(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for GeneralizedDavidsonSolMgr.
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)
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
Abstract class definition for Anasazi output stream.
bool isParameter(const std::string &name) const
Solver Manager for GeneralizedDavidson.
A status test for testing the norm of the eigenvectors residuals.
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::OrthoManager class.
Implementation of a block Generalized Davidson eigensolver.
int getNumIters() const
Get the iteration count for the most recent call to solve()