Belos Version of the Day
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00047 #ifndef BELOS_MATORTHOMANAGER_HPP
00048 #define BELOS_MATORTHOMANAGER_HPP
00049 
00065 #include "BelosConfigDefs.hpp"
00066 #include "BelosTypes.hpp"
00067 #include "BelosOrthoManager.hpp"
00068 #include "BelosMultiVecTraits.hpp"
00069 #include "BelosOperatorTraits.hpp"
00070 
00071 namespace Belos {
00072 
00073   template <class ScalarType, class MV, class OP>
00074   class MatOrthoManager : public OrthoManager<ScalarType,MV> {
00075   protected:
00076     Teuchos::RCP<const OP> _Op;
00077     bool _hasOp;
00078 
00079   public:
00081 
00082 
00083     MatOrthoManager(Teuchos::RCP<const OP> Op = Teuchos::null) : _Op(Op), _hasOp(Op!=Teuchos::null) {};
00084 
00086     virtual ~MatOrthoManager() {};
00088 
00090 
00091 
00093     void setOp( Teuchos::RCP<const OP> Op ) { 
00094       _Op = Op; 
00095       _hasOp = (_Op != Teuchos::null);
00096     };
00097 
00099     Teuchos::RCP<const OP> getOp() const { return _Op; } 
00100 
00102 
00103 
00105 
00106 
00111     void innerProd( const MV& X, const MV& Y, 
00112                                   Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const {
00113       typedef Teuchos::ScalarTraits<ScalarType> SCT;
00114       typedef MultiVecTraits<ScalarType,MV>     MVT;
00115       typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00116 
00117       Teuchos::RCP<const MV> P,Q;
00118       Teuchos::RCP<MV> R;
00119 
00120       if (_hasOp) {
00121         // attempt to minimize the amount of work in applying 
00122         if ( MVT::GetNumberVecs(X) < MVT::GetNumberVecs(Y) ) {
00123           R = MVT::Clone(X,MVT::GetNumberVecs(X));
00124           OPT::Apply(*_Op,X,*R);
00125           P = R;
00126           Q = Teuchos::rcp( &Y, false );
00127         }
00128         else {
00129           P = Teuchos::rcp( &X, false );
00130           R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
00131           OPT::Apply(*_Op,Y,*R);
00132           Q = R;
00133         }
00134       }
00135       else {
00136         P = Teuchos::rcp( &X, false );
00137         Q = Teuchos::rcp( &Y, false );
00138       }
00139       
00140       MVT::MvTransMv(SCT::one(),*P,*Q,Z);
00141     }
00142 
00149     void innerProd( const MV& X, const MV& Y, Teuchos::RCP<const MV> MY, 
00150                             Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const {
00151       typedef Teuchos::ScalarTraits<ScalarType> SCT;
00152       typedef MultiVecTraits<ScalarType,MV>     MVT;
00153       typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00154 
00155       Teuchos::RCP<MV> P,Q;
00156 
00157       if ( MY == Teuchos::null ) {
00158         innerProd(X,Y,Z);
00159       }
00160       else if ( _hasOp ) {
00161         // the user has done the matrix std::vector for us
00162         MVT::MvTransMv(SCT::one(),X,*MY,Z);
00163       }
00164       else {
00165         // there is no matrix std::vector
00166         MVT::MvTransMv(SCT::one(),X,Y,Z);
00167       }
00168     }
00169 
00172     void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >& normvec ) const {
00173       norm(X,Teuchos::null,normvec);
00174     }
00175 
00193     void 
00194     norm (const MV& X, 
00195     Teuchos::RCP<const MV> MX, 
00196     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec) const 
00197     {
00198       typedef Teuchos::ScalarTraits<ScalarType> SCT;
00199       typedef MultiVecTraits<ScalarType,MV>     MVT;
00200       typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00201 
00202       const int numColsX = MVT::GetNumberVecs(X);
00203       if (!_hasOp) {
00204   // X == MX, since the operator M is the identity.
00205         MX = Teuchos::rcp(&X, false);
00206       }
00207       else if (MX.is_null()) {
00208   // The caller didn't give us a previously computed MX, so
00209   // apply the operator.  We assign to MX only after applying
00210   // the operator, so that if the application fails, MX won't be
00211   // modified.
00212         Teuchos::RCP<MV> R = MVT::Clone(X, numColsX);
00213         OPT::Apply(*_Op,X,*R);
00214         MX = R;
00215       } 
00216       else {
00217   // The caller gave us a previously computed MX.  Make sure
00218   // that it has at least as many columns as X.
00219   const int numColsMX = MVT::GetNumberVecs(*MX);
00220   TEST_FOR_EXCEPTION(numColsMX < numColsX, std::invalid_argument,
00221          "MatOrthoManager::norm(X, MX, normvec): "
00222          "MX has fewer columns than X: "
00223          "MX has " << numColsMX << " columns, "
00224          "and X has " << numColsX << " columns.");
00225       }
00226 
00227       // Make sure that normvec has enough entries to hold the norms
00228       // of all the columns of X.  std::vector<T>::size_type is
00229       // unsigned, so do the appropriate cast to avoid signed/unsigned
00230       // comparisons that trigger compiler warnings.
00231       if (normvec.size() < static_cast<size_t>(numColsX))
00232   normvec.resize (numColsX);
00233 
00234       Teuchos::SerialDenseMatrix<int,ScalarType> z(1,1);
00235       Teuchos::RCP<const MV> Xi, MXi;
00236       std::vector<int> ind(1);
00237       for (int i = 0; i < numColsX; ++i) {
00238         ind[0] = i;
00239         Xi = MVT::CloneView(X,ind);
00240         MXi = MVT::CloneView(*MX,ind);
00241         MVT::MvTransMv(SCT::one(),*Xi,*MXi,z);
00242         normvec[i] = SCT::magnitude( SCT::squareroot( z(0,0) ) );
00243       }
00244     }
00245 
00246 
00268     virtual void project ( MV &X, Teuchos::RCP<MV> MX, 
00269                            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00270                            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
00271 
00272 
00273     
00276     virtual void project ( MV &X, 
00277                            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00278                            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00279       project(X,Teuchos::null,C,Q);
00280     }
00281 
00303     virtual int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00304                             Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const = 0;
00305 
00306 
00309     virtual int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00310       return normalize(X,Teuchos::null,B);
00311     }
00312 
00313 
00348     virtual int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00349                                       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00350                                       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00351                                       Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const = 0;
00352 
00355     virtual int projectAndNormalize ( MV &X, 
00356                                       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00357                                       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00358                                       Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const {
00359       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00360     }
00361 
00363 
00365 
00366 
00369     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00370     orthonormError(const MV &X) const {
00371       return orthonormError(X,Teuchos::null);
00372     }
00373 
00377     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00378     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const = 0;
00379 
00382     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00383     orthogError(const MV &X1, const MV &X2) const {
00384       return orthogError(X1,Teuchos::null,X2);
00385     }
00386 
00391     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00392     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const = 0;
00393 
00395 
00396   };
00397 
00398 } // end of Belos namespace
00399 
00400 
00401 #endif
00402 
00403 // end of file BelosMatOrthoManager.hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines