42 #ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 43 #define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 49 #include "Teuchos_ParameterList.hpp" 50 #include "Teuchos_RCPDecl.hpp" 91 template <
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;
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,
177 Teuchos::ParameterList &pl )
180 MagnitudeType::this_class_is_missing_a_specialization();
191 template <
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.");
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>
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;
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...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
Solves eigenvalue problem using generalized Davidson method.
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.
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...
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.
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()