Anasazi  Version of the Day
AnasaziMatOrthoManager.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 
34 #ifndef ANASAZI_MATORTHOMANAGER_HPP
35 #define ANASAZI_MATORTHOMANAGER_HPP
36 
54 #include "AnasaziConfigDefs.hpp"
55 #include "AnasaziTypes.hpp"
56 #include "AnasaziOrthoManager.hpp"
59 
60 namespace Anasazi {
61 
62  template <class ScalarType, class MV, class OP>
63  class MatOrthoManager : public OrthoManager<ScalarType,MV> {
64  public:
66 
67  MatOrthoManager(Teuchos::RCP<const OP> Op = Teuchos::null);
69 
71  virtual ~MatOrthoManager() {};
73 
75 
76 
78  virtual void setOp( Teuchos::RCP<const OP> Op );
79 
81  virtual Teuchos::RCP<const OP> getOp() const;
82 
84 
89  int getOpCounter() const;
90 
92 
94  void resetOpCounter();
95 
97 
99 
100 
110  void innerProdMat(
111  const MV& X, const MV& Y,
113  Teuchos::RCP<const MV> MX = Teuchos::null,
114  Teuchos::RCP<const MV> MY = Teuchos::null
115  ) const;
116 
125  void normMat(
126  const MV& X,
127  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
128  Teuchos::RCP<const MV> MX = Teuchos::null
129  ) const;
130 
141  virtual void projectMat (
142  MV &X,
145  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
146  Teuchos::RCP<MV> MX = Teuchos::null,
148  = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
149  ) const = 0;
150 
163  virtual int normalizeMat (
164  MV &X,
166  Teuchos::RCP<MV> MX = Teuchos::null
167  ) const = 0;
168 
169 
184  virtual int projectAndNormalizeMat (
185  MV &X,
188  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
190  Teuchos::RCP<MV> MX = Teuchos::null,
192  = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
193  ) const = 0;
194 
200  orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const = 0;
201 
208  const MV &X,
209  const MV &Y,
210  Teuchos::RCP<const MV> MX = Teuchos::null,
211  Teuchos::RCP<const MV> MY = Teuchos::null
212  ) const = 0;
213 
215 
217 
218 
226  void innerProd( const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const;
227 
235  void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const;
236 
244  void project (
245  MV &X,
248  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null))
249  ) const;
250 
258  int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null) const;
259 
267  int projectAndNormalize (
268  MV &X,
271  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
273  ) const;
274 
283  orthonormError(const MV &X) const;
284 
293  orthogError(const MV &X1, const MV &X2) const;
294 
296 
297  protected:
299  bool _hasOp;
300  mutable int _OpCounter;
301 
302  };
303 
304  template <class ScalarType, class MV, class OP>
306  : _Op(Op), _hasOp(Op!=Teuchos::null), _OpCounter(0) {}
307 
308  template <class ScalarType, class MV, class OP>
310  {
311  _Op = Op;
312  _hasOp = (_Op != Teuchos::null);
313  }
314 
315  template <class ScalarType, class MV, class OP>
317  {
318  return _Op;
319  }
320 
321  template <class ScalarType, class MV, class OP>
323  {
324  return _OpCounter;
325  }
326 
327  template <class ScalarType, class MV, class OP>
329  {
330  _OpCounter = 0;
331  }
332 
333  template <class ScalarType, class MV, class OP>
335  const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const
336  {
338  typedef MultiVecTraits<ScalarType,MV> MVT;
340 
343 
344  if (_hasOp) {
345  // attempt to minimize the amount of work in applying
346  if ( MVT::GetNumberVecs(X) < MVT::GetNumberVecs(Y) ) {
347  R = MVT::Clone(X,MVT::GetNumberVecs(X));
348  OPT::Apply(*_Op,X,*R);
349  _OpCounter += MVT::GetNumberVecs(X);
350  P = R;
351  Q = Teuchos::rcpFromRef(Y);
352  }
353  else {
354  P = Teuchos::rcpFromRef(X);
355  R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
356  OPT::Apply(*_Op,Y,*R);
357  _OpCounter += MVT::GetNumberVecs(Y);
358  Q = R;
359  }
360  }
361  else {
362  P = Teuchos::rcpFromRef(X);
363  Q = Teuchos::rcpFromRef(Y);
364  }
365 
366  MVT::MvTransMv(SCT::one(),*P,*Q,Z);
367  }
368 
369  template <class ScalarType, class MV, class OP>
372  {
373  (void) MX;
375  typedef MultiVecTraits<ScalarType,MV> MVT;
376  // typedef OperatorTraits<ScalarType,MV,OP> OPT; // unused
377 
378  Teuchos::RCP<MV> P,Q;
379 
380  if ( MY == Teuchos::null ) {
381  innerProd(X,Y,Z);
382  }
383  else if ( _hasOp ) {
384  // the user has done the matrix vector for us
385  MVT::MvTransMv(SCT::one(),X,*MY,Z);
386  }
387  else {
388  // there is no matrix vector
389  MVT::MvTransMv(SCT::one(),X,Y,Z);
390  }
391 #ifdef TEUCHOS_DEBUG
392  for (int j=0; j<Z.numCols(); j++) {
393  for (int i=0; i<Z.numRows(); i++) {
394  TEUCHOS_TEST_FOR_EXCEPTION(SCT::isnaninf(Z(i,j)), std::logic_error,
395  "Anasazi::MatOrthoManager::innerProdMat(): detected NaN/inf.");
396  }
397  }
398 #endif
399  }
400 
401  template <class ScalarType, class MV, class OP>
403  const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const
404  {
405  this->normMat(X,normvec);
406  }
407 
408  template <class ScalarType, class MV, class OP>
410  const MV& X,
411  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
412  Teuchos::RCP<const MV> MX) const
413  {
416  typedef MultiVecTraits<ScalarType,MV> MVT;
418 
419  int nvecs = MVT::GetNumberVecs(X);
420 
421  // Make sure that normvec has enough entries to hold the norms
422  // of all the columns of X. std::vector<T>::size_type is
423  // unsigned, so do the appropriate cast to avoid signed/unsigned
424  // comparisons that trigger compiler warnings.
425  if (normvec.size() < static_cast<size_t>(nvecs))
426  normvec.resize (nvecs);
427 
428  if (!_hasOp) {
429  // X == MX, since the operator M is the identity.
430  MX = Teuchos::rcp(&X, false);
431  MVT::MvNorm(X, normvec);
432  }
433  else {
434  // The caller didn't give us a previously computed MX, so
435  // apply the operator. We assign to MX only after applying
436  // the operator, so that if the application fails, MX won't be
437  // modified.
438  if(MX == Teuchos::null) {
439  Teuchos::RCP<MV> tempVec = MVT::Clone(X,nvecs);
440  OPT::Apply(*_Op,X,*tempVec);
441  _OpCounter += nvecs;
442  MX = tempVec;
443  }
444  else {
445  // The caller gave us a previously computed MX. Make sure
446  // that it has at least as many columns as X.
447  const int numColsMX = MVT::GetNumberVecs(*MX);
448  TEUCHOS_TEST_FOR_EXCEPTION(numColsMX < nvecs, std::invalid_argument,
449  "MatOrthoManager::norm(X, MX, normvec): "
450  "MX has fewer columns than X: "
451  "MX has " << numColsMX << " columns, "
452  "and X has " << nvecs << " columns.");
453  }
454 
455  std::vector<ScalarType> dotvec(nvecs);
456  MVT::MvDot(X,*MX,dotvec);
457  for (int i=0; i<nvecs; i++) {
458  normvec[i] = MT::squareroot( SCT::magnitude(dotvec[i]) );
459  }
460  }
461  }
462 
463  template <class ScalarType, class MV, class OP>
465  MV &X,
468  ) const
469  {
470  this->projectMat(X,Q,C);
471  }
472 
473  template <class ScalarType, class MV, class OP>
476  {
477  return this->normalizeMat(X,B);
478  }
479 
480  template <class ScalarType, class MV, class OP>
482  MV &X,
486  ) const
487  {
488  return this->projectAndNormalizeMat(X,Q,C,B);
489  }
490 
491  template <class ScalarType, class MV, class OP>
494  {
495  return this->orthonormErrorMat(X,Teuchos::null);
496  }
497 
498  template <class ScalarType, class MV, class OP>
500  MatOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, const MV &X2) const
501  {
502  return this->orthogErrorMat(X1,X2);
503  }
504 
505 } // end of Anasazi namespace
506 
507 
508 #endif
509 
510 // end of file AnasaziMatOrthoManager.hpp
void project(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null))) const
Implements the interface OrthoManager::project().
virtual void setOp(Teuchos::RCP< const OP > Op)
Set operator used for inner product.
Declaration of basic traits for the multivector type.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void norm(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec) const
Implements the interface OrthoManager::norm().
Virtual base class which defines basic traits for the operator type.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogError(const MV &X1, const MV &X2) const
Implements the interface OrthoManager::orthogError().
virtual void projectMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const =0
Provides matrix-based projection method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null) const
Implements the interface OrthoManager::normalize().
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
OrdinalType numRows() const
virtual Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormErrorMat(const MV &X, Teuchos::RCP< const MV > MX=Teuchos::null) const =0
This method computes the error in orthonormality of a multivector.
MatOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null)
Default constructor.
Templated virtual class for providing orthogonalization/orthonormalization methods.
void normMat(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, Teuchos::RCP< const MV > MX=Teuchos::null) const
Provides the norm induced by the matrix-based inner product.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
virtual ~MatOrthoManager()
Destructor.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
void innerProdMat(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
Provides a matrix-based inner product.
virtual Teuchos::RCP< const OP > getOp() const
Get operator used for inner product.
void resetOpCounter()
Reset the operator counter to zero.
int projectAndNormalize(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null) const
Implements the interface OrthoManager::projectAndNormalize().
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
Implements the interface OrthoManager::orthonormError().
virtual int normalizeMat(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null) const =0
Provides matrix-based orthonormalization method.
int getOpCounter() const
Retrieve operator counter.
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
OrdinalType numCols() const
virtual Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogErrorMat(const MV &X, const MV &Y, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const =0
This method computes the error in orthogonality of two multivectors.
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Implements the interface OrthoManager::innerProd().
virtual int projectAndNormalizeMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const =0
Provides matrix-based projection/orthonormalization method.