43#ifndef ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
44#define ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
68#include "Teuchos_TimeMonitor.hpp"
69#include "Teuchos_FancyOStream.hpp"
104template<
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_;
179template<
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_);
240template<
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;
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();
Basic implementation of the Anasazi::SortManager class.
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Implementation of the locally-optimal block preconditioned conjugate gradient (LOBPCG) method.
Abstract class definition for Anasazi Output Managers.
Abstract class definition for Anasazi output stream.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Pure virtual base class which describes the basic interface for a solver manager.
Class which provides internal utilities for the Anasazi solvers.
Status test for forming logical combinations of other status tests.
Status test for testing the number of iterations.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals.
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
This class defines the interface required by an eigensolver and status test class to compute solution...
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration,...
Traits class which defines basic operations on multivectors.
Output managers remove the need for the eigensolver to know any information about the required output...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
The Anasazi::SimpleLOBPCGSolMgr provides a simple solver manager over the LOBPCG eigensolver.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
SimpleLOBPCGSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for SimpleLOBPCGSolMgr.
virtual ~SimpleLOBPCGSolMgr()
Destructor.
int getNumIters() const
Get the iteration count for the most recent call to solve().
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Anasazi's templated, static class providing utilities for the solvers.
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...
Status test for forming logical combinations of other status tests.
A status test for testing the number of iterations.
A special StatusTest for printing other status tests.
A status test for testing the norm of the eigenvectors residuals.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ReturnType
Enumerated type used to pass back information from a solver manager.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
int numVecs
The number of computed eigenpairs.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Output managers remove the need for the eigensolver to know any information about the required output...