29 #ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 30 #define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 77 template <
class ScalarType,
class MV,
class OP>
152 template <
class MagnitudeType,
class MV,
class OP>
153 class GeneralizedDavidsonSolMgr<std::complex<MagnitudeType>,MV,OP>
157 typedef std::complex<MagnitudeType> ScalarType;
159 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
163 MagnitudeType::this_class_is_missing_a_specialization();
174 template <
class ScalarType,
class MV,
class OP>
183 d_problem->getOperator() == Teuchos::null, std::invalid_argument,
"A operator not supplied on Eigenproblem." );
184 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getInitVec() == Teuchos::null, std::invalid_argument,
"No vector to clone from on Eigenproblem." );
185 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV() <= 0, std::invalid_argument,
"Number of requested eigenvalues must be positive.");
187 if( !pl.
isType<
int>(
"Block Size") )
189 pl.
set<
int>(
"Block Size",1);
192 if( !pl.
isType<
int>(
"Maximum Subspace Dimension") )
194 pl.
set<
int>(
"Maximum Subspace Dimension",3*problem->getNEV()*pl.
get<
int>(
"Block Size"));
197 if( !pl.
isType<
int>(
"Print Number of Ritz Values") )
199 int numToPrint = std::max( pl.
get<
int>(
"Block Size"), d_problem->getNEV() );
200 pl.
set<
int>(
"Print Number of Ritz Values",numToPrint);
204 MagnitudeType tol = pl.
get<MagnitudeType>(
"Convergence Tolerance",
MT::eps() );
206 std::invalid_argument,
"Convergence Tolerance must be greater than zero." );
209 if( pl.
isType<
int>(
"Maximum Restarts") )
211 d_maxRestarts = pl.
get<
int>(
"Maximum Restarts");
220 d_restartDim = pl.
get<
int>(
"Restart Dimension",d_problem->getNEV());
222 std::invalid_argument,
"Restart Dimension must be at least NEV" );
225 std::string initType;
226 if( pl.
isType<std::string>(
"Initial Guess") )
228 initType = pl.
get<std::string>(
"Initial Guess");
230 "Initial Guess type must be 'User' or 'Random'." );
239 if( pl.
isType<std::string>(
"Which") )
241 which = pl.
get<std::string>(
"Which");
243 std::invalid_argument,
244 "Which must be one of LM,SM,LR,SR,LI,SI." );
255 std::string ortho = pl.
get<std::string>(
"Orthogonalization",
"SVQB");
257 "Anasazi::GeneralizedDavidsonSolMgr::constructor: Invalid orthogonalization type" );
263 else if( ortho==
"SVQB" )
267 else if( ortho==
"ICGS" )
273 bool scaleRes =
false;
274 bool failOnNaN =
false;
277 RES_2NORM,scaleRes,failOnNaN) );
281 int verbosity = pl.
get<
int>(
"Verbosity",
Errors);
283 d_outputMan->setVerbosity( verbosity );
286 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Building solver" << std::endl;
293 template <
class ScalarType,
class MV,
class OP>
298 d_problem->setSolution(sol);
300 d_solver->initialize();
308 if( d_tester->getStatus() ==
Passed )
312 if( restarts == d_maxRestarts )
316 d_solver->sortProblem( d_restartDim );
318 getRestartState( state );
319 d_solver->initialize( state );
325 d_solver->currentStatus(d_outputMan->stream(
FinalSummary));
328 sol.
numVecs = d_tester->howMany();
331 std::vector<int> whichVecs = d_tester->whichVecs();
332 std::vector<int> origIndex = d_solver->getRitzIndex();
336 for(
int i=0; i<sol.
numVecs; ++i )
338 if( origIndex[ whichVecs[i] ] == 1 )
340 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]+1 ) == whichVecs.end() )
342 whichVecs.push_back( whichVecs[i]+1 );
346 else if( origIndex[ whichVecs[i] ] == -1 )
348 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]-1 ) == whichVecs.end() )
350 whichVecs.push_back( whichVecs[i]-1 );
356 if( d_outputMan->isVerbosity(
Debug) )
358 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: " 359 << sol.
numVecs <<
" eigenpairs converged" << std::endl;
363 std::vector< Value<ScalarType> > origVals = d_solver->getRitzValues();
364 std::vector<MagnitudeType> realParts;
365 std::vector<MagnitudeType> imagParts;
366 for(
int i=0; i<sol.
numVecs; ++i )
368 realParts.push_back( origVals[whichVecs[i]].realpart );
369 imagParts.push_back( origVals[whichVecs[i]].imagpart );
372 std::vector<int> permVec(sol.
numVecs);
373 d_sortMan->sort( realParts, imagParts, Teuchos::rcpFromRef(permVec), sol.
numVecs );
376 std::vector<int> newWhich;
377 for(
int i=0; i<sol.
numVecs; ++i )
378 newWhich.push_back( whichVecs[permVec[i]] );
382 for(
int i=0; i<sol.
numVecs; ++i )
394 sol.
index = origIndex;
396 sol.
Evals = d_solver->getRitzValues();
406 for(
int i=0; i<sol.
numVecs; ++i )
408 sol.
index[i] = origIndex[ newWhich[i] ];
409 sol.
Evals[i] = origVals[ newWhich[i] ];
412 sol.
Evecs = MVT::CloneCopy( *(d_solver->getRitzVectors()), newWhich );
414 d_problem->setSolution(sol);
417 if( sol.
numVecs < d_problem->getNEV() )
426 template <
class ScalarType,
class MV,
class OP>
431 "Anasazi::GeneralizedDavidsonSolMgr: State dimension at restart is smaller than Restart Dimension" );
433 std::vector<int> ritzIndex = d_solver->getRitzIndex();
436 int restartDim = d_restartDim;
437 if( ritzIndex[d_restartDim-1]==1 )
440 d_outputMan->stream(
Debug) <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Restarting with " 441 << restartDim <<
" vectors" << std::endl;
454 std::vector<int> allIndices(state.
curDim);
455 for(
int i=0; i<state.
curDim; ++i )
461 std::vector<int> restartIndices(restartDim);
462 for(
int i=0; i<restartDim; ++i )
463 restartIndices[i] = i;
466 RCP<MV> V_restart = MVT::CloneViewNonConst( *state.
V, restartIndices );
469 RCP<MV> restartVecs = MVT::Clone(*state.
V,restartDim);
472 MVT::MvTimesMatAddMv(ST::one(),*V_orig,Z_wanted,ST::zero(),*restartVecs);
473 MVT::SetBlock(*restartVecs,restartIndices,*V_restart);
476 if( d_outputMan->isVerbosity(
Debug) )
478 MagnitudeType orthErr = d_orthoMan->orthonormError(*V_restart);
479 std::stringstream os;
480 os <<
" >> Anasazi::GeneralizedDavidsonSolMgr: Error in V^T V == I after restart : " << orthErr << std::endl;
481 d_outputMan->print(
Debug,os.str());
485 RCP<MV> AV_restart = MVT::CloneViewNonConst( *state.
AV, restartIndices );
488 MVT::MvTimesMatAddMv(ST::one(),*AV_orig,Z_wanted,ST::zero(),*restartVecs);
489 MVT::SetBlock(*restartVecs,restartIndices,*AV_restart);
497 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
501 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
503 if( d_problem->getM() != Teuchos::null )
507 RCP<MV> BV_restart = MVT::CloneViewNonConst( *state.
BV, restartIndices );
509 MVT::MvTimesMatAddMv(ST::one(),*BV_orig,Z_wanted,ST::zero(),*restartVecs);
510 MVT::SetBlock(*restartVecs,restartIndices,*BV_restart);
516 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
520 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
524 state.
Q->putScalar( ST::zero() );
525 state.
Z->putScalar( ST::zero() );
526 for(
int ii=0; ii<restartDim; ii++ )
528 (*state.
Q)(ii,ii)= ST::one();
529 (*state.
Z)(ii,ii)= ST::one();
533 state.
curDim = restartDim;
538 #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...
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.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Anasazi's basic output manager for sending information of select verbosity levels to the appropriate ...
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.
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.
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()