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 
00154       Teuchos::RCP<MV> P,Q;
00155 
00156       if ( MY == Teuchos::null ) {
00157         innerProd(X,Y,Z);
00158       }
00159       else if ( _hasOp ) {
00160         // the user has done the matrix vector for us
00161         MVT::MvTransMv(SCT::one(),X,*MY,Z);
00162       }
00163       else {
00164         // there is no matrix vector
00165         MVT::MvTransMv(SCT::one(),X,Y,Z);
00166       }
00167     }
00168 
00171     void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >& normvec ) const {
00172       norm(X,Teuchos::null,normvec);
00173     }
00174 
00192     void 
00193     norm (const MV& X, 
00194     Teuchos::RCP<const MV> MX, 
00195     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec) const 
00196     {
00197       typedef Teuchos::ScalarTraits<ScalarType> SCT;
00198       typedef MultiVecTraits<ScalarType,MV>     MVT;
00199       typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00200 
00201       const int numColsX = MVT::GetNumberVecs(X);
00202       if (!_hasOp) {
00203   // X == MX, since the operator M is the identity.
00204         MX = Teuchos::rcp(&X, false);
00205       }
00206       else if (MX.is_null()) {
00207   // The caller didn't give us a previously computed MX, so
00208   // apply the operator.  We assign to MX only after applying
00209   // the operator, so that if the application fails, MX won't be
00210   // modified.
00211         Teuchos::RCP<MV> R = MVT::Clone(X, numColsX);
00212         OPT::Apply(*_Op,X,*R);
00213         MX = R;
00214       } 
00215       else {
00216   // The caller gave us a previously computed MX.  Make sure
00217   // that it has at least as many columns as X.
00218   const int numColsMX = MVT::GetNumberVecs(*MX);
00219   TEUCHOS_TEST_FOR_EXCEPTION(numColsMX < numColsX, std::invalid_argument,
00220          "MatOrthoManager::norm(X, MX, normvec): "
00221          "MX has fewer columns than X: "
00222          "MX has " << numColsMX << " columns, "
00223          "and X has " << numColsX << " columns.");
00224       }
00225 
00226       // Make sure that normvec has enough entries to hold the norms
00227       // of all the columns of X.  std::vector<T>::size_type is
00228       // unsigned, so do the appropriate cast to avoid signed/unsigned
00229       // comparisons that trigger compiler warnings.
00230       if (normvec.size() < static_cast<size_t>(numColsX))
00231   normvec.resize (numColsX);
00232 
00233       Teuchos::SerialDenseMatrix<int,ScalarType> z(1,1);
00234       Teuchos::RCP<const MV> Xi, MXi;
00235       std::vector<int> ind(1);
00236       for (int i = 0; i < numColsX; ++i) {
00237         ind[0] = i;
00238         Xi = MVT::CloneView(X,ind);
00239         MXi = MVT::CloneView(*MX,ind);
00240         MVT::MvTransMv(SCT::one(),*Xi,*MXi,z);
00241         normvec[i] = SCT::magnitude( SCT::squareroot( z(0,0) ) );
00242       }
00243     }
00244 
00245 
00267     virtual void project ( MV &X, Teuchos::RCP<MV> MX, 
00268                            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00269                            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
00270 
00271 
00272     
00275     virtual void project ( MV &X, 
00276                            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00277                            Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00278       project(X,Teuchos::null,C,Q);
00279     }
00280 
00302     virtual int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00303                             Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const = 0;
00304 
00305 
00308     virtual int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00309       return normalize(X,Teuchos::null,B);
00310     }
00311 
00312 
00347     virtual int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00348                                       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00349                                       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00350                                       Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const = 0;
00351 
00354     virtual int projectAndNormalize ( MV &X, 
00355                                       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00356                                       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00357                                       Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const {
00358       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00359     }
00360 
00362 
00364 
00365 
00368     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00369     orthonormError(const MV &X) const {
00370       return orthonormError(X,Teuchos::null);
00371     }
00372 
00376     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00377     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const = 0;
00378 
00381     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00382     orthogError(const MV &X1, const MV &X2) const {
00383       return orthogError(X1,Teuchos::null,X2);
00384     }
00385 
00390     virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00391     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const = 0;
00392 
00394 
00395   };
00396 
00397 } // end of Belos namespace
00398 
00399 
00400 #endif
00401 
00402 // end of file BelosMatOrthoManager.hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines