Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AnasaziSimpleLOBPCGSolMgr.hpp
Go to the documentation of this file.
1 
2 // @HEADER
3 // *****************************************************************************
4 // Anasazi: Block Eigensolvers Package
5 //
6 // Copyright 2004 NTESS and the Anasazi contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
11 #ifndef ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
12 #define ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
13 
19 #include "AnasaziConfigDefs.hpp"
20 #include "AnasaziTypes.hpp"
21 
22 #include "AnasaziEigenproblem.hpp"
23 #include "AnasaziSolverManager.hpp"
24 
25 #include "AnasaziLOBPCG.hpp"
26 #include "AnasaziBasicSort.hpp"
32 #include "AnasaziOutputManager.hpp"
34 #include "AnasaziSolverUtils.hpp"
35 
36 #include "Teuchos_TimeMonitor.hpp"
37 #include "Teuchos_FancyOStream.hpp"
38 
46 
70 namespace Anasazi {
71 
72 template<class ScalarType, class MV, class OP>
73 class SimpleLOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
74 
75  private:
78  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
80 
81  public:
82 
84 
85 
101 
103  virtual ~SimpleLOBPCGSolMgr() {};
105 
107 
108 
110  return *problem_;
111  }
112 
113  int getNumIters() const {
114  return numIters_;
115  }
116 
118 
120 
121 
130  ReturnType solve();
132 
133  private:
136  std::string whch_;
137  MagnitudeType tol_;
138  int osProc_;
139  int verb_;
140  int blockSize_;
141  int maxIters_;
142  int numIters_;
143 };
144 
145 
147 template<class ScalarType, class MV, class OP>
150  Teuchos::ParameterList &pl ) :
151  problem_(problem),
152  whch_("LM"),
153  tol_(1e-6),
154  osProc_(0),
155  verb_(Anasazi::Errors),
156  blockSize_(0),
157  maxIters_(100),
158  numIters_(0)
159 {
160  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
161  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
162  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
163  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
164 
165  whch_ = pl.get("Which","SR");
166  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
167  AnasaziError,
168  "SimpleLOBPCGSolMgr: \"Which\" parameter must be SM, LM, SR or LR.");
169 
170  tol_ = pl.get("Convergence Tolerance",tol_);
171  TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
172  AnasaziError,
173  "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
174 
175  // Create a formatted output stream to print to.
176  // See if user requests output processor.
177  osProc_ = pl.get("Output Processor", osProc_);
178 
179  // If not passed in by user, it will be chosen based upon operator type.
180  if (pl.isParameter("Output Stream")) {
181  osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
182  }
183  else {
184  osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc_);
185  }
186 
187  // verbosity level
188  if (pl.isParameter("Verbosity")) {
189  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
190  verb_ = pl.get("Verbosity", verb_);
191  } else {
192  verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
193  }
194  }
195 
196 
197  blockSize_= pl.get("Block Size",problem_->getNEV());
198  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
199  AnasaziError,
200  "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
201 
202  maxIters_ = pl.get("Maximum Iterations",maxIters_);
203 }
204 
205 
206 
208 template<class ScalarType, class MV, class OP>
211 
212  // sort manager
214  // output manager
216  // status tests
218  if (maxIters_ > 0) {
219  max = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
220  }
221  else {
222  max = Teuchos::null;
223  }
227  alltests.push_back(norm);
228  if (max != Teuchos::null) alltests.push_back(max);
232  ));
233  // printing StatusTest
235  = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combo,1,Passed ) );
236  // orthomanager
238  = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
239  // parameter list
241  plist.set("Block Size",blockSize_);
242  plist.set("Full Ortho",true);
243 
244  // create an LOBPCG solver
246  = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
247  // add the auxillary vecs from the eigenproblem to the solver
248  if (problem_->getAuxVecs() != Teuchos::null) {
249  lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
250  }
251 
252  int numfound = 0;
253  int nev = problem_->getNEV();
256  while (numfound < nev) {
257  // reduce the strain on norm test, if we are almost done
258  if (nev - numfound < blockSize_) {
259  norm->setQuorum(nev-numfound);
260  }
261 
262  // tell the solver to iterate
263  try {
264  lobpcg_solver->iterate();
265  }
266  catch (const std::exception &e) {
267  // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
268  printer->stream(Anasazi::Errors) << "Exception: " << e.what() << std::endl;
270  sol.numVecs = 0;
271  problem_->setSolution(sol);
272  throw;
273  }
274 
275  // check the status tests
276  if (norm->getStatus() == Passed) {
277 
278  int num = norm->howMany();
279  // if num < blockSize_, it is because we are on the last iteration: num+numfound>=nev
280  TEUCHOS_TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
281  std::logic_error,
282  "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
283  std::vector<int> ind = norm->whichVecs();
284  // just grab the ones that we need
285  if (num + numfound > nev) {
286  num = nev - numfound;
287  ind.resize(num);
288  }
289 
290  // copy the converged eigenvectors
291  Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
292  // store them
293  foundvecs.push_back(newvecs);
294  // add them as auxiliary vectors
295  Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
296  auxvecs.push_back(newvecs);
297  // setAuxVecs() will reset the solver to uninitialized, without messing with numIters()
298  lobpcg_solver->setAuxVecs(auxvecs);
299 
300  // copy the converged eigenvalues
301  Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp( new std::vector<MagnitudeType>(num) );
302  std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
303  for (int i=0; i<num; i++) {
304  (*newvals)[i] = all[ind[i]].realpart;
305  }
306  foundvals.push_back(newvals);
307 
308  numfound += num;
309  }
310  else if (max != Teuchos::null && max->getStatus() == Passed) {
311 
312  int num = norm->howMany();
313  std::vector<int> ind = norm->whichVecs();
314 
315  if (num > 0) {
316  // copy the converged eigenvectors
317  Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
318  // store them
319  foundvecs.push_back(newvecs);
320  // don't bother adding them as auxiliary vectors; we have reached maxiters and are going to quit
321 
322  // copy the converged eigenvalues
323  Teuchos::RCP<std::vector<MagnitudeType> > newvals = Teuchos::rcp( new std::vector<MagnitudeType>(num) );
324  std::vector<Value<ScalarType> > all = lobpcg_solver->getRitzValues();
325  for (int i=0; i<num; i++) {
326  (*newvals)[i] = all[ind[i]].realpart;
327  }
328  foundvals.push_back(newvals);
329 
330  numfound += num;
331  }
332  break; // while(numfound < nev)
333  }
334  else {
335  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
336  }
337  } // end of while(numfound < nev)
338 
339  TEUCHOS_TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
340 
341  // create contiguous storage for all eigenvectors, eigenvalues
343  sol.numVecs = numfound;
344  if (numfound > 0) {
345  // allocate space for eigenvectors
346  sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
347  }
348  else {
349  sol.Evecs = Teuchos::null;
350  }
351  sol.Espace = sol.Evecs;
352  // allocate space for eigenvalues
353  std::vector<MagnitudeType> vals(numfound);
354  sol.Evals.resize(numfound);
355  // all real eigenvalues: set index vectors [0,...,numfound-1]
356  sol.index.resize(numfound,0);
357  // store eigenvectors, eigenvalues
358  int curttl = 0;
359  for (unsigned int i=0; i<foundvals.size(); i++) {
360  TEUCHOS_TEST_FOR_EXCEPTION((signed int)(foundvals[i]->size()) != MVT::GetNumberVecs(*foundvecs[i]), std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
361  unsigned int lclnum = foundvals[i]->size();
362  std::vector<int> lclind(lclnum);
363  for (unsigned int j=0; j<lclnum; j++) lclind[j] = curttl+j;
364  // put the eigenvectors
365  MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
366  // put the eigenvalues
367  std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
368 
369  curttl += lclnum;
370  }
371  TEUCHOS_TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
372 
373  // sort the eigenvalues and permute the eigenvectors appropriately
374  if (numfound > 0) {
375  std::vector<int> order(sol.numVecs);
376  sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
377  // store the values in the Eigensolution
378  for (int i=0; i<sol.numVecs; i++) {
379  sol.Evals[i].realpart = vals[i];
380  sol.Evals[i].imagpart = MT::zero();
381  }
382  // now permute the eigenvectors according to order
384  msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
385  }
386 
387  // print final summary
388  lobpcg_solver->currentStatus(printer->stream(FinalSummary));
389 
390  // print timing information
391 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
392  if ( printer->isVerbosity( TimingDetails ) ) {
393  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
394  }
395 #endif
396 
397  // send the solution to the eigenproblem
398  problem_->setSolution(sol);
399  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
400 
401  // get the number of iterations performed for this solve.
402  numIters_ = lobpcg_solver->getNumIters();
403 
404  // return from SolMgr::solve()
405  if (sol.numVecs < nev) return Unconverged;
406  return Converged;
407 }
408 
409 
410 
411 
412 } // end Anasazi namespace
413 
414 #endif /* ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP */
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
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...
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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...
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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...
Anasazi&#39;s templated, static class providing utilities for the solvers.
bool isParameter(const std::string &name) const
int numVecs
The number of computed eigenpairs.
Abstract class definition for Anasazi Output Managers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Abstract base class which defines the interface required by an eigensolver and status test class to c...
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
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 &quot;A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver&#39;s iterate() routine until ...
A status test for testing the number of iterations.
Status test for testing the number of iterations.
void push_back(const value_type &x)
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.
size_type size() const
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.
iterator begin()
int getNumIters() const
Get the iteration count for the most recent call to solve().
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.