Anasazi  Version of the Day
AnasaziSimpleLOBPCGSolMgr.hpp
Go to the documentation of this file.
1 
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Anasazi: Block Eigensolvers Package
6 // Copyright (2004) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // This library is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Lesser General Public License as
13 // published by the Free Software Foundation; either version 2.1 of the
14 // License, or (at your option) any later version.
15 //
16 // This library is distributed in the hope that it will be useful, but
17 // WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 // Lesser General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License along with this library; if not, write to the Free Software
23 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
24 // USA
25 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
26 //
27 // ***********************************************************************
28 // @HEADER
29 
30 #ifndef ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
31 #define ANASAZI_SIMPLE_LOBPCG_SOLMGR_HPP
32 
38 #include "AnasaziConfigDefs.hpp"
39 #include "AnasaziTypes.hpp"
40 
41 #include "AnasaziEigenproblem.hpp"
42 #include "AnasaziSolverManager.hpp"
43 
44 #include "AnasaziLOBPCG.hpp"
45 #include "AnasaziBasicSort.hpp"
52 #include "AnasaziSolverUtils.hpp"
53 
54 #include "Teuchos_TimeMonitor.hpp"
55 
63 
87 namespace Anasazi {
88 
89 template<class ScalarType, class MV, class OP>
90 class SimpleLOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
91 
92  private:
94  typedef Teuchos::ScalarTraits<ScalarType> SCT;
95  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
96  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
97 
98  public:
99 
101 
102 
113  SimpleLOBPCGSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
114  Teuchos::ParameterList &pl );
115 
117  virtual ~SimpleLOBPCGSolMgr() {};
119 
121 
122 
124  return *problem_;
125  }
126 
127  int getNumIters() const {
128  return numIters_;
129  }
130 
132 
134 
135 
144  ReturnType solve();
146 
147  private:
148  Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
149  std::string whch_;
150  MagnitudeType tol_;
151  int verb_;
152  int blockSize_;
153  int maxIters_;
154  int numIters_;
155 };
156 
157 
159 template<class ScalarType, class MV, class OP>
161  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
162  Teuchos::ParameterList &pl ) :
163  problem_(problem),
164  whch_("LM"),
165  tol_(1e-6),
166  verb_(Anasazi::Errors),
167  blockSize_(0),
168  maxIters_(100),
169  numIters_(0)
170 {
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.");
175 
176  whch_ = pl.get("Which","SR");
177  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
178  AnasaziError,
179  "SimpleLOBPCGSolMgr: \"Which\" parameter must be SM, LM, SR or LR.");
180 
181  tol_ = pl.get("Convergence Tolerance",tol_);
182  TEUCHOS_TEST_FOR_EXCEPTION(tol_ <= 0,
183  AnasaziError,
184  "SimpleLOBPCGSolMgr: \"Tolerance\" parameter must be strictly postiive.");
185 
186  // verbosity level
187  if (pl.isParameter("Verbosity")) {
188  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
189  verb_ = pl.get("Verbosity", verb_);
190  } else {
191  verb_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
192  }
193  }
194 
195 
196  blockSize_= pl.get("Block Size",problem_->getNEV());
197  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0,
198  AnasaziError,
199  "SimpleLOBPCGSolMgr: \"Block Size\" parameter must be strictly positive.");
200 
201  maxIters_ = pl.get("Maximum Iterations",maxIters_);
202 }
203 
204 
205 
207 template<class ScalarType, class MV, class OP>
210 
211  // sort manager
212  Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
213  // output manager
214  Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verb_) );
215  // status tests
216  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > max;
217  if (maxIters_ > 0) {
218  max = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
219  }
220  else {
221  max = Teuchos::null;
222  }
223  Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > norm
224  = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(tol_) );
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
229  = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>(
231  ));
232  // printing StatusTest
233  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest
234  = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combo,1,Passed ) );
235  // orthomanager
236  Teuchos::RCP<SVQBOrthoManager<ScalarType,MV,OP> > ortho
237  = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
238  // parameter list
239  Teuchos::ParameterList plist;
240  plist.set("Block Size",blockSize_);
241  plist.set("Full Ortho",true);
242 
243  // create an LOBPCG solver
244  Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
245  = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
246  // add the auxillary vecs from the eigenproblem to the solver
247  if (problem_->getAuxVecs() != Teuchos::null) {
248  lobpcg_solver->setAuxVecs( Teuchos::tuple<Teuchos::RCP<const MV> >(problem_->getAuxVecs()) );
249  }
250 
251  int numfound = 0;
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) {
256  // reduce the strain on norm test, if we are almost done
257  if (nev - numfound < blockSize_) {
258  norm->setQuorum(nev-numfound);
259  }
260 
261  // tell the solver to iterate
262  try {
263  lobpcg_solver->iterate();
264  }
265  catch (const std::exception &e) {
266  // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
267  printer->stream(Anasazi::Errors) << "Exception: " << e.what() << std::endl;
269  sol.numVecs = 0;
270  problem_->setSolution(sol);
271  throw;
272  }
273 
274  // check the status tests
275  if (norm->getStatus() == Passed) {
276 
277  int num = norm->howMany();
278  // if num < blockSize_, it is because we are on the last iteration: num+numfound>=nev
279  TEUCHOS_TEST_FOR_EXCEPTION(num < blockSize_ && num+numfound < nev,
280  std::logic_error,
281  "Anasazi::SimpleLOBPCGSolMgr::solve(): logic error.");
282  std::vector<int> ind = norm->whichVecs();
283  // just grab the ones that we need
284  if (num + numfound > nev) {
285  num = nev - numfound;
286  ind.resize(num);
287  }
288 
289  // copy the converged eigenvectors
290  Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
291  // store them
292  foundvecs.push_back(newvecs);
293  // add them as auxiliary vectors
294  Teuchos::Array<Teuchos::RCP<const MV> > auxvecs = lobpcg_solver->getAuxVecs();
295  auxvecs.push_back(newvecs);
296  // setAuxVecs() will reset the solver to uninitialized, without messing with numIters()
297  lobpcg_solver->setAuxVecs(auxvecs);
298 
299  // copy the converged eigenvalues
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;
304  }
305  foundvals.push_back(newvals);
306 
307  numfound += num;
308  }
309  else if (max != Teuchos::null && max->getStatus() == Passed) {
310 
311  int num = norm->howMany();
312  std::vector<int> ind = norm->whichVecs();
313 
314  if (num > 0) {
315  // copy the converged eigenvectors
316  Teuchos::RCP<MV> newvecs = MVT::CloneCopy(*lobpcg_solver->getRitzVectors(),ind);
317  // store them
318  foundvecs.push_back(newvecs);
319  // don't bother adding them as auxiliary vectors; we have reached maxiters and are going to quit
320 
321  // copy the converged eigenvalues
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;
326  }
327  foundvals.push_back(newvals);
328 
329  numfound += num;
330  }
331  break; // while(numfound < nev)
332  }
333  else {
334  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): solver returned without satisfy status test.");
335  }
336  } // end of while(numfound < nev)
337 
338  TEUCHOS_TEST_FOR_EXCEPTION(foundvecs.size() != foundvals.size(),std::logic_error,"Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent array sizes");
339 
340  // create contiguous storage for all eigenvectors, eigenvalues
342  sol.numVecs = numfound;
343  if (numfound > 0) {
344  // allocate space for eigenvectors
345  sol.Evecs = MVT::Clone(*problem_->getInitVec(),numfound);
346  }
347  else {
348  sol.Evecs = Teuchos::null;
349  }
350  sol.Espace = sol.Evecs;
351  // allocate space for eigenvalues
352  std::vector<MagnitudeType> vals(numfound);
353  sol.Evals.resize(numfound);
354  // all real eigenvalues: set index vectors [0,...,numfound-1]
355  sol.index.resize(numfound,0);
356  // store eigenvectors, eigenvalues
357  int curttl = 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;
363  // put the eigenvectors
364  MVT::SetBlock(*foundvecs[i],lclind,*sol.Evecs);
365  // put the eigenvalues
366  std::copy( foundvals[i]->begin(), foundvals[i]->end(), vals.begin()+curttl );
367 
368  curttl += lclnum;
369  }
370  TEUCHOS_TEST_FOR_EXCEPTION( curttl != sol.numVecs, std::logic_error, "Anasazi::SimpleLOBPCGSolMgr::solve(): inconsistent sizes");
371 
372  // sort the eigenvalues and permute the eigenvectors appropriately
373  if (numfound > 0) {
374  std::vector<int> order(sol.numVecs);
375  sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
376  // store the values in the Eigensolution
377  for (int i=0; i<sol.numVecs; i++) {
378  sol.Evals[i].realpart = vals[i];
379  sol.Evals[i].imagpart = MT::zero();
380  }
381  // now permute the eigenvectors according to order
383  msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
384  }
385 
386  // print final summary
387  lobpcg_solver->currentStatus(printer->stream(FinalSummary));
388 
389  // print timing information
390 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
391  if ( printer->isVerbosity( TimingDetails ) ) {
392  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
393  }
394 #endif
395 
396  // send the solution to the eigenproblem
397  problem_->setSolution(sol);
398  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
399 
400  // get the number of iterations performed for this solve.
401  numIters_ = lobpcg_solver->getNumIters();
402 
403  // return from SolMgr::solve()
404  if (sol.numVecs < nev) return Unconverged;
405  return Converged;
406 }
407 
408 
409 
410 
411 } // end Anasazi namespace
412 
413 #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...
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&#39;s templated, static class providing utilities for the solvers.
int numVecs
The number of computed eigenpairs.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Anasazi&#39;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.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
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.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver&#39;s iterate() routine until ...
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
A status test for testing the number of iterations.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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.