Anasazi Version of the Day
AnasaziTpetraAdapter.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 
00029 #ifndef ANASAZI_TPETRA_ADAPTER_HPP
00030 #define ANASAZI_TPETRA_ADAPTER_HPP
00031 
00032 #include <Kokkos_NodeTrace.hpp>
00033 
00043 
00044 #include <Tpetra_MultiVector.hpp>
00045 #include <Tpetra_Operator.hpp>
00046 #include <Teuchos_Assert.hpp>
00047 #include <Teuchos_ScalarTraits.hpp>
00048 #include <Teuchos_Array.hpp>
00049 #include <Teuchos_DefaultSerialComm.hpp>
00050 
00051 #include <AnasaziConfigDefs.hpp>
00052 #include <AnasaziTypes.hpp>
00053 #include <AnasaziMultiVecTraits.hpp>
00054 #include <AnasaziOperatorTraits.hpp>
00055 #include <Kokkos_NodeAPIConfigDefs.hpp>
00056 
00057 #ifdef HAVE_ANASAZI_TSQR
00058 #  include <Tpetra_TsqrAdaptor.hpp>
00059 #endif // HAVE_ANASAZI_TSQR
00060 
00061 
00062 namespace Anasazi {
00063 
00065   //
00066   // Implementation of the Anasazi::MultiVecTraits for Tpetra::MultiVector.
00067   //
00069 
00074   template<class Scalar, class LO, class GO, class Node>
00075   class MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node> >
00076   {
00077   public:
00078 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
00079     static Teuchos::RCP<Teuchos::Time> mvTimesMatAddMvTimer_, mvTransMvTimer_;
00080 #endif
00081 
00082     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > Clone( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const int numvecs )
00083     { 
00084       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>(mv.getMap(),numvecs));
00085     }
00086 
00087     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00088     {
00089       KOKKOS_NODE_TRACE("Anasazi::MVT::CloneCopy(MV)")
00090       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>( mv ) ); 
00091     }
00092 
00093     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00094     { 
00095       KOKKOS_NODE_TRACE("Anasazi::MVT::CloneCopy(MV,ind)")
00096       TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00097           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): numvecs must be greater than zero.");
00098 #ifdef HAVE_TPETRA_DEBUG
00099       TEUCHOS_TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::runtime_error,
00100           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be >= zero.");
00101       TEUCHOS_TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::runtime_error,
00102           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be < mv.getNumVectors().");
00103 #endif
00104       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00105         if (index[j] != index[j-1]+1) {
00106           // not contiguous; short circuit
00107           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00108           return mv.subCopy(stinds());
00109         }
00110       }
00111       // contiguous
00112       return mv.subCopy(Teuchos::Range1D(index.front(),index.back()));
00113     }
00114 
00115 
00116     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > 
00117     CloneCopy (const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, 
00118          const Teuchos::Range1D& index)
00119     { 
00120       KOKKOS_NODE_TRACE("Anasazi::MVT::CloneCopy(MV,ind)")
00121       const bool validRange = index.size() > 0 && 
00122   index.lbound() >= 0 && 
00123   index.ubound() < mv.getNumVectors();
00124       if (! validRange)
00125   {
00126     std::ostringstream os;
00127     os << "Anasazi::MultiVecTraits<Scalar, Tpetra::MultiVector<...> >::"
00128       "CloneCopy(mv,index=[" << index.lbound() << ", " << index.ubound() 
00129        << "]): ";
00130     TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00131            os.str() << "Empty index range is not allowed.");
00132     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00133            os.str() << "Index range includes negative "
00134            "index/ices, which is not allowed.");
00135     // Range1D bounds are signed; size_t is unsigned.
00136     TEUCHOS_TEST_FOR_EXCEPTION((size_t) index.ubound() >= mv.getNumVectors(), 
00137            std::invalid_argument, 
00138            os.str() << "Index range exceeds number of vectors " 
00139            << mv.getNumVectors() << " in the input multivector.");
00140     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, 
00141            os.str() << "Should never get here!");
00142   }
00143       return mv.subCopy (index);
00144     }
00145 
00146 
00147     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneViewNonConst( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00148     {
00149       TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00150           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): numvecs must be greater than zero.");
00151 #ifdef HAVE_TPETRA_DEBUG
00152       TEUCHOS_TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00153           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): indices must be >= zero.");
00154       TEUCHOS_TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00155           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): indices must be < mv.getNumVectors().");
00156 #endif
00157       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00158         if (index[j] != index[j-1]+1) {
00159           // not contiguous; short circuit
00160           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00161           return mv.subViewNonConst(stinds);
00162         }
00163       }
00164       // contiguous
00165       return mv.subViewNonConst(Teuchos::Range1D(index.front(),index.back()));
00166     }
00167 
00168 
00169     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > 
00170     CloneViewNonConst (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, 
00171            const Teuchos::Range1D& index)
00172     {
00173       // FIXME (mfh 11 Jan 2011) Should check for overflowing int!
00174       const int numCols = static_cast<int> (mv.getNumVectors());
00175       const bool validRange = index.size() > 0 && 
00176   index.lbound() >= 0 && index.ubound() < numCols;
00177       if (! validRange)
00178   {
00179     std::ostringstream os;
00180     os << "Anasazi::MultiVecTraits<Scalar, Tpetra::MultiVector<...> >::"
00181       "CloneViewNonConst(mv,index=[" << index.lbound() << ", " 
00182        << index.ubound() << "]): ";
00183     TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00184            os.str() << "Empty index range is not allowed.");
00185     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00186            os.str() << "Index range includes negative "
00187            "index/ices, which is not allowed.");
00188     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numCols, std::invalid_argument, 
00189            os.str() << "Index range exceeds number of "
00190            "vectors " << numCols << " in the input "
00191            "multivector.");
00192     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, 
00193            os.str() << "Should never get here!");
00194   }
00195       return mv.subViewNonConst (index);
00196     }
00197 
00198 
00199     static Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneView(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00200     {
00201       TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00202           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
00203 #ifdef HAVE_TPETRA_DEBUG
00204       TEUCHOS_TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00205           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be >= zero.");
00206       TEUCHOS_TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00207           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be < mv.getNumVectors().");
00208 #endif
00209       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00210         if (index[j] != index[j-1]+1) {
00211           // not contiguous; short circuit
00212           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00213           return mv.subView(stinds);
00214         }
00215       }
00216       // contiguous
00217       return mv.subView(Teuchos::Range1D(index.front(),index.back()));
00218     }
00219 
00220     static Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > 
00221     CloneView (const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, 
00222          const Teuchos::Range1D& index)
00223     {
00224       // FIXME (mfh 11 Jan 2011) Should check for overflowing int!
00225       const int numCols = static_cast<int> (mv.getNumVectors());
00226       const bool validRange = index.size() > 0 && 
00227   index.lbound() >= 0 && index.ubound() < numCols;
00228       if (! validRange)
00229   {
00230     std::ostringstream os;
00231     os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<...> >::"
00232       "CloneView(mv, index=[" << index.lbound() << ", " 
00233        << index.ubound() << "]): ";
00234     TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
00235            os.str() << "Empty index range is not allowed.");
00236     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00237            os.str() << "Index range includes negative "
00238            "index/ices, which is not allowed.");
00239     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numCols, std::invalid_argument, 
00240            os.str() << "Index range exceeds number of "
00241            "vectors " << numCols << " in the input "
00242            "multivector.");
00243     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, 
00244            os.str() << "Should never get here!");
00245   }
00246       return mv.subView (index);
00247     }
00248 
00249     static int GetVecLength( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00250     { return mv.getGlobalLength(); }
00251 
00252     static int GetNumberVecs( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00253     { return mv.getNumVectors(); }
00254 
00255     static void MvTimesMatAddMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, 
00256                                  const Teuchos::SerialDenseMatrix<int,Scalar>& B, 
00257                                  Scalar beta, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00258     {
00259       KOKKOS_NODE_TRACE("Anasazi::MVT::MvTimesMatAddMv()")
00260 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
00261       Teuchos::TimeMonitor lcltimer(*mvTimesMatAddMvTimer_);
00262 #endif
00263       // create local map
00264       Teuchos::SerialComm<int> scomm;
00265       Tpetra::Map<LO,GO,Node> LocalMap(B.numRows(), 0, Teuchos::rcpFromRef< const Teuchos::Comm<int> >(scomm), Tpetra::LocallyReplicated, A.getMap()->getNode());
00266       // encapsulate Teuchos::SerialDenseMatrix data in ArrayView
00267       Teuchos::ArrayView<const Scalar> Bvalues(B.values(),B.stride()*B.numCols());
00268       // create locally replicated MultiVector with a copy of this data
00269       Tpetra::MultiVector<Scalar,LO,GO,Node> B_mv(Teuchos::rcpFromRef(LocalMap),Bvalues,B.stride(),B.numCols());
00270       // multiply
00271       mv.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
00272     }
00273 
00274     static void MvAddMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, Scalar beta, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00275     {
00276       mv.update(alpha,A,beta,B,Teuchos::ScalarTraits<Scalar>::zero());
00277     }
00278 
00279     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha )
00280     { mv.scale(alpha); }
00281 
00282     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<Scalar>& alphas )
00283     { mv.scale(alphas); }
00284 
00285     static void MvTransMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, Teuchos::SerialDenseMatrix<int,Scalar>& C)
00286     { 
00287       KOKKOS_NODE_TRACE("Anasazi::MVT::MvTransMv()")
00288 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
00289       Teuchos::TimeMonitor lcltimer(*mvTransMvTimer_);
00290 #endif
00291       // form alpha * A^H * B, then copy into SDM
00292       // we will create a multivector C_mv from a a local map
00293       // this map has a serial comm, the purpose being to short-circuit the MultiVector::reduce() call at the end of MultiVector::multiply()
00294       // otherwise, the reduced multivector data would be copied back to the GPU, only to turn around and have to get it back here.
00295       // this saves us a round trip for this data.
00296       const int numRowsC = C.numRows(),
00297                 numColsC = C.numCols(),
00298                 strideC  = C.stride();
00299       Teuchos::SerialComm<int> scomm;
00300       // create local map with serial comm
00301       Tpetra::Map<LO,GO,Node> LocalMap(numRowsC, 0, Teuchos::rcpFromRef< const Teuchos::Comm<int> >(scomm), Tpetra::LocallyReplicated, A.getMap()->getNode());
00302       // create local multivector to hold the result
00303       const bool INIT_TO_ZERO = true;
00304       Tpetra::MultiVector<Scalar,LO,GO,Node> C_mv(Teuchos::rcpFromRef(LocalMap),numColsC, INIT_TO_ZERO);
00305       // multiply result into local multivector
00306       C_mv.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,alpha,A,B,Teuchos::ScalarTraits<Scalar>::zero());
00307       // get comm
00308       Teuchos::RCP< const Teuchos::Comm<int> > pcomm = A.getMap()->getComm();
00309       // create arrayview encapsulating the Teuchos::SerialDenseMatrix
00310       Teuchos::ArrayView<Scalar> C_view(C.values(),strideC*numColsC);
00311       if (pcomm->getSize() == 1) {
00312         // no accumulation to do; simply extract the multivector data into C
00313         // extract a copy of the result into the array view (and therefore, the SerialDenseMatrix)
00314         C_mv.get1dCopy(C_view,strideC);
00315       }  
00316       else {
00317         // get a const host view of the data in C_mv
00318         Teuchos::ArrayRCP<const Scalar> C_mv_view = C_mv.get1dView();
00319         if (strideC == numRowsC) {
00320           // sumall into C
00321           Teuchos::reduceAll<int,Scalar>(*pcomm,Teuchos::REDUCE_SUM,numColsC*numRowsC,C_mv_view.getRawPtr(),C_view.getRawPtr());
00322         }
00323         else {
00324           // sumall into temp, copy into C
00325           Teuchos::Array<Scalar> destBuff(numColsC*numRowsC);
00326           Teuchos::reduceAll<int,Scalar>(*pcomm,Teuchos::REDUCE_SUM,numColsC*numRowsC,C_mv_view.getRawPtr(),destBuff.getRawPtr());
00327           for (int j=0; j < numColsC; ++j) {
00328             for (int i=0; i < numRowsC; ++i) {
00329               C_view[strideC*j+i] = destBuff[numRowsC*j+i];
00330             }
00331           }
00332         }
00333       }
00334     }
00335 
00336     static void MvDot( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, std::vector<Scalar> &dots)
00337     {
00338       TEUCHOS_TEST_FOR_EXCEPTION(A.getNumVectors() != B.getNumVectors(),std::invalid_argument,
00339           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): A and B must have the same number of vectors.");
00340 #ifdef HAVE_TPETRA_DEBUG
00341       TEUCHOS_TEST_FOR_EXCEPTION(dots.size() < (typename std::vector<int>::size_type)A.getNumVectors(),std::invalid_argument,
00342           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): dots must have room for all dot products.");
00343 #endif
00344       Teuchos::ArrayView<Scalar> av(dots);
00345       A.dot(B,av(0,A.getNumVectors()));
00346     }
00347 
00348     static void MvNorm(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
00349     { 
00350 #ifdef HAVE_TPETRA_DEBUG
00351       TEUCHOS_TEST_FOR_EXCEPTION(normvec.size() < (typename std::vector<int>::size_type)mv.getNumVectors(),std::invalid_argument,
00352           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvNorm(mv,normvec): normvec must have room for all norms.");
00353 #endif
00354       Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> av(normvec);
00355       mv.norm2(av(0,mv.getNumVectors()));
00356     }
00357 
00358     static void SetBlock( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const std::vector<int>& index, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00359     {
00360       KOKKOS_NODE_TRACE("Anasazi::MVT::SetBlock()")
00361 #ifdef HAVE_TPETRA_DEBUG
00362       TEUCHOS_TEST_FOR_EXCEPTION((typename std::vector<int>::size_type)A.getNumVectors() < index.size(),std::invalid_argument,
00363           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::SetBlock(A,index,mv): index must be the same size as A.");
00364 #endif
00365       Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > mvsub = CloneViewNonConst(mv,index);
00366       if ((typename std::vector<int>::size_type)A.getNumVectors() > index.size()) {
00367         Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > Asub = A.subView(Teuchos::Range1D(0,index.size()-1));
00368         (*mvsub) = (*Asub);
00369       }
00370       else {
00371         (*mvsub) = A;
00372       }
00373       mvsub = Teuchos::null;
00374     }
00375 
00376     static void
00377     SetBlock (const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, 
00378         const Teuchos::Range1D& index, 
00379         Tpetra::MultiVector<Scalar,LO,GO,Node>& mv)
00380     {
00381       KOKKOS_NODE_TRACE("Anasazi::MVT::SetBlock()")
00382 
00383       // Range1D bounds are signed; size_t is unsigned.
00384       // Assignment of Tpetra::MultiVector is a deep copy.
00385 
00386       // Tpetra::MultiVector::getNumVectors() returns size_t.  It's
00387       // fair to assume that the number of vectors won't overflow int,
00388       // since the typical use case of multivectors involves few
00389       // columns, but it's friendly to check just in case.
00390       const size_t maxInt = static_cast<size_t> (Teuchos::OrdinalTraits<int>::max());
00391       const bool overflow = maxInt < A.getNumVectors() && maxInt < mv.getNumVectors();
00392       if (overflow)
00393   {
00394     std::ostringstream os;
00395     os << "Anasazi::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
00396       "> >::SetBlock(A, index=[" << index.lbound() << ", " 
00397        << index.ubound() << "], mv): ";
00398     TEUCHOS_TEST_FOR_EXCEPTION(maxInt < A.getNumVectors(), std::range_error,
00399            os.str() << "Number of columns in the input multi"
00400            "vector 'A' (a size_t) overflows int.");
00401     TEUCHOS_TEST_FOR_EXCEPTION(maxInt < mv.getNumVectors(), std::range_error,
00402            os.str() << "Number of columns in the output multi"
00403            "vector 'mv' (a size_t) overflows int.");
00404   }
00405       // We've already validated the static casts above.
00406       const int numColsA = static_cast<int> (A.getNumVectors());
00407       const int numColsMv = static_cast<int> (mv.getNumVectors());
00408       // 'index' indexes into mv; it's the index set of the target.
00409       const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
00410       // We can't take more columns out of A than A has.
00411       const bool validSource = index.size() <= numColsA;
00412 
00413       if (! validIndex || ! validSource)
00414   {
00415     std::ostringstream os;
00416     os << "Anasazi::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
00417       "> >::SetBlock(A, index=[" << index.lbound() << ", " 
00418        << index.ubound() << "], mv): ";
00419     TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
00420            os.str() << "Range lower bound must be nonnegative.");
00421     TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
00422            os.str() << "Range upper bound must be less than "
00423            "the number of columns " << numColsA << " in the "
00424            "'mv' output argument.");
00425     TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
00426            os.str() << "Range must have no more elements than"
00427            " the number of columns " << numColsA << " in the "
00428            "'A' input argument.");
00429     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00430   }
00431       typedef Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > MV_ptr;
00432       typedef Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > const_MV_ptr;
00433 
00434       // View of the relevant column(s) of the target multivector mv.
00435       // We avoid view creation overhead by only creating a view if
00436       // the index range is different than [0, (# columns in mv) - 1].
00437       MV_ptr mv_view;
00438       if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
00439   mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
00440       else
00441   mv_view = CloneViewNonConst (mv, index);
00442 
00443       // View of the relevant column(s) of the source multivector A.
00444       // If A has fewer columns than mv_view, then create a view of
00445       // the first index.size() columns of A.
00446       const_MV_ptr A_view;
00447       if (index.size() == numColsA)
00448   A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
00449       else
00450   A_view = CloneView (A, Teuchos::Range1D(0, index.size()-1));
00451 
00452       // Assignment of Tpetra::MultiVector objects via operator=()
00453       // assumes that both arguments have compatible Maps.  If
00454       // HAVE_TPETRA_DEBUG is defined at compile time, operator=()
00455       // will throw an std::runtime_error if the Maps are
00456       // incompatible.
00457       *mv_view = *A_view; 
00458     }
00459 
00460     static void
00461     Assign (const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, 
00462       Tpetra::MultiVector<Scalar,LO,GO,Node>& mv)
00463     {
00464       KOKKOS_NODE_TRACE("Anasazi::MVT::Assign()")
00465 
00466       // Range1D bounds are signed; size_t is unsigned.
00467       // Assignment of Tpetra::MultiVector is a deep copy.
00468 
00469       // Tpetra::MultiVector::getNumVectors() returns size_t.  It's
00470       // fair to assume that the number of vectors won't overflow int,
00471       // since the typical use case of multivectors involves few
00472       // columns, but it's friendly to check just in case.
00473       const size_t maxInt = static_cast<size_t> (Teuchos::OrdinalTraits<int>::max());
00474       const bool overflow = maxInt < A.getNumVectors() && maxInt < mv.getNumVectors();
00475       if (overflow)
00476   {
00477     std::ostringstream os;
00478     os << "Anasazi::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
00479       "> >::Assign(A, mv): ";
00480     TEUCHOS_TEST_FOR_EXCEPTION(maxInt < A.getNumVectors(), std::range_error,
00481            os.str() << "Number of columns in the input multi"
00482            "vector 'A' (a size_t) overflows int.");
00483     TEUCHOS_TEST_FOR_EXCEPTION(maxInt < mv.getNumVectors(), std::range_error,
00484            os.str() << "Number of columns in the output multi"
00485            "vector 'mv' (a size_t) overflows int.");
00486     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00487   }
00488       // We've already validated the static casts above.
00489       const int numColsA = static_cast<int> (A.getNumVectors());
00490       const int numColsMv = static_cast<int> (mv.getNumVectors());
00491       if (numColsA > numColsMv)
00492   {
00493     std::ostringstream os;
00494     os << "Anasazi::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
00495       "> >::Assign(A, mv): ";
00496     TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
00497            os.str() << "Input multivector 'A' has " 
00498            << numColsA << " columns, but output multivector "
00499            "'mv' has only " << numColsMv << " columns.");
00500     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
00501   }
00502       // Assignment of Tpetra::MultiVector objects via operator=()
00503       // assumes that both arguments have compatible Maps.  If
00504       // HAVE_TPETRA_DEBUG is defined at compile time, operator=()
00505       // will throw an std::runtime_error if the Maps are
00506       // incompatible.
00507       if (numColsA == numColsMv)
00508   mv = A;
00509       else
00510   {
00511     Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > mv_view = 
00512       CloneViewNonConst (mv, Teuchos::Range1D(0, numColsA-1));
00513     *mv_view = A;
00514   }
00515     }
00516 
00517     static void MvRandom( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00518     { 
00519       KOKKOS_NODE_TRACE("Anasazi::MVT::randomize()")
00520       mv.randomize(); 
00521     }
00522 
00523     static void MvInit( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero() )
00524     { mv.putScalar(alpha); }
00525 
00526     static void MvPrint( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::ostream& os )
00527     { mv.print(os); }
00528 
00529 #ifdef HAVE_ANASAZI_TSQR
00530 
00531 
00532 
00533     typedef Tpetra::TsqrAdaptor< Tpetra::MultiVector< Scalar, LO, GO, Node > > tsqr_adaptor_type;
00534 #endif // HAVE_ANASAZI_TSQR
00535   };        
00536 
00538   //
00539   // Implementation of the Anasazi::OperatorTraits for Tpetra::Operator.
00540   //
00542 
00543   template <class Scalar, class LO, class GO, class Node> 
00544   class OperatorTraits < Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node>, Tpetra::Operator<Scalar,LO,GO,Node> >
00545   {
00546   public:
00547     static void Apply ( const Tpetra::Operator<Scalar,LO,GO,Node> & Op, 
00548                         const Tpetra::MultiVector<Scalar,LO,GO,Node> & X,
00549                               Tpetra::MultiVector<Scalar,LO,GO,Node> & Y)
00550     { 
00551       Op.apply(X,Y,Teuchos::NO_TRANS);
00552     }
00553   };
00554 
00555 } // end of Anasazi namespace 
00556 
00557 #endif 
00558 // end of file ANASAZI_TPETRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends