Anasazi  Version of the Day
AnasaziBasicEigenproblem.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_BASIC_EIGENPROBLEM_H
30 #define ANASAZI_BASIC_EIGENPROBLEM_H
31 
36 #include "AnasaziEigenproblem.hpp"
39 
45 namespace Anasazi {
46 
47  template<class ScalarType, class MV, class OP>
48  class BasicEigenproblem : public virtual Eigenproblem<ScalarType, MV, OP> {
49 
50  public:
51 
53 
54 
57 
59  BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<MV>& InitVec );
60 
62  BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<const OP>& B, const Teuchos::RCP<MV>& InitVec );
63 
66 
68  virtual ~BasicEigenproblem() {};
70 
72 
73 
80  void setOperator( const Teuchos::RCP<const OP>& Op ) { _Op = Op; _isSet=false; };
81 
84  void setA( const Teuchos::RCP<const OP>& A ) { _AOp = A; _isSet=false; };
85 
88  void setM( const Teuchos::RCP<const OP>& M ) { _MOp = M; _isSet=false; };
89 
92  void setPrec( const Teuchos::RCP<const OP>& Prec ) { _Prec = Prec; _isSet=false; };
93 
101  void setInitVec( const Teuchos::RCP<MV>& InitVec ) { _InitVec = InitVec; _isSet=false; };
102 
108  void setAuxVecs( const Teuchos::RCP<const MV>& AuxVecs ) { _AuxVecs = AuxVecs; _isSet=false; };
109 
111  void setNEV( int nev ){ _nev = nev; _isSet=false; };
112 
114 
117  void setHermitian( bool isSym ){ _isSym = isSym; _isSet=false; };
118 
134  bool setProblem();
135 
143 
145 
147 
148 
150  Teuchos::RCP<const OP> getOperator() const { return( _Op ); };
151 
153  Teuchos::RCP<const OP> getA() const { return( _AOp ); };
154 
156  Teuchos::RCP<const OP> getM() const { return( _MOp ); };
157 
159  Teuchos::RCP<const OP> getPrec() const { return( _Prec ); };
160 
162  Teuchos::RCP<const MV> getInitVec() const { return( _InitVec ); };
163 
165  Teuchos::RCP<const MV> getAuxVecs() const { return( _AuxVecs ); };
166 
168  int getNEV() const { return( _nev ); }
169 
171  bool isHermitian() const { return( _isSym ); }
172 
174  bool isProblemSet() const { return( _isSet ); }
175 
181  const Eigensolution<ScalarType,MV> & getSolution() const { return(_sol); }
182 
184 
186 
187 
195  Teuchos::RCP<const MV> computeCurrResVec() const;
196 
198 
199  protected:
200 
202  Teuchos::RCP<const OP> _AOp;
203 
205  Teuchos::RCP<const OP> _MOp;
206 
208  Teuchos::RCP<const OP> _Op;
209 
211  Teuchos::RCP<const OP> _Prec;
212 
214  Teuchos::RCP<MV> _InitVec;
215 
217  Teuchos::RCP<const MV> _AuxVecs;
218 
220  int _nev;
221 
223 
226  bool _isSym;
227 
229  bool _isSet;
230 
235 
238  };
239 
240 
241  //=============================================================================
242  // Implementations (Constructors / Destructors)
243  //=============================================================================
244  template <class ScalarType, class MV, class OP>
246  _nev(0),
247  _isSym(false),
248  _isSet(false)
249  {
250  }
251 
252 
253  //=============================================================================
254  template <class ScalarType, class MV, class OP>
255  BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<MV>& InitVec ) :
256  _Op(Op),
257  _InitVec(InitVec),
258  _nev(0),
259  _isSym(false),
260  _isSet(false)
261  {
262  }
263 
264 
265  //=============================================================================
266  template <class ScalarType, class MV, class OP>
267  BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<const OP>& M,
268  const Teuchos::RCP<MV>& InitVec ) :
269  _MOp(M),
270  _Op(Op),
271  _InitVec(InitVec),
272  _nev(0),
273  _isSym(false),
274  _isSet(false)
275  {
276  }
277 
278 
279  //=============================================================================
280  template <class ScalarType, class MV, class OP>
282  _AOp(Problem._AOp),
283  _MOp(Problem._MOp),
284  _Op(Problem._Op),
285  _Prec(Problem._Prec),
286  _InitVec(Problem._InitVec),
287  _nev(Problem._nev),
288  _isSym(Problem._isSym),
289  _isSet(Problem._isSet),
290  _sol(Problem._sol)
291  {
292  }
293 
294 
295  //=============================================================================
296  // SetProblem (sanity check method)
297  //=============================================================================
298  template <class ScalarType, class MV, class OP>
300  {
301  //----------------------------------------------------------------
302  // Sanity Checks
303  //----------------------------------------------------------------
304  // If there is no operator, then we can't proceed.
305  if ( !_AOp.get() && !_Op.get() ) { return false; }
306 
307  // If there is no initial vector, then we don't have anything to clone workspace from.
308  if ( !_InitVec.get() ) { return false; }
309 
310  // If we don't need any eigenvalues, we don't need to continue.
311  if (_nev == 0) { return false; }
312 
313  // If there is an A, but no operator, we can set them equal.
314  if (_AOp.get() && !_Op.get()) { _Op = _AOp; }
315 
316  // Clear the storage from any previous call to setSolution()
318  _sol = emptysol;
319 
320  // mark the problem as set and return no-error
321  _isSet=true;
322  return true;
323  }
324 
325  //=============================================================================
326  // Computes the residual vector
327  //=============================================================================
328  template <class ScalarType, class MV, class OP>
330  {
331  using Teuchos::RCP;
332 
333  TEUCHOS_TEST_FOR_EXCEPTION(!isHermitian(), std::invalid_argument,
334  "BasicEigenproblem::computeCurrResVec: This method is not currently "
335  "implemented for non-Hermitian problems. Sorry for any inconvenience.");
336 
337  const Eigensolution<ScalarType,MV> sol = getSolution();
338 
339  if(sol.numVecs <= 0)
340  return Teuchos::null;
341 
342  // Copy the eigenvalues
343  RCP<MV> X = sol.Evecs;
344  std::vector<ScalarType> Lambda(sol.numVecs);
345  for(int i = 0; i < sol.numVecs; i++)
346  Lambda[i] = sol.Evals[i].realpart;
347 
348  // Compute AX
349  RCP<MV> AX = MVT::Clone(*X,sol.numVecs);
350  if(_AOp != Teuchos::null)
351  OPT::Apply(*_AOp,*X,*AX);
352  else
353  OPT::Apply(*_Op,*X,*AX);
354 
355  // Compute MX if necessary
356  RCP<MV> MX;
357  if(_MOp != Teuchos::null)
358  {
359  MX = MVT::Clone(*X,sol.numVecs);
360  OPT::Apply(*_MOp,*X,*MX);
361  }
362  else
363  {
364  MX = MVT::CloneCopy(*X);
365  }
366 
367  // Compute R = AX - MX \Lambda
368  RCP<MV> R = MVT::Clone(*X,sol.numVecs);
369  MVT::MvScale(*MX,Lambda);
370  MVT::MvAddMv(1.0,*AX,-1.0,*MX,*R);
371 
372  return R;
373  }
374 
375 } // end Anasazi namespace
376 #endif
377 
378 // end AnasaziBasicEigenproblem.hpp
bool setProblem()
Specify that this eigenproblem is fully defined.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
MultiVecTraits< ScalarType, MV > MVT
Type-definition for the MultiVecTraits class corresponding to the MV type.
void setInitVec(const Teuchos::RCP< MV > &InitVec)
Set the initial guess.
This class defines the interface required by an eigensolver and status test class to compute solution...
bool isHermitian() const
Get the symmetry information for this eigenproblem.
Declaration of basic traits for the multivector type.
OperatorTraits< ScalarType, MV, OP > OPT
Type-definition for the OperatorTraits class corresponding to the OP type.
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
int getNEV() const
Get the number of eigenvalues (NEV) that are required by this eigenproblem.
Teuchos::RCP< const MV > _AuxVecs
Reference-counted pointer for the auxiliary vector of the eigenproblem .
Teuchos::RCP< const OP > getA() const
Get a pointer to the operator A of the eigenproblem .
bool _isSym
Symmetry of the eigenvalue problem.
int _nev
Number of eigenvalues requested.
Teuchos::RCP< MV > _InitVec
Reference-counted pointer for the initial vector of the eigenproblem .
Teuchos::RCP< const OP > _Op
Reference-counted pointer for the operator of the eigenproblem .
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
void setAuxVecs(const Teuchos::RCP< const MV > &AuxVecs)
Set auxiliary vectors.
int numVecs
The number of computed eigenpairs.
Teuchos::RCP< const OP > getOperator() const
Get a pointer to the operator for which eigenvalues will be computed.
Teuchos::RCP< const OP > getM() const
Get a pointer to the operator M of the eigenproblem .
Teuchos::RCP< const MV > computeCurrResVec() const
Returns the residual vector corresponding to the computed solution.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
const Eigensolution< ScalarType, MV > & getSolution() const
Get the solution to the eigenproblem.
void setSolution(const Eigensolution< ScalarType, MV > &sol)
Set the solution to the eigenproblem.
virtual ~BasicEigenproblem()
Destructor.
Teuchos::RCP< const OP > _MOp
Reference-counted pointer for M of the eigenproblem .
Teuchos::RCP< const MV > getInitVec() const
Get a pointer to the initial vector.
Struct for storing an eigenproblem solution.
void setPrec(const Teuchos::RCP< const OP > &Prec)
Set the preconditioner for this eigenvalue problem .
bool isProblemSet() const
If the problem has been set, this method will return true.
void setHermitian(bool isSym)
Specify the symmetry of this eigenproblem.
This provides a basic implementation for defining standard or generalized eigenvalue problems...
Teuchos::RCP< const MV > getAuxVecs() const
Get a pointer to the auxiliary vector.
Teuchos::RCP< const OP > getPrec() const
Get a pointer to the preconditioner of the eigenproblem .
Teuchos::RCP< const OP > _Prec
Reference-counted pointer for the preconditioner of the eigenproblem .
void setM(const Teuchos::RCP< const OP > &M)
Set the operator M of the eigenvalue problem .
BasicEigenproblem()
Empty constructor - allows Anasazi::BasicEigenproblem to be described at a later time through "Set Me...
void setA(const Teuchos::RCP< const OP > &A)
Set the operator A of the eigenvalue problem .
Teuchos::RCP< const OP > _AOp
Reference-counted pointer for A of the eigenproblem .
void setNEV(int nev)
Specify the number of eigenvalues (NEV) that are requested.
Eigensolution< ScalarType, MV > _sol
Solution to problem.
void setOperator(const Teuchos::RCP< const OP > &Op)
Set the operator for which eigenvalues will be computed.