BelosEpetraAdapter.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 
00029 #ifndef BELOS_EPETRA_ADAPTER_HPP
00030 #define BELOS_EPETRA_ADAPTER_HPP
00031 
00036 #include "Epetra_MultiVector.h"
00037 #include "Epetra_Operator.h"
00038 #include "Epetra_Map.h"
00039 #include "Epetra_LocalMap.h"
00040 #include "Teuchos_SerialDenseMatrix.hpp"
00041 
00042 #include "BelosConfigDefs.hpp"
00043 #include "BelosMultiVec.hpp"
00044 #include "BelosOperator.hpp"
00045 #include "BelosTypes.hpp"
00046 
00047 namespace Belos {
00048  
00050 
00051 
00055   class EpetraMultiVecFailure : public BelosError {public:
00056     EpetraMultiVecFailure(const std::string& what_arg) : BelosError(what_arg)
00057     {}};
00058 
00062   class EpetraOpFailure : public BelosError {public:
00063     EpetraOpFailure(const std::string& what_arg) : BelosError(what_arg)
00064     {}};
00065 
00067  
00068   //--------template class BelosEpetraMultiVec-------------------------------------
00069   class EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector {
00070   public:
00071     // constructors
00072     EpetraMultiVec(const Epetra_BlockMap& Map, double * array, const int numvecs, const int stride=0);
00073     EpetraMultiVec(const Epetra_BlockMap& Map, const int numvecs, bool zeroOut=true);
00074     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00075     EpetraMultiVec& operator=(const EpetraMultiVec& pv) { Epetra_MultiVector::operator=(pv); return *this; }
00076     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00077     ~EpetraMultiVec();
00078     //
00079     //  member functions inherited from Belos::MultiVec
00080     //
00081     //  the following is a virtual copy constructor returning
00082     //  a pointer to the pure virtual class. std::vector values are
00083     //  not copied; instead a new MultiVec is created containing
00084     //  a non-zero amount of columns.  
00085     //
00086     MultiVec<double> * Clone ( const int numvecs ) const;
00087     //
00088     //  the following is a virtual copy constructor returning
00089     //  a pointer to the pure virtual class. std::vector values are
00090     //  copied and a new stand-alone MultiVector is created.
00091     //  (deep copy).
00092     //
00093     MultiVec<double> * CloneCopy () const;
00094     //
00095     //  the following is a virtual copy constructor returning
00096     //  a pointer to the pure virtual class. std::vector values are
00097     //  copied and a new stand-alone MultiVector is created
00098     //  where only selected columns are chosen.  (deep copy).
00099     //
00100     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00101     //
00102     //  the following is a virtual view constructor returning
00103     //  a pointer to the pure virtual class. std::vector values are 
00104     //  shared and hence no memory is allocated for the columns.
00105     //
00106     MultiVec<double> * CloneView ( const std::vector<int>& index );
00107     //
00108     //  this routine sets a subblock of the multivector, which
00109     //  need not be contiguous, and is given by the indices.
00110     //
00111     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00112     //
00113     int GetNumberVecs () const { return NumVectors(); }
00114     int GetVecLength () const { return GlobalLength(); }
00115     //
00116     // *this <- alpha * A * B + beta * (*this)
00117     //
00118     void MvTimesMatAddMv ( const double alpha, const MultiVec<double>& A, 
00119          const Teuchos::SerialDenseMatrix<int,double>& B, const double beta );
00120     //
00121     // *this <- alpha * A + beta * B
00122     //
00123     void MvAddMv ( const double alpha, const MultiVec<double>& A, const double beta,
00124        const MultiVec<double>& B);
00125 
00128     void MvScale ( double alpha ) { int ret = this->Scale( alpha ); assert(ret == 0); }
00129 
00132     void MvScale ( const std::vector<double>& alpha );
00133     //
00134     // B <- alpha * A^T * (*this)
00135     //
00136     void MvTransMv ( const double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B ) const;
00137     //
00138     // b[i] = A[i]^T * this[i]
00139     // 
00140     void MvDot ( const MultiVec<double>& A, std::vector<double>& b ) const;
00141     //
00142     // alpha[i] = norm of i-th column of (*this)
00143     //  
00144     void MvNorm ( std::vector<double>& normvec, NormType norm_type = TwoNorm ) const;
00145     //
00146     // random vectors in i-th column of (*this)
00147     //
00148     void MvRandom() { 
00149       TEST_FOR_EXCEPTION( Random()!=0, EpetraMultiVecFailure, 
00150         "Belos::EpetraMultiVec::MvRandom() call to Random() returned a nonzero value."); }
00151     //
00152     // initializes each element of (*this) with alpha
00153     //
00154     void MvInit ( const double alpha ) { 
00155       TEST_FOR_EXCEPTION( PutScalar(alpha)!=0, EpetraMultiVecFailure, 
00156         "Belos::EpetraMultiVec::MvInit() call to PutScalar() returned a nonzero value."); }
00157     //
00158     // print (*this)
00159     //
00160     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00161   private:
00162   };
00163   
00164   
00166   //--------template class BelosEpetraOp---------------------
00167   
00168   class EpetraOp : public virtual Operator<double> {
00169   public:
00170     EpetraOp( const Teuchos::RCP<Epetra_Operator> &Op );
00171     ~EpetraOp() {};
00172     void Apply ( const MultiVec<double>& x, MultiVec<double>& y, ETrans trans=NOTRANS ) const;
00173   private:
00174     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00175   };
00176   
00178   //--------template class BelosEpetraPrecOp---------------------
00179   
00180   class EpetraPrecOp : public virtual Operator<double>, public virtual Epetra_Operator {
00181   public:
00183     EpetraPrecOp( const Teuchos::RCP<Epetra_Operator> &Op );
00184 
00186     virtual ~EpetraPrecOp() {};
00187 
00189     void Apply ( const MultiVec<double>& x, MultiVec<double>& y, ETrans trans=NOTRANS ) const;
00190 
00192     int Apply( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const;
00193 
00195     int ApplyInverse( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const;
00196 
00198     const char* Label() const { return "Epetra_Operator applying A^{-1} as A"; };
00199 
00201     bool UseTranspose() const { return (false); };
00202 
00204     int SetUseTranspose(bool UseTranspose) { return 0; };
00205 
00207     bool HasNormInf() const { return (false); };
00208 
00210     double NormInf() const  { return (-1.0); };
00211 
00213     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00214 
00216     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00217 
00219     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00220     
00221   private:
00222 
00223     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00224 
00225   };
00226 
00227   
00229   //
00230   // Implementation of the Belos::MultiVecTraits for Epetra::MultiVector.
00231   //
00233 
00234   template<>
00235   class MultiVecTraits<double, Epetra_MultiVector>
00236   {
00237   public:
00239     static Teuchos::RCP<Epetra_MultiVector> Clone( const Epetra_MultiVector& mv, const int numvecs )
00240     { return Teuchos::rcp( new Epetra_MultiVector(mv.Map(), numvecs, false) ); }
00242     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv )
00243     { return Teuchos::rcp( new Epetra_MultiVector( mv ) ); }
00245     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv, const std::vector<int>& index )
00246     { 
00247       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00248       return Teuchos::rcp( new Epetra_MultiVector(Copy, mv, &tmp_index[0], index.size()) ); 
00249     }
00251     static Teuchos::RCP<Epetra_MultiVector> CloneView( Epetra_MultiVector& mv, const std::vector<int>& index )
00252     { 
00253       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00254       return Teuchos::rcp( new Epetra_MultiVector(View, mv, &tmp_index[0], index.size()) ); 
00255     }
00257     static Teuchos::RCP<const Epetra_MultiVector> CloneView( const Epetra_MultiVector& mv, const std::vector<int>& index )
00258     { 
00259       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00260       return Teuchos::rcp( new Epetra_MultiVector(View, mv, &tmp_index[0], index.size()) ); 
00261     }
00263     static int GetVecLength( const Epetra_MultiVector& mv )
00264     { return mv.GlobalLength(); }
00266     static int GetNumberVecs( const Epetra_MultiVector& mv )
00267     { return mv.NumVectors(); }
00269     static void MvTimesMatAddMv( const double alpha, const Epetra_MultiVector& A, 
00270          const Teuchos::SerialDenseMatrix<int,double>& B, 
00271          const double beta, Epetra_MultiVector& mv )
00272     { 
00273       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00274       Epetra_MultiVector B_Pvec(Copy, LocalMap, B.values(), B.stride(), B.numCols());
00275       
00276       int info = mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta );
00277       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00278        "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvTimesMatAddMv call to Multiply() returned a nonzero value.");
00279     }
00281     static void MvAddMv( const double alpha, const Epetra_MultiVector& A, const double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
00282     { 
00283       int info = mv.Update( alpha, A, beta, B, 0.0 );
00284       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00285        "Belos::MultiVecTraits<double,Epetra_MultiVec>::MvAddMv call to Update() returned a nonzero value.");
00286     }
00287 
00290     static void MvScale ( Epetra_MultiVector& mv, double alpha )
00291     { int ret = mv.Scale( alpha );
00292       assert( ret == 0 );
00293     }
00294 
00297     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
00298     {
00299       // Check to make sure the vector is as long as the multivector has columns.
00300       int numvecs = mv.NumVectors();
00301       assert( (int)alpha.size() == numvecs );
00302 
00303       int ret = 0;
00304       std::vector<int> tmp_index( 1, 0 );
00305       for (int i=0; i<numvecs; i++) {
00306         Epetra_MultiVector temp_vec(::View, mv, &tmp_index[0], 1);
00307         ret = temp_vec.Scale( alpha[i] );
00308         assert (ret == 0);
00309         tmp_index[0]++;
00310       }
00311     }
00312 
00314     static void MvTransMv( const double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B )
00315     { 
00316       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00317       Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00318       
00319       int info = B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 );
00320       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00321        "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvTransMv call to Multiply() returned a nonzero value.");
00322     }
00324     static void MvDot( const Epetra_MultiVector& mv, const Epetra_MultiVector& A, std::vector<double>& b )
00325     {
00326       int info = mv.Dot( A, &b[0] );
00327       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00328        "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvDot call to Dot() returned a nonzero value.");   
00329     }
00331     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double>& normvec, NormType type = TwoNorm )
00332     { 
00333       if ((int)normvec.size() >= mv.NumVectors()) {
00334         int info = 0;
00335   switch( type ) {
00336   case ( OneNorm ) :
00337     info = mv.Norm1(&normvec[0]);
00338     break;
00339   case ( TwoNorm ) :
00340     info = mv.Norm2(&normvec[0]);
00341     break;
00342   case ( InfNorm ) :  
00343     info = mv.NormInf(&normvec[0]);
00344     break;
00345   default:
00346     break;
00347   }
00348         TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00349          "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvNorm call to Norm() returned a nonzero value.");
00350       }
00351     }
00353     static void SetBlock( const Epetra_MultiVector& A, const std::vector<int>& index, Epetra_MultiVector& mv )
00354     { 
00355       // Extract the "numvecs" columns of mv indicated by the index std::vector.
00356       int numvecs = index.size(), info = 0;
00357       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00358       Epetra_MultiVector temp_vec(View, mv, &tmp_index[0], numvecs);
00359       
00360       if ( A.NumVectors() != numvecs ) {
00361         std::vector<int> index2( numvecs );
00362         for(int i=0; i<numvecs; i++)
00363     index2[i] = i;
00364         Epetra_MultiVector A_vec(View, A, &index2[0], numvecs);      
00365         info = temp_vec.Update( 1.0, A_vec, 0.0, A_vec, 0.0 );
00366       }
00367       else {
00368         info = temp_vec.Update( 1.0, A, 0.0, A, 0.0 );
00369       }
00370       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00371        "Belos::MultiVecTraits<double,Epetra_MultiVector>::SetBlock call to Update() returned a nonzero value.");
00372     }
00374     static void MvRandom( Epetra_MultiVector& mv )
00375     { TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure, 
00376         "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvRandom() call to Random() returned a nonzero value.");
00377     }
00379     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00380     { TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure, 
00381         "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvInit() call to PutScalar() returned a nonzero value.");
00382     }
00384     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
00385     { os << mv << std::endl; }
00386     
00387   };        
00388   
00390   //
00391   // Implementation of the Belos::OperatorTraits for Epetra::Operator.
00392   //
00394   
00395   template <> 
00396   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
00397   {
00398   public:
00399     
00401     static void Apply ( const Epetra_Operator& Op, 
00402                 const Epetra_MultiVector& x, 
00403       Epetra_MultiVector& y,
00404       ETrans trans=NOTRANS )
00405     { 
00406       int info = 0;
00407       if ( trans )
00408   const_cast<Epetra_Operator &>(Op).SetUseTranspose( true );
00409       info = Op.Apply( x, y );
00410       if ( trans )
00411   const_cast<Epetra_Operator &>(Op).SetUseTranspose( false );      
00412       TEST_FOR_EXCEPTION(info!=0, EpetraOpFailure, 
00413        "Belos::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply call to Apply() returned a nonzero value.");
00414     }
00415     
00416   };
00417   
00418 } // end of Belos namespace 
00419 
00420 #endif 
00421 // end of file BELOS_EPETRA_ADAPTER_HPP

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7