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 "Anasaziepetra_DLLExportMacro.h"
00038 #include "AnasaziTypes.hpp"
00039 #include "AnasaziMultiVec.hpp"
00040 #include "AnasaziOperator.hpp"
00041 
00042 #include "Teuchos_TestForException.hpp"
00043 #include "Teuchos_SerialDenseMatrix.hpp"
00044 #include "Epetra_MultiVector.h"
00045 #include "Epetra_Vector.h"
00046 #include "Epetra_Operator.h"
00047 #include "Epetra_Map.h"
00048 #include "Epetra_LocalMap.h"
00049 
00050 namespace Anasazi {
00051 
00053 
00054 
00058   class EpetraMultiVecFailure : public AnasaziError {public:
00059     EpetraMultiVecFailure(const std::string& what_arg) : AnasaziError(what_arg)
00060     {}};
00061 
00065   class EpetraOpFailure : public AnasaziError {public:
00066     EpetraOpFailure(const std::string& what_arg) : AnasaziError(what_arg)
00067     {}};
00068 
00070   
00072   //
00073   //--------template class AnasaziEpetraMultiVec-----------------
00074   //
00076   
00083   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector {
00084   public:
00086 
00087 
00089 
00094     EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs);
00095 
00097     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00098     
00100 
00108     EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride=0);
00109 
00111 
00117     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00118 
00120     virtual ~EpetraMultiVec() {};
00121 
00123 
00125 
00126 
00131     MultiVec<double> * Clone ( const int numvecs ) const;
00132 
00138     MultiVec<double> * CloneCopy () const;
00139 
00147     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00148     
00156     MultiVec<double> * CloneView ( const std::vector<int>& index );
00157 
00159 
00161 
00162 
00164     int GetNumberVecs () const { return NumVectors(); }
00165 
00167     int GetVecLength () const { return GlobalLength(); }
00168 
00170 
00172 
00173 
00175     void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00176                            const Teuchos::SerialDenseMatrix<int,double>& B, 
00177                            double beta );
00178 
00181     void MvAddMv ( double alpha, const MultiVec<double>& A, 
00182                    double beta, const MultiVec<double>& B);
00183 
00186     void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00187 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00188         , ConjType conj = Anasazi::CONJ
00189 #endif
00190         ) const;
00191   
00194     void MvDot ( const MultiVec<double>& A, std::vector<double> &b
00195 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00196         , ConjType conj = Anasazi::CONJ
00197 #endif
00198         ) const;
00199 
00202     void MvScale ( double alpha ) { 
00203       TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
00204           "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00205     }
00206     
00209     void MvScale ( const std::vector<double>& alpha );
00210 
00212 
00213 
00214     
00218     void MvNorm ( std::vector<double> & normvec ) const {
00219       if (((int)normvec.size() >= GetNumberVecs()) ) {
00220         TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00221             "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
00222       }
00223     };
00225     
00227 
00228 
00233     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00234 
00237     void MvRandom() { 
00238       TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
00239           "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00240     };
00241 
00244     void MvInit ( double alpha ) { 
00245       TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
00246           "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00247     };
00248     
00250 
00251 
00252 
00254     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00256 
00257   private:
00258   };
00259   //-------------------------------------------------------------
00260   
00262   //
00263   //--------template class AnasaziEpetraOp---------------------
00264   //
00266   
00273   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraOp : public virtual Operator<double> {
00274   public:
00276 
00277     
00279     EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
00280     
00282     ~EpetraOp();
00284     
00286 
00287     
00291     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00293     
00294   private:
00295 //use pragmas to disable some false-positive warnings for windows 
00296 // sharedlibs export
00297 #ifdef _MSC_VER
00298 #pragma warning(push)
00299 #pragma warning(disable:4251)
00300 #endif
00301     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00302 #ifdef _MSC_VER
00303 #pragma warning(pop)
00304 #endif
00305   };
00306   //-------------------------------------------------------------
00307 
00309   //
00310   //--------template class AnasaziEpetraGenOp--------------------
00311   //
00313   
00327   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00328   public:
00330 
00333     EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp, 
00334                 const Teuchos::RCP<Epetra_Operator> &MOp,
00335                 bool isAInverse = true );
00336 
00338     ~EpetraGenOp();
00339     
00341 
00343     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00344 
00346 
00348     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00349 
00351 
00353     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00354 
00356     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00357     
00359     bool UseTranspose() const { return (false); };
00360 
00362     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00363     
00365     bool HasNormInf() const { return (false); };
00366     
00368     double NormInf() const  { return (-1.0); };
00369     
00371     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00372 
00374     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00375 
00377     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00378 
00379   private:
00380     bool isAInverse;
00381 
00382 //use pragmas to disable some false-positive warnings for windows 
00383 // sharedlibs export
00384 #ifdef _MSC_VER
00385 #pragma warning(push)
00386 #pragma warning(disable:4251)
00387 #endif
00388     Teuchos::RCP<Epetra_Operator> Epetra_AOp;
00389     Teuchos::RCP<Epetra_Operator> Epetra_MOp;
00390 #ifdef _MSC_VER
00391 #pragma warning(pop)
00392 #endif
00393   };
00394   
00396   //
00397   //--------template class AnasaziEpetraSymOp--------------------
00398   //
00400 
00413   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00414   public:
00416 
00418     EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
00419 
00421     ~EpetraSymOp();
00422     
00424 
00426     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00427 
00429 
00431     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00432 
00434 
00437     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00438 
00440     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00441     
00443     bool UseTranspose() const { return (false); };
00444 
00446     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00447     
00449     bool HasNormInf() const { return (false); };
00450     
00452     double NormInf() const  { return (-1.0); };
00453     
00455     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00456 
00458     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00459 
00461     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00462 
00463   private:
00464 
00465 //use pragmas to disable false-positive warnings in generating windows sharedlib exports
00466 #ifdef _MSC_VER
00467 #pragma warning(push)
00468 #pragma warning(disable:4251)
00469 #endif
00470     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00471 #ifdef _MSC_VER
00472 #pragma warning(pop)
00473 #endif
00474 
00475     bool isTrans_;
00476   };
00477 
00478 
00480   //
00481   //--------template class AnasaziEpetraSymMVOp---------------------
00482   //
00484 
00497   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymMVOp : public virtual Operator<double> {
00498   public:
00500 
00502     EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00503                   bool isTrans = false );
00504     
00506     ~EpetraSymMVOp() {};
00507     
00509 
00511     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00512 
00513   private:
00514 
00515 //use pragmas to disable some false-positive warnings for windows 
00516 // sharedlibs export
00517 #ifdef _MSC_VER
00518 #pragma warning(push)
00519 #pragma warning(disable:4251)
00520 #endif
00521     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00522     Teuchos::RCP<const Epetra_Map> MV_localmap;
00523     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00524 #ifdef _MSC_VER
00525 #pragma warning(pop)
00526 #endif
00527 
00528     bool isTrans_;
00529   };
00530 
00532   //
00533   //--------template class AnasaziEpetraWSymMVOp---------------------
00534   //
00536 
00549   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraWSymMVOp : public virtual Operator<double> {
00550   public:
00552     EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00553                    const Teuchos::RCP<Epetra_Operator> &OP );
00554     
00556     ~EpetraWSymMVOp() {};
00557     
00559 
00561     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00562 
00563   private:
00564 //use pragmas to disable some false-positive warnings for windows 
00565 // sharedlibs export
00566 #ifdef _MSC_VER
00567 #pragma warning(push)
00568 #pragma warning(disable:4251)
00569 #endif
00570     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00571     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00572     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00573     Teuchos::RCP<const Epetra_Map> MV_localmap;
00574     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00575 #ifdef _MSC_VER
00576 #pragma warning(pop)
00577 #endif
00578   };
00579 
00581   //
00582   //--------template class AnasaziEpetraW2SymMVOp---------------------
00583   //
00585 
00598   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraW2SymMVOp : public virtual Operator<double> {
00599   public:
00601     EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00602                    const Teuchos::RCP<Epetra_Operator> &OP );
00603     
00605     ~EpetraW2SymMVOp() {};
00606     
00608 
00610     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00611 
00612   private:
00613 //use pragmas to disable some false-positive warnings for windows 
00614 // sharedlibs export
00615 #ifdef _MSC_VER
00616 #pragma warning(push)
00617 #pragma warning(disable:4251)
00618 #endif
00619     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00620     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00621     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00622     Teuchos::RCP<const Epetra_Map> MV_localmap;
00623     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00624 #ifdef _MSC_VER
00625 #pragma warning(pop)
00626 #endif
00627   };
00628 
00629   
00631   //
00632   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00633   //
00635 
00646   template<>
00647   class MultiVecTraits<double, Epetra_MultiVector>
00648   {
00649   public:
00650 
00652 
00653 
00658     static Teuchos::RCP<Epetra_MultiVector> Clone( const Epetra_MultiVector& mv, const int numvecs )
00659     { 
00660 #ifdef TEUCHOS_DEBUG
00661       TEST_FOR_EXCEPTION(numvecs <= 0,std::invalid_argument,
00662           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,numvecs): numvecs must be greater than zero.");
00663 #endif
00664       return Teuchos::rcp( new Epetra_MultiVector(mv.Map(), numvecs) ); 
00665     }
00666 
00671     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv )
00672     { return Teuchos::rcp( new Epetra_MultiVector( mv ) ); }
00673 
00679     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv, const std::vector<int>& index )
00680     { 
00681 #ifdef TEUCHOS_DEBUG
00682       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00683           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): numvecs must be greater than zero.");
00684       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00685           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): indices must be >= zero.");
00686       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00687           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): indices must be < mv.NumVectors().");
00688 #endif
00689       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00690       return Teuchos::rcp( new Epetra_MultiVector(::Copy, mv, &tmp_index[0], index.size()) ); 
00691     }
00692 
00698     static Teuchos::RCP<Epetra_MultiVector> CloneView( Epetra_MultiVector& mv, const std::vector<int>& index )
00699     { 
00700 #ifdef TEUCHOS_DEBUG
00701       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00702           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
00703       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00704           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView(mv,index): indices must be >= zero.");
00705       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00706           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView(mv,index): indices must be < mv.NumVectors().");
00707 #endif
00708       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00709       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00710     }
00711 
00717     static Teuchos::RCP<const Epetra_MultiVector> CloneView( const Epetra_MultiVector& mv, const std::vector<int>& index )
00718     { 
00719 #ifdef TEUCHOS_DEBUG
00720       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00721           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): numvecs must be greater than zero.");
00722       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00723           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): indices must be >= zero.");
00724       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00725           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): indices must be < mv.NumVectors().");
00726 #endif
00727       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00728       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00729     }
00730 
00732 
00734 
00735 
00737     static int GetVecLength( const Epetra_MultiVector& mv )
00738     { return mv.GlobalLength(); }
00739 
00741     static int GetNumberVecs( const Epetra_MultiVector& mv )
00742     { return mv.NumVectors(); }
00744 
00746 
00747 
00750     static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A, 
00751                                  const Teuchos::SerialDenseMatrix<int,double>& B, 
00752                                  double beta, Epetra_MultiVector& mv )
00753     { 
00754       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00755       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00756 
00757       TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
00758           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00759     }
00760 
00763     static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
00764     { 
00765       // epetra mv.Update(alpha,A,beta,B,gamma) will check 
00766       //   alpha == 0.0 
00767       // and 
00768       //   beta == 0.0 
00769       // and based on this will call 
00770       //   mv.Update(beta,B,gamma) 
00771       // or 
00772       //   mv.Update(alpha,A,gamma)
00773       //
00774       // mv.Update(alpha,A,gamma) 
00775       // will then check for one of 
00776       //   gamma == 0
00777       // or 
00778       //   gamma == 1
00779       // or 
00780       //   alpha == 1 
00781       // in that order. however, it will not look for the combination
00782       //   alpha == 1 and gamma = 0
00783       // which is a common use case when we wish to assign 
00784       //   mv = A   (in which case alpha == 1, beta == gamma == 0)
00785       // or 
00786       //   mv = B   (in which case beta == 1, alpha == gamma == 0)
00787       // 
00788       // therefore, we will check for these use cases ourselves
00789       if (beta == 0.0) {
00790         if (alpha == 1.0) {
00791           // assign
00792           mv = A;
00793         }
00794         else {
00795           // single update
00796           TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
00797               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
00798         }
00799       }
00800       else if (alpha == 0.0) {
00801         if (beta == 1.0) {
00802           // assign 
00803           mv = B;
00804         }
00805         else {
00806           // single update
00807           TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
00808               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
00809         }
00810       }
00811       else {
00812         // double update
00813         TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
00814             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
00815       }
00816     }
00817 
00820     static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
00821 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00822                           , ConjType conj = Anasazi::CONJ
00823 #endif
00824                         )
00825     { 
00826       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00827       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00828       
00829       TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
00830           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00831     }
00832     
00835     static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
00836 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00837                       , ConjType conj = Anasazi::CONJ
00838 #endif
00839                       )
00840     {
00841 #ifdef TEUCHOS_DEBUG
00842       TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
00843           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
00844       TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
00845           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
00846 #endif
00847       TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
00848           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");     
00849     }
00850 
00852 
00853 
00854 
00858     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
00859     { 
00860 #ifdef TEUCHOS_DEBUG
00861       TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
00862           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
00863 #endif
00864       TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00865           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value."); 
00866     }
00867 
00869     
00871 
00872 
00874     static void SetBlock( const Epetra_MultiVector& A, const std::vector<int>& index, Epetra_MultiVector& mv )
00875     { 
00876       TEST_FOR_EXCEPTION((unsigned int)A.NumVectors() != index.size(),std::invalid_argument,
00877           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::SetBlock(A,index,mv): index must be the same size as A");
00878       Teuchos::RCP<Epetra_MultiVector> mv_sub = CloneView(mv,index);
00879       *mv_sub = A;
00880     }
00881 
00884     static void MvScale ( Epetra_MultiVector& mv, double alpha ) 
00885     { 
00886       TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
00887           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value."); 
00888     }
00889 
00892     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
00893     { 
00894       // Check to make sure the vector is as long as the multivector has columns.
00895       int numvecs = mv.NumVectors();
00896 #ifdef TEUCHOS_DEBUG
00897       TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
00898                           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
00899 #endif
00900       for (int i=0; i<numvecs; i++) {
00901         TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i])!=0, EpetraMultiVecFailure,
00902             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00903       }
00904     }
00905 
00908     static void MvRandom( Epetra_MultiVector& mv )
00909     { 
00910       TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
00911           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00912     }
00913 
00916     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00917     { 
00918       TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
00919           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00920     }
00921 
00923 
00925 
00926 
00929     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
00930     { os << mv << std::endl; }
00931 
00933   };        
00934 
00936   //
00937   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
00938   //
00940 
00952   template <> 
00953   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
00954   {
00955   public:
00956     
00960     static void Apply ( const Epetra_Operator& Op, 
00961                         const Epetra_MultiVector& x, 
00962                         Epetra_MultiVector& y )
00963     { 
00964 #ifdef TEUCHOS_DEBUG
00965       TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
00966           "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
00967 #endif
00968       int ret = Op.Apply(x,y);
00969       TEST_FOR_EXCEPTION(ret != 0, OperatorError, 
00970           "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
00971     }
00972     
00973   };
00974   
00975 } // end of Anasazi namespace 
00976 
00977 #endif 
00978 // end of file ANASAZI_EPETRA_ADAPTER_HPP

Generated on Tue Jul 13 09:22:47 2010 for Anasazi by  doxygen 1.4.7