AnasaziEpetraAdapter.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers 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 
00033 #ifndef ANASAZI_EPETRA_ADAPTER_HPP
00034 #define ANASAZI_EPETRA_ADAPTER_HPP
00035 
00036 #include "AnasaziConfigDefs.hpp"
00037 #include "AnasaziTypes.hpp"
00038 #include "AnasaziMultiVec.hpp"
00039 #include "AnasaziOperator.hpp"
00040 
00041 #include "Teuchos_SerialDenseMatrix.hpp"
00042 #include "Epetra_MultiVector.h"
00043 #include "Epetra_Operator.h"
00044 #include "Epetra_Map.h"
00045 #include "Epetra_LocalMap.h"
00046 
00047 namespace Anasazi {
00048   
00050   //
00051   //--------template class AnasaziEpetraMultiVec-----------------
00052   //
00054   
00061   class EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector {
00062   public:
00064 
00065 
00067 
00072     EpetraMultiVec(const Epetra_BlockMap& Map, const int numvecs);
00073 
00075     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00076     
00078 
00086     EpetraMultiVec(const Epetra_BlockMap& Map, double * array, const int numvecs, const int stride=0);
00087 
00089 
00095     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00096 
00098     virtual ~EpetraMultiVec() {};
00099 
00101 
00103 
00104 
00109     MultiVec<double> * Clone ( const int numvecs ) const;
00110 
00116     MultiVec<double> * CloneCopy () const;
00117 
00125     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00126     
00134     MultiVec<double> * CloneView ( const std::vector<int>& index );
00135 
00137 
00139 
00140 
00142     int GetNumberVecs () const { return NumVectors(); }
00143 
00145     int GetVecLength () const { return GlobalLength(); }
00146 
00148 
00150 
00151 
00153     void MvTimesMatAddMv ( const double alpha, const MultiVec<double>& A, 
00154          const Teuchos::SerialDenseMatrix<int,double>& B, const double beta );
00155 
00158     void MvAddMv ( const double alpha, const MultiVec<double>& A, const double beta,
00159        const MultiVec<double>& B);
00160 
00163     void MvTransMv ( const double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00164 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00165          , ConjType conj = Anasazi::CONJ
00166 #endif
00167          ) const;
00168   
00171     void MvDot ( const MultiVec<double>& A, std::vector<double>* b
00172 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00173      , ConjType conj = Anasazi::CONJ
00174 #endif
00175      ) const;
00176 
00179     void MvScale ( const double alpha ) { int ret = this->Scale( alpha ); assert(ret == 0); }
00180     
00183     void MvScale ( const std::vector<double>& alpha );
00184 
00186 
00187 
00188     
00192     void MvNorm ( std::vector<double>* normvec ) const {
00193       if ((normvec!=NULL) && ((int)normvec->size() >= GetNumberVecs()) ) {
00194   int ret = Norm2(&(*normvec)[0]);
00195   assert( ret == 0 );
00196       }
00197     };
00199 
00201 
00202 
00207     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00208 
00211     void MvRandom() { int ret = Random(); assert( ret == 0 ); };
00212 
00215     void MvInit ( const double alpha ) { int ret = PutScalar( alpha ); assert( ret == 0 ); };
00216 
00218 
00219 
00220 
00222     void MvPrint( ostream& os ) const { os << *this << endl; };
00224 
00225   private:
00226   };
00227   //-------------------------------------------------------------
00228   
00230   //
00231   //--------template class AnasaziEpetraOp---------------------
00232   //
00234   
00241   class EpetraOp : public virtual Operator<double> {
00242   public:
00244 
00245     
00247     EpetraOp(const Teuchos::RefCountPtr<Epetra_Operator> &Op );
00248     
00250     ~EpetraOp();
00252     
00254 
00255     
00259     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00261     
00262   private:
00263     Teuchos::RefCountPtr<Epetra_Operator> Epetra_Op;
00264   };
00265   //-------------------------------------------------------------
00266 
00268   //
00269   //--------template class AnasaziEpetraGenOp--------------------
00270   //
00272   
00286   class EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00287   public:
00289 
00292     EpetraGenOp(const Teuchos::RefCountPtr<Epetra_Operator> &AOp, 
00293                 const Teuchos::RefCountPtr<Epetra_Operator> &MOp,
00294     bool isAInverse = true );
00295 
00297     ~EpetraGenOp();
00298     
00300 
00302     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00303 
00305 
00307     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00308 
00310 
00312     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00313 
00315     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00316     
00318     bool UseTranspose() const { return (false); };
00319 
00321     int SetUseTranspose(bool UseTranspose) { return 0; };
00322     
00324     bool HasNormInf() const { return (false); };
00325     
00327     double NormInf() const  { return (-1.0); };
00328     
00330     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00331 
00333     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00334 
00336     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00337 
00338   private:
00339     bool isAInverse;
00340     Teuchos::RefCountPtr<Epetra_Operator> Epetra_AOp;
00341     Teuchos::RefCountPtr<Epetra_Operator> Epetra_MOp;
00342   };
00343   
00345   //
00346   //--------template class AnasaziEpetraSymOp--------------------
00347   //
00349 
00362   class EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00363   public:
00365 
00367     EpetraSymOp(const Teuchos::RefCountPtr<Epetra_Operator> &Op, const bool isTrans = false );
00368 
00370     ~EpetraSymOp();
00371     
00373 
00375     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00376 
00378 
00380     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00381 
00383 
00386     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00387 
00389     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00390     
00392     bool UseTranspose() const { return (false); };
00393 
00395     int SetUseTranspose(bool UseTranspose) { return 0; };
00396     
00398     bool HasNormInf() const { return (false); };
00399     
00401     double NormInf() const  { return (-1.0); };
00402     
00404     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00405 
00407     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00408 
00410     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00411 
00412   private:
00413     Teuchos::RefCountPtr<Epetra_Operator> Epetra_Op;
00414     bool isTrans_;
00415   };
00416 
00417 
00419   //
00420   //--------template class AnasaziEpetraSymMVOp---------------------
00421   //
00423 
00436   class EpetraSymMVOp : public virtual Operator<double> {
00437   public:
00439 
00441     EpetraSymMVOp(const Teuchos::RefCountPtr<Epetra_MultiVector> &MV, 
00442       const bool isTrans = false );
00443     
00445     ~EpetraSymMVOp() {};
00446     
00448 
00450     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00451 
00452   private:
00453     Teuchos::RefCountPtr<Epetra_MultiVector> Epetra_MV;
00454     Teuchos::RefCountPtr<const Epetra_Map> MV_localmap;
00455     Teuchos::RefCountPtr<const Epetra_BlockMap> MV_blockmap;
00456     bool isTrans_;
00457   };
00458 
00459   
00461   //
00462   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00463   //
00465 
00476   template<>
00477   class MultiVecTraits<double, Epetra_MultiVector>
00478   {
00479   public:
00480 
00482 
00483 
00488     static Teuchos::RefCountPtr<Epetra_MultiVector> Clone( const Epetra_MultiVector& mv, const int numvecs )
00489     { return Teuchos::rcp( new Epetra_MultiVector(mv.Map(), numvecs) ); }
00490 
00495     static Teuchos::RefCountPtr<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv )
00496     { return Teuchos::rcp( new Epetra_MultiVector( mv ) ); }
00497 
00503     static Teuchos::RefCountPtr<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv, const std::vector<int>& index )
00504     { 
00505       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00506       return Teuchos::rcp( new Epetra_MultiVector(::Copy, mv, &tmp_index[0], index.size()) ); 
00507     }
00508 
00514     static Teuchos::RefCountPtr<Epetra_MultiVector> CloneView( Epetra_MultiVector& mv, const std::vector<int>& index )
00515     { 
00516       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00517       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00518     }
00519 
00525     static Teuchos::RefCountPtr<const Epetra_MultiVector> CloneView( const Epetra_MultiVector& mv, const std::vector<int>& index )
00526     { 
00527       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00528       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00529     }
00530 
00532 
00534 
00535 
00537     static int GetVecLength( const Epetra_MultiVector& mv )
00538     { return mv.GlobalLength(); }
00539 
00541     static int GetNumberVecs( const Epetra_MultiVector& mv )
00542     { return mv.NumVectors(); }
00544 
00546 
00547 
00550     static void MvTimesMatAddMv( const double alpha, const Epetra_MultiVector& A, 
00551          const Teuchos::SerialDenseMatrix<int,double>& B, 
00552          const double beta, Epetra_MultiVector& mv )
00553     { 
00554       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00555       Epetra_MultiVector B_Pvec(::Copy, LocalMap, B.values(), B.stride(), B.numCols());
00556 
00557       int ret = mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta );
00558       assert( ret == 0 );   
00559     }
00560 
00563     static void MvAddMv( const double alpha, const Epetra_MultiVector& A, const double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
00564     { 
00565       int ret = mv.Update( alpha, A, beta, B, 0.0 );
00566       assert( ret == 0 );
00567     }
00568 
00571     static void MvTransMv( const double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
00572 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00573          , ConjType conj = Anasazi::CONJ
00574 #endif
00575          )
00576     { 
00577       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00578       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00579       
00580       int ret = B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 );
00581       assert( ret == 0 );
00582     }
00583     
00586     static void MvDot( const Epetra_MultiVector& mv, const Epetra_MultiVector& A, std::vector<double>* b
00587 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00588            , ConjType conj = Anasazi::CONJ
00589 #endif
00590            )
00591     {
00592       int ret = mv.Dot( A, &(*b)[0] );
00593       assert( ret == 0 );
00594     }
00595 
00597 
00598 
00599 
00603     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double>* normvec )
00604     { 
00605       int ret = mv.Norm2(&(*normvec)[0]);
00606       assert( ret == 0 );
00607     }
00608 
00610 
00612 
00613 
00615     static void SetBlock( const Epetra_MultiVector& A, const std::vector<int>& index, Epetra_MultiVector& mv )
00616     { 
00617       // Extract the "numvecs" columns of mv indicated by the index vector.
00618       int numvecs = index.size();
00619       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00620       Epetra_MultiVector temp_vec(::View, mv, &tmp_index[0], numvecs);
00621 
00622       if ( A.NumVectors() != numvecs ) {
00623         std::vector<int> index2( numvecs );
00624         for(int i=0; i<numvecs; i++)
00625     index2[i] = i;
00626         Epetra_MultiVector A_vec(::View, A, &index2[0], numvecs);      
00627         int ret = temp_vec.Update( 1.0, A_vec, 0.0, A_vec, 0.0 );
00628   assert( ret == 0 );
00629       }
00630       else {
00631         int ret = temp_vec.Update( 1.0, A, 0.0, A, 0.0 );
00632   assert( ret == 0 );
00633       }
00634     }
00635 
00638     void MvScale ( const double alpha, Epetra_MultiVector& mv ) 
00639     { int ret = mv.Scale( alpha ); 
00640       assert( ret == 0 );
00641     }
00642     
00645     void MvScale ( const std::vector<double>& alpha, Epetra_MultiVector& mv )
00646     { 
00647       // Check to make sure the vector is as long as the multivector has columns.
00648       int numvecs = mv.NumVectors();
00649       assert( (int)alpha.size() == numvecs );
00650 
00651       int ret = 0;
00652       std::vector<int> tmp_index( 1, 0 );
00653       for (int i=0; i<numvecs; i++) {
00654   Epetra_MultiVector temp_vec(::View, mv, &tmp_index[0], 1);
00655         ret = temp_vec.Scale( alpha[i] );
00656   assert (ret == 0);
00657   tmp_index[0]++;
00658       }
00659     }
00660 
00663     static void MvRandom( Epetra_MultiVector& mv )
00664     { int ret = mv.Random(); assert( ret == 0 ); }
00665 
00668     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00669     { int ret = mv.PutScalar(alpha); assert( ret == 0 ); }
00670 
00672 
00674 
00675 
00678     static void MvPrint( const Epetra_MultiVector& mv, ostream& os )
00679     { os << mv << endl; }
00680 
00682   };        
00683 
00685   //
00686   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
00687   //
00689 
00701   template <> 
00702   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
00703   {
00704   public:
00705     
00709     static void Apply ( const Epetra_Operator& Op, 
00710             const Epetra_MultiVector& x, 
00711             Epetra_MultiVector& y )
00712     { 
00713       TEST_FOR_EXCEPTION( Op.Apply( x, y ) != 0, OperatorError, "Error in Epetra_Operator::Apply()!" );
00714     }
00715     
00716   };
00717   
00718 } // end of Anasazi namespace 
00719 
00720 #endif 
00721 // end of file ANASAZI_EPETRA_ADAPTER_HPP

Generated on Thu Sep 18 12:31:36 2008 for Anasazi by doxygen 1.3.9.1