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 
00055   class EpetraMultiVecFailure : public AnasaziError {public:
00056     EpetraMultiVecFailure(const std::string& what_arg) : AnasaziError(what_arg)
00057     {}};
00058 
00062   class EpetraOpFailure : public AnasaziError {public:
00063     EpetraOpFailure(const std::string& what_arg) : AnasaziError(what_arg)
00064     {}};
00065 
00067   
00069   //
00070   //--------template class AnasaziEpetraMultiVec-----------------
00071   //
00073   
00080   class EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector {
00081   public:
00083 
00084 
00086 
00091     EpetraMultiVec(const Epetra_BlockMap& Map, const int numvecs);
00092 
00094     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00095     
00097 
00105     EpetraMultiVec(const Epetra_BlockMap& Map, double * array, const int numvecs, const int stride=0);
00106 
00108 
00114     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00115 
00117     virtual ~EpetraMultiVec() {};
00118 
00120 
00122 
00123 
00128     MultiVec<double> * Clone ( const int numvecs ) const;
00129 
00135     MultiVec<double> * CloneCopy () const;
00136 
00144     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00145     
00153     MultiVec<double> * CloneView ( const std::vector<int>& index );
00154 
00156 
00158 
00159 
00161     int GetNumberVecs () const { return NumVectors(); }
00162 
00164     int GetVecLength () const { return GlobalLength(); }
00165 
00167 
00169 
00170 
00172     void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00173                            const Teuchos::SerialDenseMatrix<int,double>& B, 
00174                            double beta );
00175 
00178     void MvAddMv ( double alpha, const MultiVec<double>& A, 
00179                    double beta, const MultiVec<double>& B);
00180 
00183     void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00184 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00185         , ConjType conj = Anasazi::CONJ
00186 #endif
00187         ) const;
00188   
00191     void MvDot ( const MultiVec<double>& A, std::vector<double>* b
00192 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00193         , ConjType conj = Anasazi::CONJ
00194 #endif
00195         ) const;
00196 
00199     void MvScale ( double alpha ) { 
00200       TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
00201           "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00202     }
00203     
00206     void MvScale ( const std::vector<double>& alpha );
00207 
00209 
00210 
00211     
00215     void MvNorm ( std::vector<double>* normvec ) const {
00216       if ((normvec!=NULL) && ((int)normvec->size() >= GetNumberVecs()) ) {
00217         TEST_FOR_EXCEPTION( this->Norm2(&(*normvec)[0])!=0, EpetraMultiVecFailure,
00218             "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
00219       }
00220     };
00222     
00224 
00225 
00230     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00231 
00234     void MvRandom() { 
00235       TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
00236           "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00237     };
00238 
00241     void MvInit ( double alpha ) { 
00242       TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
00243           "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00244     };
00245     
00247 
00248 
00249 
00251     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00253 
00254   private:
00255   };
00256   //-------------------------------------------------------------
00257   
00259   //
00260   //--------template class AnasaziEpetraOp---------------------
00261   //
00263   
00270   class EpetraOp : public virtual Operator<double> {
00271   public:
00273 
00274     
00276     EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
00277     
00279     ~EpetraOp();
00281     
00283 
00284     
00288     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00290     
00291   private:
00292     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00293   };
00294   //-------------------------------------------------------------
00295 
00297   //
00298   //--------template class AnasaziEpetraGenOp--------------------
00299   //
00301   
00315   class EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00316   public:
00318 
00321     EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp, 
00322                 const Teuchos::RCP<Epetra_Operator> &MOp,
00323                 bool isAInverse = true );
00324 
00326     ~EpetraGenOp();
00327     
00329 
00331     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00332 
00334 
00336     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00337 
00339 
00341     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00342 
00344     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00345     
00347     bool UseTranspose() const { return (false); };
00348 
00350     int SetUseTranspose(bool UseTranspose) { return 0; };
00351     
00353     bool HasNormInf() const { return (false); };
00354     
00356     double NormInf() const  { return (-1.0); };
00357     
00359     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00360 
00362     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00363 
00365     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00366 
00367   private:
00368     bool isAInverse;
00369     Teuchos::RCP<Epetra_Operator> Epetra_AOp;
00370     Teuchos::RCP<Epetra_Operator> Epetra_MOp;
00371   };
00372   
00374   //
00375   //--------template class AnasaziEpetraSymOp--------------------
00376   //
00378 
00391   class EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00392   public:
00394 
00396     EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
00397 
00399     ~EpetraSymOp();
00400     
00402 
00404     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00405 
00407 
00409     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00410 
00412 
00415     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00416 
00418     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00419     
00421     bool UseTranspose() const { return (false); };
00422 
00424     int SetUseTranspose(bool UseTranspose) { return 0; };
00425     
00427     bool HasNormInf() const { return (false); };
00428     
00430     double NormInf() const  { return (-1.0); };
00431     
00433     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00434 
00436     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00437 
00439     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00440 
00441   private:
00442     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00443     bool isTrans_;
00444   };
00445 
00446 
00448   //
00449   //--------template class AnasaziEpetraSymMVOp---------------------
00450   //
00452 
00465   class EpetraSymMVOp : public virtual Operator<double> {
00466   public:
00468 
00470     EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00471                   bool isTrans = false );
00472     
00474     ~EpetraSymMVOp() {};
00475     
00477 
00479     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00480 
00481   private:
00482     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00483     Teuchos::RCP<const Epetra_Map> MV_localmap;
00484     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00485     bool isTrans_;
00486   };
00487 
00489   //
00490   //--------template class AnasaziEpetraWSymMVOp---------------------
00491   //
00493 
00506   class EpetraWSymMVOp : public virtual Operator<double> {
00507   public:
00509 
00511     EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00512                    const Teuchos::RCP<Epetra_Operator> &OP );
00513     
00515     ~EpetraWSymMVOp() {};
00516     
00518 
00520     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00521 
00522   private:
00523     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00524     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00525     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00526     Teuchos::RCP<const Epetra_Map> MV_localmap;
00527     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00528   };
00529 
00530   
00531   
00533   //
00534   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00535   //
00537 
00548   template<>
00549   class MultiVecTraits<double, Epetra_MultiVector>
00550   {
00551   public:
00552 
00554 
00555 
00560     static Teuchos::RCP<Epetra_MultiVector> Clone( const Epetra_MultiVector& mv, const int numvecs )
00561     { return Teuchos::rcp( new Epetra_MultiVector(mv.Map(), numvecs) ); }
00562 
00567     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv )
00568     { return Teuchos::rcp( new Epetra_MultiVector( mv ) ); }
00569 
00575     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv, const std::vector<int>& index )
00576     { 
00577       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00578       return Teuchos::rcp( new Epetra_MultiVector(::Copy, mv, &tmp_index[0], index.size()) ); 
00579     }
00580 
00586     static Teuchos::RCP<Epetra_MultiVector> CloneView( Epetra_MultiVector& mv, const std::vector<int>& index )
00587     { 
00588       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00589       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00590     }
00591 
00597     static Teuchos::RCP<const Epetra_MultiVector> CloneView( const Epetra_MultiVector& mv, const std::vector<int>& index )
00598     { 
00599       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00600       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00601     }
00602 
00604 
00606 
00607 
00609     static int GetVecLength( const Epetra_MultiVector& mv )
00610     { return mv.GlobalLength(); }
00611 
00613     static int GetNumberVecs( const Epetra_MultiVector& mv )
00614     { return mv.NumVectors(); }
00616 
00618 
00619 
00622     static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A, 
00623                                  const Teuchos::SerialDenseMatrix<int,double>& B, 
00624                                  double beta, Epetra_MultiVector& mv )
00625     { 
00626       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00627       Epetra_MultiVector B_Pvec(::Copy, LocalMap, B.values(), B.stride(), B.numCols());
00628 
00629       TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
00630           "MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00631     }
00632 
00635     static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
00636     { 
00637       TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
00638           "MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update() returned a nonzero value.");
00639     }
00640 
00643     static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
00644 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00645                           , ConjType conj = Anasazi::CONJ
00646 #endif
00647                         )
00648     { 
00649       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00650       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00651       
00652       TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
00653           "MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00654     }
00655     
00658     static void MvDot( const Epetra_MultiVector& mv, const Epetra_MultiVector& A, std::vector<double>* b
00659 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00660                       , ConjType conj = Anasazi::CONJ
00661 #endif
00662                       )
00663     {
00664       TEST_FOR_EXCEPTION( mv.Dot( A, &(*b)[0] )!=0, EpetraMultiVecFailure,
00665           "MultiVecTraits<double, Epetra_MultiVector>::MvDot call to Epetra_MultiVector::Dot() returned a nonzero value.");     
00666     }
00667 
00669 
00670 
00671 
00675     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double>* normvec )
00676     { 
00677       TEST_FOR_EXCEPTION( mv.Norm2(&(*normvec)[0])!=0, EpetraMultiVecFailure,
00678           "MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value."); 
00679     }
00680 
00682     
00684 
00685 
00687     static void SetBlock( const Epetra_MultiVector& A, const std::vector<int>& index, Epetra_MultiVector& mv )
00688     { 
00689       // Extract the "numvecs" columns of mv indicated by the index vector.
00690       int numvecs = index.size();
00691       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00692       Epetra_MultiVector temp_vec(::View, mv, &tmp_index[0], numvecs);
00693 
00694       if ( A.NumVectors() != numvecs ) {
00695         std::vector<int> index2( numvecs );
00696         for(int i=0; i<numvecs; i++)
00697           index2[i] = i;
00698         Epetra_MultiVector A_vec(::View, A, &index2[0], numvecs);      
00699         TEST_FOR_EXCEPTION( temp_vec.Update( 1.0, A_vec, 0.0, A_vec, 0.0 )!=0, EpetraMultiVecFailure,
00700             "MultiVecTraits<double, Epetra_MultiVector>::SetBlock call to Epetra_MultiVector::Update() returned a nonzero value."); 
00701       }
00702       else {
00703         TEST_FOR_EXCEPTION( temp_vec.Update( 1.0, A, 0.0, A, 0.0 )!=0, EpetraMultiVecFailure,
00704             "MultiVecTraits<double, Epetra_MultiVector>::SetBlock call to Epetra_MultiVector::Update() returned a nonzero value.");
00705       }
00706     }
00707     
00710     static void MvScale ( Epetra_MultiVector& mv, double alpha ) 
00711     { 
00712       TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
00713           "MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value."); 
00714     }
00715     
00718     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
00719     { 
00720       // Check to make sure the vector is as long as the multivector has columns.
00721       int numvecs = mv.NumVectors();
00722       TEST_FOR_EXCEPTION( (int)alpha.size() != numvecs, std::invalid_argument,
00723                           "MultiVecTraits<double, Epetra_MultiVector>::MvScale(MV mv,vector alpha) alpha argument size was inconsistent with number of vectors in mv.")
00724 
00725       std::vector<int> tmp_index( 1, 0 );
00726       for (int i=0; i<numvecs; i++) {
00727         Epetra_MultiVector temp_vec(::View, mv, &tmp_index[0], 1);
00728         TEST_FOR_EXCEPTION( temp_vec.Scale( alpha[i] )!=0, EpetraMultiVecFailure,
00729             "MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00730         tmp_index[0]++;
00731       }
00732     }
00733 
00736     static void MvRandom( Epetra_MultiVector& mv )
00737     { 
00738       TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
00739           "MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00740     }
00741 
00744     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00745     { 
00746       TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
00747           "MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00748     }
00749     
00751 
00753 
00754 
00757     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
00758     { os << mv << std::endl; }
00759 
00761   };        
00762 
00764   //
00765   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
00766   //
00768 
00780   template <> 
00781   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
00782   {
00783   public:
00784     
00788     static void Apply ( const Epetra_Operator& Op, 
00789                         const Epetra_MultiVector& x, 
00790                         Epetra_MultiVector& y )
00791     { 
00792       TEST_FOR_EXCEPTION( Op.Apply( x, y ) != 0, OperatorError, "Error in Epetra_Operator::Apply()!" );
00793     }
00794     
00795   };
00796   
00797 } // end of Anasazi namespace 
00798 
00799 #endif 
00800 // end of file ANASAZI_EPETRA_ADAPTER_HPP

Generated on Tue Oct 20 12:45:18 2009 for Anasazi by doxygen 1.4.7