AnasaziMatOrthoManager.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
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       // attempt to minimize the amount of work in applying 
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       // the user has done the matrix vector for us
00385       MVT::MvTransMv(SCT::one(),X,*MY,Z);
00386     }
00387     else {
00388       // there is no matrix vector
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 } // end of Anasazi namespace
00484 
00485 
00486 #endif
00487 
00488 // end of file AnasaziMatOrthoManager.hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 09:56:58 2011 for Anasazi by  doxygen 1.6.3