30 #ifndef ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP 31 #define ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP 54 #include "Teuchos_TimeMonitor.hpp" 89 template<
class ScalarType,
class MV,
class OP>
94 typedef Teuchos::ScalarTraits<ScalarType> SCT;
95 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
96 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
114 Teuchos::ParameterList &pl );
148 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
159 template<
class ScalarType,
class MV,
class OP>
162 Teuchos::ParameterList &pl ) :
171 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
172 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
173 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument,
"Problem not symmetric.");
174 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
176 whch_ = pl.get(
"Which",
"SR");
177 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR",
179 "SimpleLOBPCGSolMgr: \"Which\" parameter must be SM, LM, SR or LR.");
181 tol_ = pl.get(
"Convergence Tolerance",tol_);
182 TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
184 "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
187 if (pl.isParameter(
"Verbosity")) {
188 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
189 verb_ = pl.get(
"Verbosity", verb_);
191 verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
196 blockSize_= pl.get(
"Block Size",problem_->getNEV());
197 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
199 "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
201 maxIters_ = pl.get(
"Maximum Iterations",maxIters_);
207 template<
class ScalarType,
class MV,
class OP>
216 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > max;
223 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > norm
225 Teuchos::Array< Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
226 alltests.push_back(norm);
227 if (max != Teuchos::null) alltests.push_back(max);
228 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combo
233 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
236 Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho
239 Teuchos::ParameterList plist;
240 plist.set(
"Block Size",blockSize_);
241 plist.set(
"Full Ortho",
true);
244 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
247 if (problem_->getAuxVecs() != Teuchos::null) {
248 lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
252 int nev = problem_->getNEV();
253 Teuchos::Array< Teuchos::RCP<MV> > foundvecs;
254 Teuchos::Array< Teuchos::RCP< std::vector<MagnitudeType> > > foundvals;
255 while (numfound < nev) {
257 if (nev - numfound < blockSize_) {
258 norm->setQuorum(nev-numfound);
263 lobpcg_solver->iterate();
265 catch (
const std::exception &e) {
267 printer->stream(
Anasazi::Errors) <<
"Exception: " << e.what() << std::endl;
270 problem_->setSolution(sol);
275 if (norm->getStatus() ==
Passed) {
277 int num = norm->howMany();
279 TEUCHOS_TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
281 "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
282 std::vector<int> ind = norm->whichVecs();
284 if (num + numfound > nev) {
285 num = nev - numfound;
290 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
292 foundvecs.push_back(newvecs);
294 Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
295 auxvecs.push_back(newvecs);
297 lobpcg_solver->setAuxVecs(auxvecs);
300 Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp(
new std::vector<MagnitudeType>(num) );
301 std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
302 for (
int i=0; i<num; i++) {
303 (*newvals)[i] = all[ind[i]].realpart;
305 foundvals.push_back(newvals);
309 else if (max != Teuchos::null && max->getStatus() ==
Passed) {
311 int num = norm->howMany();
312 std::vector<int> ind = norm->whichVecs();
316 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
318 foundvecs.push_back(newvecs);
322 Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp(
new std::vector<MagnitudeType>(num) );
323 std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
324 for (
int i=0; i<num; i++) {
325 (*newvals)[i] = all[ind[i]].realpart;
327 foundvals.push_back(newvals);
334 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
338 TEUCHOS_TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
345 sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
348 sol.Evecs = Teuchos::null;
350 sol.Espace = sol.Evecs;
352 std::vector<MagnitudeType> vals(numfound);
353 sol.Evals.resize(numfound);
355 sol.index.resize(numfound,0);
358 for (
unsigned int i=0; i<foundvals.size(); i++) {
359 TEUCHOS_TEST_FOR_EXCEPTION((
signed int)(foundvals[i]->size()) != MVT::GetNumberVecs(*foundvecs[i]), std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
360 unsigned int lclnum = foundvals[i]->size();
361 std::vector<int> lclind(lclnum);
362 for (
unsigned int j=0; j<lclnum; j++) lclind[j] = curttl+j;
364 MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
366 std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
370 TEUCHOS_TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
374 std::vector<int> order(sol.numVecs);
375 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
377 for (
int i=0; i<sol.numVecs; i++) {
378 sol.Evals[i].realpart = vals[i];
379 sol.Evals[i].imagpart = MT::zero();
387 lobpcg_solver->currentStatus(printer->stream(
FinalSummary));
390 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 392 Teuchos::TimeMonitor::summarize( printer->stream(
TimingDetails ) );
397 problem_->setSolution(sol);
398 printer->stream(
Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
401 numIters_ = lobpcg_solver->getNumIters();
Pure virtual base class which describes the basic interface for a solver manager. ...
A special StatusTest for printing other status tests.
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...
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration...
Status test for forming logical combinations of other status tests.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Basic implementation of the Anasazi::SortManager class.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
An exception class parent to all Anasazi exceptions.
Implementation of the locally-optimal block preconditioned conjugate gradient (LOBPCG) method...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi's templated, static class providing utilities for the solvers.
virtual ~SimpleLOBPCGSolMgr()
Destructor.
int numVecs
The number of computed eigenpairs.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
int getNumIters() const
Get the iteration count for the most recent call to solve().
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.
Traits class which defines basic operations on multivectors.
The Anasazi::SimpleLOBPCGSolMgr provides a simple solver manager over the LOBPCG eigensolver.
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.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
A status test for testing the number of iterations.
Status test for testing the number of iterations.
static void permuteVectors(const int n, const std::vector< int > &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *resids=0)
Permute the vectors in a multivector according to the permutation vector perm, and optionally the res...
Special StatusTest for printing status tests.
Status test for forming logical combinations of other status tests.
Types and exceptions used within Anasazi solvers and interfaces.
A status test for testing the norm of the eigenvectors residuals.
SimpleLOBPCGSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for SimpleLOBPCGSolMgr.
Class which provides internal utilities for the Anasazi solvers.