00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00034 #ifndef ANASAZI_MATORTHOMANAGER_HPP
00035 #define ANASAZI_MATORTHOMANAGER_HPP
00036
00054 #include "AnasaziConfigDefs.hpp"
00055 #include "AnasaziTypes.hpp"
00056 #include "AnasaziOrthoManager.hpp"
00057 #include "AnasaziMultiVecTraits.hpp"
00058 #include "AnasaziOperatorTraits.hpp"
00059
00060 namespace Anasazi {
00061
00062 template <class ScalarType, class MV, class OP>
00063 class MatOrthoManager : public OrthoManager<ScalarType,MV> {
00064 public:
00066
00067
00068 MatOrthoManager(Teuchos::RCP<const OP> Op = Teuchos::null);
00069
00071 virtual ~MatOrthoManager() {};
00073
00075
00076
00078 void setOp( Teuchos::RCP<const OP> Op );
00079
00081 Teuchos::RCP<const OP> getOp() const;
00082
00084
00089 int getOpCounter() const;
00090
00092
00094 void resetOpCounter();
00095
00097
00099
00100
00110 void innerProdMat(
00111 const MV& X, const MV& Y,
00112 Teuchos::SerialDenseMatrix<int,ScalarType>& Z,
00113 Teuchos::RCP<const MV> MX = Teuchos::null,
00114 Teuchos::RCP<const MV> MY = Teuchos::null
00115 ) const;
00116
00125 void normMat(
00126 const MV& X,
00127 std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
00128 Teuchos::RCP<const MV> MX = Teuchos::null
00129 ) const;
00130
00141 virtual void projectMat (
00142 MV &X,
00143 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00144 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00145 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00146 Teuchos::RCP<MV> MX = Teuchos::null,
00147 Teuchos::Array<Teuchos::RCP<const MV> > MQ
00148 = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00149 ) const = 0;
00150
00163 virtual int normalizeMat (
00164 MV &X,
00165 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00166 Teuchos::RCP<MV> MX = Teuchos::null
00167 ) const = 0;
00168
00169
00184 virtual int projectAndNormalizeMat (
00185 MV &X,
00186 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00187 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00188 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00189 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00190 Teuchos::RCP<MV> MX = Teuchos::null,
00191 Teuchos::Array<Teuchos::RCP<const MV> > MQ
00192 = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00193 ) const = 0;
00194
00199 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00200 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const = 0;
00201
00206 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00207 orthogErrorMat(
00208 const MV &X,
00209 const MV &Y,
00210 Teuchos::RCP<const MV> MX = Teuchos::null,
00211 Teuchos::RCP<const MV> MY = Teuchos::null
00212 ) const = 0;
00213
00215
00217
00218
00226 void innerProd( const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const;
00227
00235 void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const;
00236
00244 void project (
00245 MV &X,
00246 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00247 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00248 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null))
00249 ) const;
00250
00258 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null) const;
00259
00267 int projectAndNormalize (
00268 MV &X,
00269 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00270 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00271 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00272 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null
00273 ) const;
00274
00282 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00283 orthonormError(const MV &X) const;
00284
00292 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00293 orthogError(const MV &X1, const MV &X2) const;
00294
00296
00297 protected:
00298 Teuchos::RCP<const OP> _Op;
00299 bool _hasOp;
00300 mutable int _OpCounter;
00301
00302 };
00303
00304 template <class ScalarType, class MV, class OP>
00305 MatOrthoManager<ScalarType,MV,OP>::MatOrthoManager(Teuchos::RCP<const OP> Op)
00306 : _Op(Op), _hasOp(Op!=Teuchos::null), _OpCounter(0) {}
00307
00308 template <class ScalarType, class MV, class OP>
00309 void MatOrthoManager<ScalarType,MV,OP>::setOp( Teuchos::RCP<const OP> Op )
00310 {
00311 _Op = Op;
00312 _hasOp = (_Op != Teuchos::null);
00313 }
00314
00315 template <class ScalarType, class MV, class OP>
00316 Teuchos::RCP<const OP> MatOrthoManager<ScalarType,MV,OP>::getOp() const
00317 {
00318 return _Op;
00319 }
00320
00321 template <class ScalarType, class MV, class OP>
00322 int MatOrthoManager<ScalarType,MV,OP>::getOpCounter() const
00323 {
00324 return _OpCounter;
00325 }
00326
00327 template <class ScalarType, class MV, class OP>
00328 void MatOrthoManager<ScalarType,MV,OP>::resetOpCounter()
00329 {
00330 _OpCounter = 0;
00331 }
00332
00333 template <class ScalarType, class MV, class OP>
00334 void MatOrthoManager<ScalarType,MV,OP>::innerProd(
00335 const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const
00336 {
00337 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00338 typedef MultiVecTraits<ScalarType,MV> MVT;
00339 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00340
00341 Teuchos::RCP<const MV> P,Q;
00342 Teuchos::RCP<MV> R;
00343
00344 if (_hasOp) {
00345
00346 if ( MVT::GetNumberVecs(X) < MVT::GetNumberVecs(Y) ) {
00347 R = MVT::Clone(X,MVT::GetNumberVecs(X));
00348 OPT::Apply(*_Op,X,*R);
00349 _OpCounter += MVT::GetNumberVecs(X);
00350 P = R;
00351 Q = Teuchos::rcpFromRef(Y);
00352 }
00353 else {
00354 P = Teuchos::rcpFromRef(X);
00355 R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
00356 OPT::Apply(*_Op,Y,*R);
00357 _OpCounter += MVT::GetNumberVecs(Y);
00358 Q = R;
00359 }
00360 }
00361 else {
00362 P = Teuchos::rcpFromRef(X);
00363 Q = Teuchos::rcpFromRef(Y);
00364 }
00365
00366 MVT::MvTransMv(SCT::one(),*P,*Q,Z);
00367 }
00368
00369 template <class ScalarType, class MV, class OP>
00370 void MatOrthoManager<ScalarType,MV,OP>::innerProdMat(
00371 const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z, Teuchos::RCP<const MV> MX, Teuchos::RCP<const MV> MY) const
00372 {
00373 (void)MX;
00374 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00375 typedef MultiVecTraits<ScalarType,MV> MVT;
00376 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00377
00378 Teuchos::RCP<MV> P,Q;
00379
00380 if ( MY == Teuchos::null ) {
00381 innerProd(X,Y,Z);
00382 }
00383 else if ( _hasOp ) {
00384
00385 MVT::MvTransMv(SCT::one(),X,*MY,Z);
00386 }
00387 else {
00388
00389 MVT::MvTransMv(SCT::one(),X,Y,Z);
00390 }
00391 #ifdef TEUCHOS_DEBUG
00392 for (int j=0; j<Z.numCols(); j++) {
00393 for (int i=0; i<Z.numRows(); i++) {
00394 TEST_FOR_EXCEPTION(SCT::isnaninf(Z(i,j)), std::logic_error,
00395 "Anasazi::MatOrthoManager::innerProdMat(): detected NaN/inf.");
00396 }
00397 }
00398 #endif
00399 }
00400
00401 template <class ScalarType, class MV, class OP>
00402 void MatOrthoManager<ScalarType,MV,OP>::norm(
00403 const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const
00404 {
00405 this->normMat(X,normvec);
00406 }
00407
00408 template <class ScalarType, class MV, class OP>
00409 void MatOrthoManager<ScalarType,MV,OP>::normMat(
00410 const MV& X,
00411 std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
00412 Teuchos::RCP<const MV> MX) const
00413 {
00414 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00415 typedef Teuchos::ScalarTraits<typename SCT::magnitudeType> MT;
00416 typedef MultiVecTraits<ScalarType,MV> MVT;
00417 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00418
00419 if (!_hasOp) {
00420 MX = Teuchos::rcpFromRef(X);
00421 }
00422 else if (MX == Teuchos::null) {
00423 Teuchos::RCP<MV> R = MVT::Clone(X,MVT::GetNumberVecs(X));
00424 OPT::Apply(*_Op,X,*R);
00425 _OpCounter += MVT::GetNumberVecs(X);
00426 MX = R;
00427 }
00428
00429 Teuchos::SerialDenseMatrix<int,ScalarType> z(1,1);
00430 Teuchos::RCP<const MV> Xi, MXi;
00431 std::vector<int> ind(1);
00432 for (int i=0; i<MVT::GetNumberVecs(X); i++) {
00433 ind[0] = i;
00434 Xi = MVT::CloneView(X,ind);
00435 MXi = MVT::CloneView(*MX,ind);
00436 MVT::MvTransMv(SCT::one(),*Xi,*MXi,z);
00437 normvec[i] = MT::squareroot( SCT::magnitude(z(0,0)) );
00438 }
00439 }
00440
00441 template <class ScalarType, class MV, class OP>
00442 void MatOrthoManager<ScalarType,MV,OP>::project (
00443 MV &X,
00444 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00445 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00446 ) const
00447 {
00448 this->projectMat(X,Q,C);
00449 }
00450
00451 template <class ScalarType, class MV, class OP>
00452 int MatOrthoManager<ScalarType,MV,OP>::normalize (
00453 MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const
00454 {
00455 return this->normalizeMat(X,B);
00456 }
00457
00458 template <class ScalarType, class MV, class OP>
00459 int MatOrthoManager<ScalarType,MV,OP>::projectAndNormalize (
00460 MV &X,
00461 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00462 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00463 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B
00464 ) const
00465 {
00466 return this->projectAndNormalizeMat(X,Q,C,B);
00467 }
00468
00469 template <class ScalarType, class MV, class OP>
00470 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00471 MatOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X) const
00472 {
00473 return this->orthonormErrorMat(X,Teuchos::null);
00474 }
00475
00476 template <class ScalarType, class MV, class OP>
00477 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00478 MatOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, const MV &X2) const
00479 {
00480 return this->orthogErrorMat(X1,X2);
00481 }
00482
00483 }
00484
00485
00486 #endif
00487
00488