Anasazi Version of the Day
Loading...
Searching...
No Matches
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"
72
73namespace Anasazi {
74
75 template <class ScalarType, class MV, class OP>
76 class MatOrthoManager : public OrthoManager<ScalarType,MV> {
77 public:
79
80
81 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
108
110
112
113
124 const MV& X, const MV& Y,
125 Teuchos::SerialDenseMatrix<int,ScalarType>& Z,
126 Teuchos::RCP<const MV> MX = Teuchos::null,
127 Teuchos::RCP<const MV> MY = Teuchos::null
128 ) const;
129
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,
156 Teuchos::Array<Teuchos::RCP<const MV> > Q,
157 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
158 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
159 Teuchos::RCP<MV> MX = Teuchos::null,
160 Teuchos::Array<Teuchos::RCP<const MV> > MQ
161 = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
162 ) const = 0;
163
176 virtual int normalizeMat (
177 MV &X,
178 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
179 Teuchos::RCP<MV> MX = Teuchos::null
180 ) const = 0;
181
182
198 MV &X,
199 Teuchos::Array<Teuchos::RCP<const MV> > Q,
200 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
201 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
202 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
203 Teuchos::RCP<MV> MX = Teuchos::null,
204 Teuchos::Array<Teuchos::RCP<const MV> > MQ
205 = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
206 ) const = 0;
207
212 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
213 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const = 0;
214
219 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
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,
259 Teuchos::Array<Teuchos::RCP<const MV> > Q,
260 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
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
281 MV &X,
282 Teuchos::Array<Teuchos::RCP<const MV> > Q,
283 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
284 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
285 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null
286 ) const;
287
295 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
296 orthonormError(const MV &X) const;
297
305 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
306 orthogError(const MV &X1, const MV &X2) const;
307
309
310 protected:
311 Teuchos::RCP<const OP> _Op;
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>
322 void MatOrthoManager<ScalarType,MV,OP>::setOp( Teuchos::RCP<const OP> Op )
323 {
324 _Op = Op;
325 _hasOp = (_Op != Teuchos::null);
326 }
327
328 template <class ScalarType, class MV, class OP>
329 Teuchos::RCP<const OP> MatOrthoManager<ScalarType,MV,OP>::getOp() const
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 {
350 typedef Teuchos::ScalarTraits<ScalarType> SCT;
353
354 Teuchos::RCP<const MV> P,Q;
355 Teuchos::RCP<MV> R;
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>
384 const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z, Teuchos::RCP<const MV> MX, Teuchos::RCP<const MV> MY) const
385 {
386 (void) MX;
387 typedef Teuchos::ScalarTraits<ScalarType> SCT;
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 {
427 typedef Teuchos::ScalarTraits<ScalarType> SCT;
428 typedef Teuchos::ScalarTraits<typename SCT::magnitudeType> MT;
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,
479 Teuchos::Array<Teuchos::RCP<const MV> > Q,
480 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
481 ) const
482 {
483 this->projectMat(X,Q,C);
484 }
485
486 template <class ScalarType, class MV, class OP>
488 MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const
489 {
490 return this->normalizeMat(X,B);
491 }
492
493 template <class ScalarType, class MV, class OP>
495 MV &X,
496 Teuchos::Array<Teuchos::RCP<const MV> > Q,
497 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
498 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B
499 ) const
500 {
501 return this->projectAndNormalizeMat(X,Q,C,B);
502 }
503
504 template <class ScalarType, class MV, class OP>
505 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
507 {
508 return this->orthonormErrorMat(X,Teuchos::null);
509 }
510
511 template <class ScalarType, class MV, class OP>
512 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
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
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
MatOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null)
Default constructor.
virtual Teuchos::RCP< const OP > getOp() const
Get operator used for inner product.
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.
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Implements the interface OrthoManager::innerProd().
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().
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().
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.
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.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
Implements the interface OrthoManager::orthonormError().
void norm(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec) const
Implements the interface OrthoManager::norm().
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 void setOp(Teuchos::RCP< const OP > Op)
Set operator used for inner product.
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 ~MatOrthoManager()
Destructor.
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.
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 orthogError(const MV &X1, const MV &X2) const
Implements the interface OrthoManager::orthogError().
int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null) const
Implements the interface OrthoManager::normalize().
int getOpCounter() const
Retrieve operator counter.
void resetOpCounter()
Reset the operator counter to zero.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.