Anasazi  Version of the Day
AnasaziTraceMinDavidsonSolMgr.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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
48 #ifndef ANASAZI_TRACEMIN_DAVIDSON_SOLMGR_HPP
49 #define ANASAZI_TRACEMIN_DAVIDSON_SOLMGR_HPP
50 
51 #include "AnasaziConfigDefs.hpp"
54 
73 namespace Anasazi {
74 namespace Experimental {
75 
109 template<class ScalarType, class MV, class OP>
110 class TraceMinDavidsonSolMgr : public TraceMinBaseSolMgr<ScalarType,MV,OP> {
111 
112  private:
115  typedef Teuchos::ScalarTraits<ScalarType> SCT;
116  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
117  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
118 
119  public:
120 
122 
123 
148  TraceMinDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
149  Teuchos::ParameterList &pl );
151 
152  private:
153  int maxRestarts_;
154 
155  // Returns true if the subspace is full
156  bool needToRestart(const Teuchos::RCP< TraceMinBase<ScalarType,MV,OP> > solver)
157  {
158  return (solver->getCurSubspaceDim() == solver->getMaxSubspaceDim());
159  };
160 
161  // Performs a restart by reinitializing TraceMinDavidson with the most significant part of the basis
162  bool performRestart(int &numRestarts, Teuchos::RCP< TraceMinBase<ScalarType,MV,OP> > solver);
163 
164  // Returns a TraceMinDavidson solver
165  Teuchos::RCP< TraceMinBase<ScalarType,MV,OP> > createSolver(
166  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
167  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &outputtest,
168  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
169  Teuchos::ParameterList &plist)
170  {
171  return Teuchos::rcp( new TraceMinDavidson<ScalarType,MV,OP>(this->problem_,sorter,this->printer_,outputtest,ortho,plist) );
172  };
173 };
174 
175 
176 //---------------------------------------------------------------------------//
177 // Prevent instantiation on complex scalar type
178 // FIXME: this really is just a current flaw in the implementation, TraceMin
179 // *should* work for Hermitian matrices
180 //---------------------------------------------------------------------------//
181 template <class MagnitudeType, class MV, class OP>
182 class TraceMinDavidsonSolMgr<std::complex<MagnitudeType>,MV,OP>
183 {
184  public:
185 
186  typedef std::complex<MagnitudeType> ScalarType;
188  const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
189  Teuchos::ParameterList &pl )
190  {
191  // Provide a compile error when attempting to instantiate on complex type
192  MagnitudeType::this_class_is_missing_a_specialization();
193  }
194 };
195 
197 // Basic constructor for TraceMinDavidsonSolMgr
198 template<class ScalarType, class MV, class OP>
199 TraceMinDavidsonSolMgr<ScalarType,MV,OP>::TraceMinDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, Teuchos::ParameterList &pl ) :
200  TraceMinBaseSolMgr<ScalarType,MV,OP>(problem,pl)
201 {
202  // TODO: Come back tot these exceptions and make the descriptions better.
203  maxRestarts_ = pl.get("Maximum Restarts", 50);
204  TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts_ <= 0, std::invalid_argument,
205  "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Maximum Restarts\" must be strictly positive.");
206 
207  this->useHarmonic_ = pl.get("Use Harmonic Ritz Values", false);
208 
209  TEUCHOS_TEST_FOR_EXCEPTION(this->useHarmonic_ && problem->getM() != Teuchos::null, std::invalid_argument, "Anasazi::TraceMinDavidsonSolMgr::constructor(): Harmonic Ritz values do not currently work with generalized eigenvalue problems. Please disable \"Use Harmonic Ritz Values\".");
210 
211  // block size: default is 1
212  this->blockSize_ = pl.get("Block Size", 1);
213  TEUCHOS_TEST_FOR_EXCEPTION(this->blockSize_ <= 0, std::invalid_argument,
214  "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Block Size\" must be strictly positive.");
215 
216  this->numRestartBlocks_ = (int)std::ceil(this->problem_->getNEV() / this->blockSize_);
217  this->numRestartBlocks_ = pl.get("Num Restart Blocks", this->numRestartBlocks_);
218  TEUCHOS_TEST_FOR_EXCEPTION(this->numRestartBlocks_ <= 0, std::invalid_argument,
219  "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Num Restart Blocks\" must be strictly positive.");
220 
221  this->numBlocks_ = pl.get("Num Blocks", 3*this->numRestartBlocks_);
222  TEUCHOS_TEST_FOR_EXCEPTION(this->numBlocks_ <= 1, std::invalid_argument,
223  "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Num Blocks\" must be greater than 1. If you only wish to use one block, please use TraceMinSolMgr instead of TraceMinDavidsonSolMgr.");
224 
225  TEUCHOS_TEST_FOR_EXCEPTION(this->numRestartBlocks_ >= this->numBlocks_, std::invalid_argument,
226  "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Num Blocks\" must be strictly greater than \"Num Restart Blocks\".");
227 
228  std::stringstream ss;
229  ss << "Anasazi::TraceMinDavidsonSolMgr::constructor(): Potentially impossible orthogonality requests. Reduce basis size (" << static_cast<ptrdiff_t>(this->numBlocks_)*this->blockSize_ << ") or locking size (" << this->maxLocked_ << ") because " << static_cast<ptrdiff_t>(this->numBlocks_) << "*" << this->blockSize_ << " + " << this->maxLocked_ << " > " << MVT::GetGlobalLength(*this->problem_->getInitVec()) << ".";
230  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(this->numBlocks_)*this->blockSize_ + this->maxLocked_ > MVT::GetGlobalLength(*this->problem_->getInitVec()),
231  std::invalid_argument, ss.str());
232 
233  TEUCHOS_TEST_FOR_EXCEPTION(this->maxLocked_ + this->blockSize_ < this->problem_->getNEV(), std::invalid_argument,
234  "Anasazi::TraceMinDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
235 }
236 
237 
239 // Performs a restart by reinitializing TraceMinDavidson with the most significant part of the basis
240 template <class ScalarType, class MV, class OP>
242 {
243 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
244  Teuchos::TimeMonitor restimer(*this->_timerRestarting);
245 #endif
246 
247  if ( numRestarts >= maxRestarts_ ) {
248  return false; // break from while(1){tm_solver->iterate()}
249  }
250  numRestarts++;
251 
252  this->printer_->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
253 
254  TraceMinBaseState<ScalarType,MV> oldstate = solver->getState();
255  TraceMinBaseState<ScalarType,MV> newstate;
256  int newdim = this->numRestartBlocks_*this->blockSize_;
257  std::vector<int> indToCopy(newdim);
258  for(int i=0; i<newdim; i++) indToCopy[i] = i;
259 
260  // Copy the relevant parts of the old state to the new one
261  // This may involve computing parts of X
262  if(this->useHarmonic_)
263  {
264  newstate.V = MVT::CloneView(*solver->getRitzVectors(),indToCopy);
265  newstate.curDim = newdim;
266 
267  }
268  else
269  {
270  this->copyPartOfState (oldstate, newstate, indToCopy);
271  }
272 
273  // send the new state to the solver
274  newstate.NEV = oldstate.NEV;
275  solver->initialize(newstate);
276 
277  return true;
278 }
279 
280 
281 }} // end Anasazi namespace
282 
283 #endif /* ANASAZI_TraceMinDavidson_SOLMGR_HPP */
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
This class defines the interface required by an eigensolver and status test class to compute solution...
Virtual base class which defines basic traits for the operator type.
TraceMinDavidsonSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for TraceMinDavidsonSolMgr.
The Anasazi::TraceMinBaseSolMgr provides an abstract base class for the TraceMin series of solver man...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
Traits class which defines basic operations on multivectors.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
The Anasazi::TraceMinBaseSolMgr provides an abstract base class for the TraceMin series of solver man...
The Anasazi::TraceMinDavidsonSolMgr provides a flexible solver manager over the TraceMinDavidson eige...
This class implements a TraceMin-Davidson iteration for solving symmetric generalized eigenvalue prob...
Implementation of the TraceMin-Davidson method.
This is an abstract base class for the trace minimization eigensolvers.
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi&#39;s solvers.