Anasazi  Version of the Day
AnasaziTraceMin.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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
33 #ifndef ANASAZI_TRACEMIN_HPP
34 #define ANASAZI_TRACEMIN_HPP
35 
36 #include "AnasaziTypes.hpp"
37 #include "AnasaziBasicSort.hpp"
38 #include "AnasaziTraceMinBase.hpp"
39 
40 #include "Epetra_Operator.h"
41 
42 #include "AnasaziEigensolver.hpp"
45 #include "Teuchos_ScalarTraits.hpp"
46 
48 #include "AnasaziSolverUtils.hpp"
49 
50 #include "Teuchos_LAPACK.hpp"
51 #include "Teuchos_BLAS.hpp"
52 #include "Teuchos_SerialDenseMatrix.hpp"
53 #include "Teuchos_SerialDenseSolver.hpp"
54 #include "Teuchos_ParameterList.hpp"
55 #include "Teuchos_TimeMonitor.hpp"
56 
57 // TODO: TraceMin performs some unnecessary computations upon restarting. Fix it!
58 
59 namespace Anasazi {
60 namespace Experimental {
112  template <class ScalarType, class MV, class OP>
113  class TraceMin : public TraceMinBase<ScalarType,MV,OP> {
114  public:
116 
117 
158  TraceMin( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
159  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
160  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
161  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
162  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
163  Teuchos::ParameterList &params
164  );
165 
166  private:
167  //
168  // Convenience typedefs
169  //
173  typedef Teuchos::ScalarTraits<ScalarType> SCT;
174  typedef typename SCT::magnitudeType MagnitudeType;
175  const MagnitudeType ONE;
176  const MagnitudeType ZERO;
177  const MagnitudeType NANVAL;
178 
179  // TraceMin specific methods
180  void addToBasis(const Teuchos::RCP<const MV> Delta);
181 
182  void harmonicAddToBasis(const Teuchos::RCP<const MV> Delta);
183  };
184 
187  //
188  // Implementations
189  //
192 
193 
195  // Constructor
196  template <class ScalarType, class MV, class OP>
198  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
199  const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
200  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
201  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
202  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
203  Teuchos::ParameterList &params
204  ) :
205  TraceMinBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params),
206  ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
207  ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
208  NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan())
209  {
210  }
211 
212 
213  template <class ScalarType, class MV, class OP>
214  void TraceMin<ScalarType,MV,OP>::addToBasis(const Teuchos::RCP<const MV> Delta)
215  {
216  MVT::MvAddMv(ONE,*this->X_,-ONE,*Delta,*this->V_);
217 
218  if(this->hasM_)
219  {
220 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
221  Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
222 #endif
223  this->count_ApplyM_+= this->blockSize_;
224 
225  OPT::Apply(*this->MOp_,*this->V_,*this->MV_);
226  }
227 
228  int rank;
229  {
230 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
231  Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
232 #endif
233 
234  if(this->numAuxVecs_ > 0)
235  {
236  rank = this->orthman_->projectAndNormalizeMat(*this->V_,this->auxVecs_,
237  Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
238  Teuchos::null,this->MV_,this->MauxVecs_);
239  }
240  else
241  {
242  rank = this->orthman_->normalizeMat(*this->V_,Teuchos::null,this->MV_);
243  }
244  }
245 
246  // FIXME (mfh 07 Oct 2014) This variable is currently unused, but
247  // it would make sense to use it to check whether the block is
248  // rank deficient.
249  (void) rank;
250 
251  if(this->Op_ != Teuchos::null)
252  {
253 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
254  Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
255 #endif
256  this->count_ApplyOp_+= this->blockSize_;
257  OPT::Apply(*this->Op_,*this->V_,*this->KV_);
258  }
259  }
260 
261 
262 
263  template <class ScalarType, class MV, class OP>
264  void TraceMin<ScalarType,MV,OP>::harmonicAddToBasis(const Teuchos::RCP<const MV> Delta)
265  {
266  // V = X - Delta
267  MVT::MvAddMv(ONE,*this->X_,-ONE,*Delta,*this->V_);
268 
269  // Project out auxVecs
270  if(this->numAuxVecs_ > 0)
271  {
272 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
273  Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
274 #endif
275  this->orthman_->projectMat(*this->V_,this->auxVecs_);
276  }
277 
278  // Compute KV
279  if(this->Op_ != Teuchos::null)
280  {
281 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
282  Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
283 #endif
284  this->count_ApplyOp_+= this->blockSize_;
285 
286  OPT::Apply(*this->Op_,*this->V_,*this->KV_);
287  }
288 
289  // Normalize lclKV
290  RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(this->blockSize_,this->blockSize_));
291  int rank = this->orthman_->normalizeMat(*this->KV_,gamma);
292  // FIXME (mfh 18 Feb 2015) It would make sense to check the rank.
293  (void) rank;
294 
295  // lclV = lclV/gamma
296  Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
297  SDsolver.setMatrix(gamma);
298  SDsolver.invert();
299  RCP<MV> tempMV = MVT::CloneCopy(*this->V_);
300  MVT::MvTimesMatAddMv(ONE,*tempMV,*gamma,ZERO,*this->V_);
301 
302  // Compute MV
303  if(this->hasM_)
304  {
305 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
306  Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
307 #endif
308  this->count_ApplyM_+= this->blockSize_;
309 
310  OPT::Apply(*this->MOp_,*this->V_,*this->MV_);
311  }
312  }
313 
314 }} // End of namespace Anasazi
315 
316 #endif
317 
318 // End of file AnasaziTraceMin.hpp
Abstract base class for trace minimization eigensolvers.
This class implements a TraceMIN iteration, a preconditioned iteration for solving linear symmetric p...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
This class defines the interface required by an eigensolver and status test class to compute solution...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
TraceMin(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
TraceMin constructor with eigenproblem, solver utilities, and parameter list of solver options...
Basic implementation of the Anasazi::SortManager class.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
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.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
Output managers remove the need for the eigensolver to know any information about the required output...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Types and exceptions used within Anasazi solvers and interfaces.
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.
Class which provides internal utilities for the Anasazi solvers.