43 #ifndef ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP 44 #define ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP 68 #include "Teuchos_TimeMonitor.hpp" 69 #include "Teuchos_FancyOStream.hpp" 104 template<
class ScalarType,
class MV,
class OP>
109 typedef Teuchos::ScalarTraits<ScalarType> SCT;
110 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
111 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
132 Teuchos::ParameterList &pl );
166 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
167 Teuchos::RCP<Teuchos::FancyOStream> osp_;
179 template<
class ScalarType,
class MV,
class OP>
182 Teuchos::ParameterList &pl ) :
192 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
193 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument,
"Problem not set.");
194 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument,
"Problem not symmetric.");
195 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument,
"Problem does not contain initial vectors to clone from.");
197 whch_ = pl.get(
"Which",
"SR");
198 TEUCHOS_TEST_FOR_EXCEPTION(whch_ !=
"SM" && whch_ !=
"LM" && whch_ !=
"SR" && whch_ !=
"LR",
200 "SimpleLOBPCGSolMgr: \"Which\" parameter must be SM, LM, SR or LR.");
202 tol_ = pl.get(
"Convergence Tolerance",tol_);
203 TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
205 "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
209 osProc_ = pl.get(
"Output Processor", osProc_);
212 if (pl.isParameter(
"Output Stream")) {
213 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,
"Output Stream");
220 if (pl.isParameter(
"Verbosity")) {
221 if (Teuchos::isParameterType<int>(pl,
"Verbosity")) {
222 verb_ = pl.get(
"Verbosity", verb_);
224 verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,
"Verbosity");
229 blockSize_= pl.get(
"Block Size",problem_->getNEV());
230 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
232 "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
234 maxIters_ = pl.get(
"Maximum Iterations",maxIters_);
240 template<
class ScalarType,
class MV,
class OP>
249 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > max;
256 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > norm
258 Teuchos::Array< Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
259 alltests.push_back(norm);
260 if (max != Teuchos::null) alltests.push_back(max);
261 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combo
266 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
269 Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho
272 Teuchos::ParameterList plist;
273 plist.set(
"Block Size",blockSize_);
274 plist.set(
"Full Ortho",
true);
277 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
280 if (problem_->getAuxVecs() != Teuchos::null) {
281 lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
285 int nev = problem_->getNEV();
286 Teuchos::Array< Teuchos::RCP<MV> > foundvecs;
287 Teuchos::Array< Teuchos::RCP< std::vector<MagnitudeType> > > foundvals;
288 while (numfound < nev) {
290 if (nev - numfound < blockSize_) {
291 norm->setQuorum(nev-numfound);
296 lobpcg_solver->iterate();
298 catch (
const std::exception &e) {
300 printer->stream(
Anasazi::Errors) <<
"Exception: " << e.what() << std::endl;
303 problem_->setSolution(sol);
308 if (norm->getStatus() ==
Passed) {
310 int num = norm->howMany();
312 TEUCHOS_TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
314 "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
315 std::vector<int> ind = norm->whichVecs();
317 if (num + numfound > nev) {
318 num = nev - numfound;
323 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
325 foundvecs.push_back(newvecs);
327 Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
328 auxvecs.push_back(newvecs);
330 lobpcg_solver->setAuxVecs(auxvecs);
333 Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp(
new std::vector<MagnitudeType>(num) );
334 std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
335 for (
int i=0; i<num; i++) {
336 (*newvals)[i] = all[ind[i]].realpart;
338 foundvals.push_back(newvals);
342 else if (max != Teuchos::null && max->getStatus() ==
Passed) {
344 int num = norm->howMany();
345 std::vector<int> ind = norm->whichVecs();
349 Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
351 foundvecs.push_back(newvecs);
355 Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp(
new std::vector<MagnitudeType>(num) );
356 std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
357 for (
int i=0; i<num; i++) {
358 (*newvals)[i] = all[ind[i]].realpart;
360 foundvals.push_back(newvals);
367 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
371 TEUCHOS_TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
378 sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
381 sol.Evecs = Teuchos::null;
383 sol.Espace = sol.Evecs;
385 std::vector<MagnitudeType> vals(numfound);
386 sol.Evals.resize(numfound);
388 sol.index.resize(numfound,0);
391 for (
unsigned int i=0; i<foundvals.size(); i++) {
392 TEUCHOS_TEST_FOR_EXCEPTION((
signed int)(foundvals[i]->size()) != MVT::GetNumberVecs(*foundvecs[i]), std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
393 unsigned int lclnum = foundvals[i]->size();
394 std::vector<int> lclind(lclnum);
395 for (
unsigned int j=0; j<lclnum; j++) lclind[j] = curttl+j;
397 MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
399 std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
403 TEUCHOS_TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error,
"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
407 std::vector<int> order(sol.numVecs);
408 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
410 for (
int i=0; i<sol.numVecs; i++) {
411 sol.Evals[i].realpart = vals[i];
412 sol.Evals[i].imagpart = MT::zero();
420 lobpcg_solver->currentStatus(printer->stream(
FinalSummary));
423 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 425 Teuchos::TimeMonitor::summarize( printer->stream(
TimingDetails ) );
430 problem_->setSolution(sol);
431 printer->stream(
Debug) <<
"Returning " << sol.numVecs <<
" eigenpairs to eigenproblem." << std::endl;
434 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.
Output managers remove the need for the eigensolver to know any information about the required output...
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
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.
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.
Abstract class definition for Anasazi output stream.
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.