Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 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 
47 #ifndef ANASAZI_MATORTHOMANAGER_HPP
48 #define ANASAZI_MATORTHOMANAGER_HPP
49 
67 #include "AnasaziConfigDefs.hpp"
68 #include "AnasaziTypes.hpp"
69 #include "AnasaziOrthoManager.hpp"
72 
73 namespace Anasazi {
74 
75  template <class ScalarType, class MV, class OP>
76  class MatOrthoManager : public OrthoManager<ScalarType,MV> {
77  public:
79 
80  MatOrthoManager(Teuchos::RCP<const OP> Op = Teuchos::null);
82 
84  virtual ~MatOrthoManager() {};
86 
88 
89 
91  virtual void setOp( Teuchos::RCP<const OP> Op );
92 
94  virtual Teuchos::RCP<const OP> getOp() const;
95 
97 
102  int getOpCounter() const;
103 
105 
107  void resetOpCounter();
108 
110 
112 
113 
123  void innerProdMat(
124  const MV& X, const MV& Y,
126  Teuchos::RCP<const MV> MX = Teuchos::null,
127  Teuchos::RCP<const MV> MY = Teuchos::null
128  ) const;
129 
138  void normMat(
139  const MV& X,
140  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
141  Teuchos::RCP<const MV> MX = Teuchos::null
142  ) const;
143 
154  virtual void projectMat (
155  MV &X,
158  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
159  Teuchos::RCP<MV> MX = Teuchos::null,
161  = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
162  ) const = 0;
163 
176  virtual int normalizeMat (
177  MV &X,
179  Teuchos::RCP<MV> MX = Teuchos::null
180  ) const = 0;
181 
182 
197  virtual int projectAndNormalizeMat (
198  MV &X,
201  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
203  Teuchos::RCP<MV> MX = Teuchos::null,
205  = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
206  ) const = 0;
207 
213  orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const = 0;
214 
221  const MV &X,
222  const MV &Y,
223  Teuchos::RCP<const MV> MX = Teuchos::null,
224  Teuchos::RCP<const MV> MY = Teuchos::null
225  ) const = 0;
226 
228 
230 
231 
239  void innerProd( const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const;
240 
248  void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const;
249 
257  void project (
258  MV &X,
261  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null))
262  ) const;
263 
271  int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null) const;
272 
280  int projectAndNormalize (
281  MV &X,
284  = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
286  ) const;
287 
296  orthonormError(const MV &X) const;
297 
306  orthogError(const MV &X1, const MV &X2) const;
307 
309 
310  protected:
312  bool _hasOp;
313  mutable int _OpCounter;
314 
315  };
316 
317  template <class ScalarType, class MV, class OP>
319  : _Op(Op), _hasOp(Op!=Teuchos::null), _OpCounter(0) {}
320 
321  template <class ScalarType, class MV, class OP>
323  {
324  _Op = Op;
325  _hasOp = (_Op != Teuchos::null);
326  }
327 
328  template <class ScalarType, class MV, class OP>
330  {
331  return _Op;
332  }
333 
334  template <class ScalarType, class MV, class OP>
336  {
337  return _OpCounter;
338  }
339 
340  template <class ScalarType, class MV, class OP>
342  {
343  _OpCounter = 0;
344  }
345 
346  template <class ScalarType, class MV, class OP>
348  const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const
349  {
351  typedef MultiVecTraits<ScalarType,MV> MVT;
353 
356 
357  if (_hasOp) {
358  // attempt to minimize the amount of work in applying
359  if ( MVT::GetNumberVecs(X) < MVT::GetNumberVecs(Y) ) {
360  R = MVT::Clone(X,MVT::GetNumberVecs(X));
361  OPT::Apply(*_Op,X,*R);
362  _OpCounter += MVT::GetNumberVecs(X);
363  P = R;
364  Q = Teuchos::rcpFromRef(Y);
365  }
366  else {
367  P = Teuchos::rcpFromRef(X);
368  R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
369  OPT::Apply(*_Op,Y,*R);
370  _OpCounter += MVT::GetNumberVecs(Y);
371  Q = R;
372  }
373  }
374  else {
375  P = Teuchos::rcpFromRef(X);
376  Q = Teuchos::rcpFromRef(Y);
377  }
378 
379  MVT::MvTransMv(SCT::one(),*P,*Q,Z);
380  }
381 
382  template <class ScalarType, class MV, class OP>
385  {
386  (void) MX;
388  typedef MultiVecTraits<ScalarType,MV> MVT;
389  // typedef OperatorTraits<ScalarType,MV,OP> OPT; // unused
390 
391  Teuchos::RCP<MV> P,Q;
392 
393  if ( MY == Teuchos::null ) {
394  innerProd(X,Y,Z);
395  }
396  else if ( _hasOp ) {
397  // the user has done the matrix vector for us
398  MVT::MvTransMv(SCT::one(),X,*MY,Z);
399  }
400  else {
401  // there is no matrix vector
402  MVT::MvTransMv(SCT::one(),X,Y,Z);
403  }
404 #ifdef TEUCHOS_DEBUG
405  for (int j=0; j<Z.numCols(); j++) {
406  for (int i=0; i<Z.numRows(); i++) {
407  TEUCHOS_TEST_FOR_EXCEPTION(SCT::isnaninf(Z(i,j)), std::logic_error,
408  "Anasazi::MatOrthoManager::innerProdMat(): detected NaN/inf.");
409  }
410  }
411 #endif
412  }
413 
414  template <class ScalarType, class MV, class OP>
416  const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const
417  {
418  this->normMat(X,normvec);
419  }
420 
421  template <class ScalarType, class MV, class OP>
423  const MV& X,
424  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
425  Teuchos::RCP<const MV> MX) const
426  {
429  typedef MultiVecTraits<ScalarType,MV> MVT;
431 
432  int nvecs = MVT::GetNumberVecs(X);
433 
434  // Make sure that normvec has enough entries to hold the norms
435  // of all the columns of X. std::vector<T>::size_type is
436  // unsigned, so do the appropriate cast to avoid signed/unsigned
437  // comparisons that trigger compiler warnings.
438  if (normvec.size() < static_cast<size_t>(nvecs))
439  normvec.resize (nvecs);
440 
441  if (!_hasOp) {
442  // X == MX, since the operator M is the identity.
443  MX = Teuchos::rcp(&X, false);
444  MVT::MvNorm(X, normvec);
445  }
446  else {
447  // The caller didn't give us a previously computed MX, so
448  // apply the operator. We assign to MX only after applying
449  // the operator, so that if the application fails, MX won't be
450  // modified.
451  if(MX == Teuchos::null) {
452  Teuchos::RCP<MV> tempVec = MVT::Clone(X,nvecs);
453  OPT::Apply(*_Op,X,*tempVec);
454  _OpCounter += nvecs;
455  MX = tempVec;
456  }
457  else {
458  // The caller gave us a previously computed MX. Make sure
459  // that it has at least as many columns as X.
460  const int numColsMX = MVT::GetNumberVecs(*MX);
461  TEUCHOS_TEST_FOR_EXCEPTION(numColsMX < nvecs, std::invalid_argument,
462  "MatOrthoManager::norm(X, MX, normvec): "
463  "MX has fewer columns than X: "
464  "MX has " << numColsMX << " columns, "
465  "and X has " << nvecs << " columns.");
466  }
467 
468  std::vector<ScalarType> dotvec(nvecs);
469  MVT::MvDot(X,*MX,dotvec);
470  for (int i=0; i<nvecs; i++) {
471  normvec[i] = MT::squareroot( SCT::magnitude(dotvec[i]) );
472  }
473  }
474  }
475 
476  template <class ScalarType, class MV, class OP>
478  MV &X,
481  ) const
482  {
483  this->projectMat(X,Q,C);
484  }
485 
486  template <class ScalarType, class MV, class OP>
489  {
490  return this->normalizeMat(X,B);
491  }
492 
493  template <class ScalarType, class MV, class OP>
495  MV &X,
499  ) const
500  {
501  return this->projectAndNormalizeMat(X,Q,C,B);
502  }
503 
504  template <class ScalarType, class MV, class OP>
507  {
508  return this->orthonormErrorMat(X,Teuchos::null);
509  }
510 
511  template <class ScalarType, class MV, class OP>
513  MatOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, const MV &X2) const
514  {
515  return this->orthogErrorMat(X1,X2);
516  }
517 
518 } // end of Anasazi namespace
519 
520 
521 #endif
522 
523 // end of file AnasaziMatOrthoManager.hpp
void norm(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec) const
Implements the interface OrthoManager::norm().
virtual void setOp(Teuchos::RCP< const OP > Op)
Set operator used for inner product.
Declaration of basic traits for the multivector type.
int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null) const
Implements the interface OrthoManager::normalize().
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
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.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
Implements the interface OrthoManager::orthonormError().
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.
virtual Teuchos::RCP< const OP > getOp() const
Get operator used for inner product.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Implements the interface OrthoManager::innerProd().
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
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 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().
int getOpCounter() const
Retrieve operator counter.
OrdinalType numCols() const
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogError(const MV &X1, const MV &X2) const
Implements the interface OrthoManager::orthogError().
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.
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi&#39;s templated virtual class for providing routines for orthogonalization and orthonormalization...
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().
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.
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.
OrdinalType numRows() const
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.