Anasazi Version of the Day
Loading...
Searching...
No Matches
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 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
42#ifndef ANASAZI_BASIC_EIGENPROBLEM_H
43#define ANASAZI_BASIC_EIGENPROBLEM_H
44
52
58namespace Anasazi {
59
60 template<class ScalarType, class MV, class OP>
61 class BasicEigenproblem : public virtual Eigenproblem<ScalarType, MV, OP> {
62
63 public:
64
66
67
70
72 BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<MV>& InitVec );
73
75 BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<const OP>& B, const Teuchos::RCP<MV>& InitVec );
76
79
81 virtual ~BasicEigenproblem() {};
83
85
86
93 void setOperator( const Teuchos::RCP<const OP>& Op ) { _Op = Op; _isSet=false; };
94
97 void setA( const Teuchos::RCP<const OP>& A ) { _AOp = A; _isSet=false; };
98
101 void setM( const Teuchos::RCP<const OP>& M ) { _MOp = M; _isSet=false; };
102
105 void setPrec( const Teuchos::RCP<const OP>& Prec ) { _Prec = Prec; _isSet=false; };
106
114 void setInitVec( const Teuchos::RCP<MV>& InitVec ) { _InitVec = InitVec; _isSet=false; };
115
121 void setAuxVecs( const Teuchos::RCP<const MV>& AuxVecs ) { _AuxVecs = AuxVecs; _isSet=false; };
122
124 void setNEV( int nev ){ _nev = nev; _isSet=false; };
125
127
130 void setHermitian( bool isSym ){ _isSym = isSym; _isSet=false; };
131
147 bool setProblem();
148
156
158
160
161
163 Teuchos::RCP<const OP> getOperator() const { return( _Op ); };
164
166 Teuchos::RCP<const OP> getA() const { return( _AOp ); };
167
169 Teuchos::RCP<const OP> getM() const { return( _MOp ); };
170
172 Teuchos::RCP<const OP> getPrec() const { return( _Prec ); };
173
175 Teuchos::RCP<const MV> getInitVec() const { return( _InitVec ); };
176
178 Teuchos::RCP<const MV> getAuxVecs() const { return( _AuxVecs ); };
179
181 int getNEV() const { return( _nev ); }
182
184 bool isHermitian() const { return( _isSym ); }
185
187 bool isProblemSet() const { return( _isSet ); }
188
194 const Eigensolution<ScalarType,MV> & getSolution() const { return(_sol); }
195
197
199
200
208 Teuchos::RCP<const MV> computeCurrResVec() const;
209
211
212 protected:
213
215 Teuchos::RCP<const OP> _AOp;
216
218 Teuchos::RCP<const OP> _MOp;
219
221 Teuchos::RCP<const OP> _Op;
222
224 Teuchos::RCP<const OP> _Prec;
225
227 Teuchos::RCP<MV> _InitVec;
228
230 Teuchos::RCP<const MV> _AuxVecs;
231
233 int _nev;
234
236
239 bool _isSym;
240
242 bool _isSet;
243
248
251 };
252
253
254 //=============================================================================
255 // Implementations (Constructors / Destructors)
256 //=============================================================================
257 template <class ScalarType, class MV, class OP>
259 _nev(0),
260 _isSym(false),
261 _isSet(false)
262 {
263 }
264
265
266 //=============================================================================
267 template <class ScalarType, class MV, class OP>
268 BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<MV>& InitVec ) :
269 _Op(Op),
270 _InitVec(InitVec),
271 _nev(0),
272 _isSym(false),
273 _isSet(false)
274 {
275 }
276
277
278 //=============================================================================
279 template <class ScalarType, class MV, class OP>
280 BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<const OP>& M,
281 const Teuchos::RCP<MV>& InitVec ) :
282 _MOp(M),
283 _Op(Op),
284 _InitVec(InitVec),
285 _nev(0),
286 _isSym(false),
287 _isSet(false)
288 {
289 }
290
291
292 //=============================================================================
293 template <class ScalarType, class MV, class OP>
295 _AOp(Problem._AOp),
296 _MOp(Problem._MOp),
297 _Op(Problem._Op),
298 _Prec(Problem._Prec),
299 _InitVec(Problem._InitVec),
300 _nev(Problem._nev),
301 _isSym(Problem._isSym),
302 _isSet(Problem._isSet),
303 _sol(Problem._sol)
304 {
305 }
306
307
308 //=============================================================================
309 // SetProblem (sanity check method)
310 //=============================================================================
311 template <class ScalarType, class MV, class OP>
313 {
314 //----------------------------------------------------------------
315 // Sanity Checks
316 //----------------------------------------------------------------
317 // If there is no operator, then we can't proceed.
318 if ( !_AOp.get() && !_Op.get() ) { return false; }
319
320 // If there is no initial vector, then we don't have anything to clone workspace from.
321 if ( !_InitVec.get() ) { return false; }
322
323 // If we don't need any eigenvalues, we don't need to continue.
324 if (_nev == 0) { return false; }
325
326 // If there is an A, but no operator, we can set them equal.
327 if (_AOp.get() && !_Op.get()) { _Op = _AOp; }
328
329 // Clear the storage from any previous call to setSolution()
331 _sol = emptysol;
332
333 // mark the problem as set and return no-error
334 _isSet=true;
335 return true;
336 }
337
338 //=============================================================================
339 // Computes the residual vector
340 //=============================================================================
341 template <class ScalarType, class MV, class OP>
343 {
344 using Teuchos::RCP;
345
346 TEUCHOS_TEST_FOR_EXCEPTION(!isHermitian(), std::invalid_argument,
347 "BasicEigenproblem::computeCurrResVec: This method is not currently "
348 "implemented for non-Hermitian problems. Sorry for any inconvenience.");
349
350 const Eigensolution<ScalarType,MV> sol = getSolution();
351
352 if(sol.numVecs <= 0)
353 return Teuchos::null;
354
355 // Copy the eigenvalues
356 RCP<MV> X = sol.Evecs;
357 std::vector<ScalarType> Lambda(sol.numVecs);
358 for(int i = 0; i < sol.numVecs; i++)
359 Lambda[i] = sol.Evals[i].realpart;
360
361 // Compute AX
362 RCP<MV> AX = MVT::Clone(*X,sol.numVecs);
363 if(_AOp != Teuchos::null)
364 OPT::Apply(*_AOp,*X,*AX);
365 else
366 OPT::Apply(*_Op,*X,*AX);
367
368 // Compute MX if necessary
369 RCP<MV> MX;
370 if(_MOp != Teuchos::null)
371 {
372 MX = MVT::Clone(*X,sol.numVecs);
373 OPT::Apply(*_MOp,*X,*MX);
374 }
375 else
376 {
377 MX = MVT::CloneCopy(*X);
378 }
379
380 // Compute R = AX - MX \Lambda
381 RCP<MV> R = MVT::Clone(*X,sol.numVecs);
382 MVT::MvScale(*MX,Lambda);
383 MVT::MvAddMv(1.0,*AX,-1.0,*MX,*R);
384
385 return R;
386 }
387
388} // end Anasazi namespace
389#endif
390
391// end AnasaziBasicEigenproblem.hpp
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
This provides a basic implementation for defining standard or generalized eigenvalue problems.
void setInitVec(const Teuchos::RCP< MV > &InitVec)
Set the initial guess.
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 setNEV(int nev)
Specify the number of eigenvalues (NEV) that are requested.
Teuchos::RCP< const MV > getInitVec() const
Get a pointer to the initial vector.
void setPrec(const Teuchos::RCP< const OP > &Prec)
Set the preconditioner for this eigenvalue problem .
void setHermitian(bool isSym)
Specify the symmetry of this eigenproblem.
void setA(const Teuchos::RCP< const OP > &A)
Set the operator A of the eigenvalue problem .
Eigensolution< ScalarType, MV > _sol
Solution to problem.
bool setProblem()
Specify that this eigenproblem is fully defined.
void setOperator(const Teuchos::RCP< const OP > &Op)
Set the operator for which eigenvalues will be computed.
Teuchos::RCP< const MV > _AuxVecs
Reference-counted pointer for the auxiliary vector of the eigenproblem .
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 .
MultiVecTraits< ScalarType, MV > MVT
Type-definition for the MultiVecTraits class corresponding to the MV type.
Teuchos::RCP< const OP > getA() const
Get a pointer to the operator A of the eigenproblem .
Teuchos::RCP< const OP > _MOp
Reference-counted pointer for M of the eigenproblem .
Teuchos::RCP< const MV > getAuxVecs() const
Get a pointer to the auxiliary vector.
int _nev
Number of eigenvalues requested.
void setM(const Teuchos::RCP< const OP > &M)
Set the operator M of the eigenvalue problem .
void setSolution(const Eigensolution< ScalarType, MV > &sol)
Set the solution to the eigenproblem.
Teuchos::RCP< const MV > computeCurrResVec() const
Returns the residual vector corresponding to the computed solution.
Teuchos::RCP< MV > _InitVec
Reference-counted pointer for the initial vector of the eigenproblem .
Teuchos::RCP< const OP > _AOp
Reference-counted pointer for A of the eigenproblem .
bool isHermitian() const
Get the symmetry information for this eigenproblem.
const Eigensolution< ScalarType, MV > & getSolution() const
Get the solution to the eigenproblem.
bool _isSym
Symmetry of the eigenvalue problem.
bool isProblemSet() const
If the problem has been set, this method will return true.
OperatorTraits< ScalarType, MV, OP > OPT
Type-definition for the OperatorTraits class corresponding to the OP type.
void setAuxVecs(const Teuchos::RCP< const MV > &AuxVecs)
Set auxiliary vectors.
Teuchos::RCP< const OP > _Op
Reference-counted pointer for the operator of the eigenproblem .
BasicEigenproblem()
Empty constructor - allows Anasazi::BasicEigenproblem to be described at a later time through "Set Me...
int getNEV() const
Get the number of eigenvalues (NEV) that are required by this eigenproblem.
This class defines the interface required by an eigensolver and status test class to compute solution...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.