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