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_TestForException.hpp"
00042 #include "Teuchos_SerialDenseMatrix.hpp"
00043 #include "Epetra_MultiVector.h"
00044 #include "Epetra_Vector.h"
00045 #include "Epetra_Operator.h"
00046 #include "Epetra_Map.h"
00047 #include "Epetra_LocalMap.h"
00048 
00049 namespace Anasazi {
00050 
00052 
00053 
00057   class EpetraMultiVecFailure : public AnasaziError {public:
00058     EpetraMultiVecFailure(const std::string& what_arg) : AnasaziError(what_arg)
00059     {}};
00060 
00064   class EpetraOpFailure : public AnasaziError {public:
00065     EpetraOpFailure(const std::string& what_arg) : AnasaziError(what_arg)
00066     {}};
00067 
00069   
00071   //
00072   //--------template class AnasaziEpetraMultiVec-----------------
00073   //
00075   
00082   class EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector {
00083   public:
00085 
00086 
00088 
00093     EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs);
00094 
00096     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00097     
00099 
00107     EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride=0);
00108 
00110 
00116     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00117 
00119     virtual ~EpetraMultiVec() {};
00120 
00122 
00124 
00125 
00130     MultiVec<double> * Clone ( const int numvecs ) const;
00131 
00137     MultiVec<double> * CloneCopy () const;
00138 
00146     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00147     
00155     MultiVec<double> * CloneView ( const std::vector<int>& index );
00156 
00158 
00160 
00161 
00163     int GetNumberVecs () const { return NumVectors(); }
00164 
00166     int GetVecLength () const { return GlobalLength(); }
00167 
00169 
00171 
00172 
00174     void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00175                            const Teuchos::SerialDenseMatrix<int,double>& B, 
00176                            double beta );
00177 
00180     void MvAddMv ( double alpha, const MultiVec<double>& A, 
00181                    double beta, const MultiVec<double>& B);
00182 
00185     void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00186 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00187         , ConjType conj = Anasazi::CONJ
00188 #endif
00189         ) const;
00190   
00193     void MvDot ( const MultiVec<double>& A, std::vector<double> &b
00194 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00195         , ConjType conj = Anasazi::CONJ
00196 #endif
00197         ) const;
00198 
00201     void MvScale ( double alpha ) { 
00202       TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
00203           "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00204     }
00205     
00208     void MvScale ( const std::vector<double>& alpha );
00209 
00211 
00212 
00213     
00217     void MvNorm ( std::vector<double> & normvec ) const {
00218       if (((int)normvec.size() >= GetNumberVecs()) ) {
00219         TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00220             "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
00221       }
00222     };
00224     
00226 
00227 
00232     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00233 
00236     void MvRandom() { 
00237       TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
00238           "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00239     };
00240 
00243     void MvInit ( double alpha ) { 
00244       TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
00245           "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00246     };
00247     
00249 
00250 
00251 
00253     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00255 
00256   private:
00257   };
00258   //-------------------------------------------------------------
00259   
00261   //
00262   //--------template class AnasaziEpetraOp---------------------
00263   //
00265   
00272   class EpetraOp : public virtual Operator<double> {
00273   public:
00275 
00276     
00278     EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
00279     
00281     ~EpetraOp();
00283     
00285 
00286     
00290     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00292     
00293   private:
00294     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00295   };
00296   //-------------------------------------------------------------
00297 
00299   //
00300   //--------template class AnasaziEpetraGenOp--------------------
00301   //
00303   
00317   class EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00318   public:
00320 
00323     EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp, 
00324                 const Teuchos::RCP<Epetra_Operator> &MOp,
00325                 bool isAInverse = true );
00326 
00328     ~EpetraGenOp();
00329     
00331 
00333     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00334 
00336 
00338     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00339 
00341 
00343     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00344 
00346     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00347     
00349     bool UseTranspose() const { return (false); };
00350 
00352     int SetUseTranspose(bool UseTranspose_in) { return 0; };
00353     
00355     bool HasNormInf() const { return (false); };
00356     
00358     double NormInf() const  { return (-1.0); };
00359     
00361     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00362 
00364     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00365 
00367     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00368 
00369   private:
00370     bool isAInverse;
00371     Teuchos::RCP<Epetra_Operator> Epetra_AOp;
00372     Teuchos::RCP<Epetra_Operator> Epetra_MOp;
00373   };
00374   
00376   //
00377   //--------template class AnasaziEpetraSymOp--------------------
00378   //
00380 
00393   class EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00394   public:
00396 
00398     EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
00399 
00401     ~EpetraSymOp();
00402     
00404 
00406     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00407 
00409 
00411     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00412 
00414 
00417     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00418 
00420     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00421     
00423     bool UseTranspose() const { return (false); };
00424 
00426     int SetUseTranspose(bool UseTranspose_in) { return 0; };
00427     
00429     bool HasNormInf() const { return (false); };
00430     
00432     double NormInf() const  { return (-1.0); };
00433     
00435     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00436 
00438     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00439 
00441     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00442 
00443   private:
00444     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00445     bool isTrans_;
00446   };
00447 
00448 
00450   //
00451   //--------template class AnasaziEpetraSymMVOp---------------------
00452   //
00454 
00467   class EpetraSymMVOp : public virtual Operator<double> {
00468   public:
00470 
00472     EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00473                   bool isTrans = false );
00474     
00476     ~EpetraSymMVOp() {};
00477     
00479 
00481     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00482 
00483   private:
00484     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00485     Teuchos::RCP<const Epetra_Map> MV_localmap;
00486     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00487     bool isTrans_;
00488   };
00489 
00491   //
00492   //--------template class AnasaziEpetraWSymMVOp---------------------
00493   //
00495 
00508   class EpetraWSymMVOp : public virtual Operator<double> {
00509   public:
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 
00531   //
00532   //--------template class AnasaziEpetraW2SymMVOp---------------------
00533   //
00535 
00548   class EpetraW2SymMVOp : public virtual Operator<double> {
00549   public:
00551     EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00552                    const Teuchos::RCP<Epetra_Operator> &OP );
00553     
00555     ~EpetraW2SymMVOp() {};
00556     
00558 
00560     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00561 
00562   private:
00563     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00564     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00565     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00566     Teuchos::RCP<const Epetra_Map> MV_localmap;
00567     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00568   };
00569 
00570   
00572   //
00573   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00574   //
00576 
00587   template<>
00588   class MultiVecTraits<double, Epetra_MultiVector>
00589   {
00590   public:
00591 
00593 
00594 
00599     static Teuchos::RCP<Epetra_MultiVector> Clone( const Epetra_MultiVector& mv, const int numvecs )
00600     { 
00601 #ifdef TEUCHOS_DEBUG
00602       TEST_FOR_EXCEPTION(numvecs <= 0,std::invalid_argument,
00603           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,numvecs): numvecs must be greater than zero.");
00604 #endif
00605       return Teuchos::rcp( new Epetra_MultiVector(mv.Map(), numvecs) ); 
00606     }
00607 
00612     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv )
00613     { return Teuchos::rcp( new Epetra_MultiVector( mv ) ); }
00614 
00620     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv, const std::vector<int>& index )
00621     { 
00622 #ifdef TEUCHOS_DEBUG
00623       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00624           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): numvecs must be greater than zero.");
00625       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00626           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): indices must be >= zero.");
00627       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00628           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): indices must be < mv.NumVectors().");
00629 #endif
00630       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00631       return Teuchos::rcp( new Epetra_MultiVector(::Copy, mv, &tmp_index[0], index.size()) ); 
00632     }
00633 
00639     static Teuchos::RCP<Epetra_MultiVector> CloneView( Epetra_MultiVector& mv, const std::vector<int>& index )
00640     { 
00641 #ifdef TEUCHOS_DEBUG
00642       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00643           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
00644       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00645           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView(mv,index): indices must be >= zero.");
00646       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00647           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView(mv,index): indices must be < mv.NumVectors().");
00648 #endif
00649       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00650       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00651     }
00652 
00658     static Teuchos::RCP<const Epetra_MultiVector> CloneView( const Epetra_MultiVector& mv, const std::vector<int>& index )
00659     { 
00660 #ifdef TEUCHOS_DEBUG
00661       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00662           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): numvecs must be greater than zero.");
00663       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00664           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): indices must be >= zero.");
00665       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00666           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): indices must be < mv.NumVectors().");
00667 #endif
00668       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00669       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00670     }
00671 
00673 
00675 
00676 
00678     static int GetVecLength( const Epetra_MultiVector& mv )
00679     { return mv.GlobalLength(); }
00680 
00682     static int GetNumberVecs( const Epetra_MultiVector& mv )
00683     { return mv.NumVectors(); }
00685 
00687 
00688 
00691     static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A, 
00692                                  const Teuchos::SerialDenseMatrix<int,double>& B, 
00693                                  double beta, Epetra_MultiVector& mv )
00694     { 
00695       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00696       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00697 
00698       TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
00699           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00700     }
00701 
00704     static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
00705     { 
00706       // epetra mv.Update(alpha,A,beta,B,gamma) will check 
00707       //   alpha == 0.0 
00708       // and 
00709       //   beta == 0.0 
00710       // and based on this will call 
00711       //   mv.Update(beta,B,gamma) 
00712       // or 
00713       //   mv.Update(alpha,A,gamma)
00714       //
00715       // mv.Update(alpha,A,gamma) 
00716       // will then check for one of 
00717       //   gamma == 0
00718       // or 
00719       //   gamma == 1
00720       // or 
00721       //   alpha == 1 
00722       // in that order. however, it will not look for the combination
00723       //   alpha == 1 and gamma = 0
00724       // which is a common use case when we wish to assign 
00725       //   mv = A   (in which case alpha == 1, beta == gamma == 0)
00726       // or 
00727       //   mv = B   (in which case beta == 1, alpha == gamma == 0)
00728       // 
00729       // therefore, we will check for these use cases ourselves
00730       if (beta == 0.0) {
00731         if (alpha == 1.0) {
00732           // assign
00733           mv = A;
00734         }
00735         else {
00736           // single update
00737           TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
00738               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
00739         }
00740       }
00741       else if (alpha == 0.0) {
00742         if (beta == 1.0) {
00743           // assign 
00744           mv = B;
00745         }
00746         else {
00747           // single update
00748           TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
00749               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
00750         }
00751       }
00752       else {
00753         // double update
00754         TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
00755             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
00756       }
00757     }
00758 
00761     static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
00762 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00763                           , ConjType conj = Anasazi::CONJ
00764 #endif
00765                         )
00766     { 
00767       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00768       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00769       
00770       TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
00771           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00772     }
00773     
00776     static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
00777 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00778                       , ConjType conj = Anasazi::CONJ
00779 #endif
00780                       )
00781     {
00782 #ifdef TEUCHOS_DEBUG
00783       TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
00784           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
00785       TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
00786           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
00787 #endif
00788       TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
00789           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");     
00790     }
00791 
00793 
00794 
00795 
00799     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
00800     { 
00801 #ifdef TEUCHOS_DEBUG
00802       TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
00803           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
00804 #endif
00805       TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00806           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value."); 
00807     }
00808 
00810     
00812 
00813 
00815     static void SetBlock( const Epetra_MultiVector& A, const std::vector<int>& index, Epetra_MultiVector& mv )
00816     { 
00817       TEST_FOR_EXCEPTION((unsigned int)A.NumVectors() != index.size(),std::invalid_argument,
00818           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::SetBlock(A,index,mv): index must be the same size as A");
00819       Teuchos::RCP<Epetra_MultiVector> mv_sub = CloneView(mv,index);
00820       *mv_sub = A;
00821     }
00822 
00825     static void MvScale ( Epetra_MultiVector& mv, double alpha ) 
00826     { 
00827       TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
00828           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value."); 
00829     }
00830 
00833     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
00834     { 
00835       // Check to make sure the vector is as long as the multivector has columns.
00836       int numvecs = mv.NumVectors();
00837 #ifdef TEUCHOS_DEBUG
00838       TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
00839                           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
00840 #endif
00841       for (int i=0; i<numvecs; i++) {
00842         TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i]), EpetraMultiVecFailure,
00843             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00844       }
00845     }
00846 
00849     static void MvRandom( Epetra_MultiVector& mv )
00850     { 
00851       TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
00852           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00853     }
00854 
00857     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00858     { 
00859       TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
00860           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00861     }
00862 
00864 
00866 
00867 
00870     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
00871     { os << mv << std::endl; }
00872 
00874   };        
00875 
00877   //
00878   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
00879   //
00881 
00893   template <> 
00894   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
00895   {
00896   public:
00897     
00901     static void Apply ( const Epetra_Operator& Op, 
00902                         const Epetra_MultiVector& x, 
00903                         Epetra_MultiVector& y )
00904     { 
00905 #ifdef TEUCHOS_DEBUG
00906       TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
00907           "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
00908 #endif
00909       int ret = Op.Apply(x,y);
00910       TEST_FOR_EXCEPTION(ret != 0, OperatorError, 
00911           "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
00912     }
00913     
00914   };
00915   
00916 } // end of Anasazi namespace 
00917 
00918 #endif 
00919 // end of file ANASAZI_EPETRA_ADAPTER_HPP

Generated on Wed May 12 21:40:22 2010 for Anasazi by  doxygen 1.4.7