BelosMatOrthoManager.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers 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 BELOS_MATORTHOMANAGER_HPP
00035 #define BELOS_MATORTHOMANAGER_HPP
00036 
00052 #include "BelosConfigDefs.hpp"
00053 #include "BelosTypes.hpp"
00054 #include "BelosOrthoManager.hpp"
00055 #include "BelosMultiVecTraits.hpp"
00056 #include "BelosOperatorTraits.hpp"
00057 
00058 namespace Belos {
00059 
00060   template <class ScalarType, class MV, class OP>
00061   class MatOrthoManager : public OrthoManager<ScalarType,MV> {
00062   protected:
00063     Teuchos::RCP<const OP> _Op;
00064     bool _hasOp;
00065 
00066   public:
00068 
00069 
00070     MatOrthoManager(Teuchos::RCP<const OP> Op = Teuchos::null) : _Op(Op), _hasOp(Op!=Teuchos::null) {};
00071 
00073     virtual ~MatOrthoManager() {};
00075 
00077 
00078 
00080     void setOp( Teuchos::RCP<const OP> Op ) { 
00081       _Op = Op; 
00082       _hasOp = (_Op != Teuchos::null);
00083     };
00084 
00086     Teuchos::RCP<const OP> getOp() const { return _Op; } 
00087 
00089 
00090 
00092 
00093 
00098     void innerProd( const MV& X, const MV& Y, 
00099                                   Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const {
00100       typedef Teuchos::ScalarTraits<ScalarType> SCT;
00101       typedef MultiVecTraits<ScalarType,MV>     MVT;
00102       typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00103 
00104       Teuchos::RCP<const MV> P,Q;
00105       Teuchos::RCP<MV> R;
00106 
00107       if (_hasOp) {
00108         // attempt to minimize the amount of work in applying 
00109         if ( MVT::GetNumberVecs(X) < MVT::GetNumberVecs(Y) ) {
00110           R = MVT::Clone(X,MVT::GetNumberVecs(X));
00111           OPT::Apply(*_Op,X,*R);
00112           P = R;
00113           Q = Teuchos::rcp( &Y, false );
00114         }
00115         else {
00116           P = Teuchos::rcp( &X, false );
00117           R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
00118           OPT::Apply(*_Op,Y,*R);
00119           Q = R;
00120         }
00121       }
00122       else {
00123         P = Teuchos::rcp( &X, false );
00124         Q = Teuchos::rcp( &Y, false );
00125       }
00126       
00127       MVT::MvTransMv(SCT::one(),*P,*Q,Z);
00128     }
00129 
00136     void innerProd( const MV& X, const MV& Y, Teuchos::RCP<const MV> MY, 
00137                             Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const {
00138       typedef Teuchos::ScalarTraits<ScalarType> SCT;
00139       typedef MultiVecTraits<ScalarType,MV>     MVT;
00140       typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00141 
00142       Teuchos::RCP<MV> P,Q;
00143 
00144       if ( MY == Teuchos::null ) {
00145         innerProd(X,Y,Z);
00146       }
00147       else if ( _hasOp ) {
00148         // the user has done the matrix std::vector for us
00149         MVT::MvTransMv(SCT::one(),X,*MY,Z);
00150       }
00151       else {
00152         // there is no matrix std::vector
00153         MVT::MvTransMv(SCT::one(),X,Y,Z);
00154       }
00155     }
00156 
00159     void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > normvec ) const {
00160       norm(X,Teuchos::null,normvec);
00161     }
00162 
00166     void norm( const MV& X, Teuchos::RCP<const MV> MX, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > normvec ) const {
00167 
00168       typedef Teuchos::ScalarTraits<ScalarType> SCT;
00169       typedef MultiVecTraits<ScalarType,MV>     MVT;
00170       typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00171       
00172       if (!_hasOp) {
00173         MX = Teuchos::rcp(&X,false);
00174       }
00175       else if (MX == Teuchos::null) {
00176         Teuchos::RCP<MV> R = MVT::Clone(X,MVT::GetNumberVecs(X));
00177         OPT::Apply(*_Op,X,*R);
00178         MX = R;
00179       }
00180 
00181       Teuchos::SerialDenseMatrix<int,ScalarType> z(1,1);
00182       Teuchos::RCP<const MV> Xi, MXi;
00183       std::vector<int> ind(1);
00184       for (int i=0; i<MVT::GetNumberVecs(X); i++) {
00185         ind[0] = i;
00186         Xi = MVT::CloneView(X,ind);
00187         MXi = MVT::CloneView(*MX,ind);
00188         MVT::MvTransMv(SCT::one(),*Xi,*MXi,z);
00189         normvec[i] = SCT::magnitude( SCT::squareroot( z(0,0) ) );
00190       }
00191     }
00192 
00193 
00215     virtual void project ( MV &X, Teuchos::RCP<MV> MX, 
00216                            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00217                            Teuchos::Array<Teuchos::RCP<const MV> > Q) const = 0;
00218 
00219 
00220     
00223     virtual void project ( MV &X, 
00224                            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00225                            Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00226       project(X,Teuchos::null,C,Q);
00227     }
00228 
00250     virtual int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00251                             Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const = 0;
00252 
00253 
00256     virtual int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00257       return normalize(X,Teuchos::null,B);
00258     }
00259 
00260 
00295     virtual int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00296                                       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00297                                       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00298                                       Teuchos::Array<Teuchos::RCP<const MV> > Q ) const = 0;
00299 
00302     virtual int projectAndNormalize ( MV &X, 
00303                                       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00304                                       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00305                                       Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00306       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00307     }
00308 
00310 
00312 
00313 
00316     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00317     orthonormError(const MV &X) const {
00318       return orthonormError(X,Teuchos::null);
00319     }
00320 
00324     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00325     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const = 0;
00326 
00329     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00330     orthogError(const MV &X1, const MV &X2) const {
00331       return orthogError(X1,Teuchos::null,X2);
00332     }
00333 
00338     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00339     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const = 0;
00340 
00342 
00343   };
00344 
00345 } // end of Belos namespace
00346 
00347 
00348 #endif
00349 
00350 // end of file BelosMatOrthoManager.hpp

Generated on Wed May 12 21:45:51 2010 for Belos by  doxygen 1.4.7