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_Assert.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 #if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
00051 #  include <Tpetra_ConfigDefs.hpp> // HAVE_TPETRA_EPETRA 
00052 #  if defined(HAVE_TPETRA_EPETRA)
00053 #    include <Epetra_TsqrAdaptor.hpp>
00054 #  endif // defined(HAVE_TPETRA_EPETRA)
00055 #endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
00056 
00057 namespace Anasazi {
00058 
00060 
00061 
00065   class EpetraMultiVecFailure : public AnasaziError {public:
00066     EpetraMultiVecFailure(const std::string& what_arg) : AnasaziError(what_arg)
00067     {}};
00068 
00072   class EpetraOpFailure : public AnasaziError {public:
00073     EpetraOpFailure(const std::string& what_arg) : AnasaziError(what_arg)
00074     {}};
00075 
00077 
00079 
00080 
00084   class EpetraMultiVecAccessor {
00085  
00086   public:
00088     virtual Epetra_MultiVector* GetEpetraMultiVec() { return 0; }
00089 
00091     virtual const Epetra_MultiVector* GetEpetraMultiVec() const { return 0; }
00092   };
00093 
00095   
00097   //
00098   //--------template class AnasaziEpetraMultiVec-----------------
00099   //
00101   
00108   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraMultiVec : public MultiVec<double>,  public Epetra_MultiVector, public EpetraMultiVecAccessor {
00109   public:
00111 
00112 
00114 
00119     EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs);
00120 
00122     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00123     
00125 
00133     EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride=0);
00134 
00136 
00142     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00143 
00145     virtual ~EpetraMultiVec() {};
00146 
00148 
00150 
00151 
00156     MultiVec<double> * Clone ( const int numvecs ) const;
00157 
00163     MultiVec<double> * CloneCopy () const;
00164 
00172     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00173     
00181     MultiVec<double> * CloneViewNonConst ( const std::vector<int>& index );
00182 
00190     const MultiVec<double> * CloneView ( const std::vector<int>& index ) const;
00191 
00193 
00196     int GetNumberVecs () const { return NumVectors(); }
00197 
00199     int GetVecLength () const { return GlobalLength(); }
00200 
00202 
00204 
00205 
00207     void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00208                            const Teuchos::SerialDenseMatrix<int,double>& B, 
00209                            double beta );
00210 
00213     void MvAddMv ( double alpha, const MultiVec<double>& A, 
00214                    double beta, const MultiVec<double>& B);
00215 
00218     void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00219 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00220         , ConjType conj = Anasazi::CONJ
00221 #endif
00222         ) const;
00223   
00226     void MvDot ( const MultiVec<double>& A, std::vector<double> &b
00227 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00228         , ConjType conj = Anasazi::CONJ
00229 #endif
00230         ) const;
00231 
00234     void MvScale ( double alpha ) { 
00235       TEUCHOS_TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
00236           "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00237     }
00238     
00241     void MvScale ( const std::vector<double>& alpha );
00242 
00244 
00245 
00246     
00250     void MvNorm ( std::vector<double> & normvec ) const {
00251       if (((int)normvec.size() >= GetNumberVecs()) ) {
00252         TEUCHOS_TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00253             "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
00254       }
00255     };
00257     
00259 
00260 
00265     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00266 
00269     void MvRandom() { 
00270       TEUCHOS_TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
00271           "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00272     };
00273 
00276     void MvInit ( double alpha ) { 
00277       TEUCHOS_TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
00278           "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00279     };
00280    
00282 
00283    
00285     Epetra_MultiVector* GetEpetraMultiVec() { return this; };
00286  
00288     const Epetra_MultiVector* GetEpetraMultiVec() const { return this; };
00289 
00291  
00293 
00294 
00295 
00297     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00299 
00300   private:
00301   };
00302   //-------------------------------------------------------------
00303   
00305   //
00306   //--------template class AnasaziEpetraOp---------------------
00307   //
00309   
00316   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraOp : public virtual Operator<double> {
00317   public:
00319 
00320     
00322     EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
00323     
00325     ~EpetraOp();
00327     
00329 
00330     
00334     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00336     
00337   private:
00338 //use pragmas to disable some false-positive warnings for windows 
00339 // sharedlibs export
00340 #ifdef _MSC_VER
00341 #pragma warning(push)
00342 #pragma warning(disable:4251)
00343 #endif
00344     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00345 #ifdef _MSC_VER
00346 #pragma warning(pop)
00347 #endif
00348   };
00349   //-------------------------------------------------------------
00350 
00352   //
00353   //--------template class AnasaziEpetraGenOp--------------------
00354   //
00356   
00370   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00371   public:
00373 
00376     EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp, 
00377                 const Teuchos::RCP<Epetra_Operator> &MOp,
00378                 bool isAInverse = true );
00379 
00381     ~EpetraGenOp();
00382     
00384 
00386     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00387 
00389 
00391     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00392 
00394 
00396     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00397 
00399     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00400     
00402     bool UseTranspose() const { return (false); };
00403 
00405     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00406     
00408     bool HasNormInf() const { return (false); };
00409     
00411     double NormInf() const  { return (-1.0); };
00412     
00414     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00415 
00417     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00418 
00420     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00421 
00422   private:
00423     bool isAInverse;
00424 
00425 //use pragmas to disable some false-positive warnings for windows 
00426 // sharedlibs export
00427 #ifdef _MSC_VER
00428 #pragma warning(push)
00429 #pragma warning(disable:4251)
00430 #endif
00431     Teuchos::RCP<Epetra_Operator> Epetra_AOp;
00432     Teuchos::RCP<Epetra_Operator> Epetra_MOp;
00433 #ifdef _MSC_VER
00434 #pragma warning(pop)
00435 #endif
00436   };
00437   
00439   //
00440   //--------template class AnasaziEpetraSymOp--------------------
00441   //
00443 
00456   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00457   public:
00459 
00461     EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
00462 
00464     ~EpetraSymOp();
00465     
00467 
00469     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00470 
00472 
00474     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00475 
00477 
00480     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00481 
00483     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00484     
00486     bool UseTranspose() const { return (false); };
00487 
00489     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00490     
00492     bool HasNormInf() const { return (false); };
00493     
00495     double NormInf() const  { return (-1.0); };
00496     
00498     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00499 
00501     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00502 
00504     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00505 
00506   private:
00507 
00508 //use pragmas to disable false-positive warnings in generating windows sharedlib exports
00509 #ifdef _MSC_VER
00510 #pragma warning(push)
00511 #pragma warning(disable:4251)
00512 #endif
00513     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00514 #ifdef _MSC_VER
00515 #pragma warning(pop)
00516 #endif
00517 
00518     bool isTrans_;
00519   };
00520 
00521 
00523   //
00524   //--------template class AnasaziEpetraSymMVOp---------------------
00525   //
00527 
00540   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymMVOp : public virtual Operator<double> {
00541   public:
00543 
00545     EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00546                   bool isTrans = false );
00547     
00549     ~EpetraSymMVOp() {};
00550     
00552 
00554     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00555 
00556   private:
00557 
00558 //use pragmas to disable some false-positive warnings for windows 
00559 // sharedlibs export
00560 #ifdef _MSC_VER
00561 #pragma warning(push)
00562 #pragma warning(disable:4251)
00563 #endif
00564     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00565     Teuchos::RCP<const Epetra_Map> MV_localmap;
00566     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00567 #ifdef _MSC_VER
00568 #pragma warning(pop)
00569 #endif
00570 
00571     bool isTrans_;
00572   };
00573 
00575   //
00576   //--------template class AnasaziEpetraWSymMVOp---------------------
00577   //
00579 
00592   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraWSymMVOp : public virtual Operator<double> {
00593   public:
00595     EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00596                    const Teuchos::RCP<Epetra_Operator> &OP );
00597     
00599     ~EpetraWSymMVOp() {};
00600     
00602 
00604     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00605 
00606   private:
00607 //use pragmas to disable some false-positive warnings for windows 
00608 // sharedlibs export
00609 #ifdef _MSC_VER
00610 #pragma warning(push)
00611 #pragma warning(disable:4251)
00612 #endif
00613     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00614     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00615     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00616     Teuchos::RCP<const Epetra_Map> MV_localmap;
00617     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00618 #ifdef _MSC_VER
00619 #pragma warning(pop)
00620 #endif
00621   };
00622 
00624   //
00625   //--------template class AnasaziEpetraW2SymMVOp---------------------
00626   //
00628 
00641   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraW2SymMVOp : public virtual Operator<double> {
00642   public:
00644     EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00645                    const Teuchos::RCP<Epetra_Operator> &OP );
00646     
00648     ~EpetraW2SymMVOp() {};
00649     
00651 
00653     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00654 
00655   private:
00656 //use pragmas to disable some false-positive warnings for windows 
00657 // sharedlibs export
00658 #ifdef _MSC_VER
00659 #pragma warning(push)
00660 #pragma warning(disable:4251)
00661 #endif
00662     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00663     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00664     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00665     Teuchos::RCP<const Epetra_Map> MV_localmap;
00666     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00667 #ifdef _MSC_VER
00668 #pragma warning(pop)
00669 #endif
00670   };
00671 
00672   
00674   //
00675   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00676   //
00678 
00689   template<>
00690   class MultiVecTraits<double, Epetra_MultiVector>
00691   {
00692   public:
00693 
00695 
00696 
00701     static Teuchos::RCP<Epetra_MultiVector> 
00702     Clone (const Epetra_MultiVector& mv, const int outNumVecs)
00703     { 
00704       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs <= 0, std::invalid_argument,
00705        "Belos::MultiVecTraits<double, Epetra_MultiVector>::"
00706        "Clone(mv, outNumVecs = " << outNumVecs << "): "
00707        "outNumVecs must be positive.");
00708       // FIXME (mfh 13 Jan 2011) Anasazi currently lets Epetra fill in
00709       // the entries of the returned multivector with zeros, but Belos
00710       // does not.  We retain this different behavior for now, but the
00711       // two versions will need to be reconciled.
00712       return Teuchos::rcp (new Epetra_MultiVector (mv.Map(), outNumVecs)); 
00713     }
00714 
00719     static Teuchos::RCP<Epetra_MultiVector> 
00720     CloneCopy (const Epetra_MultiVector& mv)
00721     { 
00722       return Teuchos::rcp (new Epetra_MultiVector (mv)); 
00723     }
00724 
00730     static Teuchos::RCP<Epetra_MultiVector> 
00731     CloneCopy (const Epetra_MultiVector& mv, const std::vector<int>& index)
00732     { 
00733       const int inNumVecs = GetNumberVecs (mv);
00734       const int outNumVecs = index.size();
00735 
00736       // Simple, inexpensive tests of the index vector.
00737       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00738        "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00739        "CloneCopy(mv, index = {}): At least one vector must be"
00740        " cloned from mv.");
00741       if (outNumVecs > inNumVecs)
00742   {
00743     std::ostringstream os;
00744     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00745       "CloneCopy(mv, index = {";
00746     for (int k = 0; k < outNumVecs - 1; ++k)
00747       os << index[k] << ", ";
00748     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00749        << " indices to copy, but only " << inNumVecs << " columns of mv.";
00750     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00751   }
00752 #ifdef TEUCHOS_DEBUG
00753       // In debug mode, we perform more expensive tests of the index
00754       // vector, to ensure all the elements are in range.
00755       // Dereferencing the iterator is valid because index has length
00756       // > 0.
00757       const int minIndex = *std::min_element (index.begin(), index.end());
00758       const int maxIndex = *std::max_element (index.begin(), index.end());
00759 
00760       if (minIndex < 0)
00761   {
00762     std::ostringstream os;
00763     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00764       "CloneCopy(mv, index = {";
00765     for (int k = 0; k < outNumVecs - 1; ++k)
00766       os << index[k] << ", ";
00767     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00768       "the smallest index " << minIndex << " is negative.";
00769     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00770   }
00771       if (maxIndex >= inNumVecs)
00772   {
00773     std::ostringstream os;
00774     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00775       "CloneCopy(mv, index = {";
00776     for (int k = 0; k < outNumVecs - 1; ++k)
00777       os << index[k] << ", ";
00778     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00779       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00780        << maxIndex << " is out of bounds.";
00781     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00782   }
00783 #endif // TEUCHOS_DEBUG
00784       // Cast to nonconst, because Epetra_MultiVector's constructor
00785       // wants a nonconst int array argument.  It doesn't actually
00786       // change the entries of the array.
00787       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00788       return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, &tmpind[0], index.size())); 
00789       // return Teuchos::rcp (new Epetra_MultiVector (::Copy, mv, &tmpind[0], index.size())); 
00790     }
00791 
00792     static Teuchos::RCP<Epetra_MultiVector> 
00793     CloneCopy (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00794     { 
00795       const int inNumVecs = GetNumberVecs (mv);
00796       const int outNumVecs = index.size();
00797       const bool validRange = outNumVecs > 0 && index.lbound() >= 0 && 
00798   index.ubound() < inNumVecs;
00799       if (! validRange)
00800   {
00801     std::ostringstream os;
00802     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,"
00803       "index=[" << index.lbound() << ", " << index.ubound() << "]): ";
00804     TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00805            os.str() << "Column index range must be nonempty.");
00806     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00807            os.str() << "Column index range must be nonnegative.");
00808     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= inNumVecs, std::invalid_argument,
00809            os.str() << "Column index range must not exceed "
00810            "number of vectors " << inNumVecs << " in the "
00811            "input multivector.");
00812   }
00813       return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, index.lbound(), index.size()));
00814     }
00815 
00821     static Teuchos::RCP<Epetra_MultiVector> 
00822     CloneViewNonConst (Epetra_MultiVector& mv, const std::vector<int>& index)
00823     { 
00824       const int inNumVecs = GetNumberVecs (mv);
00825       const int outNumVecs = index.size();
00826 
00827       // Simple, inexpensive tests of the index vector.
00828       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00829        "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00830        "CloneViewNonConst(mv, index = {}): The output view "
00831        "must have at least one column.");
00832       if (outNumVecs > inNumVecs)
00833   {
00834     std::ostringstream os;
00835     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00836       "CloneViewNonConst(mv, index = {";
00837     for (int k = 0; k < outNumVecs - 1; ++k)
00838       os << index[k] << ", ";
00839     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00840        << " indices to view, but only " << inNumVecs << " columns of mv.";
00841     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00842   }
00843 #ifdef TEUCHOS_DEBUG
00844       // In debug mode, we perform more expensive tests of the index
00845       // vector, to ensure all the elements are in range.
00846       // Dereferencing the iterator is valid because index has length
00847       // > 0.
00848       const int minIndex = *std::min_element (index.begin(), index.end());
00849       const int maxIndex = *std::max_element (index.begin(), index.end());
00850 
00851       if (minIndex < 0)
00852   {
00853     std::ostringstream os;
00854     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00855       "CloneViewNonConst(mv, index = {";
00856     for (int k = 0; k < outNumVecs - 1; ++k)
00857       os << index[k] << ", ";
00858     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00859       "the smallest index " << minIndex << " is negative.";
00860     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00861   }
00862       if (maxIndex >= inNumVecs)
00863   {
00864     std::ostringstream os;
00865     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00866       "CloneViewNonConst(mv, index = {";
00867     for (int k = 0; k < outNumVecs - 1; ++k)
00868       os << index[k] << ", ";
00869     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00870       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00871        << maxIndex << " is out of bounds.";
00872     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00873   }
00874 #endif // TEUCHOS_DEBUG
00875       // Cast to nonconst, because Epetra_MultiVector's constructor
00876       // wants a nonconst int array argument.  It doesn't actually
00877       // change the entries of the array.
00878       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00879       return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
00880     }
00881 
00882     static Teuchos::RCP<Epetra_MultiVector> 
00883     CloneViewNonConst (Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00884     { 
00885       const bool validRange = index.size() > 0 && 
00886   index.lbound() >= 0 && 
00887   index.ubound() < mv.NumVectors();
00888       if (! validRange)
00889   {
00890     std::ostringstream os;
00891     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
00892       "NonConst(mv,index=[" << index.lbound() << ", " << index.ubound() 
00893        << "]): ";
00894     TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00895            os.str() << "Column index range must be nonempty.");
00896     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00897            os.str() << "Column index range must be nonnegative.");
00898     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(), 
00899            std::invalid_argument,
00900            os.str() << "Column index range must not exceed "
00901            "number of vectors " << mv.NumVectors() << " in "
00902            "the input multivector.");
00903   }
00904       return Teuchos::rcp (new Epetra_MultiVector (View, mv, index.lbound(), index.size()));
00905     }
00906 
00912     static Teuchos::RCP<const Epetra_MultiVector> 
00913     CloneView (const Epetra_MultiVector& mv, const std::vector<int>& index)
00914     { 
00915       const int inNumVecs = GetNumberVecs (mv);
00916       const int outNumVecs = index.size();
00917 
00918       // Simple, inexpensive tests of the index vector.
00919       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00920        "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00921        "CloneView(mv, index = {}): The output view "
00922        "must have at least one column.");
00923       if (outNumVecs > inNumVecs)
00924   {
00925     std::ostringstream os;
00926     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00927       "CloneView(mv, index = {";
00928     for (int k = 0; k < outNumVecs - 1; ++k)
00929       os << index[k] << ", ";
00930     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00931        << " indices to view, but only " << inNumVecs << " columns of mv.";
00932     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00933   }
00934 #ifdef TEUCHOS_DEBUG
00935       // In debug mode, we perform more expensive tests of the index
00936       // vector, to ensure all the elements are in range.
00937       // Dereferencing the iterator is valid because index has length
00938       // > 0.
00939       const int minIndex = *std::min_element (index.begin(), index.end());
00940       const int maxIndex = *std::max_element (index.begin(), index.end());
00941 
00942       if (minIndex < 0)
00943   {
00944     std::ostringstream os;
00945     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00946       "CloneView(mv, index = {";
00947     for (int k = 0; k < outNumVecs - 1; ++k)
00948       os << index[k] << ", ";
00949     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00950       "the smallest index " << minIndex << " is negative.";
00951     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00952   }
00953       if (maxIndex >= inNumVecs)
00954   {
00955     std::ostringstream os;
00956     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00957       "CloneView(mv, index = {";
00958     for (int k = 0; k < outNumVecs - 1; ++k)
00959       os << index[k] << ", ";
00960     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00961       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00962        << maxIndex << " is out of bounds.";
00963     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00964   }
00965 #endif // TEUCHOS_DEBUG
00966       // Cast to nonconst, because Epetra_MultiVector's constructor
00967       // wants a nonconst int array argument.  It doesn't actually
00968       // change the entries of the array.
00969       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00970       return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
00971     }
00972 
00973     static Teuchos::RCP<Epetra_MultiVector> 
00974     CloneView (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00975     { 
00976       const bool validRange = index.size() > 0 && 
00977   index.lbound() >= 0 && 
00978   index.ubound() < mv.NumVectors();
00979       if (! validRange)
00980   {
00981     std::ostringstream os;
00982     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
00983       "(mv,index=[" << index.lbound() << ", " << index.ubound() 
00984        << "]): ";
00985     TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00986            os.str() << "Column index range must be nonempty.");
00987     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00988            os.str() << "Column index range must be nonnegative.");
00989     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(), 
00990            std::invalid_argument,
00991            os.str() << "Column index range must not exceed "
00992            "number of vectors " << mv.NumVectors() << " in "
00993            "the input multivector.");
00994   }
00995       return Teuchos::rcp (new Epetra_MultiVector(View, mv, index.lbound(), index.size()));
00996     }
00997 
00999 
01001 
01002 
01004     static int GetVecLength( const Epetra_MultiVector& mv )
01005     { return mv.GlobalLength(); }
01006 
01008     static int GetNumberVecs( const Epetra_MultiVector& mv )
01009     { return mv.NumVectors(); }
01010 
01011     static bool HasConstantStride( const Epetra_MultiVector& mv )
01012     { return mv.ConstantStride(); }
01014 
01016 
01017 
01020     static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A, 
01021                                  const Teuchos::SerialDenseMatrix<int,double>& B, 
01022                                  double beta, Epetra_MultiVector& mv )
01023     { 
01024       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
01025       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
01026 
01027       TEUCHOS_TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
01028           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Multiply() returned a nonzero value.");
01029     }
01030 
01033     static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
01034     { 
01035       // epetra mv.Update(alpha,A,beta,B,gamma) will check 
01036       //   alpha == 0.0 
01037       // and 
01038       //   beta == 0.0 
01039       // and based on this will call 
01040       //   mv.Update(beta,B,gamma) 
01041       // or 
01042       //   mv.Update(alpha,A,gamma)
01043       //
01044       // mv.Update(alpha,A,gamma) 
01045       // will then check for one of 
01046       //   gamma == 0
01047       // or 
01048       //   gamma == 1
01049       // or 
01050       //   alpha == 1 
01051       // in that order. however, it will not look for the combination
01052       //   alpha == 1 and gamma = 0
01053       // which is a common use case when we wish to assign 
01054       //   mv = A   (in which case alpha == 1, beta == gamma == 0)
01055       // or 
01056       //   mv = B   (in which case beta == 1, alpha == gamma == 0)
01057       // 
01058       // therefore, we will check for these use cases ourselves
01059       if (beta == 0.0) {
01060         if (alpha == 1.0) {
01061           // assign
01062           mv = A;
01063         }
01064         else {
01065           // single update
01066           TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
01067               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
01068         }
01069       }
01070       else if (alpha == 0.0) {
01071         if (beta == 1.0) {
01072           // assign 
01073           mv = B;
01074         }
01075         else {
01076           // single update
01077           TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
01078               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
01079         }
01080       }
01081       else {
01082         // double update
01083         TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
01084             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
01085       }
01086     }
01087 
01090     static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
01091 #ifdef HAVE_ANASAZI_EXPERIMENTAL
01092                           , ConjType conj = Anasazi::CONJ
01093 #endif
01094                         )
01095     { 
01096       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
01097       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
01098       
01099       TEUCHOS_TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
01100           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
01101     }
01102     
01105     static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
01106 #ifdef HAVE_ANASAZI_EXPERIMENTAL
01107                       , ConjType conj = Anasazi::CONJ
01108 #endif
01109                       )
01110     {
01111 #ifdef TEUCHOS_DEBUG
01112       TEUCHOS_TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
01113           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
01114       TEUCHOS_TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
01115           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
01116 #endif
01117       TEUCHOS_TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
01118           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");     
01119     }
01120 
01122 
01123 
01124 
01128     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
01129     { 
01130 #ifdef TEUCHOS_DEBUG
01131       TEUCHOS_TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
01132           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
01133 #endif
01134       TEUCHOS_TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
01135           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value."); 
01136     }
01137 
01139     
01141 
01142 
01144     static void 
01145     SetBlock (const Epetra_MultiVector& A, 
01146         const std::vector<int>& index, 
01147         Epetra_MultiVector& mv)
01148     { 
01149       const int inNumVecs = GetNumberVecs (A);
01150       const int outNumVecs = index.size();
01151 
01152       // FIXME (mfh 13 Jan 2011) Belos allows A to have more columns
01153       // than index.size(), in which case we just take the first
01154       // index.size() columns of A.  Anasazi requires that A have the
01155       // same number of columns as index.size().  Changing Anasazi's
01156       // behavior should not break existing Anasazi solvers, but the
01157       // tests need to be done.
01158       if (inNumVecs != outNumVecs)
01159   {
01160     std::ostringstream os;
01161     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
01162       "SetBlock(A, mv, index = {";
01163     if (outNumVecs > 0)
01164       {
01165         for (int k = 0; k < outNumVecs - 1; ++k)
01166     os << index[k] << ", ";
01167         os << index[outNumVecs-1];
01168       }
01169     os << "}): A has only " << inNumVecs << " columns, but there are "
01170        << outNumVecs << " indices in the index vector.";
01171     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
01172   }
01173       // Make a view of the columns of mv indicated by the index std::vector.
01174       Teuchos::RCP<Epetra_MultiVector> mv_view = CloneViewNonConst (mv, index);
01175 
01176       // View of columns [0, outNumVecs-1] of the source multivector A.
01177       // If A has fewer columns than mv_view, then create a view of
01178       // the first outNumVecs columns of A.
01179       Teuchos::RCP<const Epetra_MultiVector> A_view;
01180       if (outNumVecs == inNumVecs)
01181   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
01182       else
01183   A_view = CloneView (A, Teuchos::Range1D(0, outNumVecs - 1));
01184 
01185       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01186       // copies the data directly, ignoring the underlying
01187       // Epetra_Map(s).  If A and mv don't have the same data
01188       // distribution (Epetra_Map), this may result in incorrect or
01189       // undefined behavior.  Epetra_MultiVector::Update() also
01190       // ignores the Epetra_Maps, so we might as well just use the
01191       // (perhaps slightly cheaper) Assign() method via operator=().
01192       *mv_view = *A_view;
01193     }
01194 
01195     static void 
01196     SetBlock (const Epetra_MultiVector& A, 
01197         const Teuchos::Range1D& index, 
01198         Epetra_MultiVector& mv)
01199     { 
01200       const int numColsA = A.NumVectors();
01201       const int numColsMv = mv.NumVectors();
01202       // 'index' indexes into mv; it's the index set of the target.
01203       const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
01204       // We can't take more columns out of A than A has.
01205       const bool validSource = index.size() <= numColsA;
01206 
01207       if (! validIndex || ! validSource)
01208   {
01209     std::ostringstream os;
01210     os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::SetBlock"
01211       "(A, index=[" << index.lbound() << ", " << index.ubound() << "], "
01212       "mv): ";
01213     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
01214            os.str() << "Range lower bound must be nonnegative.");
01215     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
01216            os.str() << "Range upper bound must be less than "
01217            "the number of columns " << numColsA << " in the "
01218            "'mv' output argument.");
01219     TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
01220            os.str() << "Range must have no more elements than"
01221            " the number of columns " << numColsA << " in the "
01222            "'A' input argument.");
01223     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
01224   }
01225 
01226       // View of columns [index.lbound(), index.ubound()] of the
01227       // target multivector mv.  We avoid view creation overhead by
01228       // only creating a view if the index range is different than [0,
01229       // (# columns in mv) - 1].
01230       Teuchos::RCP<Epetra_MultiVector> mv_view;
01231       if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
01232   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
01233       else
01234   mv_view = CloneViewNonConst (mv, index);
01235 
01236       // View of columns [0, index.size()-1] of the source multivector
01237       // A.  If A has fewer columns than mv_view, then create a view
01238       // of the first index.size() columns of A.
01239       Teuchos::RCP<const Epetra_MultiVector> A_view;
01240       if (index.size() == numColsA)
01241   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
01242       else
01243   A_view = CloneView (A, Teuchos::Range1D(0, index.size()-1));
01244 
01245       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01246       // copies the data directly, ignoring the underlying
01247       // Epetra_Map(s).  If A and mv don't have the same data
01248       // distribution (Epetra_Map), this may result in incorrect or
01249       // undefined behavior.  Epetra_MultiVector::Update() also
01250       // ignores the Epetra_Maps, so we might as well just use the
01251       // (perhaps slightly cheaper) Assign() method via operator=().
01252       *mv_view = *A_view; 
01253     }
01254 
01255     static void 
01256     Assign (const Epetra_MultiVector& A, 
01257       Epetra_MultiVector& mv)
01258     {
01259       const int numColsA = GetNumberVecs (A);
01260       const int numColsMv = GetNumberVecs (mv);
01261       if (numColsA > numColsMv)
01262   {
01263     std::ostringstream os;
01264     os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::Assign"
01265       "(A, mv): ";
01266     TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
01267            os.str() << "Input multivector 'A' has " 
01268            << numColsA << " columns, but output multivector "
01269            "'mv' has only " << numColsMv << " columns.");
01270     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
01271   }
01272       // View of the first [0, numColsA-1] columns of mv.
01273       Teuchos::RCP<Epetra_MultiVector> mv_view;
01274       if (numColsMv == numColsA)
01275   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
01276       else // numColsMv > numColsA
01277   mv_view = CloneView (mv, Teuchos::Range1D(0, numColsA - 1));
01278       
01279       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01280       // copies the data directly, ignoring the underlying
01281       // Epetra_Map(s).  If A and mv don't have the same data
01282       // distribution (Epetra_Map), this may result in incorrect or
01283       // undefined behavior.  Epetra_MultiVector::Update() also
01284       // ignores the Epetra_Maps, so we might as well just use the
01285       // (perhaps slightly cheaper) Assign() method via operator=().
01286       *mv_view = A;
01287     }
01288 
01291     static void MvScale ( Epetra_MultiVector& mv, double alpha ) 
01292     { 
01293       TEUCHOS_TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
01294           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value."); 
01295     }
01296 
01299     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
01300     { 
01301       // Check to make sure the vector is as long as the multivector has columns.
01302       int numvecs = mv.NumVectors();
01303 #ifdef TEUCHOS_DEBUG
01304       TEUCHOS_TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
01305                           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
01306 #endif
01307       for (int i=0; i<numvecs; i++) {
01308         TEUCHOS_TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i])!=0, EpetraMultiVecFailure,
01309             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
01310       }
01311     }
01312 
01315     static void MvRandom( Epetra_MultiVector& mv )
01316     { 
01317       TEUCHOS_TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
01318           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
01319     }
01320 
01323     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
01324     { 
01325       TEUCHOS_TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
01326           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
01327     }
01328 
01330 
01332 
01333 
01336     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
01337     { os << mv << std::endl; }
01338 
01340 
01341 #if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
01342 #  if defined(HAVE_TPETRA_EPETRA)
01343 
01344 
01345 
01346 
01347 
01348     typedef Epetra::TsqrAdaptor tsqr_adaptor_type;
01349 #  endif // defined(HAVE_TPETRA_EPETRA)
01350 #endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
01351   };        
01352 
01354   //
01355   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
01356   //
01358 
01370   template <> 
01371   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
01372   {
01373   public:
01374     
01378     static void Apply ( const Epetra_Operator& Op, 
01379                         const Epetra_MultiVector& x, 
01380                         Epetra_MultiVector& y )
01381     { 
01382 #ifdef TEUCHOS_DEBUG
01383       TEUCHOS_TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
01384           "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
01385 #endif
01386       int ret = Op.Apply(x,y);
01387       TEUCHOS_TEST_FOR_EXCEPTION(ret != 0, OperatorError, 
01388           "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
01389     }
01390     
01391   };
01392   
01393 } // end of Anasazi namespace 
01394 
01395 #endif 
01396 // end of file ANASAZI_EPETRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends