47#ifndef ANASAZI_MATORTHOMANAGER_HPP
48#define ANASAZI_MATORTHOMANAGER_HPP
75 template <
class ScalarType,
class MV,
class OP>
91 virtual void setOp( Teuchos::RCP<const OP> Op );
94 virtual Teuchos::RCP<const OP>
getOp()
const;
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
140 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
141 Teuchos::RCP<const MV> MX = Teuchos::null
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))
178 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
179 Teuchos::RCP<MV> MX = Teuchos::null
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))
212 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
219 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
223 Teuchos::RCP<const MV> MX = Teuchos::null,
224 Teuchos::RCP<const MV> MY = Teuchos::null
239 void innerProd(
const MV& X,
const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z )
const;
248 void norm(
const MV& X, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec )
const;
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))
271 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null)
const;
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
295 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
305 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
311 Teuchos::RCP<const OP> _Op;
313 mutable int _OpCounter;
317 template <
class ScalarType,
class MV,
class OP>
319 : _Op(Op), _hasOp(Op!=Teuchos::null), _OpCounter(0) {}
321 template <
class ScalarType,
class MV,
class OP>
325 _hasOp = (_Op != Teuchos::null);
328 template <
class ScalarType,
class MV,
class OP>
334 template <
class ScalarType,
class MV,
class OP>
340 template <
class ScalarType,
class MV,
class OP>
346 template <
class ScalarType,
class MV,
class OP>
348 const MV& X,
const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z )
const
350 typedef Teuchos::ScalarTraits<ScalarType> SCT;
351 typedef MultiVecTraits<ScalarType,MV> MVT;
352 typedef OperatorTraits<ScalarType,MV,OP> OPT;
354 Teuchos::RCP<const MV> P,Q;
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);
364 Q = Teuchos::rcpFromRef(Y);
367 P = Teuchos::rcpFromRef(X);
368 R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
369 OPT::Apply(*_Op,Y,*R);
370 _OpCounter += MVT::GetNumberVecs(Y);
375 P = Teuchos::rcpFromRef(X);
376 Q = Teuchos::rcpFromRef(Y);
379 MVT::MvTransMv(SCT::one(),*P,*Q,Z);
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
387 typedef Teuchos::ScalarTraits<ScalarType> SCT;
388 typedef MultiVecTraits<ScalarType,MV> MVT;
391 Teuchos::RCP<MV> P,Q;
393 if ( MY == Teuchos::null ) {
398 MVT::MvTransMv(SCT::one(),X,*MY,Z);
402 MVT::MvTransMv(SCT::one(),X,Y,Z);
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.");
414 template <
class ScalarType,
class MV,
class OP>
416 const MV& X, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec )
const
418 this->normMat(X,normvec);
421 template <
class ScalarType,
class MV,
class OP>
424 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
425 Teuchos::RCP<const MV> MX)
const
427 typedef Teuchos::ScalarTraits<ScalarType> SCT;
428 typedef Teuchos::ScalarTraits<typename SCT::magnitudeType> MT;
429 typedef MultiVecTraits<ScalarType,MV> MVT;
430 typedef OperatorTraits<ScalarType,MV,OP> OPT;
432 int nvecs = MVT::GetNumberVecs(X);
438 if (normvec.size() <
static_cast<size_t>(nvecs))
439 normvec.resize (nvecs);
443 MX = Teuchos::rcp(&X,
false);
444 MVT::MvNorm(X, normvec);
451 if(MX == Teuchos::null) {
452 Teuchos::RCP<MV> tempVec = MVT::Clone(X,nvecs);
453 OPT::Apply(*_Op,X,*tempVec);
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.");
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]) );
476 template <
class ScalarType,
class MV,
class OP>
479 Teuchos::Array<Teuchos::RCP<const MV> > Q,
480 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
483 this->projectMat(X,Q,C);
486 template <
class ScalarType,
class MV,
class OP>
488 MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const
490 return this->normalizeMat(X,B);
493 template <
class ScalarType,
class MV,
class OP>
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
501 return this->projectAndNormalizeMat(X,Q,C,B);
504 template <
class ScalarType,
class MV,
class OP>
505 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
508 return this->orthonormErrorMat(X,Teuchos::null);
511 template <
class ScalarType,
class MV,
class OP>
512 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
515 return this->orthogErrorMat(X1,X2);
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.
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.