10 #ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
11 #define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
59 template <
class ScalarType,
class MV,
class OP>
137 template <
class MagnitudeType,
class MV,
class OP>
138 class GeneralizedDavidsonSolMgr<std::complex<MagnitudeType>,MV,OP>
142 typedef std::complex<MagnitudeType> ScalarType;
144 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
148 MagnitudeType::this_class_is_missing_a_specialization();
159 template <
class ScalarType,
class MV,
class OP>
168 d_problem->getOperator() == Teuchos::null, std::invalid_argument,
"A operator not supplied on Eigenproblem." );
169 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getInitVec() == Teuchos::null, std::invalid_argument,
"No vector to clone from on Eigenproblem." );
170 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV() <= 0, std::invalid_argument,
"Number of requested eigenvalues must be positive.");
172 if( !pl.
isType<
int>(
"Block Size") )
174 pl.
set<
int>(
"Block Size",1);
177 if( !pl.
isType<
int>(
"Maximum Subspace Dimension") )
179 pl.
set<
int>(
"Maximum Subspace Dimension",3*problem->getNEV()*pl.
get<
int>(
"Block Size"));
182 if( !pl.
isType<
int>(
"Print Number of Ritz Values") )
184 int numToPrint = std::max( pl.
get<
int>(
"Block Size"), d_problem->getNEV() );
185 pl.
set<
int>(
"Print Number of Ritz Values",numToPrint);
189 MagnitudeType tol = pl.
get<MagnitudeType>(
"Convergence Tolerance",
MT::eps() );
191 std::invalid_argument,
"Convergence Tolerance must be greater than zero." );
194 if( pl.
isType<
int>(
"Maximum Restarts") )
196 d_maxRestarts = pl.
get<
int>(
"Maximum Restarts");
197 TEUCHOS_TEST_FOR_EXCEPTION( d_maxRestarts < 0, std::invalid_argument,
"Maximum Restarts must be non-negative" );
205 d_restartDim = pl.
get<
int>(
"Restart Dimension",d_problem->getNEV());
206 TEUCHOS_TEST_FOR_EXCEPTION( d_restartDim < d_problem->getNEV(),
207 std::invalid_argument,
"Restart Dimension must be at least NEV" );
210 std::string initType;
211 if( pl.
isType<std::string>(
"Initial Guess") )
213 initType = pl.
get<std::string>(
"Initial Guess");
214 TEUCHOS_TEST_FOR_EXCEPTION( initType!=
"User" && initType!=
"Random", std::invalid_argument,
215 "Initial Guess type must be 'User' or 'Random'." );
224 if( pl.
isType<std::string>(
"Which") )
226 which = pl.
get<std::string>(
"Which");
227 TEUCHOS_TEST_FOR_EXCEPTION( which!=
"LM" && which!=
"SM" && which!=
"LR" && which!=
"SR" && which!=
"LI" && which!=
"SI",
228 std::invalid_argument,
229 "Which must be one of LM,SM,LR,SR,LI,SI." );
240 std::string ortho = pl.
get<std::string>(
"Orthogonalization",
"SVQB");
241 TEUCHOS_TEST_FOR_EXCEPTION( ortho!=
"DGKS" && ortho!=
"SVQB" && ortho!=
"ICGS", std::invalid_argument,
242 "Anasazi::GeneralizedDavidsonSolMgr::constructor: Invalid orthogonalization type" );
248 else if( ortho==
"SVQB" )
252 else if( ortho==
"ICGS" )
258 bool scaleRes =
false;
259 bool failOnNaN =
false;
262 RES_2NORM,scaleRes,failOnNaN) );
269 int osProc = pl.
get(
"Output Processor", 0);
275 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
284 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
285 verbosity = pl.
get(
"Verbosity", verbosity);
287 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
294 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Building solver" << std::endl;
297 TEUCHOS_TEST_FOR_EXCEPTION(d_solver->getMaxSubspaceDim() < d_restartDim, std::invalid_argument,
298 "The maximum size of the subspace dimension (" << d_solver->getMaxSubspaceDim() <<
") must "
299 "not be smaller than the size of the restart space (" << d_restartDim <<
"). "
300 "Please adjust \"Restart Dimension\" and/or \"Maximum Subspace Dimension\" parameters.");
307 template <
class ScalarType,
class MV,
class OP>
312 d_problem->setSolution(sol);
314 d_solver->initialize();
322 if( d_tester->getStatus() ==
Passed )
326 if( restarts == d_maxRestarts )
330 d_solver->sortProblem( d_restartDim );
332 getRestartState( state );
333 d_solver->initialize( state );
339 d_solver->currentStatus(d_outputMan->stream(
FinalSummary));
342 sol.
numVecs = d_tester->howMany();
345 std::vector<int> whichVecs = d_tester->whichVecs();
346 std::vector<int> origIndex = d_solver->getRitzIndex();
350 for(
int i=0; i<sol.
numVecs; ++i )
352 if( origIndex[ whichVecs[i] ] == 1 )
354 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]+1 ) == whichVecs.end() )
356 whichVecs.push_back( whichVecs[i]+1 );
360 else if( origIndex[ whichVecs[i] ] == -1 )
362 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]-1 ) == whichVecs.end() )
364 whichVecs.push_back( whichVecs[i]-1 );
370 if( d_outputMan->isVerbosity(
Debug) )
372 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: "
373 << sol.
numVecs <<
" eigenpairs converged" << std::endl;
377 std::vector< Value<ScalarType> > origVals = d_solver->getRitzValues();
378 std::vector<MagnitudeType> realParts;
379 std::vector<MagnitudeType> imagParts;
380 for(
int i=0; i<sol.
numVecs; ++i )
382 realParts.push_back( origVals[whichVecs[i]].realpart );
383 imagParts.push_back( origVals[whichVecs[i]].imagpart );
386 std::vector<int> permVec(sol.
numVecs);
387 d_sortMan->sort( realParts, imagParts, Teuchos::rcpFromRef(permVec), sol.
numVecs );
390 std::vector<int> newWhich;
391 for(
int i=0; i<sol.
numVecs; ++i )
392 newWhich.push_back( whichVecs[permVec[i]] );
396 for(
int i=0; i<sol.
numVecs; ++i )
408 sol.
index = origIndex;
410 sol.
Evals = d_solver->getRitzValues();
420 for(
int i=0; i<sol.
numVecs; ++i )
422 sol.
index[i] = origIndex[ newWhich[i] ];
423 sol.
Evals[i] = origVals[ newWhich[i] ];
426 sol.
Evecs = MVT::CloneCopy( *(d_solver->getRitzVectors()), newWhich );
428 d_problem->setSolution(sol);
431 if( sol.
numVecs < d_problem->getNEV() )
440 template <
class ScalarType,
class MV,
class OP>
445 "Anasazi::GeneralizedDavidsonSolMgr: State dimension at restart is smaller than Restart Dimension" );
447 std::vector<int> ritzIndex = d_solver->getRitzIndex();
450 int restartDim = d_restartDim;
451 if( ritzIndex[d_restartDim-1]==1 )
454 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Restarting with "
455 << restartDim <<
" vectors" << std::endl;
468 std::vector<int> allIndices(state.
curDim);
469 for(
int i=0; i<state.
curDim; ++i )
475 std::vector<int> restartIndices(restartDim);
476 for(
int i=0; i<restartDim; ++i )
477 restartIndices[i] = i;
480 RCP<MV> V_restart = MVT::CloneViewNonConst( *state.
V, restartIndices );
483 RCP<MV> restartVecs = MVT::Clone(*state.
V,restartDim);
486 MVT::MvTimesMatAddMv(ST::one(),*V_orig,Z_wanted,ST::zero(),*restartVecs);
487 MVT::SetBlock(*restartVecs,restartIndices,*V_restart);
490 if( d_outputMan->isVerbosity(
Debug) )
492 MagnitudeType orthErr = d_orthoMan->orthonormError(*V_restart);
493 std::stringstream os;
494 os <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Error in V^T V == I after restart : " << orthErr << std::endl;
495 d_outputMan->print(
Debug,os.str());
499 RCP<MV> AV_restart = MVT::CloneViewNonConst( *state.
AV, restartIndices );
502 MVT::MvTimesMatAddMv(ST::one(),*AV_orig,Z_wanted,ST::zero(),*restartVecs);
503 MVT::SetBlock(*restartVecs,restartIndices,*AV_restart);
511 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
515 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
517 if( d_problem->getM() != Teuchos::null )
521 RCP<MV> BV_restart = MVT::CloneViewNonConst( *state.
BV, restartIndices );
523 MVT::MvTimesMatAddMv(ST::one(),*BV_orig,Z_wanted,ST::zero(),*restartVecs);
524 MVT::SetBlock(*restartVecs,restartIndices,*BV_restart);
530 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
534 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
538 state.
Q->putScalar( ST::zero() );
539 state.
Z->putScalar( ST::zero() );
540 for(
int ii=0; ii<restartDim; ii++ )
542 (*state.
Q)(ii,ii)= ST::one();
543 (*state.
Z)(ii,ii)= ST::one();
547 state.
curDim = restartDim;
552 #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)
#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.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Basic implementation of the Anasazi::SortManager class.
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...
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.
bool isParameter(const std::string &name) const
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...
int getNumIters() const
Get the iteration count for the most recent call to solve()
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
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.
bool isType(const std::string &name) const
Abstract class definition for Anasazi output stream.
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.