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 #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   //--------template class AnasaziEpetraMultiVec-----------------
00081   //
00083   
00090   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector {
00091   public:
00093 
00094 
00096 
00101     EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs);
00102 
00104     EpetraMultiVec(const Epetra_MultiVector & P_vec);
00105     
00107 
00115     EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride=0);
00116 
00118 
00124     EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
00125 
00127     virtual ~EpetraMultiVec() {};
00128 
00130 
00132 
00133 
00138     MultiVec<double> * Clone ( const int numvecs ) const;
00139 
00145     MultiVec<double> * CloneCopy () const;
00146 
00154     MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
00155     
00163     MultiVec<double> * CloneViewNonConst ( const std::vector<int>& index );
00164 
00172     const MultiVec<double> * CloneView ( const std::vector<int>& index ) const;
00173 
00175 
00177 
00178 
00180     int GetNumberVecs () const { return NumVectors(); }
00181 
00183     int GetVecLength () const { return GlobalLength(); }
00184 
00186 
00188 
00189 
00191     void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
00192                            const Teuchos::SerialDenseMatrix<int,double>& B, 
00193                            double beta );
00194 
00197     void MvAddMv ( double alpha, const MultiVec<double>& A, 
00198                    double beta, const MultiVec<double>& B);
00199 
00202     void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B 
00203 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00204         , ConjType conj = Anasazi::CONJ
00205 #endif
00206         ) const;
00207   
00210     void MvDot ( const MultiVec<double>& A, std::vector<double> &b
00211 #ifdef HAVE_ANASAZI_EXPERIMENTAL
00212         , ConjType conj = Anasazi::CONJ
00213 #endif
00214         ) const;
00215 
00218     void MvScale ( double alpha ) { 
00219       TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
00220           "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
00221     }
00222     
00225     void MvScale ( const std::vector<double>& alpha );
00226 
00228 
00229 
00230     
00234     void MvNorm ( std::vector<double> & normvec ) const {
00235       if (((int)normvec.size() >= GetNumberVecs()) ) {
00236         TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
00237             "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
00238       }
00239     };
00241     
00243 
00244 
00249     void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
00250 
00253     void MvRandom() { 
00254       TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
00255           "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
00256     };
00257 
00260     void MvInit ( double alpha ) { 
00261       TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
00262           "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
00263     };
00264     
00266 
00267 
00268 
00270     void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
00272 
00273   private:
00274   };
00275   //-------------------------------------------------------------
00276   
00278   //
00279   //--------template class AnasaziEpetraOp---------------------
00280   //
00282   
00289   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraOp : public virtual Operator<double> {
00290   public:
00292 
00293     
00295     EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
00296     
00298     ~EpetraOp();
00300     
00302 
00303     
00307     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
00309     
00310   private:
00311 //use pragmas to disable some false-positive warnings for windows 
00312 // sharedlibs export
00313 #ifdef _MSC_VER
00314 #pragma warning(push)
00315 #pragma warning(disable:4251)
00316 #endif
00317     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00318 #ifdef _MSC_VER
00319 #pragma warning(pop)
00320 #endif
00321   };
00322   //-------------------------------------------------------------
00323 
00325   //
00326   //--------template class AnasaziEpetraGenOp--------------------
00327   //
00329   
00343   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
00344   public:
00346 
00349     EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp, 
00350                 const Teuchos::RCP<Epetra_Operator> &MOp,
00351                 bool isAInverse = true );
00352 
00354     ~EpetraGenOp();
00355     
00357 
00359     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00360 
00362 
00364     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00365 
00367 
00369     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00370 
00372     const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
00373     
00375     bool UseTranspose() const { return (false); };
00376 
00378     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00379     
00381     bool HasNormInf() const { return (false); };
00382     
00384     double NormInf() const  { return (-1.0); };
00385     
00387     const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
00388 
00390     const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
00391 
00393     const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
00394 
00395   private:
00396     bool isAInverse;
00397 
00398 //use pragmas to disable some false-positive warnings for windows 
00399 // sharedlibs export
00400 #ifdef _MSC_VER
00401 #pragma warning(push)
00402 #pragma warning(disable:4251)
00403 #endif
00404     Teuchos::RCP<Epetra_Operator> Epetra_AOp;
00405     Teuchos::RCP<Epetra_Operator> Epetra_MOp;
00406 #ifdef _MSC_VER
00407 #pragma warning(pop)
00408 #endif
00409   };
00410   
00412   //
00413   //--------template class AnasaziEpetraSymOp--------------------
00414   //
00416 
00429   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
00430   public:
00432 
00434     EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
00435 
00437     ~EpetraSymOp();
00438     
00440 
00442     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00443 
00445 
00447     int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00448 
00450 
00453     int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
00454 
00456     const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
00457     
00459     bool UseTranspose() const { return (false); };
00460 
00462     int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
00463     
00465     bool HasNormInf() const { return (false); };
00466     
00468     double NormInf() const  { return (-1.0); };
00469     
00471     const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
00472 
00474     const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
00475 
00477     const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
00478 
00479   private:
00480 
00481 //use pragmas to disable false-positive warnings in generating windows sharedlib exports
00482 #ifdef _MSC_VER
00483 #pragma warning(push)
00484 #pragma warning(disable:4251)
00485 #endif
00486     Teuchos::RCP<Epetra_Operator> Epetra_Op;
00487 #ifdef _MSC_VER
00488 #pragma warning(pop)
00489 #endif
00490 
00491     bool isTrans_;
00492   };
00493 
00494 
00496   //
00497   //--------template class AnasaziEpetraSymMVOp---------------------
00498   //
00500 
00513   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymMVOp : public virtual Operator<double> {
00514   public:
00516 
00518     EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00519                   bool isTrans = false );
00520     
00522     ~EpetraSymMVOp() {};
00523     
00525 
00527     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00528 
00529   private:
00530 
00531 //use pragmas to disable some false-positive warnings for windows 
00532 // sharedlibs export
00533 #ifdef _MSC_VER
00534 #pragma warning(push)
00535 #pragma warning(disable:4251)
00536 #endif
00537     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00538     Teuchos::RCP<const Epetra_Map> MV_localmap;
00539     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00540 #ifdef _MSC_VER
00541 #pragma warning(pop)
00542 #endif
00543 
00544     bool isTrans_;
00545   };
00546 
00548   //
00549   //--------template class AnasaziEpetraWSymMVOp---------------------
00550   //
00552 
00565   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraWSymMVOp : public virtual Operator<double> {
00566   public:
00568     EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00569                    const Teuchos::RCP<Epetra_Operator> &OP );
00570     
00572     ~EpetraWSymMVOp() {};
00573     
00575 
00577     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00578 
00579   private:
00580 //use pragmas to disable some false-positive warnings for windows 
00581 // sharedlibs export
00582 #ifdef _MSC_VER
00583 #pragma warning(push)
00584 #pragma warning(disable:4251)
00585 #endif
00586     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00587     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00588     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00589     Teuchos::RCP<const Epetra_Map> MV_localmap;
00590     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00591 #ifdef _MSC_VER
00592 #pragma warning(pop)
00593 #endif
00594   };
00595 
00597   //
00598   //--------template class AnasaziEpetraW2SymMVOp---------------------
00599   //
00601 
00614   class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraW2SymMVOp : public virtual Operator<double> {
00615   public:
00617     EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV, 
00618                    const Teuchos::RCP<Epetra_Operator> &OP );
00619     
00621     ~EpetraW2SymMVOp() {};
00622     
00624 
00626     void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const; 
00627 
00628   private:
00629 //use pragmas to disable some false-positive warnings for windows 
00630 // sharedlibs export
00631 #ifdef _MSC_VER
00632 #pragma warning(push)
00633 #pragma warning(disable:4251)
00634 #endif
00635     Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
00636     Teuchos::RCP<Epetra_Operator> Epetra_OP;
00637     Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
00638     Teuchos::RCP<const Epetra_Map> MV_localmap;
00639     Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
00640 #ifdef _MSC_VER
00641 #pragma warning(pop)
00642 #endif
00643   };
00644 
00645   
00647   //
00648   // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
00649   //
00651 
00662   template<>
00663   class MultiVecTraits<double, Epetra_MultiVector>
00664   {
00665   public:
00666 
00668 
00669 
00674     static Teuchos::RCP<Epetra_MultiVector> 
00675     Clone (const Epetra_MultiVector& mv, const int outNumVecs)
00676     { 
00677       TEST_FOR_EXCEPTION(outNumVecs <= 0, std::invalid_argument,
00678        "Belos::MultiVecTraits<double, Epetra_MultiVector>::"
00679        "Clone(mv, outNumVecs = " << outNumVecs << "): "
00680        "outNumVecs must be positive.");
00681       // FIXME (mfh 13 Jan 2011) Anasazi currently lets Epetra fill in
00682       // the entries of the returned multivector with zeros, but Belos
00683       // does not.  We retain this different behavior for now, but the
00684       // two versions will need to be reconciled.
00685       return Teuchos::rcp (new Epetra_MultiVector (mv.Map(), outNumVecs)); 
00686     }
00687 
00692     static Teuchos::RCP<Epetra_MultiVector> 
00693     CloneCopy (const Epetra_MultiVector& mv)
00694     { 
00695       return Teuchos::rcp (new Epetra_MultiVector (mv)); 
00696     }
00697 
00703     static Teuchos::RCP<Epetra_MultiVector> 
00704     CloneCopy (const Epetra_MultiVector& mv, const std::vector<int>& index)
00705     { 
00706       const int inNumVecs = GetNumberVecs (mv);
00707       const int outNumVecs = index.size();
00708 
00709       // Simple, inexpensive tests of the index vector.
00710       TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00711        "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00712        "CloneCopy(mv, index = {}): At least one vector must be"
00713        " cloned from mv.");
00714       if (outNumVecs > inNumVecs)
00715   {
00716     std::ostringstream os;
00717     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00718       "CloneCopy(mv, index = {";
00719     for (int k = 0; k < outNumVecs - 1; ++k)
00720       os << index[k] << ", ";
00721     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00722        << " indices to copy, but only " << inNumVecs << " columns of mv.";
00723     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00724   }
00725 #ifdef TEUCHOS_DEBUG
00726       // In debug mode, we perform more expensive tests of the index
00727       // vector, to ensure all the elements are in range.
00728       // Dereferencing the iterator is valid because index has length
00729       // > 0.
00730       const int minIndex = *std::min_element (index.begin(), index.end());
00731       const int maxIndex = *std::max_element (index.begin(), index.end());
00732 
00733       if (minIndex < 0)
00734   {
00735     std::ostringstream os;
00736     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00737       "CloneCopy(mv, index = {";
00738     for (int k = 0; k < outNumVecs - 1; ++k)
00739       os << index[k] << ", ";
00740     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00741       "the smallest index " << minIndex << " is negative.";
00742     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00743   }
00744       if (maxIndex >= inNumVecs)
00745   {
00746     std::ostringstream os;
00747     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00748       "CloneCopy(mv, index = {";
00749     for (int k = 0; k < outNumVecs - 1; ++k)
00750       os << index[k] << ", ";
00751     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00752       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00753        << maxIndex << " is out of bounds.";
00754     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00755   }
00756 #endif // TEUCHOS_DEBUG
00757       // Cast to nonconst, because Epetra_MultiVector's constructor
00758       // wants a nonconst int array argument.  It doesn't actually
00759       // change the entries of the array.
00760       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00761       return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, &tmpind[0], index.size())); 
00762       // return Teuchos::rcp (new Epetra_MultiVector (::Copy, mv, &tmpind[0], index.size())); 
00763     }
00764 
00765     static Teuchos::RCP<Epetra_MultiVector> 
00766     CloneCopy (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00767     { 
00768       const int inNumVecs = GetNumberVecs (mv);
00769       const int outNumVecs = index.size();
00770       const bool validRange = outNumVecs > 0 && index.lbound() >= 0 && 
00771   index.ubound() < inNumVecs;
00772       if (! validRange)
00773   {
00774     std::ostringstream os;
00775     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,"
00776       "index=[" << index.lbound() << ", " << index.ubound() << "]): ";
00777     TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00778            os.str() << "Column index range must be nonempty.");
00779     TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00780            os.str() << "Column index range must be nonnegative.");
00781     TEST_FOR_EXCEPTION(index.ubound() >= inNumVecs, std::invalid_argument,
00782            os.str() << "Column index range must not exceed "
00783            "number of vectors " << inNumVecs << " in the "
00784            "input multivector.");
00785   }
00786       return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, index.lbound(), index.size()));
00787     }
00788 
00794     static Teuchos::RCP<Epetra_MultiVector> 
00795     CloneViewNonConst (Epetra_MultiVector& mv, const std::vector<int>& index)
00796     { 
00797       const int inNumVecs = GetNumberVecs (mv);
00798       const int outNumVecs = index.size();
00799 
00800       // Simple, inexpensive tests of the index vector.
00801       TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00802        "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00803        "CloneViewNonConst(mv, index = {}): The output view "
00804        "must have at least one column.");
00805       if (outNumVecs > inNumVecs)
00806   {
00807     std::ostringstream os;
00808     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00809       "CloneViewNonConst(mv, index = {";
00810     for (int k = 0; k < outNumVecs - 1; ++k)
00811       os << index[k] << ", ";
00812     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00813        << " indices to view, but only " << inNumVecs << " columns of mv.";
00814     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00815   }
00816 #ifdef TEUCHOS_DEBUG
00817       // In debug mode, we perform more expensive tests of the index
00818       // vector, to ensure all the elements are in range.
00819       // Dereferencing the iterator is valid because index has length
00820       // > 0.
00821       const int minIndex = *std::min_element (index.begin(), index.end());
00822       const int maxIndex = *std::max_element (index.begin(), index.end());
00823 
00824       if (minIndex < 0)
00825   {
00826     std::ostringstream os;
00827     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00828       "CloneViewNonConst(mv, index = {";
00829     for (int k = 0; k < outNumVecs - 1; ++k)
00830       os << index[k] << ", ";
00831     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00832       "the smallest index " << minIndex << " is negative.";
00833     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00834   }
00835       if (maxIndex >= inNumVecs)
00836   {
00837     std::ostringstream os;
00838     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
00839       "CloneViewNonConst(mv, index = {";
00840     for (int k = 0; k < outNumVecs - 1; ++k)
00841       os << index[k] << ", ";
00842     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00843       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00844        << maxIndex << " is out of bounds.";
00845     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00846   }
00847 #endif // TEUCHOS_DEBUG
00848       // Cast to nonconst, because Epetra_MultiVector's constructor
00849       // wants a nonconst int array argument.  It doesn't actually
00850       // change the entries of the array.
00851       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00852       return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
00853     }
00854 
00855     static Teuchos::RCP<Epetra_MultiVector> 
00856     CloneViewNonConst (Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00857     { 
00858       const bool validRange = index.size() > 0 && 
00859   index.lbound() >= 0 && 
00860   index.ubound() < mv.NumVectors();
00861       if (! validRange)
00862   {
00863     std::ostringstream os;
00864     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
00865       "NonConst(mv,index=[" << index.lbound() << ", " << index.ubound() 
00866        << "]): ";
00867     TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00868            os.str() << "Column index range must be nonempty.");
00869     TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00870            os.str() << "Column index range must be nonnegative.");
00871     TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(), 
00872            std::invalid_argument,
00873            os.str() << "Column index range must not exceed "
00874            "number of vectors " << mv.NumVectors() << " in "
00875            "the input multivector.");
00876   }
00877       return Teuchos::rcp (new Epetra_MultiVector (View, mv, index.lbound(), index.size()));
00878     }
00879 
00885     static Teuchos::RCP<const Epetra_MultiVector> 
00886     CloneView (const Epetra_MultiVector& mv, const std::vector<int>& index)
00887     { 
00888       const int inNumVecs = GetNumberVecs (mv);
00889       const int outNumVecs = index.size();
00890 
00891       // Simple, inexpensive tests of the index vector.
00892       TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
00893        "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00894        "CloneView(mv, index = {}): The output view "
00895        "must have at least one column.");
00896       if (outNumVecs > inNumVecs)
00897   {
00898     std::ostringstream os;
00899     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00900       "CloneView(mv, index = {";
00901     for (int k = 0; k < outNumVecs - 1; ++k)
00902       os << index[k] << ", ";
00903     os << index[outNumVecs-1] << "}): There are " << outNumVecs 
00904        << " indices to view, but only " << inNumVecs << " columns of mv.";
00905     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00906   }
00907 #ifdef TEUCHOS_DEBUG
00908       // In debug mode, we perform more expensive tests of the index
00909       // vector, to ensure all the elements are in range.
00910       // Dereferencing the iterator is valid because index has length
00911       // > 0.
00912       const int minIndex = *std::min_element (index.begin(), index.end());
00913       const int maxIndex = *std::max_element (index.begin(), index.end());
00914 
00915       if (minIndex < 0)
00916   {
00917     std::ostringstream os;
00918     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00919       "CloneView(mv, index = {";
00920     for (int k = 0; k < outNumVecs - 1; ++k)
00921       os << index[k] << ", ";
00922     os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
00923       "the smallest index " << minIndex << " is negative.";
00924     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00925   }
00926       if (maxIndex >= inNumVecs)
00927   {
00928     std::ostringstream os;
00929     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
00930       "CloneView(mv, index = {";
00931     for (int k = 0; k < outNumVecs - 1; ++k)
00932       os << index[k] << ", ";
00933     os << index[outNumVecs-1] << "}): Indices must be strictly less than "
00934       "the number of vectors " << inNumVecs << " in mv; the largest index " 
00935        << maxIndex << " is out of bounds.";
00936     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00937   }
00938 #endif // TEUCHOS_DEBUG
00939       // Cast to nonconst, because Epetra_MultiVector's constructor
00940       // wants a nonconst int array argument.  It doesn't actually
00941       // change the entries of the array.
00942       std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
00943       return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
00944     }
00945 
00946     static Teuchos::RCP<Epetra_MultiVector> 
00947     CloneView (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
00948     { 
00949       const bool validRange = index.size() > 0 && 
00950   index.lbound() >= 0 && 
00951   index.ubound() < mv.NumVectors();
00952       if (! validRange)
00953   {
00954     std::ostringstream os;
00955     os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
00956       "(mv,index=[" << index.lbound() << ", " << index.ubound() 
00957        << "]): ";
00958     TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00959            os.str() << "Column index range must be nonempty.");
00960     TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00961            os.str() << "Column index range must be nonnegative.");
00962     TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(), 
00963            std::invalid_argument,
00964            os.str() << "Column index range must not exceed "
00965            "number of vectors " << mv.NumVectors() << " in "
00966            "the input multivector.");
00967   }
00968       return Teuchos::rcp (new Epetra_MultiVector(View, mv, index.lbound(), index.size()));
00969     }
00970 
00972 
00974 
00975 
00977     static int GetVecLength( const Epetra_MultiVector& mv )
00978     { return mv.GlobalLength(); }
00979 
00981     static int GetNumberVecs( const Epetra_MultiVector& mv )
00982     { return mv.NumVectors(); }
00983 
00984     static bool HasConstantStride( const Epetra_MultiVector& mv )
00985     { return mv.ConstantStride(); }
00987 
00989 
00990 
00993     static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A, 
00994                                  const Teuchos::SerialDenseMatrix<int,double>& B, 
00995                                  double beta, Epetra_MultiVector& mv )
00996     { 
00997       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
00998       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
00999 
01000       TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
01001           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Multiply() returned a nonzero value.");
01002     }
01003 
01006     static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
01007     { 
01008       // epetra mv.Update(alpha,A,beta,B,gamma) will check 
01009       //   alpha == 0.0 
01010       // and 
01011       //   beta == 0.0 
01012       // and based on this will call 
01013       //   mv.Update(beta,B,gamma) 
01014       // or 
01015       //   mv.Update(alpha,A,gamma)
01016       //
01017       // mv.Update(alpha,A,gamma) 
01018       // will then check for one of 
01019       //   gamma == 0
01020       // or 
01021       //   gamma == 1
01022       // or 
01023       //   alpha == 1 
01024       // in that order. however, it will not look for the combination
01025       //   alpha == 1 and gamma = 0
01026       // which is a common use case when we wish to assign 
01027       //   mv = A   (in which case alpha == 1, beta == gamma == 0)
01028       // or 
01029       //   mv = B   (in which case beta == 1, alpha == gamma == 0)
01030       // 
01031       // therefore, we will check for these use cases ourselves
01032       if (beta == 0.0) {
01033         if (alpha == 1.0) {
01034           // assign
01035           mv = A;
01036         }
01037         else {
01038           // single update
01039           TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
01040               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
01041         }
01042       }
01043       else if (alpha == 0.0) {
01044         if (beta == 1.0) {
01045           // assign 
01046           mv = B;
01047         }
01048         else {
01049           // single update
01050           TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
01051               "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
01052         }
01053       }
01054       else {
01055         // double update
01056         TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
01057             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
01058       }
01059     }
01060 
01063     static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
01064 #ifdef HAVE_ANASAZI_EXPERIMENTAL
01065                           , ConjType conj = Anasazi::CONJ
01066 #endif
01067                         )
01068     { 
01069       Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
01070       Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
01071       
01072       TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
01073           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
01074     }
01075     
01078     static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
01079 #ifdef HAVE_ANASAZI_EXPERIMENTAL
01080                       , ConjType conj = Anasazi::CONJ
01081 #endif
01082                       )
01083     {
01084 #ifdef TEUCHOS_DEBUG
01085       TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
01086           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
01087       TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
01088           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
01089 #endif
01090       TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
01091           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");     
01092     }
01093 
01095 
01096 
01097 
01101     static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
01102     { 
01103 #ifdef TEUCHOS_DEBUG
01104       TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
01105           "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
01106 #endif
01107       TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
01108           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value."); 
01109     }
01110 
01112     
01114 
01115 
01117     static void 
01118     SetBlock (const Epetra_MultiVector& A, 
01119         const std::vector<int>& index, 
01120         Epetra_MultiVector& mv)
01121     { 
01122       const int inNumVecs = GetNumberVecs (A);
01123       const int outNumVecs = index.size();
01124 
01125       // FIXME (mfh 13 Jan 2011) Belos allows A to have more columns
01126       // than index.size(), in which case we just take the first
01127       // index.size() columns of A.  Anasazi requires that A have the
01128       // same number of columns as index.size().  Changing Anasazi's
01129       // behavior should not break existing Anasazi solvers, but the
01130       // tests need to be done.
01131       if (inNumVecs != outNumVecs)
01132   {
01133     std::ostringstream os;
01134     os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
01135       "SetBlock(A, mv, index = {";
01136     if (outNumVecs > 0)
01137       {
01138         for (int k = 0; k < outNumVecs - 1; ++k)
01139     os << index[k] << ", ";
01140         os << index[outNumVecs-1];
01141       }
01142     os << "}): A has only " << inNumVecs << " columns, but there are "
01143        << outNumVecs << " indices in the index vector.";
01144     TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
01145   }
01146       // Make a view of the columns of mv indicated by the index std::vector.
01147       Teuchos::RCP<Epetra_MultiVector> mv_view = CloneViewNonConst (mv, index);
01148 
01149       // View of columns [0, outNumVecs-1] of the source multivector A.
01150       // If A has fewer columns than mv_view, then create a view of
01151       // the first outNumVecs columns of A.
01152       Teuchos::RCP<const Epetra_MultiVector> A_view;
01153       if (outNumVecs == inNumVecs)
01154   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
01155       else
01156   A_view = CloneView (A, Teuchos::Range1D(0, outNumVecs - 1));
01157 
01158       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01159       // copies the data directly, ignoring the underlying
01160       // Epetra_Map(s).  If A and mv don't have the same data
01161       // distribution (Epetra_Map), this may result in incorrect or
01162       // undefined behavior.  Epetra_MultiVector::Update() also
01163       // ignores the Epetra_Maps, so we might as well just use the
01164       // (perhaps slightly cheaper) Assign() method via operator=().
01165       *mv_view = *A_view;
01166     }
01167 
01168     static void 
01169     SetBlock (const Epetra_MultiVector& A, 
01170         const Teuchos::Range1D& index, 
01171         Epetra_MultiVector& mv)
01172     { 
01173       const int numColsA = A.NumVectors();
01174       const int numColsMv = mv.NumVectors();
01175       // 'index' indexes into mv; it's the index set of the target.
01176       const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
01177       // We can't take more columns out of A than A has.
01178       const bool validSource = index.size() <= numColsA;
01179 
01180       if (! validIndex || ! validSource)
01181   {
01182     std::ostringstream os;
01183     os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::SetBlock"
01184       "(A, index=[" << index.lbound() << ", " << index.ubound() << "], "
01185       "mv): ";
01186     TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
01187            os.str() << "Range lower bound must be nonnegative.");
01188     TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
01189            os.str() << "Range upper bound must be less than "
01190            "the number of columns " << numColsA << " in the "
01191            "'mv' output argument.");
01192     TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
01193            os.str() << "Range must have no more elements than"
01194            " the number of columns " << numColsA << " in the "
01195            "'A' input argument.");
01196     TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
01197   }
01198 
01199       // View of columns [index.lbound(), index.ubound()] of the
01200       // target multivector mv.  We avoid view creation overhead by
01201       // only creating a view if the index range is different than [0,
01202       // (# columns in mv) - 1].
01203       Teuchos::RCP<Epetra_MultiVector> mv_view;
01204       if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
01205   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
01206       else
01207   mv_view = CloneViewNonConst (mv, index);
01208 
01209       // View of columns [0, index.size()-1] of the source multivector
01210       // A.  If A has fewer columns than mv_view, then create a view
01211       // of the first index.size() columns of A.
01212       Teuchos::RCP<const Epetra_MultiVector> A_view;
01213       if (index.size() == numColsA)
01214   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
01215       else
01216   A_view = CloneView (A, Teuchos::Range1D(0, index.size()-1));
01217 
01218       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01219       // copies the data directly, ignoring the underlying
01220       // Epetra_Map(s).  If A and mv don't have the same data
01221       // distribution (Epetra_Map), this may result in incorrect or
01222       // undefined behavior.  Epetra_MultiVector::Update() also
01223       // ignores the Epetra_Maps, so we might as well just use the
01224       // (perhaps slightly cheaper) Assign() method via operator=().
01225       *mv_view = *A_view; 
01226     }
01227 
01228     static void 
01229     Assign (const Epetra_MultiVector& A, 
01230       Epetra_MultiVector& mv)
01231     {
01232       const int numColsA = GetNumberVecs (A);
01233       const int numColsMv = GetNumberVecs (mv);
01234       if (numColsA > numColsMv)
01235   {
01236     std::ostringstream os;
01237     os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::Assign"
01238       "(A, mv): ";
01239     TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
01240            os.str() << "Input multivector 'A' has " 
01241            << numColsA << " columns, but output multivector "
01242            "'mv' has only " << numColsMv << " columns.");
01243     TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
01244   }
01245       // View of the first [0, numColsA-1] columns of mv.
01246       Teuchos::RCP<Epetra_MultiVector> mv_view;
01247       if (numColsMv == numColsA)
01248   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
01249       else // numColsMv > numColsA
01250   mv_view = CloneView (mv, Teuchos::Range1D(0, numColsA - 1));
01251       
01252       // Assignment calls Epetra_MultiVector::Assign(), which deeply
01253       // copies the data directly, ignoring the underlying
01254       // Epetra_Map(s).  If A and mv don't have the same data
01255       // distribution (Epetra_Map), this may result in incorrect or
01256       // undefined behavior.  Epetra_MultiVector::Update() also
01257       // ignores the Epetra_Maps, so we might as well just use the
01258       // (perhaps slightly cheaper) Assign() method via operator=().
01259       *mv_view = A;
01260     }
01261 
01264     static void MvScale ( Epetra_MultiVector& mv, double alpha ) 
01265     { 
01266       TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
01267           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value."); 
01268     }
01269 
01272     static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
01273     { 
01274       // Check to make sure the vector is as long as the multivector has columns.
01275       int numvecs = mv.NumVectors();
01276 #ifdef TEUCHOS_DEBUG
01277       TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
01278                           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
01279 #endif
01280       for (int i=0; i<numvecs; i++) {
01281         TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i])!=0, EpetraMultiVecFailure,
01282             "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
01283       }
01284     }
01285 
01288     static void MvRandom( Epetra_MultiVector& mv )
01289     { 
01290       TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
01291           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
01292     }
01293 
01296     static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
01297     { 
01298       TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
01299           "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
01300     }
01301 
01303 
01305 
01306 
01309     static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
01310     { os << mv << std::endl; }
01311 
01313 
01314 #if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
01315 #  if defined(HAVE_TPETRA_EPETRA)
01316 
01317 
01318 
01319 
01320 
01321     typedef Epetra::TsqrAdaptor tsqr_adaptor_type;
01322 #  endif // defined(HAVE_TPETRA_EPETRA)
01323 #endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
01324   };        
01325 
01327   //
01328   // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
01329   //
01331 
01343   template <> 
01344   class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
01345   {
01346   public:
01347     
01351     static void Apply ( const Epetra_Operator& Op, 
01352                         const Epetra_MultiVector& x, 
01353                         Epetra_MultiVector& y )
01354     { 
01355 #ifdef TEUCHOS_DEBUG
01356       TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
01357           "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
01358 #endif
01359       int ret = Op.Apply(x,y);
01360       TEST_FOR_EXCEPTION(ret != 0, OperatorError, 
01361           "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
01362     }
01363     
01364   };
01365   
01366 } // end of Anasazi namespace 
01367 
01368 #endif 
01369 // end of file ANASAZI_EPETRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends