Anasazi Version of the Day
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> * CloneViewNonConst ( const std::vector<int>& index );
00157 
00165     const MultiVec<double> * CloneView ( const std::vector<int>& index ) const;
00166 
00168 
00170 
00171 
00173     int GetNumberVecs () const { return NumVectors(); }
00174 
00176     int GetVecLength () const { return GlobalLength(); }
00177 
00179 
00181 
00182 
00184     void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00185                            const Teuchos::SerialDenseMatrix<int,double>& B, 
00186                            double beta );
00187 
00190     void MvAddMv ( double alpha, const MultiVec<double>& A, 
00191                    double beta, const MultiVec<double>& B);
00192 
00195     void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00196 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00197         , ConjType conj = Anasazi::CONJ
00198 #endif
00199         ) const;
00200   
00203     void MvDot ( const MultiVec<double>& A, std::vector<double> &b
00204 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00205         , ConjType conj = Anasazi::CONJ
00206 #endif
00207         ) const;
00208 
00211     void MvScale ( double alpha ) { 
00212       TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
00213           "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00214     }
00215     
00218     void MvScale ( const std::vector<double>& alpha );
00219 
00221 
00222 
00223     
00227     void MvNorm ( std::vector<double> & normvec ) const {
00228       if (((int)normvec.size() >= GetNumberVecs()) ) {
00229         TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00230             "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
00231       }
00232     };
00234     
00236 
00237 
00242     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00243 
00246     void MvRandom() { 
00247       TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
00248           "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00249     };
00250 
00253     void MvInit ( double alpha ) { 
00254       TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
00255           "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00256     };
00257     
00259 
00260 
00261 
00263     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00265 
00266   private:
00267   };
00268   //-------------------------------------------------------------
00269   
00271   //
00272   //--------template class AnasaziEpetraOp---------------------
00273   //
00275   
00282   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraOp : public virtual Operator<double> {
00283   public:
00285 
00286     
00288     EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
00289     
00291     ~EpetraOp();
00293     
00295 
00296     
00300     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00302     
00303   private:
00304 //use pragmas to disable some false-positive warnings for windows 
00305 // sharedlibs export
00306 #ifdef _MSC_VER
00307 #pragma warning(push)
00308 #pragma warning(disable:4251)
00309 #endif
00310     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00311 #ifdef _MSC_VER
00312 #pragma warning(pop)
00313 #endif
00314   };
00315   //-------------------------------------------------------------
00316 
00318   //
00319   //--------template class AnasaziEpetraGenOp--------------------
00320   //
00322   
00336   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00337   public:
00339 
00342     EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp, 
00343                 const Teuchos::RCP<Epetra_Operator> &MOp,
00344                 bool isAInverse = true );
00345 
00347     ~EpetraGenOp();
00348     
00350 
00352     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00353 
00355 
00357     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00358 
00360 
00362     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00363 
00365     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00366     
00368     bool UseTranspose() const { return (false); };
00369 
00371     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00372     
00374     bool HasNormInf() const { return (false); };
00375     
00377     double NormInf() const  { return (-1.0); };
00378     
00380     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00381 
00383     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00384 
00386     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00387 
00388   private:
00389     bool isAInverse;
00390 
00391 //use pragmas to disable some false-positive warnings for windows 
00392 // sharedlibs export
00393 #ifdef _MSC_VER
00394 #pragma warning(push)
00395 #pragma warning(disable:4251)
00396 #endif
00397     Teuchos::RCP<Epetra_Operator> Epetra_AOp;
00398     Teuchos::RCP<Epetra_Operator> Epetra_MOp;
00399 #ifdef _MSC_VER
00400 #pragma warning(pop)
00401 #endif
00402   };
00403   
00405   //
00406   //--------template class AnasaziEpetraSymOp--------------------
00407   //
00409 
00422   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00423   public:
00425 
00427     EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
00428 
00430     ~EpetraSymOp();
00431     
00433 
00435     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00436 
00438 
00440     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00441 
00443 
00446     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00447 
00449     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00450     
00452     bool UseTranspose() const { return (false); };
00453 
00455     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00456     
00458     bool HasNormInf() const { return (false); };
00459     
00461     double NormInf() const  { return (-1.0); };
00462     
00464     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00465 
00467     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00468 
00470     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00471 
00472   private:
00473 
00474 //use pragmas to disable false-positive warnings in generating windows sharedlib exports
00475 #ifdef _MSC_VER
00476 #pragma warning(push)
00477 #pragma warning(disable:4251)
00478 #endif
00479     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00480 #ifdef _MSC_VER
00481 #pragma warning(pop)
00482 #endif
00483 
00484     bool isTrans_;
00485   };
00486 
00487 
00489   //
00490   //--------template class AnasaziEpetraSymMVOp---------------------
00491   //
00493 
00506   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymMVOp : public virtual Operator<double> {
00507   public:
00509 
00511     EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00512                   bool isTrans = false );
00513     
00515     ~EpetraSymMVOp() {};
00516     
00518 
00520     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00521 
00522   private:
00523 
00524 //use pragmas to disable some false-positive warnings for windows 
00525 // sharedlibs export
00526 #ifdef _MSC_VER
00527 #pragma warning(push)
00528 #pragma warning(disable:4251)
00529 #endif
00530     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00531     Teuchos::RCP<const Epetra_Map> MV_localmap;
00532     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00533 #ifdef _MSC_VER
00534 #pragma warning(pop)
00535 #endif
00536 
00537     bool isTrans_;
00538   };
00539 
00541   //
00542   //--------template class AnasaziEpetraWSymMVOp---------------------
00543   //
00545 
00558   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraWSymMVOp : public virtual Operator<double> {
00559   public:
00561     EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00562                    const Teuchos::RCP<Epetra_Operator> &OP );
00563     
00565     ~EpetraWSymMVOp() {};
00566     
00568 
00570     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00571 
00572   private:
00573 //use pragmas to disable some false-positive warnings for windows 
00574 // sharedlibs export
00575 #ifdef _MSC_VER
00576 #pragma warning(push)
00577 #pragma warning(disable:4251)
00578 #endif
00579     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00580     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00581     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00582     Teuchos::RCP<const Epetra_Map> MV_localmap;
00583     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00584 #ifdef _MSC_VER
00585 #pragma warning(pop)
00586 #endif
00587   };
00588 
00590   //
00591   //--------template class AnasaziEpetraW2SymMVOp---------------------
00592   //
00594 
00607   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraW2SymMVOp : public virtual Operator<double> {
00608   public:
00610     EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00611                    const Teuchos::RCP<Epetra_Operator> &OP );
00612     
00614     ~EpetraW2SymMVOp() {};
00615     
00617 
00619     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00620 
00621   private:
00622 //use pragmas to disable some false-positive warnings for windows 
00623 // sharedlibs export
00624 #ifdef _MSC_VER
00625 #pragma warning(push)
00626 #pragma warning(disable:4251)
00627 #endif
00628     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00629     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00630     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00631     Teuchos::RCP<const Epetra_Map> MV_localmap;
00632     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00633 #ifdef _MSC_VER
00634 #pragma warning(pop)
00635 #endif
00636   };
00637 
00638   
00640   //
00641   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00642   //
00644 
00655   template<>
00656   class MultiVecTraits<double, Epetra_MultiVector>
00657   {
00658   public:
00659 
00661 
00662 
00667     static Teuchos::RCP<Epetra_MultiVector> Clone( const Epetra_MultiVector& mv, const int numvecs )
00668     { 
00669 #ifdef TEUCHOS_DEBUG
00670       TEST_FOR_EXCEPTION(numvecs <= 0,std::invalid_argument,
00671           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,numvecs): numvecs must be greater than zero.");
00672 #endif
00673       return Teuchos::rcp( new Epetra_MultiVector(mv.Map(), numvecs) ); 
00674     }
00675 
00680     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv )
00681     { return Teuchos::rcp( new Epetra_MultiVector( mv ) ); }
00682 
00688     static Teuchos::RCP<Epetra_MultiVector> CloneCopy( const Epetra_MultiVector& mv, const std::vector<int>& index )
00689     { 
00690 #ifdef TEUCHOS_DEBUG
00691       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00692           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): numvecs must be greater than zero.");
00693       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00694           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): indices must be >= zero.");
00695       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00696           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,index): indices must be < mv.NumVectors().");
00697 #endif
00698       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00699       return Teuchos::rcp( new Epetra_MultiVector(::Copy, mv, &tmp_index[0], index.size()) ); 
00700     }
00701 
00707     static Teuchos::RCP<Epetra_MultiVector> CloneViewNonConst( Epetra_MultiVector& mv, const std::vector<int>& index )
00708     { 
00709 #ifdef TEUCHOS_DEBUG
00710       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00711           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneViewNonConst(mv,index): numvecs must be greater than zero.");
00712       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00713           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneViewNonConst(mv,index): indices must be >= zero.");
00714       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00715           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneViewNonConst(mv,index): indices must be < mv.NumVectors().");
00716 #endif
00717       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00718       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00719     }
00720 
00726     static Teuchos::RCP<const Epetra_MultiVector> CloneView( const Epetra_MultiVector& mv, const std::vector<int>& index )
00727     { 
00728 #ifdef TEUCHOS_DEBUG
00729       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00730           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): numvecs must be greater than zero.");
00731       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00732           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): indices must be >= zero.");
00733       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.NumVectors(), std::invalid_argument,
00734           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(const mv,index): indices must be < mv.NumVectors().");
00735 #endif
00736       std::vector<int>& tmp_index = const_cast<std::vector<int> &>( index );
00737       return Teuchos::rcp( new Epetra_MultiVector(::View, mv, &tmp_index[0], index.size()) ); 
00738     }
00739 
00741 
00743 
00744 
00746     static int GetVecLength( const Epetra_MultiVector& mv )
00747     { return mv.GlobalLength(); }
00748 
00750     static int GetNumberVecs( const Epetra_MultiVector& mv )
00751     { return mv.NumVectors(); }
00753 
00755 
00756 
00759     static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A, 
00760                                  const Teuchos::SerialDenseMatrix<int,double>& B, 
00761                                  double beta, Epetra_MultiVector& mv )
00762     { 
00763       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00764       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00765 
00766       TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
00767           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00768     }
00769 
00772     static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
00773     { 
00774       // epetra mv.Update(alpha,A,beta,B,gamma) will check 
00775       //   alpha == 0.0 
00776       // and 
00777       //   beta == 0.0 
00778       // and based on this will call 
00779       //   mv.Update(beta,B,gamma) 
00780       // or 
00781       //   mv.Update(alpha,A,gamma)
00782       //
00783       // mv.Update(alpha,A,gamma) 
00784       // will then check for one of 
00785       //   gamma == 0
00786       // or 
00787       //   gamma == 1
00788       // or 
00789       //   alpha == 1 
00790       // in that order. however, it will not look for the combination
00791       //   alpha == 1 and gamma = 0
00792       // which is a common use case when we wish to assign 
00793       //   mv = A   (in which case alpha == 1, beta == gamma == 0)
00794       // or 
00795       //   mv = B   (in which case beta == 1, alpha == gamma == 0)
00796       // 
00797       // therefore, we will check for these use cases ourselves
00798       if (beta == 0.0) {
00799         if (alpha == 1.0) {
00800           // assign
00801           mv = A;
00802         }
00803         else {
00804           // single update
00805           TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
00806               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
00807         }
00808       }
00809       else if (alpha == 0.0) {
00810         if (beta == 1.0) {
00811           // assign 
00812           mv = B;
00813         }
00814         else {
00815           // single update
00816           TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
00817               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
00818         }
00819       }
00820       else {
00821         // double update
00822         TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
00823             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
00824       }
00825     }
00826 
00829     static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
00830 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00831                           , ConjType conj = Anasazi::CONJ
00832 #endif
00833                         )
00834     { 
00835       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00836       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00837       
00838       TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
00839           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
00840     }
00841     
00844     static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
00845 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00846                       , ConjType conj = Anasazi::CONJ
00847 #endif
00848                       )
00849     {
00850 #ifdef TEUCHOS_DEBUG
00851       TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
00852           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
00853       TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
00854           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
00855 #endif
00856       TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
00857           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");     
00858     }
00859 
00861 
00862 
00863 
00867     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
00868     { 
00869 #ifdef TEUCHOS_DEBUG
00870       TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
00871           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
00872 #endif
00873       TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00874           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value."); 
00875     }
00876 
00878     
00880 
00881 
00883     static void SetBlock( const Epetra_MultiVector& A, const std::vector<int>& index, Epetra_MultiVector& mv )
00884     { 
00885       TEST_FOR_EXCEPTION((unsigned int)A.NumVectors() != index.size(),std::invalid_argument,
00886           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::SetBlock(A,index,mv): index must be the same size as A");
00887       Teuchos::RCP<Epetra_MultiVector> mv_sub = CloneViewNonConst(mv,index);
00888       *mv_sub = A;
00889     }
00890 
00893     static void MvScale ( Epetra_MultiVector& mv, double alpha ) 
00894     { 
00895       TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
00896           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value."); 
00897     }
00898 
00901     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
00902     { 
00903       // Check to make sure the vector is as long as the multivector has columns.
00904       int numvecs = mv.NumVectors();
00905 #ifdef TEUCHOS_DEBUG
00906       TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
00907                           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
00908 #endif
00909       for (int i=0; i<numvecs; i++) {
00910         TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i])!=0, EpetraMultiVecFailure,
00911             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00912       }
00913     }
00914 
00917     static void MvRandom( Epetra_MultiVector& mv )
00918     { 
00919       TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
00920           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00921     }
00922 
00925     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
00926     { 
00927       TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
00928           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00929     }
00930 
00932 
00934 
00935 
00938     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
00939     { os << mv << std::endl; }
00940 
00942   };        
00943 
00945   //
00946   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
00947   //
00949 
00961   template <> 
00962   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
00963   {
00964   public:
00965     
00969     static void Apply ( const Epetra_Operator& Op, 
00970                         const Epetra_MultiVector& x, 
00971                         Epetra_MultiVector& y )
00972     { 
00973 #ifdef TEUCHOS_DEBUG
00974       TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
00975           "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
00976 #endif
00977       int ret = Op.Apply(x,y);
00978       TEST_FOR_EXCEPTION(ret != 0, OperatorError, 
00979           "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
00980     }
00981     
00982   };
00983   
00984 } // end of Anasazi namespace 
00985 
00986 #endif 
00987 // end of file ANASAZI_EPETRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends