Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AnasaziGeneralizedDavidsonSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Anasazi: Block Eigensolvers Package
4 //
5 // Copyright 2004 NTESS and the Anasazi contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
11 #define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
12 
18 #include "Teuchos_RCPDecl.hpp"
19 
20 #include "AnasaziConfigDefs.hpp"
21 #include "AnasaziTypes.hpp"
22 #include "AnasaziEigenproblem.hpp"
23 #include "AnasaziSolverManager.hpp"
27 #include "AnasaziOutputManager.hpp"
29 #include "AnasaziBasicSort.hpp"
33 
34 using Teuchos::RCP;
35 
39 namespace Anasazi {
40 
59 template <class ScalarType, class MV, class OP>
60 class GeneralizedDavidsonSolMgr : public SolverManager<ScalarType,MV,OP>
61 {
62  public:
63 
97 
101  const Eigenproblem<ScalarType,MV,OP> & getProblem() const { return *d_problem; }
102 
106  int getNumIters() const { return d_solver->getNumIters(); }
107 
112  ReturnType solve();
113 
114  private:
115 
116  void getRestartState( GeneralizedDavidsonState<ScalarType,MV> &state );
117 
118  typedef MultiVecTraits<ScalarType,MV> MVT;
120  typedef typename ST::magnitudeType MagnitudeType;
122 
125  RCP< OutputManager<ScalarType> > d_outputMan;
129  int d_maxRestarts;
130  int d_restartDim;
131 
132 }; // class GeneralizedDavidsonSolMgr
133 
134 //---------------------------------------------------------------------------//
135 // Prevent instantiation on complex scalar type
136 //---------------------------------------------------------------------------//
137 template <class MagnitudeType, class MV, class OP>
138 class GeneralizedDavidsonSolMgr<std::complex<MagnitudeType>,MV,OP>
139 {
140  public:
141 
142  typedef std::complex<MagnitudeType> ScalarType;
144  const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
146  {
147  // Provide a compile error when attempting to instantiate on complex type
148  MagnitudeType::this_class_is_missing_a_specialization();
149  }
150 };
151 
152 //---------------------------------------------------------------------------//
153 // Start member definitions
154 //---------------------------------------------------------------------------//
155 
156 //---------------------------------------------------------------------------//
157 // Constructor
158 //---------------------------------------------------------------------------//
159 template <class ScalarType, class MV, class OP>
161  const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
163  : d_problem(problem)
164 {
165  TEUCHOS_TEST_FOR_EXCEPTION( d_problem == Teuchos::null, std::invalid_argument, "Problem not given to solver manager." );
166  TEUCHOS_TEST_FOR_EXCEPTION( !d_problem->isProblemSet(), std::invalid_argument, "Problem not set." );
167  TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getA() == Teuchos::null &&
168  d_problem->getOperator() == Teuchos::null, std::invalid_argument, "A operator not supplied on Eigenproblem." );
169  TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getInitVec() == Teuchos::null, std::invalid_argument, "No vector to clone from on Eigenproblem." );
170  TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV() <= 0, std::invalid_argument, "Number of requested eigenvalues must be positive.");
171 
172  if( !pl.isType<int>("Block Size") )
173  {
174  pl.set<int>("Block Size",1);
175  }
176 
177  if( !pl.isType<int>("Maximum Subspace Dimension") )
178  {
179  pl.set<int>("Maximum Subspace Dimension",3*problem->getNEV()*pl.get<int>("Block Size"));
180  }
181 
182  if( !pl.isType<int>("Print Number of Ritz Values") )
183  {
184  int numToPrint = std::max( pl.get<int>("Block Size"), d_problem->getNEV() );
185  pl.set<int>("Print Number of Ritz Values",numToPrint);
186  }
187 
188  // Get convergence info
189  MagnitudeType tol = pl.get<MagnitudeType>("Convergence Tolerance", MT::eps() );
190  TEUCHOS_TEST_FOR_EXCEPTION( pl.get<MagnitudeType>("Convergence Tolerance") <= MT::zero(),
191  std::invalid_argument, "Convergence Tolerance must be greater than zero." );
192 
193  // Get maximum restarts
194  if( pl.isType<int>("Maximum Restarts") )
195  {
196  d_maxRestarts = pl.get<int>("Maximum Restarts");
197  TEUCHOS_TEST_FOR_EXCEPTION( d_maxRestarts < 0, std::invalid_argument, "Maximum Restarts must be non-negative" );
198  }
199  else
200  {
201  d_maxRestarts = 20;
202  }
203 
204  // Get maximum restarts
205  d_restartDim = pl.get<int>("Restart Dimension",d_problem->getNEV());
206  TEUCHOS_TEST_FOR_EXCEPTION( d_restartDim < d_problem->getNEV(),
207  std::invalid_argument, "Restart Dimension must be at least NEV" );
208 
209  // Get initial guess type
210  std::string initType;
211  if( pl.isType<std::string>("Initial Guess") )
212  {
213  initType = pl.get<std::string>("Initial Guess");
214  TEUCHOS_TEST_FOR_EXCEPTION( initType!="User" && initType!="Random", std::invalid_argument,
215  "Initial Guess type must be 'User' or 'Random'." );
216  }
217  else
218  {
219  initType = "User";
220  }
221 
222  // Get sort type
223  std::string which;
224  if( pl.isType<std::string>("Which") )
225  {
226  which = pl.get<std::string>("Which");
227  TEUCHOS_TEST_FOR_EXCEPTION( which!="LM" && which!="SM" && which!="LR" && which!="SR" && which!="LI" && which!="SI",
228  std::invalid_argument,
229  "Which must be one of LM,SM,LR,SR,LI,SI." );
230  }
231  else
232  {
233  which = "LM";
234  }
235 
236  // Build sort manager (currently must be stored as pointer to derived class)
237  d_sortMan = Teuchos::rcp( new BasicSort<MagnitudeType>(which) );
238 
239  // Build orthogonalization manager
240  std::string ortho = pl.get<std::string>("Orthogonalization","SVQB");
241  TEUCHOS_TEST_FOR_EXCEPTION( ortho!="DGKS" && ortho!= "SVQB" && ortho!="ICGS", std::invalid_argument,
242  "Anasazi::GeneralizedDavidsonSolMgr::constructor: Invalid orthogonalization type" );
243 
244  if( ortho=="DGKS" )
245  {
247  }
248  else if( ortho=="SVQB" )
249  {
250  d_orthoMan = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>() );
251  }
252  else if( ortho=="ICGS" )
253  {
254  d_orthoMan = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>() );
255  }
256 
257  // Build StatusTest
258  bool scaleRes = false; // Always false, scaling the residual is handled by the solver
259  bool failOnNaN = false;
261  new StatusTestResNorm<ScalarType,MV,OP>(tol,d_problem->getNEV(),
262  RES_2NORM,scaleRes,failOnNaN) );
263  d_tester = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(resNormTest,d_sortMan,d_problem->getNEV()) );
264 
265  // Build output manager
266 
267  // Create a formatted output stream to print to.
268  // See if user requests output processor.
269  int osProc = pl.get("Output Processor", 0);
270 
271  // If not passed in by user, it will be chosen based upon operator type.
273 
274  if (pl.isParameter("Output Stream")) {
275  osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
276  }
277  else {
278  osp = OutputStreamTraits<OP>::getOutputStream (*d_problem->getOperator(), osProc);
279  }
280 
281  // verbosity
282  int verbosity = Anasazi::Errors;
283  if (pl.isParameter("Verbosity")) {
284  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
285  verbosity = pl.get("Verbosity", verbosity);
286  } else {
287  verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
288  }
289  }
290 
291  d_outputMan = Teuchos::rcp( new OutputManager<ScalarType>(verbosity,osp) );
292 
293  // Build solver
294  d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidsonSolMgr: Building solver" << std::endl;
295  d_solver = Teuchos::rcp( new GeneralizedDavidson<ScalarType,MV,OP>( problem, d_sortMan, d_outputMan, d_tester, d_orthoMan, pl ) );
296 
297  TEUCHOS_TEST_FOR_EXCEPTION(d_solver->getMaxSubspaceDim() < d_restartDim, std::invalid_argument,
298  "The maximum size of the subspace dimension (" << d_solver->getMaxSubspaceDim() << ") must "
299  "not be smaller than the size of the restart space (" << d_restartDim << "). "
300  "Please adjust \"Restart Dimension\" and/or \"Maximum Subspace Dimension\" parameters.");
301 
302 }
303 
304 //---------------------------------------------------------------------------//
305 // Solve
306 //---------------------------------------------------------------------------//
307 template <class ScalarType, class MV, class OP>
309 {
311  sol.numVecs = 0;
312  d_problem->setSolution(sol);
313 
314  d_solver->initialize();
315  int restarts = 0;
316  while( 1 )
317  {
318  // Call iterate on the solver
319  d_solver->iterate();
320 
321  // If the solver converged, we're done
322  if( d_tester->getStatus() == Passed )
323  break;
324 
325  // If we're already at maximum number of restarts, wrap it up
326  if( restarts == d_maxRestarts )
327  break;
328 
329  // We need to restart
330  d_solver->sortProblem( d_restartDim );
331  GeneralizedDavidsonState<ScalarType,MV> state = d_solver->getState();
332  getRestartState( state );
333  d_solver->initialize( state );
334  restarts++;
335  }
336 
337  // Output final state
338  if( d_outputMan->isVerbosity(FinalSummary) )
339  d_solver->currentStatus(d_outputMan->stream(FinalSummary));
340 
341  // Fill solution struct
342  sol.numVecs = d_tester->howMany();
343  if( sol.numVecs > 0 )
344  {
345  std::vector<int> whichVecs = d_tester->whichVecs();
346  std::vector<int> origIndex = d_solver->getRitzIndex();
347 
348  // Make sure no conjugate pairs are split
349  // Because these are not sorted we have to check all values
350  for( int i=0; i<sol.numVecs; ++i )
351  {
352  if( origIndex[ whichVecs[i] ] == 1 )
353  {
354  if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]+1 ) == whichVecs.end() )
355  {
356  whichVecs.push_back( whichVecs[i]+1 );
357  sol.numVecs++;
358  }
359  }
360  else if( origIndex[ whichVecs[i] ] == -1 )
361  {
362  if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]-1 ) == whichVecs.end() )
363  {
364  whichVecs.push_back( whichVecs[i]-1 );
365  sol.numVecs++;
366  }
367  }
368  }
369 
370  if( d_outputMan->isVerbosity(Debug) )
371  {
372  d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidsonSolMgr: "
373  << sol.numVecs << " eigenpairs converged" << std::endl;
374  }
375 
376  // Sort converged values
377  std::vector< Value<ScalarType> > origVals = d_solver->getRitzValues();
378  std::vector<MagnitudeType> realParts;
379  std::vector<MagnitudeType> imagParts;
380  for( int i=0; i<sol.numVecs; ++i )
381  {
382  realParts.push_back( origVals[whichVecs[i]].realpart );
383  imagParts.push_back( origVals[whichVecs[i]].imagpart );
384  }
385 
386  std::vector<int> permVec(sol.numVecs);
387  d_sortMan->sort( realParts, imagParts, Teuchos::rcpFromRef(permVec), sol.numVecs );
388 
389  // Create new which vector
390  std::vector<int> newWhich;
391  for( int i=0; i<sol.numVecs; ++i )
392  newWhich.push_back( whichVecs[permVec[i]] );
393 
394  // Check if converged vectors are ordered
395  bool ordered = true;
396  for( int i=0; i<sol.numVecs; ++i )
397  {
398  if( newWhich[i]!=i )
399  {
400  ordered = false;
401  break;
402  }
403  }
404 
405  if( ordered )
406  {
407  // Everything is ordered, pull directly from solver and resize
408  sol.index = origIndex;
409  sol.index.resize(sol.numVecs);
410  sol.Evals = d_solver->getRitzValues();
411  sol.Evals.resize(sol.numVecs);
412  }
413  else
414  {
415  // Manually copy values into sol
416 
417  sol.index.resize(sol.numVecs);
418  sol.Evals.resize(sol.numVecs);
419 
420  for( int i=0; i<sol.numVecs; ++i )
421  {
422  sol.index[i] = origIndex[ newWhich[i] ];
423  sol.Evals[i] = origVals[ newWhich[i] ];
424  }
425  }
426  sol.Evecs = MVT::CloneCopy( *(d_solver->getRitzVectors()), newWhich );
427  }
428  d_problem->setSolution(sol);
429 
430  // Return convergence status
431  if( sol.numVecs < d_problem->getNEV() )
432  return Unconverged;
433 
434  return Converged;
435 }
436 
437 //---------------------------------------------------------------------------//
438 // Update GeneralizedDavidson state for restarting
439 //---------------------------------------------------------------------------//
440 template <class ScalarType, class MV, class OP>
443 {
444  TEUCHOS_TEST_FOR_EXCEPTION( state.curDim <= d_restartDim, std::runtime_error,
445  "Anasazi::GeneralizedDavidsonSolMgr: State dimension at restart is smaller than Restart Dimension" );
446 
447  std::vector<int> ritzIndex = d_solver->getRitzIndex();
448 
449  // Don't split conjugate pair when restarting
450  int restartDim = d_restartDim;
451  if( ritzIndex[d_restartDim-1]==1 )
452  restartDim++;
453 
454  d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidsonSolMgr: Restarting with "
455  << restartDim << " vectors" << std::endl;
456 
457  // We have already sorted the problem with d_restartDim "best" values
458  // in the leading position. If we partition the Schur vectors (Z)
459  // of the projected problem as Z = [Z_wanted Z_unwanted], then the
460  // search subspace after the restart is V_restart = V*Z_wanted
461  // (same for AV,BV)
462 
463  // Get view of wanted portion of Z
466 
467  // Get indices for restart
468  std::vector<int> allIndices(state.curDim);
469  for( int i=0; i<state.curDim; ++i )
470  allIndices[i] = i;
471 
472  RCP<const MV> V_orig = MVT::CloneView( *state.V, allIndices );
473 
474  // Get indices for restart
475  std::vector<int> restartIndices(restartDim);
476  for( int i=0; i<restartDim; ++i )
477  restartIndices[i] = i;
478 
479  // Views of subspace vectors to be updated
480  RCP<MV> V_restart = MVT::CloneViewNonConst( *state.V, restartIndices );
481 
482  // Temp storage
483  RCP<MV> restartVecs = MVT::Clone(*state.V,restartDim);
484 
485  // Reset V
486  MVT::MvTimesMatAddMv(ST::one(),*V_orig,Z_wanted,ST::zero(),*restartVecs);
487  MVT::SetBlock(*restartVecs,restartIndices,*V_restart);
488 
489  // V, Z each have orthonormal columns, therefore V*Z should as well
490  if( d_outputMan->isVerbosity(Debug) )
491  {
492  MagnitudeType orthErr = d_orthoMan->orthonormError(*V_restart);
493  std::stringstream os;
494  os << " >> Anasazi::GeneralizedDavidsonSolMgr: Error in V^T V == I after restart : " << orthErr << std::endl;
495  d_outputMan->print(Debug,os.str());
496  }
497 
498  // Reset AV
499  RCP<MV> AV_restart = MVT::CloneViewNonConst( *state.AV, restartIndices );
500  RCP<const MV> AV_orig = MVT::CloneView( *state.AV, allIndices );
501 
502  MVT::MvTimesMatAddMv(ST::one(),*AV_orig,Z_wanted,ST::zero(),*restartVecs);
503  MVT::SetBlock(*restartVecs,restartIndices,*AV_restart);
504 
505  int err;
506 
507  // Update matrix projection as Z^{*}(V^{*}AV)Z
508  const Teuchos::SerialDenseMatrix<int,ScalarType> VAV_orig( Teuchos::View, *state.VAV, state.curDim, state.curDim );
509  Teuchos::SerialDenseMatrix<int,ScalarType> tmpMat(state.curDim, restartDim);
510  err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VAV_orig, Z_wanted, ST::zero() );
511  TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
512 
513  Teuchos::SerialDenseMatrix<int,ScalarType> VAV_restart( Teuchos::View, *state.VAV, restartDim, restartDim );
514  err = VAV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
515  TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
516 
517  if( d_problem->getM() != Teuchos::null )
518  {
519  // Reset BV
520  RCP<const MV> BV_orig = MVT::CloneView( *state.BV, allIndices );
521  RCP<MV> BV_restart = MVT::CloneViewNonConst( *state.BV, restartIndices );
522 
523  MVT::MvTimesMatAddMv(ST::one(),*BV_orig,Z_wanted,ST::zero(),*restartVecs);
524  MVT::SetBlock(*restartVecs,restartIndices,*BV_restart);
525 
526 
527  // Update matrix projection as Z^{*}(V^{*}BV)Z
528  const Teuchos::SerialDenseMatrix<int,ScalarType> VBV_orig( Teuchos::View, *state.VBV, state.curDim, state.curDim );
529  err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VBV_orig, Z_wanted, ST::zero() );
530  TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
531 
532  Teuchos::SerialDenseMatrix<int,ScalarType> VBV_restart( Teuchos::View, *state.VBV, restartDim, restartDim );
533  VBV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
534  TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
535  }
536 
537  // Set Q,Z to identity
538  state.Q->putScalar( ST::zero() );
539  state.Z->putScalar( ST::zero() );
540  for( int ii=0; ii<restartDim; ii++ )
541  {
542  (*state.Q)(ii,ii)= ST::one();
543  (*state.Z)(ii,ii)= ST::one();
544  }
545 
546  // Update current dimension
547  state.curDim = restartDim;
548 }
549 
550 } // namespace Anasazi
551 
552 #endif // ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
553 
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver&#39;s iterate() routine until ...
Pure virtual base class which describes the basic interface for a solver manager. ...
RCP< MV > V
Orthonormal basis for search subspace.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
This class defines the interface required by an eigensolver and status test class to compute solution...
static magnitudeType eps()
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
Solves eigenvalue problem using generalized Davidson method.
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VAV
Projection of A onto V.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
int curDim
The current subspace dimension.
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...
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 along with a set of auxiliary eigenv...
Structure to contain pointers to GeneralizedDavidson state variables.
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...
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.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Q
Left generalized Schur vectors from QZ decomposition of (VAV,VBV)
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...
int getNumIters() const
Get the iteration count for the most recent call to solve()
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
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.
GeneralizedDavidsonSolMgr(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for GeneralizedDavidsonSolMgr.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > VBV
Projection of B onto V.
RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Z
Right generalized Schur vectors from QZ decomposition of (VAV,VBV)
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
bool isType(const std::string &name) const
Abstract class definition for Anasazi output stream.
Solver Manager for GeneralizedDavidson.
A status test for testing the norm of the eigenvectors residuals.
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::OrthoManager class.
Implementation of a block Generalized Davidson eigensolver.