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_in, double * array, const int numvecs, const int stride=0);
00073     EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs, bool zeroOut=true);
00074     EpetraMultiVec(Epetra_DataAccess CV_in, 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 ) { 
00129       TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure, 
00130         "Belos::EpetraMultiVec::MvInit() call to Scale() returned a nonzero value."); }
00131 
00134     void MvScale ( const std::vector<double>& alpha );
00135     //
00136     // B <- alpha * A^T * (*this)
00137     //
00138     void MvTransMv ( const double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B ) const;
00139     //
00140     // b[i] = A[i]^T * this[i]
00141     // 
00142     void MvDot ( const MultiVec<double>& A, std::vector<double>& b ) const;
00143     //
00144     // alpha[i] = norm of i-th column of (*this)
00145     //  
00146     void MvNorm ( std::vector<double>& normvec, NormType norm_type = TwoNorm ) const;
00147     //
00148     // random vectors in i-th column of (*this)
00149     //
00150     void MvRandom() { 
00151       TEST_FOR_EXCEPTION( Random()!=0, EpetraMultiVecFailure, 
00152         "Belos::EpetraMultiVec::MvRandom() call to Random() returned a nonzero value."); }
00153     //
00154     // initializes each element of (*this) with alpha
00155     //
00156     void MvInit ( const double alpha ) { 
00157       TEST_FOR_EXCEPTION( PutScalar(alpha)!=0, EpetraMultiVecFailure, 
00158         "Belos::EpetraMultiVec::MvInit() call to PutScalar() returned a nonzero value."); }
00159     //
00160     // print (*this)
00161     //
00162     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00163   private:
00164   };
00165   
00166   
00168   //--------template class BelosEpetraOp---------------------
00169   
00170   class EpetraOp : public virtual Operator<double> {
00171   public:
00172     EpetraOp( const Teuchos::RCP<Epetra_Operator> &Op );
00173     ~EpetraOp() {};
00174     void Apply ( const MultiVec<double>& x, MultiVec<double>& y, ETrans trans=NOTRANS ) const;
00175   private:
00176     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00177   };
00178   
00180   //--------template class BelosEpetraPrecOp---------------------
00181   
00182   class EpetraPrecOp : public virtual Operator<double>, public virtual Epetra_Operator {
00183   public:
00185     EpetraPrecOp( const Teuchos::RCP<Epetra_Operator> &Op );
00186 
00188     virtual ~EpetraPrecOp() {};
00189 
00191     void Apply ( const MultiVec<double>& x, MultiVec<double>& y, ETrans trans=NOTRANS ) const;
00192 
00194     int Apply( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const;
00195 
00197     int ApplyInverse( const Epetra_MultiVector &X, Epetra_MultiVector &Y ) const;
00198 
00200     const char* Label() const { return "Epetra_Operator applying A^{-1} as A"; };
00201 
00203     bool UseTranspose() const { return (false); };
00204 
00206     int SetUseTranspose(bool UseTranspose_in) { return 0; };
00207 
00209     bool HasNormInf() const { return (false); };
00210 
00212     double NormInf() const  { return (-1.0); };
00213 
00215     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00216 
00218     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00219 
00221     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00222     
00223   private:
00224 
00225     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00226 
00227   };
00228 
00229   
00231   //
00232   // Implementation of the Belos::MultiVecTraits for Epetra::MultiVector.
00233   //
00235 
00236   template<>
00237   class MultiVecTraits<double, Epetra_MultiVector>
00238   {
00239   public:
00241     static Teuchos::RCP<Epetra_MultiVector> Clone( const Epetra_MultiVector& mv, const int numvecs )
00242     { return Teuchos::rcp( new Epetra_MultiVector(mv.Map(), numvecs, false) ); }
00244     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv )
00245     { return Teuchos::rcp( new Epetra_MultiVector( mv ) ); }
00247     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv, const std::vector<int>& index )
00248     { 
00249       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00250       return Teuchos::rcp( new Epetra_MultiVector(Copy, mv, &tmp_index[0], index.size()) ); 
00251     }
00253     static Teuchos::RCP<Epetra_MultiVector> CloneView( Epetra_MultiVector& mv, const std::vector<int>& index )
00254     { 
00255       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00256       return Teuchos::rcp( new Epetra_MultiVector(View, mv, &tmp_index[0], index.size()) ); 
00257     }
00259     static Teuchos::RCP<const Epetra_MultiVector> CloneView( const Epetra_MultiVector& mv, const std::vector<int>& index )
00260     { 
00261       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00262       return Teuchos::rcp( new Epetra_MultiVector(View, mv, &tmp_index[0], index.size()) ); 
00263     }
00265     static int GetVecLength( const Epetra_MultiVector& mv )
00266     { return mv.GlobalLength(); }
00268     static int GetNumberVecs( const Epetra_MultiVector& mv )
00269     { return mv.NumVectors(); }
00271     static void MvTimesMatAddMv( const double alpha, const Epetra_MultiVector& A, 
00272          const Teuchos::SerialDenseMatrix<int,double>& B, 
00273          const double beta, Epetra_MultiVector& mv )
00274     { 
00275       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00276       Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00277       
00278       int info = mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta );
00279       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00280        "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvTimesMatAddMv call to Multiply() returned a nonzero value.");
00281     }
00283     static void MvAddMv( const double alpha, const Epetra_MultiVector& A, const double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
00284     { 
00285       int info = mv.Update( alpha, A, beta, B, 0.0 );
00286       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00287        "Belos::MultiVecTraits<double,Epetra_MultiVec>::MvAddMv call to Update() returned a nonzero value.");
00288     }
00289 
00292     static void MvScale ( Epetra_MultiVector& mv, double alpha )
00293     { int ret = mv.Scale( alpha );
00294       TEST_FOR_EXCEPTION(ret!=0, EpetraMultiVecFailure, 
00295        "Belos::MultiVecTraits<double,Epetra_MultiVec>::MvScale call to Scale() returned a nonzero value.");
00296     }
00297 
00300     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
00301     {
00302       // Check to make sure the vector is as long as the multivector has columns.
00303       int numvecs = mv.NumVectors();
00304       TEST_FOR_EXCEPTION((int)alpha.size() != numvecs, EpetraMultiVecFailure, 
00305        "Belos::MultiVecTraits<double,Epetra_MultiVec>::MvScale scaling vector (alpha) not same size as number of input vectors (mv).");
00306 
00307       int ret = 0;
00308       std::vector<int> tmp_index( 1, 0 );
00309       for (int i=0; i<numvecs; i++) {
00310         Epetra_MultiVector temp_vec(::View, mv, &tmp_index[0], 1);
00311         ret = temp_vec.Scale( alpha[i] );
00312         TEST_FOR_EXCEPTION(ret!=0, EpetraMultiVecFailure, 
00313                      "Belos::MultiVecTraits<double,Epetra_MultiVec>::MvScale call to Scale() returned a nonzero value.");
00314         tmp_index[0]++;
00315       }
00316     }
00317 
00319     static void MvTransMv( const double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B )
00320     { 
00321       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00322       Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
00323       
00324       int info = B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 );
00325       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00326        "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvTransMv call to Multiply() returned a nonzero value.");
00327     }
00329     static void MvDot( const Epetra_MultiVector& mv, const Epetra_MultiVector& A, std::vector<double>& b )
00330     {
00331       int info = mv.Dot( A, &b[0] );
00332       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00333        "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvDot call to Dot() returned a nonzero value.");   
00334     }
00336     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double>& normvec, NormType type = TwoNorm )
00337     { 
00338       if ((int)normvec.size() >= mv.NumVectors()) {
00339         int info = 0;
00340   switch( type ) {
00341   case ( OneNorm ) :
00342     info = mv.Norm1(&normvec[0]);
00343     break;
00344   case ( TwoNorm ) :
00345     info = mv.Norm2(&normvec[0]);
00346     break;
00347   case ( InfNorm ) :  
00348     info = mv.NormInf(&normvec[0]);
00349     break;
00350   default:
00351     break;
00352   }
00353         TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00354          "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvNorm call to Norm() returned a nonzero value.");
00355       }
00356     }
00358     static void SetBlock( const Epetra_MultiVector& A, const std::vector<int>& index, Epetra_MultiVector& mv )
00359     { 
00360       // Extract the "numvecs" columns of mv indicated by the index std::vector.
00361       int numvecs = index.size(), info = 0;
00362       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00363       Epetra_MultiVector temp_vec(View, mv, &tmp_index[0], numvecs);
00364       
00365       if ( A.NumVectors() != numvecs ) {
00366         std::vector<int> index2( numvecs );
00367         for(int i=0; i<numvecs; i++)
00368     index2[i] = i;
00369         Epetra_MultiVector A_vec(View, A, &index2[0], numvecs);      
00370         info = temp_vec.Update( 1.0, A_vec, 0.0, A_vec, 0.0 );
00371       }
00372       else {
00373         info = temp_vec.Update( 1.0, A, 0.0, A, 0.0 );
00374       }
00375       TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
00376        "Belos::MultiVecTraits<double,Epetra_MultiVector>::SetBlock call to Update() returned a nonzero value.");
00377     }
00379     static void MvRandom( Epetra_MultiVector& mv )
00380     { TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure, 
00381         "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvRandom() call to Random() returned a nonzero value.");
00382     }
00384     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00385     { TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure, 
00386         "Belos::MultiVecTraits<double,Epetra_MultiVector>::MvInit() call to PutScalar() returned a nonzero value.");
00387     }
00389     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
00390     { os << mv << std::endl; }
00391     
00392   };        
00393   
00395   //
00396   // Implementation of the Belos::OperatorTraits for Epetra::Operator.
00397   //
00399   
00400   template <> 
00401   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
00402   {
00403   public:
00404     
00406     static void Apply ( const Epetra_Operator& Op, 
00407                 const Epetra_MultiVector& x, 
00408       Epetra_MultiVector& y,
00409       ETrans trans=NOTRANS )
00410     { 
00411       int info = 0;
00412       if ( trans )
00413   const_cast<Epetra_Operator &>(Op).SetUseTranspose( true );
00414       info = Op.Apply( x, y );
00415       if ( trans )
00416   const_cast<Epetra_Operator &>(Op).SetUseTranspose( false );      
00417       TEST_FOR_EXCEPTION(info!=0, EpetraOpFailure, 
00418        "Belos::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply call to Apply() returned a nonzero value.");
00419     }
00420     
00421   };
00422   
00423 } // end of Belos namespace 
00424 
00425 #endif 
00426 // end of file BELOS_EPETRA_ADAPTER_HPP

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