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 ~EpetraMultiVecAccessor() {};
00089 
00091     virtual Epetra_MultiVector* GetEpetraMultiVec() { return 0; }
00092 
00094     virtual const Epetra_MultiVector* GetEpetraMultiVec() const { return 0; }
00095   };
00096 
00098   
00100   //
00101   //--------template class AnasaziEpetraMultiVec-----------------
00102   //
00104   
00111   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraMultiVec : public MultiVec<double>,  public Epetra_MultiVector, public EpetraMultiVecAccessor {
00112   public:
00114 
00115 
00117 
00122     EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs);
00123 
00125     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00126     
00128 
00136     EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride=0);
00137 
00139 
00145     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00146 
00148     virtual ~EpetraMultiVec() {};
00149 
00151 
00153 
00154 
00159     MultiVec<double> * Clone ( const int numvecs ) const;
00160 
00166     MultiVec<double> * CloneCopy () const;
00167 
00175     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00176     
00184     MultiVec<double> * CloneViewNonConst ( const std::vector<int>& index );
00185 
00193     const MultiVec<double> * CloneView ( const std::vector<int>& index ) const;
00194 
00196 
00198     int GetVecLength () const { return GlobalLength(); }
00199 
00202     ptrdiff_t GetGlobalLength () const 
00203     {
00204        if ( Map().GlobalIndicesLongLong() )
00205           return static_cast<ptrdiff_t>( GlobalLength64() );
00206        else 
00207           return static_cast<ptrdiff_t>( GlobalLength() );
00208     }
00209 
00212     int GetNumberVecs () const { return NumVectors(); }
00213 
00215 
00217 
00218 
00220     void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00221                            const Teuchos::SerialDenseMatrix<int,double>& B, 
00222                            double beta );
00223 
00226     void MvAddMv ( double alpha, const MultiVec<double>& A, 
00227                    double beta, const MultiVec<double>& B);
00228 
00231     void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00232 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00233         , ConjType conj = Anasazi::CONJ
00234 #endif
00235         ) const;
00236   
00239     void MvDot ( const MultiVec<double>& A, std::vector<double> &b
00240 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00241         , ConjType conj = Anasazi::CONJ
00242 #endif
00243         ) const;
00244 
00247     void MvScale ( double alpha ) { 
00248       TEUCHOS_TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
00249           "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00250     }
00251     
00254     void MvScale ( const std::vector<double>& alpha );
00255 
00257 
00258 
00259     
00263     void MvNorm ( std::vector<double> & normvec ) const {
00264       if (((int)normvec.size() >= GetNumberVecs()) ) {
00265         TEUCHOS_TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00266             "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
00267       }
00268     };
00270     
00272 
00273 
00278     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00279 
00282     void MvRandom() { 
00283       TEUCHOS_TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
00284           "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00285     };
00286 
00289     void MvInit ( double alpha ) { 
00290       TEUCHOS_TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
00291           "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00292     };
00293    
00295 
00296    
00298     Epetra_MultiVector* GetEpetraMultiVec() { return this; };
00299  
00301     const Epetra_MultiVector* GetEpetraMultiVec() const { return this; };
00302 
00304  
00306 
00307 
00308 
00310     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00312 
00313   private:
00314   };
00315   //-------------------------------------------------------------
00316   
00318   //
00319   //--------template class AnasaziEpetraOp---------------------
00320   //
00322   
00329   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraOp : public virtual Operator<double> {
00330   public:
00332 
00333     
00335     EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
00336     
00338     ~EpetraOp();
00340     
00342 
00343     
00347     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00349     
00350   private:
00351 //use pragmas to disable some false-positive warnings for windows 
00352 // sharedlibs export
00353 #ifdef _MSC_VER
00354 #pragma warning(push)
00355 #pragma warning(disable:4251)
00356 #endif
00357     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00358 #ifdef _MSC_VER
00359 #pragma warning(pop)
00360 #endif
00361   };
00362   //-------------------------------------------------------------
00363 
00365   //
00366   //--------template class AnasaziEpetraGenOp--------------------
00367   //
00369   
00383   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00384   public:
00386 
00389     EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp, 
00390                 const Teuchos::RCP<Epetra_Operator> &MOp,
00391                 bool isAInverse = true );
00392 
00394     ~EpetraGenOp();
00395     
00397 
00399     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00400 
00402 
00404     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00405 
00407 
00409     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00410 
00412     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00413     
00415     bool UseTranspose() const { return (false); };
00416 
00418     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00419     
00421     bool HasNormInf() const { return (false); };
00422     
00424     double NormInf() const  { return (-1.0); };
00425     
00427     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00428 
00430     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00431 
00433     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00434 
00435   private:
00436     bool isAInverse;
00437 
00438 //use pragmas to disable some false-positive warnings for windows 
00439 // sharedlibs export
00440 #ifdef _MSC_VER
00441 #pragma warning(push)
00442 #pragma warning(disable:4251)
00443 #endif
00444     Teuchos::RCP<Epetra_Operator> Epetra_AOp;
00445     Teuchos::RCP<Epetra_Operator> Epetra_MOp;
00446 #ifdef _MSC_VER
00447 #pragma warning(pop)
00448 #endif
00449   };
00450   
00452   //
00453   //--------template class AnasaziEpetraSymOp--------------------
00454   //
00456 
00469   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00470   public:
00472 
00474     EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
00475 
00477     ~EpetraSymOp();
00478     
00480 
00482     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00483 
00485 
00487     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00488 
00490 
00493     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00494 
00496     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00497     
00499     bool UseTranspose() const { return (false); };
00500 
00502     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00503     
00505     bool HasNormInf() const { return (false); };
00506     
00508     double NormInf() const  { return (-1.0); };
00509     
00511     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00512 
00514     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00515 
00517     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00518 
00519   private:
00520 
00521 //use pragmas to disable false-positive warnings in generating windows sharedlib exports
00522 #ifdef _MSC_VER
00523 #pragma warning(push)
00524 #pragma warning(disable:4251)
00525 #endif
00526     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00527 #ifdef _MSC_VER
00528 #pragma warning(pop)
00529 #endif
00530 
00531     bool isTrans_;
00532   };
00533 
00534 
00536   //
00537   //--------template class AnasaziEpetraSymMVOp---------------------
00538   //
00540 
00553   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymMVOp : public virtual Operator<double> {
00554   public:
00556 
00558     EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00559                   bool isTrans = false );
00560     
00562     ~EpetraSymMVOp() {};
00563     
00565 
00567     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00568 
00569   private:
00570 
00571 //use pragmas to disable some false-positive warnings for windows 
00572 // sharedlibs export
00573 #ifdef _MSC_VER
00574 #pragma warning(push)
00575 #pragma warning(disable:4251)
00576 #endif
00577     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00578     Teuchos::RCP<const Epetra_Map> MV_localmap;
00579     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00580 #ifdef _MSC_VER
00581 #pragma warning(pop)
00582 #endif
00583 
00584     bool isTrans_;
00585   };
00586 
00588   //
00589   //--------template class AnasaziEpetraWSymMVOp---------------------
00590   //
00592 
00605   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraWSymMVOp : public virtual Operator<double> {
00606   public:
00608     EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00609                    const Teuchos::RCP<Epetra_Operator> &OP );
00610     
00612     ~EpetraWSymMVOp() {};
00613     
00615 
00617     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00618 
00619   private:
00620 //use pragmas to disable some false-positive warnings for windows 
00621 // sharedlibs export
00622 #ifdef _MSC_VER
00623 #pragma warning(push)
00624 #pragma warning(disable:4251)
00625 #endif
00626     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00627     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00628     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00629     Teuchos::RCP<const Epetra_Map> MV_localmap;
00630     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00631 #ifdef _MSC_VER
00632 #pragma warning(pop)
00633 #endif
00634   };
00635 
00637   //
00638   //--------template class AnasaziEpetraW2SymMVOp---------------------
00639   //
00641 
00654   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraW2SymMVOp : public virtual Operator<double> {
00655   public:
00657     EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00658                    const Teuchos::RCP<Epetra_Operator> &OP );
00659     
00661     ~EpetraW2SymMVOp() {};
00662     
00664 
00666     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00667 
00668   private:
00669 //use pragmas to disable some false-positive warnings for windows 
00670 // sharedlibs export
00671 #ifdef _MSC_VER
00672 #pragma warning(push)
00673 #pragma warning(disable:4251)
00674 #endif
00675     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00676     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00677     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00678     Teuchos::RCP<const Epetra_Map> MV_localmap;
00679     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00680 #ifdef _MSC_VER
00681 #pragma warning(pop)
00682 #endif
00683   };
00684 
00685   
00687   //
00688   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00689   //
00691 
00702   template<>
00703   class MultiVecTraits<double, Epetra_MultiVector>
00704   {
00705   public:
00706 
00708 
00709 
00714     static Teuchos::RCP<Epetra_MultiVector> 
00715     Clone (const Epetra_MultiVector& mv, const int outNumVecs)
00716     { 
00717       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs <= 0, std::invalid_argument,
00718        "Belos::MultiVecTraits<double, Epetra_MultiVector>::"
00719        "Clone(mv, outNumVecs = " << outNumVecs << "): "
00720        "outNumVecs must be positive.");
00721       // FIXME (mfh 13 Jan 2011) Anasazi currently lets Epetra fill in
00722       // the entries of the returned multivector with zeros, but Belos
00723       // does not.  We retain this different behavior for now, but the
00724       // two versions will need to be reconciled.
00725       return Teuchos::rcp (new Epetra_MultiVector (mv.Map(), outNumVecs)); 
00726     }
00727 
00732     static Teuchos::RCP<Epetra_MultiVector> 
00733     CloneCopy (const Epetra_MultiVector& mv)
00734     { 
00735       return Teuchos::rcp (new Epetra_MultiVector (mv)); 
00736     }
00737 
00743     static Teuchos::RCP<Epetra_MultiVector> 
00744     CloneCopy (const Epetra_MultiVector& mv, const std::vector<int>& index)
00745     { 
00746       const int inNumVecs = GetNumberVecs (mv);
00747       const int outNumVecs = index.size();
00748 
00749       // Simple, inexpensive tests of the index vector.
00750       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00751        "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00752        "CloneCopy(mv, index = {}): At least one vector must be"
00753        " cloned from mv.");
00754       if (outNumVecs > inNumVecs)
00755   {
00756     std::ostringstream os;
00757     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00758       "CloneCopy(mv, index = {";
00759     for (int k = 0; k < outNumVecs - 1; ++k)
00760       os << index[k] << ", ";
00761     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00762        << " indices to copy, but only " << inNumVecs << " columns of mv.";
00763     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00764   }
00765 #ifdef TEUCHOS_DEBUG
00766       // In debug mode, we perform more expensive tests of the index
00767       // vector, to ensure all the elements are in range.
00768       // Dereferencing the iterator is valid because index has length
00769       // > 0.
00770       const int minIndex = *std::min_element (index.begin(), index.end());
00771       const int maxIndex = *std::max_element (index.begin(), index.end());
00772 
00773       if (minIndex < 0)
00774   {
00775     std::ostringstream os;
00776     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00777       "CloneCopy(mv, index = {";
00778     for (int k = 0; k < outNumVecs - 1; ++k)
00779       os << index[k] << ", ";
00780     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00781       "the smallest index " << minIndex << " is negative.";
00782     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00783   }
00784       if (maxIndex >= inNumVecs)
00785   {
00786     std::ostringstream os;
00787     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00788       "CloneCopy(mv, index = {";
00789     for (int k = 0; k < outNumVecs - 1; ++k)
00790       os << index[k] << ", ";
00791     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00792       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00793        << maxIndex << " is out of bounds.";
00794     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00795   }
00796 #endif // TEUCHOS_DEBUG
00797       // Cast to nonconst, because Epetra_MultiVector's constructor
00798       // wants a nonconst int array argument.  It doesn't actually
00799       // change the entries of the array.
00800       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00801       return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, &tmpind[0], index.size())); 
00802       // return Teuchos::rcp (new Epetra_MultiVector (::Copy, mv, &tmpind[0], index.size())); 
00803     }
00804 
00805     static Teuchos::RCP<Epetra_MultiVector> 
00806     CloneCopy (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00807     { 
00808       const int inNumVecs = GetNumberVecs (mv);
00809       const int outNumVecs = index.size();
00810       const bool validRange = outNumVecs > 0 && index.lbound() >= 0 && 
00811   index.ubound() < inNumVecs;
00812       if (! validRange)
00813   {
00814     std::ostringstream os;
00815     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,"
00816       "index=[" << index.lbound() << ", " << index.ubound() << "]): ";
00817     TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00818            os.str() << "Column index range must be nonempty.");
00819     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00820            os.str() << "Column index range must be nonnegative.");
00821     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= inNumVecs, std::invalid_argument,
00822            os.str() << "Column index range must not exceed "
00823            "number of vectors " << inNumVecs << " in the "
00824            "input multivector.");
00825   }
00826       return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, index.lbound(), index.size()));
00827     }
00828 
00834     static Teuchos::RCP<Epetra_MultiVector> 
00835     CloneViewNonConst (Epetra_MultiVector& mv, const std::vector<int>& index)
00836     { 
00837       const int inNumVecs = GetNumberVecs (mv);
00838       const int outNumVecs = index.size();
00839 
00840       // Simple, inexpensive tests of the index vector.
00841       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00842        "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00843        "CloneViewNonConst(mv, index = {}): The output view "
00844        "must have at least one column.");
00845       if (outNumVecs > inNumVecs)
00846   {
00847     std::ostringstream os;
00848     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00849       "CloneViewNonConst(mv, index = {";
00850     for (int k = 0; k < outNumVecs - 1; ++k)
00851       os << index[k] << ", ";
00852     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00853        << " indices to view, but only " << inNumVecs << " columns of mv.";
00854     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00855   }
00856 #ifdef TEUCHOS_DEBUG
00857       // In debug mode, we perform more expensive tests of the index
00858       // vector, to ensure all the elements are in range.
00859       // Dereferencing the iterator is valid because index has length
00860       // > 0.
00861       const int minIndex = *std::min_element (index.begin(), index.end());
00862       const int maxIndex = *std::max_element (index.begin(), index.end());
00863 
00864       if (minIndex < 0)
00865   {
00866     std::ostringstream os;
00867     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00868       "CloneViewNonConst(mv, index = {";
00869     for (int k = 0; k < outNumVecs - 1; ++k)
00870       os << index[k] << ", ";
00871     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00872       "the smallest index " << minIndex << " is negative.";
00873     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00874   }
00875       if (maxIndex >= inNumVecs)
00876   {
00877     std::ostringstream os;
00878     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00879       "CloneViewNonConst(mv, index = {";
00880     for (int k = 0; k < outNumVecs - 1; ++k)
00881       os << index[k] << ", ";
00882     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00883       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00884        << maxIndex << " is out of bounds.";
00885     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00886   }
00887 #endif // TEUCHOS_DEBUG
00888       // Cast to nonconst, because Epetra_MultiVector's constructor
00889       // wants a nonconst int array argument.  It doesn't actually
00890       // change the entries of the array.
00891       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00892       return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
00893     }
00894 
00895     static Teuchos::RCP<Epetra_MultiVector> 
00896     CloneViewNonConst (Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00897     { 
00898       const bool validRange = index.size() > 0 && 
00899   index.lbound() >= 0 && 
00900   index.ubound() < mv.NumVectors();
00901       if (! validRange)
00902   {
00903     std::ostringstream os;
00904     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
00905       "NonConst(mv,index=[" << index.lbound() << ", " << index.ubound() 
00906        << "]): ";
00907     TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00908            os.str() << "Column index range must be nonempty.");
00909     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00910            os.str() << "Column index range must be nonnegative.");
00911     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(), 
00912            std::invalid_argument,
00913            os.str() << "Column index range must not exceed "
00914            "number of vectors " << mv.NumVectors() << " in "
00915            "the input multivector.");
00916   }
00917       return Teuchos::rcp (new Epetra_MultiVector (View, mv, index.lbound(), index.size()));
00918     }
00919 
00925     static Teuchos::RCP<const Epetra_MultiVector> 
00926     CloneView (const Epetra_MultiVector& mv, const std::vector<int>& index)
00927     { 
00928       const int inNumVecs = GetNumberVecs (mv);
00929       const int outNumVecs = index.size();
00930 
00931       // Simple, inexpensive tests of the index vector.
00932       TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00933        "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00934        "CloneView(mv, index = {}): The output view "
00935        "must have at least one column.");
00936       if (outNumVecs > inNumVecs)
00937   {
00938     std::ostringstream os;
00939     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00940       "CloneView(mv, index = {";
00941     for (int k = 0; k < outNumVecs - 1; ++k)
00942       os << index[k] << ", ";
00943     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00944        << " indices to view, but only " << inNumVecs << " columns of mv.";
00945     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00946   }
00947 #ifdef TEUCHOS_DEBUG
00948       // In debug mode, we perform more expensive tests of the index
00949       // vector, to ensure all the elements are in range.
00950       // Dereferencing the iterator is valid because index has length
00951       // > 0.
00952       const int minIndex = *std::min_element (index.begin(), index.end());
00953       const int maxIndex = *std::max_element (index.begin(), index.end());
00954 
00955       if (minIndex < 0)
00956   {
00957     std::ostringstream os;
00958     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00959       "CloneView(mv, index = {";
00960     for (int k = 0; k < outNumVecs - 1; ++k)
00961       os << index[k] << ", ";
00962     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00963       "the smallest index " << minIndex << " is negative.";
00964     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00965   }
00966       if (maxIndex >= inNumVecs)
00967   {
00968     std::ostringstream os;
00969     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00970       "CloneView(mv, index = {";
00971     for (int k = 0; k < outNumVecs - 1; ++k)
00972       os << index[k] << ", ";
00973     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00974       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00975        << maxIndex << " is out of bounds.";
00976     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00977   }
00978 #endif // TEUCHOS_DEBUG
00979       // Cast to nonconst, because Epetra_MultiVector's constructor
00980       // wants a nonconst int array argument.  It doesn't actually
00981       // change the entries of the array.
00982       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00983       return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
00984     }
00985 
00986     static Teuchos::RCP<Epetra_MultiVector> 
00987     CloneView (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00988     { 
00989       const bool validRange = index.size() > 0 && 
00990   index.lbound() >= 0 && 
00991   index.ubound() < mv.NumVectors();
00992       if (! validRange)
00993   {
00994     std::ostringstream os;
00995     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
00996       "(mv,index=[" << index.lbound() << ", " << index.ubound() 
00997        << "]): ";
00998     TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00999            os.str() << "Column index range must be nonempty.");
01000     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
01001            os.str() << "Column index range must be nonnegative.");
01002     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(), 
01003            std::invalid_argument,
01004            os.str() << "Column index range must not exceed "
01005            "number of vectors " << mv.NumVectors() << " in "
01006            "the input multivector.");
01007   }
01008       return Teuchos::rcp (new Epetra_MultiVector(View, mv, index.lbound(), index.size()));
01009     }
01010 
01012 
01014 
01015 
01017     static int GetVecLength( const Epetra_MultiVector& mv )
01018     { return mv.GlobalLength(); }
01019 
01021     static int GetNumberVecs( const Epetra_MultiVector& mv )
01022     { return mv.NumVectors(); }
01023 
01024     static bool HasConstantStride( const Epetra_MultiVector& mv )
01025     { return mv.ConstantStride(); }
01027 
01029 
01030 
01033     static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A, 
01034                                  const Teuchos::SerialDenseMatrix<int,double>& B, 
01035                                  double beta, Epetra_MultiVector& mv )
01036     { 
01037       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
01038       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
01039 
01040       TEUCHOS_TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
01041           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Multiply() returned a nonzero value.");
01042     }
01043 
01046     static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
01047     { 
01048       // epetra mv.Update(alpha,A,beta,B,gamma) will check 
01049       //   alpha == 0.0 
01050       // and 
01051       //   beta == 0.0 
01052       // and based on this will call 
01053       //   mv.Update(beta,B,gamma) 
01054       // or 
01055       //   mv.Update(alpha,A,gamma)
01056       //
01057       // mv.Update(alpha,A,gamma) 
01058       // will then check for one of 
01059       //   gamma == 0
01060       // or 
01061       //   gamma == 1
01062       // or 
01063       //   alpha == 1 
01064       // in that order. however, it will not look for the combination
01065       //   alpha == 1 and gamma = 0
01066       // which is a common use case when we wish to assign 
01067       //   mv = A   (in which case alpha == 1, beta == gamma == 0)
01068       // or 
01069       //   mv = B   (in which case beta == 1, alpha == gamma == 0)
01070       // 
01071       // therefore, we will check for these use cases ourselves
01072       if (beta == 0.0) {
01073         if (alpha == 1.0) {
01074           // assign
01075           mv = A;
01076         }
01077         else {
01078           // single update
01079           TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
01080               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
01081         }
01082       }
01083       else if (alpha == 0.0) {
01084         if (beta == 1.0) {
01085           // assign 
01086           mv = B;
01087         }
01088         else {
01089           // single update
01090           TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
01091               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
01092         }
01093       }
01094       else {
01095         // double update
01096         TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
01097             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
01098       }
01099     }
01100 
01103     static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
01104 #ifdef HAVE_ANASAZI_EXPERIMENTAL
01105                           , ConjType conj = Anasazi::CONJ
01106 #endif
01107                         )
01108     { 
01109       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
01110       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
01111       
01112       TEUCHOS_TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
01113           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
01114     }
01115     
01118     static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
01119 #ifdef HAVE_ANASAZI_EXPERIMENTAL
01120                       , ConjType conj = Anasazi::CONJ
01121 #endif
01122                       )
01123     {
01124 #ifdef TEUCHOS_DEBUG
01125       TEUCHOS_TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
01126           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
01127       TEUCHOS_TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
01128           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
01129 #endif
01130       TEUCHOS_TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
01131           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");     
01132     }
01133 
01135 
01136 
01137 
01141     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
01142     { 
01143 #ifdef TEUCHOS_DEBUG
01144       TEUCHOS_TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
01145           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
01146 #endif
01147       TEUCHOS_TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
01148           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value."); 
01149     }
01150 
01152     
01154 
01155 
01157     static void 
01158     SetBlock (const Epetra_MultiVector& A, 
01159         const std::vector<int>& index, 
01160         Epetra_MultiVector& mv)
01161     { 
01162       const int inNumVecs = GetNumberVecs (A);
01163       const int outNumVecs = index.size();
01164 
01165       // FIXME (mfh 13 Jan 2011) Belos allows A to have more columns
01166       // than index.size(), in which case we just take the first
01167       // index.size() columns of A.  Anasazi requires that A have the
01168       // same number of columns as index.size().  Changing Anasazi's
01169       // behavior should not break existing Anasazi solvers, but the
01170       // tests need to be done.
01171       if (inNumVecs != outNumVecs)
01172   {
01173     std::ostringstream os;
01174     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
01175       "SetBlock(A, mv, index = {";
01176     if (outNumVecs > 0)
01177       {
01178         for (int k = 0; k < outNumVecs - 1; ++k)
01179     os << index[k] << ", ";
01180         os << index[outNumVecs-1];
01181       }
01182     os << "}): A has only " << inNumVecs << " columns, but there are "
01183        << outNumVecs << " indices in the index vector.";
01184     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
01185   }
01186       // Make a view of the columns of mv indicated by the index std::vector.
01187       Teuchos::RCP<Epetra_MultiVector> mv_view = CloneViewNonConst (mv, index);
01188 
01189       // View of columns [0, outNumVecs-1] of the source multivector A.
01190       // If A has fewer columns than mv_view, then create a view of
01191       // the first outNumVecs columns of A.
01192       Teuchos::RCP<const Epetra_MultiVector> A_view;
01193       if (outNumVecs == inNumVecs)
01194   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
01195       else
01196   A_view = CloneView (A, Teuchos::Range1D(0, outNumVecs - 1));
01197 
01198       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01199       // copies the data directly, ignoring the underlying
01200       // Epetra_Map(s).  If A and mv don't have the same data
01201       // distribution (Epetra_Map), this may result in incorrect or
01202       // undefined behavior.  Epetra_MultiVector::Update() also
01203       // ignores the Epetra_Maps, so we might as well just use the
01204       // (perhaps slightly cheaper) Assign() method via operator=().
01205       *mv_view = *A_view;
01206     }
01207 
01208     static void 
01209     SetBlock (const Epetra_MultiVector& A, 
01210         const Teuchos::Range1D& index, 
01211         Epetra_MultiVector& mv)
01212     { 
01213       const int numColsA = A.NumVectors();
01214       const int numColsMv = mv.NumVectors();
01215       // 'index' indexes into mv; it's the index set of the target.
01216       const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
01217       // We can't take more columns out of A than A has.
01218       const bool validSource = index.size() <= numColsA;
01219 
01220       if (! validIndex || ! validSource)
01221   {
01222     std::ostringstream os;
01223     os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::SetBlock"
01224       "(A, index=[" << index.lbound() << ", " << index.ubound() << "], "
01225       "mv): ";
01226     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
01227            os.str() << "Range lower bound must be nonnegative.");
01228     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
01229            os.str() << "Range upper bound must be less than "
01230            "the number of columns " << numColsA << " in the "
01231            "'mv' output argument.");
01232     TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
01233            os.str() << "Range must have no more elements than"
01234            " the number of columns " << numColsA << " in the "
01235            "'A' input argument.");
01236     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
01237   }
01238 
01239       // View of columns [index.lbound(), index.ubound()] of the
01240       // target multivector mv.  We avoid view creation overhead by
01241       // only creating a view if the index range is different than [0,
01242       // (# columns in mv) - 1].
01243       Teuchos::RCP<Epetra_MultiVector> mv_view;
01244       if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
01245   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
01246       else
01247   mv_view = CloneViewNonConst (mv, index);
01248 
01249       // View of columns [0, index.size()-1] of the source multivector
01250       // A.  If A has fewer columns than mv_view, then create a view
01251       // of the first index.size() columns of A.
01252       Teuchos::RCP<const Epetra_MultiVector> A_view;
01253       if (index.size() == numColsA)
01254   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
01255       else
01256   A_view = CloneView (A, Teuchos::Range1D(0, index.size()-1));
01257 
01258       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01259       // copies the data directly, ignoring the underlying
01260       // Epetra_Map(s).  If A and mv don't have the same data
01261       // distribution (Epetra_Map), this may result in incorrect or
01262       // undefined behavior.  Epetra_MultiVector::Update() also
01263       // ignores the Epetra_Maps, so we might as well just use the
01264       // (perhaps slightly cheaper) Assign() method via operator=().
01265       *mv_view = *A_view; 
01266     }
01267 
01268     static void 
01269     Assign (const Epetra_MultiVector& A, 
01270       Epetra_MultiVector& mv)
01271     {
01272       const int numColsA = GetNumberVecs (A);
01273       const int numColsMv = GetNumberVecs (mv);
01274       if (numColsA > numColsMv)
01275   {
01276     std::ostringstream os;
01277     os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::Assign"
01278       "(A, mv): ";
01279     TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
01280            os.str() << "Input multivector 'A' has " 
01281            << numColsA << " columns, but output multivector "
01282            "'mv' has only " << numColsMv << " columns.");
01283     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
01284   }
01285       // View of the first [0, numColsA-1] columns of mv.
01286       Teuchos::RCP<Epetra_MultiVector> mv_view;
01287       if (numColsMv == numColsA)
01288   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
01289       else // numColsMv > numColsA
01290   mv_view = CloneView (mv, Teuchos::Range1D(0, numColsA - 1));
01291       
01292       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01293       // copies the data directly, ignoring the underlying
01294       // Epetra_Map(s).  If A and mv don't have the same data
01295       // distribution (Epetra_Map), this may result in incorrect or
01296       // undefined behavior.  Epetra_MultiVector::Update() also
01297       // ignores the Epetra_Maps, so we might as well just use the
01298       // (perhaps slightly cheaper) Assign() method via operator=().
01299       *mv_view = A;
01300     }
01301 
01304     static void MvScale ( Epetra_MultiVector& mv, double alpha ) 
01305     { 
01306       TEUCHOS_TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
01307           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value."); 
01308     }
01309 
01312     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
01313     { 
01314       // Check to make sure the vector is as long as the multivector has columns.
01315       int numvecs = mv.NumVectors();
01316 #ifdef TEUCHOS_DEBUG
01317       TEUCHOS_TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
01318                           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
01319 #endif
01320       for (int i=0; i<numvecs; i++) {
01321         TEUCHOS_TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i])!=0, EpetraMultiVecFailure,
01322             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
01323       }
01324     }
01325 
01328     static void MvRandom( Epetra_MultiVector& mv )
01329     { 
01330       TEUCHOS_TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
01331           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
01332     }
01333 
01336     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
01337     { 
01338       TEUCHOS_TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
01339           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
01340     }
01341 
01343 
01345 
01346 
01349     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
01350     { os << mv << std::endl; }
01351 
01353 
01354 #if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
01355 #  if defined(HAVE_TPETRA_EPETRA)
01356 
01357 
01358 
01359 
01360 
01361     typedef Epetra::TsqrAdaptor tsqr_adaptor_type;
01362 #  endif // defined(HAVE_TPETRA_EPETRA)
01363 #endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
01364   };        
01365 
01367   //
01368   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
01369   //
01371 
01383   template <> 
01384   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
01385   {
01386   public:
01387     
01391     static void Apply ( const Epetra_Operator& Op, 
01392                         const Epetra_MultiVector& x, 
01393                         Epetra_MultiVector& y )
01394     { 
01395 #ifdef TEUCHOS_DEBUG
01396       TEUCHOS_TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
01397           "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
01398 #endif
01399       int ret = Op.Apply(x,y);
01400       TEUCHOS_TEST_FOR_EXCEPTION(ret != 0, OperatorError, 
01401           "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
01402     }
01403     
01404   };
01405 
01406   template<>
01407   class MultiVecTraitsExt<double, Epetra_MultiVector>
01408   {
01409   public:
01411 
01412 
01415     static ptrdiff_t GetGlobalLength( const Epetra_MultiVector& mv )
01416     { 
01417       if (mv.Map().GlobalIndicesLongLong())
01418         return static_cast<ptrdiff_t>( mv.GlobalLength64() );
01419       else
01420         return static_cast<ptrdiff_t>( mv.GlobalLength() );
01421     }
01422 
01424   };
01425   
01426 } // end of Anasazi namespace 
01427 
01428 #endif 
01429 // end of file ANASAZI_EPETRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends