29 #ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 30 #define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP 36 #include "Teuchos_ParameterList.hpp" 37 #include "Teuchos_RCPDecl.hpp" 77 template <
class ScalarType,
class MV,
class OP>
111 Teuchos::ParameterList &pl );
134 typedef Teuchos::ScalarTraits<ScalarType> ST;
135 typedef typename ST::magnitudeType MagnitudeType;
136 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
138 RCP< Eigenproblem<ScalarType,MV,OP> > d_problem;
139 RCP< GeneralizedDavidson<ScalarType,MV,OP> > d_solver;
140 RCP< OutputManager<ScalarType> > d_outputMan;
141 RCP< OrthoManager<ScalarType,MV> > d_orthoMan;
142 RCP< SortManager<MagnitudeType> > d_sortMan;
143 RCP< StatusTest<ScalarType,MV,OP> > d_tester;
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,
160 Teuchos::ParameterList &pl )
163 MagnitudeType::this_class_is_missing_a_specialization();
174 template <
class ScalarType,
class MV,
class OP>
177 Teuchos::ParameterList &pl )
180 TEUCHOS_TEST_FOR_EXCEPTION( d_problem == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager." );
181 TEUCHOS_TEST_FOR_EXCEPTION( !d_problem->isProblemSet(), std::invalid_argument,
"Problem not set." );
182 TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getA() == Teuchos::null &&
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() );
205 TEUCHOS_TEST_FOR_EXCEPTION( pl.get<MagnitudeType>(
"Convergence Tolerance") <= MT::zero(),
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");
212 TEUCHOS_TEST_FOR_EXCEPTION( d_maxRestarts < 0, std::invalid_argument,
"Maximum Restarts must be non-negative" );
220 d_restartDim = pl.get<
int>(
"Restart Dimension",d_problem->getNEV());
221 TEUCHOS_TEST_FOR_EXCEPTION( d_restartDim < 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");
229 TEUCHOS_TEST_FOR_EXCEPTION( initType!=
"User" && initType!=
"Random", std::invalid_argument,
230 "Initial Guess type must be 'User' or 'Random'." );
239 if( pl.isType<std::string>(
"Which") )
241 which = pl.get<std::string>(
"Which");
242 TEUCHOS_TEST_FOR_EXCEPTION( which!=
"LM" && which!=
"SM" && which!=
"LR" && which!=
"SR" && which!=
"LI" && which!=
"SI",
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");
256 TEUCHOS_TEST_FOR_EXCEPTION( ortho!=
"DGKS" && ortho!=
"SVQB" && ortho!=
"ICGS", std::invalid_argument,
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;
275 RCP<StatusTest<ScalarType,MV,OP> > resNormTest = Teuchos::rcp(
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>
430 TEUCHOS_TEST_FOR_EXCEPTION( state.
curDim <= d_restartDim, std::runtime_error,
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;
450 const Teuchos::SerialDenseMatrix<int,ScalarType> Z_wanted =
451 Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*state.
Z,state.
curDim,restartDim);
454 std::vector<int> allIndices(state.
curDim);
455 for(
int i=0; i<state.
curDim; ++i )
458 RCP<const MV> V_orig = MVT::CloneView( *state.
V, allIndices );
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 );
486 RCP<const MV> AV_orig = MVT::CloneView( *state.
AV, allIndices );
488 MVT::MvTimesMatAddMv(ST::one(),*AV_orig,Z_wanted,ST::zero(),*restartVecs);
489 MVT::SetBlock(*restartVecs,restartIndices,*AV_restart);
494 const Teuchos::SerialDenseMatrix<int,ScalarType> VAV_orig( Teuchos::View, *state.
VAV, state.
curDim, state.
curDim );
495 Teuchos::SerialDenseMatrix<int,ScalarType> tmpMat(state.
curDim, restartDim);
496 err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VAV_orig, Z_wanted, ST::zero() );
497 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
499 Teuchos::SerialDenseMatrix<int,ScalarType> VAV_restart( Teuchos::View, *state.
VAV, restartDim, restartDim );
500 err = VAV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
501 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
503 if( d_problem->getM() != Teuchos::null )
506 RCP<const MV> BV_orig = MVT::CloneView( *state.
BV, allIndices );
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);
514 const Teuchos::SerialDenseMatrix<int,ScalarType> VBV_orig( Teuchos::View, *state.
VBV, state.
curDim, state.
curDim );
515 err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VBV_orig, Z_wanted, ST::zero() );
516 TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error,
"GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
518 Teuchos::SerialDenseMatrix<int,ScalarType> VBV_restart( Teuchos::View, *state.
VBV, restartDim, restartDim );
519 VBV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
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...
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...
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...
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...
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()