Anasazi  Version of the Day
AnasaziRTRSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_RTR_SOLMGR_HPP
30 #define ANASAZI_RTR_SOLMGR_HPP
31 
37 #include "AnasaziConfigDefs.hpp"
38 #include "AnasaziTypes.hpp"
39 
40 #include "AnasaziEigenproblem.hpp"
41 #include "AnasaziSolverManager.hpp"
42 #include "AnasaziSolverUtils.hpp"
43 
44 #include "AnasaziIRTR.hpp"
45 #include "AnasaziSIRTR.hpp"
46 #include "AnasaziBasicSort.hpp"
54 
55 #include <Teuchos_TimeMonitor.hpp>
56 #include <Teuchos_FancyOStream.hpp>
57 
67 namespace Anasazi {
68 
69 template<class ScalarType, class MV, class OP>
70 class RTRSolMgr : public SolverManager<ScalarType,MV,OP> {
71 
72  private:
75  typedef Teuchos::ScalarTraits<ScalarType> SCT;
76  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
77  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
78 
79  public:
80 
82 
83 
99  RTRSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
100  Teuchos::ParameterList &pl );
101 
103  virtual ~RTRSolMgr() {};
105 
107 
108 
111  return *problem_;
112  }
113 
119  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
120  return Teuchos::tuple(_timerSolve);
121  }
122 
124  int getNumIters() const {
125  return numIters_;
126  }
127 
128 
130 
132 
133 
142  ReturnType solve();
144 
145  private:
146  Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
147  std::string whch_;
148  std::string ortho_;
149 
150  bool skinny_;
151  MagnitudeType convtol_;
152  int maxIters_;
153  bool relconvtol_;
154  enum ResType convNorm_;
155  int numIters_;
156  int numICGS_;
157 
158  Teuchos::RCP<Teuchos::Time> _timerSolve;
159  Teuchos::RCP<BasicOutputManager<ScalarType> > printer_;
160  Teuchos::ParameterList pl_;
161 };
162 
163 
165 template<class ScalarType, class MV, class OP>
167  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
168  Teuchos::ParameterList &pl ) :
169  problem_(problem),
170  whch_("SR"),
171  skinny_(true),
172  convtol_(MT::prec()),
173  maxIters_(100),
174  relconvtol_(true),
175  numIters_(-1),
176 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
177  _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: RTRSolMgr::solve()")),
178 #endif
179  pl_(pl)
180 {
181  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
182  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
183  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
184  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
185 
186  std::string strtmp;
187 
188  whch_ = pl_.get("Which","SR");
189  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SR" && whch_ != "LR",
190  std::invalid_argument, "Anasazi::RTRSolMgr: Invalid sorting string. RTR solvers compute only LR or SR.");
191 
192  // convergence tolerance
193  convtol_ = pl_.get("Convergence Tolerance",convtol_);
194  relconvtol_ = pl_.get("Relative Convergence Tolerance",relconvtol_);
195  strtmp = pl_.get("Convergence Norm",std::string("2"));
196  if (strtmp == "2") {
197  convNorm_ = RES_2NORM;
198  }
199  else if (strtmp == "M") {
200  convNorm_ = RES_ORTH;
201  }
202  else {
203  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
204  "Anasazi::RTRSolMgr: Invalid Convergence Norm.");
205  }
206 
207 
208  // maximum number of (outer) iterations
209  maxIters_ = pl_.get("Maximum Iterations",maxIters_);
210 
211  // skinny solver or not
212  skinny_ = pl_.get("Skinny Solver",skinny_);
213 
214  // number if ICGS iterations
215  numICGS_ = pl_.get("Num ICGS",2);
216 
217  // output stream
218  std::string fntemplate = "";
219  bool allProcs = false;
220  if (pl_.isParameter("Output on all processors")) {
221  if (Teuchos::isParameterType<bool>(pl_,"Output on all processors")) {
222  allProcs = pl_.get("Output on all processors",allProcs);
223  } else {
224  allProcs = ( Teuchos::getParameter<int>(pl_,"Output on all processors") != 0 );
225  }
226  }
227  fntemplate = pl_.get("Output filename template",fntemplate);
228  int MyPID;
229 # ifdef HAVE_MPI
230  // Initialize MPI
231  int mpiStarted = 0;
232  MPI_Initialized(&mpiStarted);
233  if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
234  else MyPID=0;
235 # else
236  MyPID = 0;
237 # endif
238  if (fntemplate != "") {
239  std::ostringstream MyPIDstr;
240  MyPIDstr << MyPID;
241  // replace %d in fntemplate with MyPID
242  int pos, start=0;
243  while ( (pos = fntemplate.find("%d",start)) != -1 ) {
244  fntemplate.replace(pos,2,MyPIDstr.str());
245  start = pos+2;
246  }
247  }
248  Teuchos::RCP<ostream> osp;
249  if (fntemplate != "") {
250  osp = Teuchos::rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
251  if (!*osp) {
252  osp = Teuchos::rcpFromRef(std::cout);
253  std::cout << "Anasazi::RTRSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
254  }
255  }
256  else {
257  osp = Teuchos::rcpFromRef(std::cout);
258  }
259  // Output manager
260  int verbosity = Anasazi::Errors;
261  if (pl_.isParameter("Verbosity")) {
262  if (Teuchos::isParameterType<int>(pl_,"Verbosity")) {
263  verbosity = pl_.get("Verbosity", verbosity);
264  } else {
265  verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl_,"Verbosity");
266  }
267  }
268  if (allProcs) {
269  // print on all procs
270  printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
271  }
272  else {
273  // print only on proc 0
274  printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
275  }
276 }
277 
278 
279 // solve()
280 template<class ScalarType, class MV, class OP>
283 
284  using std::endl;
285 
286  // typedef SolverUtils<ScalarType,MV,OP> msutils; // unused
287 
288  const int nev = problem_->getNEV();
289 
290  // clear num iters
291  numIters_ = -1;
292 
293 #ifdef TEUCHOS_DEBUG
294  Teuchos::RCP<Teuchos::FancyOStream>
295  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
296  out->setShowAllFrontMatter(false).setShowProcRank(true);
297  *out << "Entering Anasazi::RTRSolMgr::solve()\n";
298 #endif
299 
301  // Sort manager
302  Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
303 
305  // Status tests
306  //
307  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > maxtest;
308  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > ordertest;
309  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > combotest;
310  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
311  // maximum iters
312  if (maxIters_ > 0) {
313  maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
314  }
315  else {
316  maxtest = Teuchos::null;
317  }
318  //
319  // convergence
320  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
321  ordertest = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
322  //
323  // combo
324  Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
325  alltests.push_back(ordertest);
326  if (maxtest != Teuchos::null) alltests.push_back(maxtest);
327  combotest = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
328  //
329  // printing StatusTest
330  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
331  if ( printer_->isVerbosity(Debug) ) {
332  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
333  }
334  else {
335  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
336  }
337 
339  // Orthomanager
340  Teuchos::RCP<ICGSOrthoManager<ScalarType,MV,OP> > ortho
341  = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM(),numICGS_) );
342 
344  // create an RTR solver
345  // leftmost or rightmost?
346  bool leftMost = true;
347  if (whch_ == "LR" || whch_ == "LM") {
348  leftMost = false;
349  }
350  pl_.set<bool>("Leftmost",leftMost);
351  Teuchos::RCP<RTRBase<ScalarType,MV,OP> > rtr_solver;
352  if (skinny_ == false) {
353  // "hefty" IRTR
354  rtr_solver = Teuchos::rcp( new IRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
355  }
356  else {
357  // "skinny" IRTR
358  rtr_solver = Teuchos::rcp( new SIRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
359  }
360  // set any auxiliary vectors defined in the problem
361  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
362  if (probauxvecs != Teuchos::null) {
363  rtr_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
364  }
365 
366  TEUCHOS_TEST_FOR_EXCEPTION(rtr_solver->getBlockSize() < problem_->getNEV(),std::logic_error,
367  "Anasazi::RTRSolMgr requires block size >= requested number of eigenvalues.");
368 
369  int numfound = 0;
370  Teuchos::RCP<MV> foundvecs;
371  std::vector<MagnitudeType> foundvals;
372 
373  // tell the solver to iterate
374  try {
375 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
376  Teuchos::TimeMonitor slvtimer(*_timerSolve);
377 #endif
378  rtr_solver->iterate();
379  numIters_ = rtr_solver->getNumIters();
380  }
381  catch (const std::exception &e) {
382  // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
383  printer_->stream(Anasazi::Errors) << "Exception: " << e.what() << endl;
385  sol.numVecs = 0;
386  problem_->setSolution(sol);
387  throw;
388  }
389 
390  // check the status tests
391  if (convtest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed))
392  {
393  int num = convtest->howMany();
394  if (num > 0) {
395  std::vector<int> ind = convtest->whichVecs();
396  // copy the converged eigenvectors
397  foundvecs = MVT::CloneCopy(*rtr_solver->getRitzVectors(),ind);
398  // copy the converged eigenvalues
399  foundvals.resize(num);
400  std::vector<Value<ScalarType> > all = rtr_solver->getRitzValues();
401  for (int i=0; i<num; i++) {
402  foundvals[i] = all[ind[i]].realpart;
403  }
404  numfound = num;
405  }
406  }
407  else {
408  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::RTRSolMgr::solve(): solver returned without satisfy status test.");
409  }
410 
411  // create contiguous storage for all eigenvectors, eigenvalues
413  sol.numVecs = numfound;
414  sol.Evecs = foundvecs;
415  sol.Espace = sol.Evecs;
416  sol.Evals.resize(sol.numVecs);
417  for (int i=0; i<sol.numVecs; i++) {
418  sol.Evals[i].realpart = foundvals[i];
419  }
420  // all real eigenvalues: set index vectors [0,...,numfound-1]
421  sol.index.resize(numfound,0);
422 
423  // print final summary
424  rtr_solver->currentStatus(printer_->stream(FinalSummary));
425 
426  // print timing information
427 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
428  if ( printer_->isVerbosity( TimingDetails ) ) {
429  Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
430  }
431 #endif
432 
433  // send the solution to the eigenproblem
434  problem_->setSolution(sol);
435  printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
436 
437  // return from SolMgr::solve()
438  if (sol.numVecs < nev) return Unconverged;
439  return Converged;
440 }
441 
442 
443 
444 
445 } // end Anasazi namespace
446 
447 #endif /* ANASAZI_RTR_SOLMGR_HPP */
Pure virtual base class which describes the basic interface for a solver manager. ...
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
A special StatusTest for printing other status tests.
This class defines the interface required by an eigensolver and status test class to compute solution...
RTRSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for RTRSolMgr.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
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...
int getNumIters() const
Get the iteration count for the most recent call to solve.
The Anasazi::RTRSolMgr provides a simple solver manager over the RTR eigensolver. For more informatio...
virtual ~RTRSolMgr()
Destructor.
Basic implementation of the Anasazi::SortManager class.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
int numVecs
The number of computed eigenpairs.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver&#39;s iterate() routine until ...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
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.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Struct for storing an eigenproblem solution.
A status test for testing the number of iterations.
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 along with a set of auxiliary eigenv...
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.
Basic implementation of the Anasazi::OrthoManager class.
Class which provides internal utilities for the Anasazi solvers.