Tpetra Matrix/Vector Services Version of the Day
Tpetra_MultiVector_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) 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 TPETRA_MULTIVECTOR_DEF_HPP
00030 #define TPETRA_MULTIVECTOR_DEF_HPP
00031 
00032 #include <Kokkos_NodeTrace.hpp>
00033 
00034 #include <Teuchos_TestForException.hpp>
00035 #include <Teuchos_as.hpp>
00036 
00037 #include "Tpetra_Vector.hpp"
00038 
00039 #ifdef DOXYGEN_USE_ONLY
00040   #include "Tpetra_MultiVector_decl.hpp"
00041 #endif
00042 
00043 namespace Tpetra {
00044 
00045   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00046   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(
00047         const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 
00048         size_t NumVectors, 
00049         bool zeroOut  /* default is true */
00050   ) 
00051   : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), 
00052     lclMV_(map->getNode()) {
00053     TEST_FOR_EXCEPTION(NumVectors < 1, std::invalid_argument,
00054         "Tpetra::MultiVector::MultiVector(): NumVectors must be strictly positive.");
00055     const size_t myLen = getLocalLength();
00056     if (myLen > 0) {
00057       Teuchos::RCP<Node> node = map->getNode();
00058       Teuchos::ArrayRCP<Scalar> data = node->template allocBuffer<Scalar>(myLen*NumVectors);
00059       MVT::initializeValues(lclMV_,myLen,NumVectors,data,myLen);
00060       if (zeroOut) {
00061         MVT::Init(lclMV_, Teuchos::ScalarTraits<Scalar>::zero());
00062       }
00063     }
00064     else {
00065       MVT::initializeValues(lclMV_,0,NumVectors,Teuchos::null,0);
00066     }
00067   }
00068 
00069   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00070   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source) 
00071   : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(source), lclMV_(MVT::getNode(source.lclMV_)) {
00072     // copy data from the source MultiVector into this multivector
00073     Teuchos::RCP<Node> node = MVT::getNode(source.lclMV_);
00074     const LocalOrdinal myLen = getLocalLength();
00075     const size_t numVecs = source.getNumVectors();
00076     if (myLen > 0) {
00077       // allocate data
00078       Teuchos::ArrayRCP<Scalar> data = node->template allocBuffer<Scalar>(myLen*numVecs);
00079       MVT::initializeValues(lclMV_,myLen,numVecs,data,myLen);
00080       // copy data
00081       {
00082         Teuchos::ArrayRCP<Scalar> dstdata = data;
00083         Teuchos::ArrayRCP<const Scalar> srcdata = MVT::getValues(source.lclMV_);
00084         for (size_t j = 0; j < numVecs; ++j) {
00085           Teuchos::ArrayRCP<const Scalar> srcj = source.getSubArrayRCP(srcdata,j);
00086           KOKKOS_NODE_TRACE("MultiVector::MultiVector(MV)")
00087           node->template copyBuffers<Scalar>(myLen,srcj,dstdata);
00088           dstdata += myLen;
00089         }
00090       }
00091     }
00092     else {
00093       MVT::initializeValues(lclMV_,0,numVecs,Teuchos::null,0);
00094     }
00095   }
00096 
00097   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00098   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(
00099                         const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 
00100                         const Teuchos::ArrayView<const Scalar> &A, size_t LDA, 
00101                         size_t NumVectors)
00102   : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()) {
00103     TEST_FOR_EXCEPTION(NumVectors < 1, std::invalid_argument,
00104         "Tpetra::MultiVector::MultiVector(A,LDA): NumVectors must be strictly positive.");
00105     const size_t myLen = getLocalLength();
00106     using Teuchos::ArrayRCP;
00107     TEST_FOR_EXCEPTION(LDA < myLen, std::runtime_error,
00108         "Tpetra::MultiVector::MultiVector(A,LDA): LDA must be large enough to accomodate the local entries.");
00109 #ifdef HAVE_TPETRA_DEBUG
00110     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(A.size()) < LDA*(NumVectors-1)+myLen, std::runtime_error,
00111         "Tpetra::MultiVector::MultiVector(A,LDA): A does not contain enough data to specify the entries in this.");
00112 #endif
00113     if (myLen > 0) {
00114       Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00115       Teuchos::ArrayRCP<Scalar> mydata = node->template allocBuffer<Scalar>(myLen*NumVectors);
00116       MVT::initializeValues(lclMV_,myLen,NumVectors,mydata,myLen);
00117       KOKKOS_NODE_TRACE("MultiVector::MultiVector(1D)")
00118       Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::WriteOnly,myLen*NumVectors,mydata);
00119       typename Teuchos::ArrayView<const Scalar>::iterator srcit = A.begin();
00120       for (size_t j = 0; j < NumVectors; ++j) {
00121         std::copy(srcit,srcit+myLen,myview);
00122         srcit += LDA;
00123         myview += myLen;
00124       }
00125       mydata = Teuchos::null;
00126       myview = Teuchos::null;
00127     }
00128     else {
00129       MVT::initializeValues(lclMV_,0,NumVectors,Teuchos::null,0);
00130     }
00131   }
00132 
00133   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00134   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 
00135                                                                    const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar> > &ArrayOfPtrs, 
00136                                                                    size_t NumVectors)
00137   : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()) {
00138     TEST_FOR_EXCEPTION(NumVectors < 1 || NumVectors != Teuchos::as<size_t>(ArrayOfPtrs.size()), std::runtime_error,
00139         "Tpetra::MultiVector::MultiVector(ArrayOfPtrs): ArrayOfPtrs.size() must be strictly positive and as large as ArrayOfPtrs.");
00140     const size_t myLen = getLocalLength();
00141     if (myLen > 0) {
00142       Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00143       Teuchos::ArrayRCP<Scalar> mydata = node->template allocBuffer<Scalar>(myLen*NumVectors);
00144       MVT::initializeValues(lclMV_,myLen,NumVectors,mydata,myLen);
00145       KOKKOS_NODE_TRACE("MultiVector::MultiVector(2D)")
00146       Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::WriteOnly,myLen*NumVectors,mydata);
00147       for (size_t j = 0; j < NumVectors; ++j) {
00148 #ifdef HAVE_TPETRA_DEBUG
00149         TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ArrayOfPtrs[j].size()) != getLocalLength(), std::runtime_error,
00150           "Tpetra::MultiVector::MultiVector(ArrayOfPtrs): ArrayOfPtrs[" << j << "].size() (== " << ArrayOfPtrs[j].size() << 
00151           ") is not equal to getLocalLength() (== " << getLocalLength());
00152 #endif
00153         typename Teuchos::ArrayView<const Scalar>::iterator src = ArrayOfPtrs[j].begin();
00154         std::copy(src,src+myLen,myview);
00155         myview += myLen;
00156       }
00157       myview = Teuchos::null;
00158       mydata = Teuchos::null;
00159     }
00160     else {
00161       MVT::initializeValues(lclMV_,0,NumVectors,Teuchos::null,0);
00162     }
00163   }
00164 
00165 
00166   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00167   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
00168               Teuchos::ArrayRCP<Scalar> data, size_t LDA, Teuchos::ArrayView<const size_t> WhichVectors)
00169   : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()), whichVectors_(WhichVectors) {
00170     const size_t myLen = getLocalLength();
00171     size_t maxVector = *std::max_element(WhichVectors.begin(), WhichVectors.end());
00172     TEST_FOR_EXCEPTION(LDA < myLen, std::runtime_error,
00173         "Tpetra::MultiVector::MultiVector(data,LDA,WhichVectors): LDA must be large enough to accomodate the local entries.");
00174 #ifdef HAVE_TPETRA_DEBUG
00175     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(data.size()) < LDA * maxVector + myLen, std::runtime_error,
00176         "Tpetra::MultiVector::MultiVector(data,LDA,WhichVectors): data does not contain enough data to specify the entries in this.");
00177 #endif
00178     if (WhichVectors.size() == 1) {
00179       // shift data so that desired vector is vector 0
00180       maxVector = 0;
00181       data += LDA*WhichVectors[0];
00182       // kill whichVectors_; we are constant stride
00183       whichVectors_.clear();
00184     }
00185     if (myLen > 0) {
00186       MVT::initializeValues(lclMV_,myLen,maxVector+1,data,LDA);
00187     }
00188     else {
00189       MVT::initializeValues(lclMV_,0,WhichVectors.size(),Teuchos::null,0);
00190     }
00191   }
00192 
00193 
00194   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00195   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MultiVector(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
00196                 Teuchos::ArrayRCP<Scalar> data, size_t LDA, size_t NumVectors)
00197   : DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map), lclMV_(map->getNode()) {
00198     TEST_FOR_EXCEPTION(NumVectors < 1, std::invalid_argument,
00199         "Tpetra::MultiVector::MultiVector(data,LDA,NumVector): NumVectors must be strictly positive.");
00200     const LocalOrdinal myLen = getLocalLength();
00201 #ifdef HAVE_TPETRA_DEBUG
00202     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(data.size()) < LDA*(NumVectors-1)+myLen, std::runtime_error,
00203         "Tpetra::MultiVector::MultiVector(data,LDA,NumVector): data does not contain enough data to specify the entries in this.");
00204 #endif
00205     MVT::initializeValues(lclMV_,myLen,NumVectors,data,LDA);
00206   }
00207 
00208 
00209   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00210   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::~MultiVector() {
00211   }
00212 
00213 
00214   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00215   bool MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isConstantStride() const {
00216     return whichVectors_.empty();
00217   }
00218 
00219 
00220   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00221   size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getLocalLength() const {
00222     return this->getMap()->getNodeNumElements();
00223   }
00224 
00225 
00226   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00227   global_size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getGlobalLength() const {
00228     return this->getMap()->getGlobalNumElements();
00229   }
00230 
00231 
00232   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00233   size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getStride() const {
00234     if (isConstantStride()) {
00235       return MVT::getStride(lclMV_);
00236     }
00237     return 0;
00238   }
00239 
00240 
00241   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00242   bool MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::checkSizes(const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node> &sourceObj) 
00243   {
00244     const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A = dynamic_cast<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(sourceObj);
00245     // objects maps have already been checked. simply check the number of vectors.
00246     bool compat = (A.getNumVectors() == this->getNumVectors());
00247     return compat;
00248   }
00249 
00250 
00251   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00252   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::copyAndPermute(
00253                           const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node> & sourceObj,
00254                           size_t numSameIDs,
00255                           const Teuchos::ArrayView<const LocalOrdinal> &permuteToLIDs,
00256                           const Teuchos::ArrayView<const LocalOrdinal> &permuteFromLIDs) {
00257     using Teuchos::ArrayView;
00258     using Teuchos::ArrayRCP;
00259     const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &sourceMV = dynamic_cast<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &>(sourceObj);
00260     typename ArrayView<const LocalOrdinal>::iterator pTo, pFrom;
00261     // any other error will be caught by Teuchos
00262     TEST_FOR_EXCEPTION(permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error,
00263         "Tpetra::MultiVector::copyAndPermute(): permuteToLIDs and permuteFromLIDs must have the same size.");
00264     // one vector at a time
00265     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00266     const size_t numCols = getNumVectors();
00267     // Get a host view of the local multivector data.
00268     for (size_t j = 0; j < numCols; ++j) {
00269       // The first numImportIDs GIDs are the same between source and target,
00270       // We can just copy them
00271       ArrayRCP<const Scalar> srcptr = sourceMV.getSubArrayRCP(sourceMV.cview_,j);
00272       ArrayRCP<      Scalar> dstptr =          getSubArrayRCP(ncview_,j);
00273       std::copy(srcptr,srcptr+numSameIDs,dstptr);
00274       // next, do permutations
00275       for (pTo = permuteToLIDs.begin(), pFrom = permuteFromLIDs.begin();
00276            pTo != permuteToLIDs.end(); ++pTo, ++pFrom) {
00277         dstptr[*pTo] = srcptr[*pFrom];
00278       }
00279     }
00280   }
00281 
00282 
00283   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00284   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::packAndPrepare(
00285           const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node> & sourceObj,
00286           const Teuchos::ArrayView<const LocalOrdinal> &exportLIDs,
00287           Teuchos::Array<Scalar> &exports,
00288           const Teuchos::ArrayView<size_t> &numExportPacketsPerLID,
00289           size_t& constantNumPackets,
00290           Distributor & /* distor */ ) {
00291     const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &sourceMV = dynamic_cast<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &>(sourceObj);
00292     using Teuchos::ArrayView;
00293     using Teuchos::ArrayRCP;
00294     /* The layout in the export for MultiVectors is as follows:
00295        exports = { all of the data from row exportLIDs.front() ; 
00296                    ....
00297                    all of the data from row exportLIDs.back() }
00298       this doesn't have the best locality, but is necessary because the data for a Packet
00299       (all data associated with an LID) is required to be contiguous */
00300     TEST_FOR_EXCEPTION(Teuchos::as<int>(numExportPacketsPerLID.size()) != exportLIDs.size(), std::runtime_error,
00301         "Tpetra::MultiVector::packAndPrepare(): size of numExportPacketsPerLID buffer should be the same as exportLIDs.");
00302     const KMV &srcData = sourceMV.lclMV_;
00303     const size_t numCols = sourceMV.getNumVectors(),
00304                   stride = MVT::getStride(srcData);
00305     constantNumPackets = numCols;
00306     exports.resize(numCols*exportLIDs.size());
00307     typename ArrayView<const LocalOrdinal>::iterator idptr;
00308     typename Teuchos::Array<Scalar>::iterator expptr;
00309     expptr = exports.begin();
00310 
00311     Teuchos::RCP<Node> node = MVT::getNode(srcData);
00312     if (sourceMV.isConstantStride()) {
00313       size_t i = 0;
00314       for (idptr = exportLIDs.begin(); idptr != exportLIDs.end(); ++idptr, ++i) {
00315         for (size_t j = 0; j < numCols; ++j) {
00316           *expptr++ = sourceMV.cview_[j*stride + (*idptr)];
00317         }
00318         //we shouldn't need to set numExportPacketsPerLID[i] since we have set
00319         //constantNumPackets to a nonzero value. But we'll set it anyway, since
00320         //I'm not sure if the API will remain the way it is.
00321         numExportPacketsPerLID[i] = numCols;
00322       }
00323     }
00324     else {
00325       size_t i = 0;
00326       for (idptr = exportLIDs.begin(); idptr != exportLIDs.end(); ++idptr, ++i) {
00327         for (size_t j = 0; j < numCols; ++j) {
00328           *expptr++ = sourceMV.cview_[sourceMV.whichVectors_[j]*stride + (*idptr)];
00329         }
00330         numExportPacketsPerLID[i] = numCols;
00331       }
00332     }
00333   }
00334 
00335 
00336   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00337   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unpackAndCombine(
00338                   const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
00339                   const Teuchos::ArrayView<const Scalar> &imports,
00340                   const Teuchos::ArrayView<size_t> &numPacketsPerLID,
00341                   size_t constantNumPackets,
00342                   Distributor & /* distor */,
00343                   CombineMode CM) {
00344     using Teuchos::ArrayView;
00345     using Teuchos::ArrayRCP;
00346     /* The layout in the export for MultiVectors is as follows:
00347        imports = { all of the data from row exportLIDs.front() ; 
00348                    ....
00349                    all of the data from row exportLIDs.back() }
00350       this doesn't have the best locality, but is necessary because the data for a Packet
00351       (all data associated with an LID) is required to be contiguous */
00352 #ifdef HAVE_TPETRA_DEBUG
00353     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(imports.size()) != getNumVectors()*importLIDs.size(), std::runtime_error,
00354         "Tpetra::MultiVector::unpackAndCombine(): sizing of imports buffer should be appropriate for the amount of data to be exported.");
00355 #endif
00356     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(constantNumPackets) == 0u, std::runtime_error,
00357         "Tpetra::MultiVector::unpackAndCombine(): 'constantNumPackets' input argument should be nonzero.");
00358 
00359     const size_t myStride = MVT::getStride(lclMV_),
00360                  numVecs  = getNumVectors();
00361     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(numPacketsPerLID.size()) != Teuchos::as<size_t>(importLIDs.size()), std::runtime_error,
00362         "Tpetra::MultiVector::unpackAndCombine(): 'numPacketsPerLID' must have same length as importLIDs.");
00363     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(numVecs) != Teuchos::as<size_t>(constantNumPackets), std::runtime_error,
00364         "Tpetra::MultiVector::unpackAndCombine(): 'constantNumPackets' must equal numVecs.");
00365 
00366     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00367     if (numVecs > 0 && importLIDs.size()) {
00368       typename ArrayView<const       Scalar>::iterator impptr;
00369       typename ArrayView<const LocalOrdinal>::iterator  idptr;
00370       impptr = imports.begin();
00371       if (CM == INSERT || CM == REPLACE) {
00372         if (isConstantStride()) {
00373           for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) {
00374             for (size_t j = 0; j < numVecs; ++j) {
00375               ncview_[myStride*j + *idptr] = *impptr++;
00376             }
00377           }
00378         }
00379         else {
00380           for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) {
00381             for (size_t j = 0; j < numVecs; ++j) {
00382               ncview_[myStride*whichVectors_[j] + *idptr] = *impptr++;
00383             }
00384           }
00385         }
00386       }
00387       else if (CM == ADD) {
00388         if (isConstantStride()) {
00389           for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) {
00390             for (size_t j = 0; j < numVecs; ++j) {
00391               ncview_[myStride*j + *idptr] += *impptr++;
00392             }
00393           }
00394         }
00395         else {
00396           for (idptr = importLIDs.begin(); idptr != importLIDs.end(); ++idptr) {
00397             for (size_t j = 0; j < numVecs; ++j) {
00398               ncview_[myStride*whichVectors_[j] + *idptr] += *impptr++;
00399             }
00400           }
00401         }
00402       }
00403       else {
00404         TEST_FOR_EXCEPTION(CM != ADD && CM != REPLACE && CM != INSERT, std::invalid_argument,
00405             "Tpetra::MultiVector::unpackAndCombine(): Invalid CombineMode: " << CM);
00406       }
00407     }
00408   }
00409 
00410 
00411   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00412   inline size_t MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNumVectors() const {
00413     size_t ret;
00414     if (isConstantStride()) {
00415       ret = MVT::getNumCols(lclMV_);
00416     }
00417     else {
00418       ret = whichVectors_.size();
00419     }
00420     return ret;
00421   }
00422 
00423 
00424   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00425   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dot(
00426       const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, 
00427       const Teuchos::ArrayView<Scalar> &dots) const 
00428   {
00429     using Teuchos::ArrayRCP;
00430     using Teuchos::arcp_const_cast;
00431     const size_t myLen   = getLocalLength(),
00432                  numVecs = getNumVectors();
00433 #ifdef HAVE_TPETRA_DEBUG
00434     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error,
00435         "Tpetra::MultiVector::dot(): MultiVectors do not have compatible Maps:" << std::endl
00436         << "this->getMap(): " << std::endl << *this->getMap() 
00437         << "A.getMap(): " << std::endl << *A.getMap() << std::endl);
00438 #else
00439     TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error,
00440         "Tpetra::MultiVector::dot(): MultiVectors do not have the same local length.");
00441 #endif
00442     TEST_FOR_EXCEPTION(A.getNumVectors() != numVecs, std::runtime_error,
00443         "Tpetra::MultiVector::dot(): MultiVectors must have the same number of vectors.");
00444     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(dots.size()) != numVecs, std::runtime_error,
00445         "Tpetra::MultiVector::dot(): dots.size() must be as large as the number of vectors in *this and A.");
00446     if (isConstantStride() && A.isConstantStride()) {
00447       MVT::Dot(lclMV_,A.lclMV_,dots);
00448     }
00449     else {
00450       KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_));
00451       ArrayRCP<Scalar> vptr = arcp_const_cast<Scalar>(MVT::getValues(lclMV_)),
00452                       avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_));
00453       for (size_t j=0; j < numVecs; ++j) {
00454         ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j),
00455                         avj = A.getSubArrayRCP(avptr,j);
00456         MVT::initializeValues(a,myLen, 1, avj, myLen);
00457         MVT::initializeValues(v,myLen, 1,  vj, myLen);
00458         dots[j] = MVT::Dot((const KMV&)v,(const KMV &)a);
00459       }
00460     }
00461     if (this->isDistributed()) {
00462       Teuchos::Array<Scalar> ldots(dots);
00463       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),ldots.getRawPtr(),dots.getRawPtr());
00464     }
00465   }
00466 
00467 
00468   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00469   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm2(
00470                   const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
00471     using Teuchos::ScalarTraits;
00472     using Teuchos::ArrayView;
00473     typedef typename ScalarTraits<Scalar>::magnitudeType Mag;
00474     const size_t numVecs = this->getNumVectors();
00475     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error,
00476         "Tpetra::MultiVector::norm2(norms): norms.size() must be as large as the number of vectors in *this.");
00477     if (isConstantStride()) {
00478       MVT::Norm2Squared(lclMV_,norms);
00479     }
00480     else {
00481       KMV v(MVT::getNode(lclMV_));
00482       Teuchos::ArrayRCP<Scalar> vi;
00483       for (size_t i=0; i < numVecs; ++i) {
00484         vi = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[i]) );
00485         MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vi, MVT::getStride(lclMV_));
00486         norms[i] = MVT::Norm2Squared(v);
00487       }
00488     }
00489     if (this->isDistributed()) {
00490       Teuchos::Array<Mag> lnorms(norms);
00491       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr());
00492     }
00493     for (typename ArrayView<Mag>::iterator n = norms.begin(); n != norms.begin()+numVecs; ++n) {
00494       (*n) = ScalarTraits<Mag>::squareroot(*n);
00495     }
00496   }
00497 
00498 
00499   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00500   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normWeighted(
00501           const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &weights,
00502           const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
00503     using Teuchos::ScalarTraits;
00504     using Teuchos::ArrayView;
00505     using Teuchos::ArrayRCP;
00506     using Teuchos::arcp_const_cast;
00507     typedef ScalarTraits<Scalar> SCT;
00508     typedef typename SCT::magnitudeType Mag;
00509     const Mag OneOverN = ScalarTraits<Mag>::one() / Teuchos::as<Mag>(getGlobalLength());
00510     bool OneW = false;
00511     const size_t numVecs = this->getNumVectors();
00512     if (weights.getNumVectors() == 1) {
00513       OneW = true;
00514     }
00515     else {
00516       TEST_FOR_EXCEPTION(weights.getNumVectors() != numVecs, std::runtime_error,
00517           "Tpetra::MultiVector::normWeighted(): MultiVector of weights must contain either one vector or the same number of vectors as this.");
00518     }
00519 #ifdef HAVE_TPETRA_DEBUG
00520     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*weights.getMap()), std::runtime_error,
00521         "Tpetra::MultiVector::normWeighted(): MultiVectors do not have compatible Maps:" << std::endl
00522         << "this->getMap(): " << std::endl << *this->getMap() 
00523         << "weights.getMap(): " << std::endl << *weights.getMap() << std::endl);
00524 #else
00525     TEST_FOR_EXCEPTION( getLocalLength() != weights.getLocalLength(), std::runtime_error,
00526         "Tpetra::MultiVector::normWeighted(): MultiVectors do not have the same local length.");
00527 #endif
00528     const size_t myLen = getLocalLength();
00529     // 
00530     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error,
00531         "Tpetra::MultiVector::normWeighted(): norms.size() must be as large as the number of vectors in *this.");
00532     if (isConstantStride() && weights.isConstantStride()) {
00533       MVT::WeightedNorm(lclMV_,weights.lclMV_,norms);
00534     }
00535     else {
00536       KMV v(MVT::getNode(lclMV_)), w(MVT::getNode(lclMV_));
00537       ArrayRCP<Scalar> vptr = arcp_const_cast<Scalar>(MVT::getValues(lclMV_)),
00538                        wptr = arcp_const_cast<Scalar>(MVT::getValues(weights.lclMV_));
00539       ArrayRCP<Scalar> wj = wptr.persistingView(0,myLen);
00540       MVT::initializeValues(w,myLen, 1, wj, myLen);
00541       for (size_t j=0; j < numVecs; ++j) {
00542         ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j);
00543         MVT::initializeValues(v,myLen, 1,  vj, myLen);
00544         if (!OneW) {
00545           wj = weights.getSubArrayRCP(wptr,j);
00546           MVT::initializeValues(w,myLen, 1, wj, myLen);
00547         }
00548         norms[j] = MVT::WeightedNorm((const KMV&)v,(const KMV &)w);
00549       }
00550     }
00551     if (this->isDistributed()) {
00552       Teuchos::Array<Mag> lnorms(norms);
00553       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr());
00554     }
00555     for (typename ArrayView<Mag>::iterator n = norms.begin(); n != norms.begin()+numVecs; ++n) {
00556       *n = ScalarTraits<Mag>::squareroot(*n * OneOverN);
00557     }
00558   }
00559 
00560 
00561   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00562   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::norm1(
00563                   const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
00564     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00565     const size_t numVecs = this->getNumVectors();
00566     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error,
00567         "Tpetra::MultiVector::norm1(norms): norms.size() must be as large as the number of vectors in *this.");
00568     if (isConstantStride()) {
00569       MVT::Norm1(lclMV_,norms);
00570     }
00571     else {
00572       KMV v(MVT::getNode(lclMV_));
00573       Teuchos::ArrayRCP<Scalar> vj;
00574       for (size_t j=0; j < numVecs; ++j) {
00575         vj = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[j]) );
00576         MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_));
00577         norms[j] = MVT::Norm1((const KMV&)v);
00578       }
00579     }
00580     if (this->isDistributed()) {
00581       Teuchos::Array<Mag> lnorms(norms);
00582       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr());
00583     }
00584   }
00585 
00586 
00587   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00588   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::normInf(
00589         const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
00590     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Mag;
00591     const size_t numVecs = this->getNumVectors();
00592     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(norms.size()) != numVecs, std::runtime_error,
00593         "Tpetra::MultiVector::normInf(norms): norms.size() must be as large as the number of vectors in *this.");
00594     if (isConstantStride()) {
00595       MVT::NormInf(lclMV_,norms);
00596     }
00597     else {
00598       KMV v(MVT::getNode(lclMV_));
00599       Teuchos::ArrayRCP<Scalar> vj;
00600       for (size_t j=0; j < numVecs; ++j) {
00601         vj = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[j]) );
00602         MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_));
00603         norms[j] = MVT::NormInf((const KMV&)v);
00604       }
00605     }
00606     if (this->isDistributed()) {
00607       Teuchos::Array<Mag> lnorms(norms);
00608       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_MAX,Teuchos::as<int>(numVecs),lnorms.getRawPtr(),norms.getRawPtr());
00609     }
00610   }
00611 
00612 
00613   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00614   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::meanValue(const Teuchos::ArrayView<Scalar> &means) const
00615   {
00616     typedef Teuchos::ScalarTraits<Scalar> SCT;
00617     using Teuchos::ArrayRCP;
00618     using Teuchos::arcp_const_cast;
00619     const size_t numVecs = getNumVectors();
00620     const size_t myLen   = getLocalLength();
00621     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(means.size()) != numVecs, std::runtime_error,
00622         "Tpetra::MultiVector::meanValue(): means.size() must be as large as the number of vectors in *this.");
00623     // compute local components of the means
00624     // sum these across all nodes
00625     if (isConstantStride()) {
00626       MVT::Sum(lclMV_,means);
00627     }
00628     else {
00629       KMV v(MVT::getNode(lclMV_));
00630       ArrayRCP<Scalar> vptr = arcp_const_cast<Scalar>(MVT::getValues(lclMV_));
00631       for (size_t j=0; j < numVecs; ++j) {
00632         ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j);
00633         MVT::initializeValues(v,myLen, 1,  vj, myLen);
00634         means[j] = MVT::Sum((const KMV &)v);
00635       }
00636     }
00637     if (this->isDistributed()) {
00638       Teuchos::Array<Scalar> lmeans(means);
00639       // only combine if we are a distributed MV
00640       Teuchos::reduceAll(*this->getMap()->getComm(),Teuchos::REDUCE_SUM,Teuchos::as<int>(numVecs),lmeans.getRawPtr(),means.getRawPtr());
00641     }
00642     const Scalar OneOverN = Teuchos::ScalarTraits<Scalar>::one() / Teuchos::as<Scalar>(getGlobalLength());
00643     for (typename Teuchos::ArrayView<Scalar>::iterator i = means.begin(); i != means.begin()+numVecs; ++i) {
00644       (*i) = (*i)*OneOverN;
00645     }
00646   }
00647 
00648 
00649   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00650   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::randomize() {
00651     if (isConstantStride()) {
00652       MVT::Random(lclMV_);
00653     }
00654     else {
00655       const size_t numVecs = this->getNumVectors();
00656       KMV v(MVT::getNode(lclMV_));
00657       Teuchos::ArrayRCP<Scalar> vj;
00658       for (size_t j=0; j < numVecs; ++j) {
00659         vj = MVT::getValuesNonConst(lclMV_,whichVectors_[j]);
00660         MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_));
00661         MVT::Random(v);
00662       }
00663     }
00664   }
00665 
00666 
00667   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00668   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::putScalar(const Scalar &alpha) {
00669     const size_t numVecs = getNumVectors();
00670     if (isConstantStride()) {
00671       MVT::Init(lclMV_,alpha);
00672     }
00673     else {
00674       KMV v(MVT::getNode(lclMV_));
00675       Teuchos::ArrayRCP<Scalar> vj;
00676       for (size_t j=0; j < numVecs; ++j) {
00677         vj = MVT::getValuesNonConst(lclMV_,whichVectors_[j]);
00678         MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_));
00679         MVT::Init(v,alpha);
00680       }
00681     }
00682   }
00683 
00684 
00685   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00686   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::scale(const Scalar &alpha) {
00687     // NOTE: can't substitute putScalar(0.0) for scale(0.0), because 
00688     //       the former will overwrite NaNs present in the MultiVector, while the 
00689     //       semantics of this call require multiplying them by 0, which IEEE requires to be NaN
00690     const size_t numVecs = getNumVectors();
00691     if (alpha == Teuchos::ScalarTraits<Scalar>::one()) {
00692       // do nothing
00693     }
00694     else if (isConstantStride()) {
00695       MVT::Scale(lclMV_,alpha);
00696     }
00697     else {
00698       KMV v(MVT::getNode(lclMV_));
00699       Teuchos::ArrayRCP<Scalar> vj;
00700       for (size_t j=0; j < numVecs; ++j) {
00701         vj = Teuchos::arcp_const_cast<Scalar>( MVT::getValues(lclMV_,whichVectors_[j]) );
00702         MVT::initializeValues(v,MVT::getNumRows(lclMV_), 1, vj, MVT::getStride(lclMV_));
00703         MVT::Scale(v,alpha);
00704       }
00705     }
00706   }
00707 
00708 
00709   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00710   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::scale(Teuchos::ArrayView<const Scalar> alphas)
00711   {
00712     using Teuchos::ArrayRCP;
00713     const size_t numVecs = this->getNumVectors();
00714     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(alphas.size()) != numVecs, std::runtime_error,
00715         "Tpetra::MultiVector::scale(alphas): alphas.size() must be as large as the number of vectors in *this.");
00716     KMV vec(MVT::getNode(lclMV_));
00717     const size_t myLen = MVT::getNumRows(lclMV_);
00718     if (myLen == 0) return;
00719     ArrayRCP<Scalar> mybuf = MVT::getValuesNonConst(lclMV_);
00720     for (size_t j = 0; j < numVecs; ++j) {
00721       if (alphas[j] == Teuchos::ScalarTraits<Scalar>::one()) {
00722         // do nothing: NaN * 1.0 == NaN, Number*1.0 == Number
00723       }
00724       else {
00725         ArrayRCP<Scalar> mybufj = getSubArrayRCP(mybuf,j);
00726         MVT::initializeValues(vec,myLen,1,mybufj,myLen);
00727         MVT::Scale(vec,alphas[j]);
00728       }
00729     }
00730   }
00731 
00732 
00733   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00734   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::scale(const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A) {
00735     using Teuchos::ArrayRCP;
00736     using Teuchos::arcp_const_cast;
00737     const size_t numVecs = getNumVectors(),
00738                  myLen   = getLocalLength();
00739 #ifdef HAVE_TPETRA_DEBUG
00740     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error,
00741         "Tpetra::MultiVector::scale(alpha,A): MultiVectors do not have compatible Maps:" << std::endl
00742         << "this->getMap(): " << std::endl << *this->getMap() 
00743         << "A.getMap(): " << std::endl << *A.getMap() << std::endl);
00744 #else
00745     TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error,
00746         "Tpetra::MultiVector::scale(alpha,A): MultiVectors do not have the same local length.");
00747 #endif
00748     TEST_FOR_EXCEPTION(A.getNumVectors() != numVecs, std::runtime_error,
00749         "Tpetra::MultiVector::scale(alpha,A): MultiVectors must have the same number of vectors.");
00750     if (isConstantStride() && A.isConstantStride()) {
00751       // set me == alpha*A
00752       MVT::Scale(lclMV_,alpha,(const KMV&)A.lclMV_);
00753     }
00754     else {
00755       KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_));
00756       ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_),
00757                       avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_));
00758       for (size_t j=0; j < numVecs; ++j) {
00759         ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j),
00760                         avj = A.getSubArrayRCP(avptr,j);
00761         MVT::initializeValues(a,myLen, 1, avj, myLen);
00762         MVT::initializeValues(v,myLen, 1,  vj, myLen);
00763         MVT::Scale(v,alpha,(const KMV &)a);
00764       }
00765     }
00766   }
00767 
00768 
00769   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00770   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::reciprocal(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A) {
00771     using Teuchos::ArrayRCP;
00772     using Teuchos::arcp_const_cast;
00773 #ifdef HAVE_TPETRA_DEBUG
00774     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error,
00775         "Tpetra::MultiVector::reciprocal(): MultiVectors do not have compatible Maps:" << std::endl
00776         << "this->getMap(): " << std::endl << *this->getMap() 
00777         << "A.getMap(): " << std::endl << *A.getMap() << std::endl);
00778 #else
00779     TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error,
00780         "Tpetra::MultiVector::reciprocal(): MultiVectors do not have the same local length.");
00781 #endif
00782     TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors(), std::runtime_error,
00783         "Tpetra::MultiVector::reciprocal(): MultiVectors must have the same number of vectors.");
00784     const size_t numVecs = getNumVectors();
00785     const size_t myLen = getLocalLength();
00786     try {
00787       if (isConstantStride() && A.isConstantStride()) {
00788         MVT::Recip(lclMV_,(const KMV&)A.lclMV_);
00789       }
00790       else {
00791         KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_));
00792         ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_),
00793                         avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_));
00794         for (size_t j=0; j < numVecs; ++j) {
00795           ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j),
00796                           avj = A.getSubArrayRCP(avptr,j);
00797           MVT::initializeValues(a,myLen, 1, avj, myLen);
00798           MVT::initializeValues(v,myLen, 1,  vj, myLen);
00799           MVT::Recip(v,(const KMV &)a);
00800         }
00801       }
00802     }
00803     catch (std::runtime_error &e) {
00804       TEST_FOR_EXCEPTION(true,std::runtime_error,
00805           "Tpetra::MultiVector::reciprocal(): caught exception from Kokkos:" << std::endl
00806           << e.what() << std::endl);
00807     }
00808   }
00809 
00810 
00811   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00812   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::abs(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A) {
00813     using Teuchos::ArrayRCP;
00814     using Teuchos::arcp_const_cast;
00815 #ifdef HAVE_TPETRA_DEBUG
00816     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error,
00817         "Tpetra::MultiVector::abs(): MultiVectors do not have compatible Maps:" << std::endl
00818         << "this->getMap(): " << std::endl << *this->getMap() 
00819         << "A.getMap(): " << std::endl << *A.getMap() << std::endl);
00820 #else
00821     TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error,
00822         "Tpetra::MultiVector::abs(): MultiVectors do not have the same local length.");
00823 #endif
00824     TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors(), std::runtime_error,
00825         "Tpetra::MultiVector::abs(): MultiVectors must have the same number of vectors.");
00826     const size_t myLen = getLocalLength();
00827     const size_t numVecs = getNumVectors();
00828     if (isConstantStride() && A.isConstantStride()) {
00829       MVT::Abs(lclMV_,(const KMV&)A.lclMV_);
00830     }
00831     else {
00832       KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_));
00833       ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_),
00834                       avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_));
00835       for (size_t j=0; j < numVecs; ++j) {
00836         ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j),
00837                         avj = A.getSubArrayRCP(avptr,j);
00838         MVT::initializeValues(a,myLen, 1, avj, myLen);
00839         MVT::initializeValues(v,myLen, 1,  vj, myLen);
00840         MVT::Abs(v,(const KMV &)a);
00841       }
00842     }
00843   }
00844 
00845 
00846   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00847   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::update(
00848                       const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, 
00849                       const Scalar &beta) {
00850     // this = beta*this + alpha*A
00851     // must support case where &this == &A
00852     // can't short circuit on alpha==0.0 or beta==0.0, because 0.0*NaN != 0.0
00853     using Teuchos::ArrayRCP;
00854     using Teuchos::arcp_const_cast;
00855 #ifdef HAVE_TPETRA_DEBUG
00856     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()), std::runtime_error,
00857         "Tpetra::MultiVector::update(): MultiVectors do not have compatible Maps:" << std::endl
00858         << "this->getMap(): " << std::endl << this->getMap() 
00859         << "A.getMap(): " << std::endl << *A.getMap() << std::endl);
00860 #else
00861     TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength(), std::runtime_error,
00862         "Tpetra::MultiVector::update(): MultiVectors do not have the same local length.");
00863 #endif
00864     const size_t myLen = getLocalLength();
00865     const size_t numVecs = getNumVectors();
00866     TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors(), std::runtime_error,
00867         "Tpetra::MultiVector::update(): MultiVectors must have the same number of vectors.");
00868     if (isConstantStride() && A.isConstantStride()) {
00869       MVT::GESUM(lclMV_,alpha,(const KMV&)A.lclMV_,beta);
00870     }
00871     else {
00872       KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_));
00873       ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_),
00874                       avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_));
00875       for (size_t j=0; j < numVecs; ++j) {
00876         ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j),
00877                         avj = A.getSubArrayRCP(avptr,j);
00878         MVT::initializeValues(a,myLen, 1, avj, myLen);
00879         MVT::initializeValues(v,myLen, 1,  vj, myLen);
00880         MVT::GESUM(v,alpha,(const KMV &)a,beta);
00881       }
00882     }
00883   }
00884 
00885 
00886   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00887   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::update(
00888                       const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, 
00889                       const Scalar &beta,  const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B, 
00890                       const Scalar &gamma) {
00891     using Teuchos::ArrayRCP;
00892     using Teuchos::arcp_const_cast;
00893     // this = alpha*A + beta*B + gamma*this
00894     // must support case where &this == &A or &this == &B
00895     // can't short circuit on alpha==0.0 or beta==0.0 or gamma==0.0, because 0.0*NaN != 0.0
00896 #ifdef HAVE_TPETRA_DEBUG
00897     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()) || !this->getMap()->isCompatible(*B.getMap()),
00898         std::runtime_error,
00899         "Tpetra::MultiVector::update(): MultiVectors do not have compatible Maps:" << std::endl
00900         << "this->getMap(): " << std::endl << *this->getMap() 
00901         << "A.getMap(): " << std::endl << *A.getMap() << std::endl
00902         << "B.getMap(): " << std::endl << *B.getMap() << std::endl);
00903 #else
00904     TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength() || getLocalLength() != B.getLocalLength(), std::runtime_error,
00905         "Tpetra::MultiVector::update(): MultiVectors do not have the same local length.");
00906 #endif
00907     TEST_FOR_EXCEPTION(A.getNumVectors() != this->getNumVectors() || B.getNumVectors() != this->getNumVectors(), std::runtime_error,
00908         "Tpetra::MultiVector::update(): MultiVectors must have the same number of vectors.");
00909     const size_t myLen = getLocalLength();
00910     const size_t numVecs = getNumVectors();
00911     if (isConstantStride() && A.isConstantStride() && B.isConstantStride()) {
00912       MVT::GESUM(lclMV_,alpha,(const KMV&)A.lclMV_,beta,(const KMV&)B.lclMV_,gamma);
00913     }
00914     else {
00915       KMV v(MVT::getNode(lclMV_)), a(MVT::getNode(lclMV_)), b(MVT::getNode(lclMV_));
00916       ArrayRCP<Scalar> vptr = MVT::getValuesNonConst(lclMV_),
00917                       avptr = arcp_const_cast<Scalar>(MVT::getValues(A.lclMV_)),
00918                       bvptr = arcp_const_cast<Scalar>(MVT::getValues(B.lclMV_));
00919       for (size_t j=0; j < numVecs; ++j) {
00920         ArrayRCP<Scalar> vj =   getSubArrayRCP( vptr,j),
00921                         avj = A.getSubArrayRCP(avptr,j),
00922                         bvj = B.getSubArrayRCP(bvptr,j);
00923         MVT::initializeValues(b,myLen, 1, bvj, myLen);
00924         MVT::initializeValues(a,myLen, 1, avj, myLen);
00925         MVT::initializeValues(v,myLen, 1,  vj, myLen);
00926         MVT::GESUM(v,alpha,(const KMV&)a,beta,(const KMV&)b,gamma);
00927       }
00928     }
00929   }
00930 
00931 
00932   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00933   Teuchos::ArrayRCP<const Scalar>
00934   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getData(size_t j) const {
00935     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00936     KOKKOS_NODE_TRACE("MultiVector::getData()")
00937     return node->template viewBuffer<Scalar>( getLocalLength(), getSubArrayRCP(MVT::getValues(lclMV_),j) );
00938   }
00939 
00940 
00941   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00942   Teuchos::ArrayRCP<Scalar>
00943   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getDataNonConst(size_t j) {
00944     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00945     KOKKOS_NODE_TRACE("MultiVector::getDataNonConst()")
00946     return node->template viewBufferNonConst<Scalar>(Kokkos::ReadWrite, getLocalLength(), getSubArrayRCP(MVT::getValuesNonConst(lclMV_),j) );
00947   }
00948 
00949 
00950   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00951   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & 
00952   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::operator=(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source) {
00953     // Check for special case of this=Source, in which case we do nothing
00954     if (this != &source) {
00955 #ifdef HAVE_TPETRA_DEBUG
00956       TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*source.getMap()), std::runtime_error,
00957           "Tpetra::MultiVector::operator=(): MultiVectors do not have compatible Maps:" << std::endl
00958           << "this->getMap(): " << std::endl << *this->getMap() 
00959           << "source.getMap(): " << std::endl << *source.getMap() << std::endl);
00960 #else
00961       TEST_FOR_EXCEPTION( getLocalLength() != source.getLocalLength(), std::runtime_error,
00962           "Tpetra::MultiVector::operator=(): MultiVectors do not have the same local length.");
00963 #endif
00964       TEST_FOR_EXCEPTION(source.getNumVectors() != getNumVectors(), std::runtime_error,
00965           "Tpetra::MultiVector::operator=(): MultiVectors must have the same number of vectors.");
00966       Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00967       const size_t numVecs = getNumVectors();
00968       if (isConstantStride() && source.isConstantStride() && getLocalLength()==getStride() && source.getLocalLength()==source.getStride()) {
00969         // we're both packed, we can copy in one call
00970         KOKKOS_NODE_TRACE("MultiVector::operator=()")
00971         node->template copyBuffers<Scalar>(getLocalLength()*numVecs, MVT::getValues(source.lclMV_), MVT::getValuesNonConst(lclMV_) );
00972       }
00973       else {
00974         for (size_t j=0; j < numVecs; ++j) {
00975           KOKKOS_NODE_TRACE("MultiVector::operator=()")
00976           node->template copyBuffers<Scalar>(getLocalLength(), source.getSubArrayRCP(MVT::getValues(source.lclMV_),j),  
00977                                                                       getSubArrayRCP(MVT::getValuesNonConst(lclMV_),j) );
00978         }
00979       }
00980     }
00981     return(*this);
00982   }
00983 
00984 
00985   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00986   Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
00987   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subCopy(const Teuchos::ArrayView<const size_t> &cols) const {
00988     TEST_FOR_EXCEPTION(cols.size() < 1, std::runtime_error,
00989         "Tpetra::MultiVector::subCopy(cols): cols must contain at least one column.");
00990     size_t numCopyVecs = cols.size();
00991     const bool zeroData = false;
00992     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
00993     Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mv; 
00994     // mv is allocated with constant stride
00995     mv = Teuchos::rcp( new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),numCopyVecs,zeroData) );
00996     // copy data from *this into mv
00997     for (size_t j=0; j<numCopyVecs; ++j) {
00998       KOKKOS_NODE_TRACE("MultiVector::subCopy()")
00999       node->template copyBuffers<Scalar>( getLocalLength(), getSubArrayRCP(MVT::getValues(lclMV_), cols[j]), 
01000                                                             MVT::getValuesNonConst(mv->lclMV_,j) );
01001     }
01002     return mv;
01003   }
01004 
01005 
01006   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01007   Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01008   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subCopy(const Teuchos::Range1D &colRng) const {
01009     TEST_FOR_EXCEPTION(colRng.size() == 0, std::runtime_error,
01010         "Tpetra::MultiVector::subCopy(Range1D): range must include at least one vector.");
01011     size_t numCopyVecs = colRng.size();
01012     const bool zeroData = false;
01013     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01014     Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mv; 
01015     // mv is allocated with constant stride
01016     mv = Teuchos::rcp( new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),numCopyVecs,zeroData) );
01017     // copy data from *this into mv
01018     for (size_t js=colRng.lbound(), jd=0; jd<numCopyVecs; ++jd, ++js) {
01019       KOKKOS_NODE_TRACE("MultiVector::subCopy()")
01020       node->template copyBuffers<Scalar>( getLocalLength(), getSubArrayRCP(MVT::getValues(lclMV_), js),
01021                                                             MVT::getValuesNonConst(mv->lclMV_,jd) );
01022     }
01023     return mv;
01024   }
01025 
01026 
01027   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01028   Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01029   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::offsetView(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &subMap, size_t offset) const 
01030   {
01031     typedef const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> CMV;
01032     TEST_FOR_EXCEPTION( subMap->getNodeNumElements() + offset > this->getLocalLength(), std::runtime_error,
01033         "Tpetra::MultiVector::offsetView(subMap,offset): sizes are not sane.\noffset == " << offset << "\nsubMap: " << subMap->description() << "\nthis->rowMap: " << this->getMap()->description());
01034     const size_t numVecs = this->getNumVectors(),
01035                 myStride = MVT::getStride(lclMV_),
01036                   newLen = subMap->getNodeNumElements();
01037     using Teuchos::ArrayRCP;
01038     // this is const, so the lclMV_ is const, so that we can only get const buffers
01039     // we will cast away the const; this is okay, because 
01040     //   a) the constructor doesn't modify the data, and 
01041     //   b) we are encapsulating in a const MV before returning
01042     ArrayRCP<const Scalar> cbuf = MVT::getValues(lclMV_);
01043     ArrayRCP<Scalar>      ncbuf = Teuchos::arcp_const_cast<Scalar>(cbuf);
01044     Teuchos::RCP<CMV> constViewMV;
01045     if (isConstantStride()) {
01046       // view goes from first entry of first vector to last entry of last vector
01047       ArrayRCP<Scalar> subdata = ncbuf.persistingView( offset, myStride * (numVecs-1) + newLen );
01048       constViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,numVecs) );
01049     }
01050     else {
01051       // use same which index, but with an offset start pointer
01052       size_t maxSubVecIndex = *std::max_element(whichVectors_.begin(), whichVectors_.end());
01053       ArrayRCP<Scalar> subdata = ncbuf.persistingView( offset, myStride * maxSubVecIndex + newLen );
01054       constViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,whichVectors_) );
01055     }
01056     return constViewMV;      
01057   }
01058 
01059 
01060   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01061   Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01062   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::offsetViewNonConst(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &subMap, size_t offset)
01063   {
01064     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
01065     TEST_FOR_EXCEPTION( subMap->getNodeNumElements() + offset > this->getLocalLength(), std::runtime_error,
01066         "Tpetra::MultiVector::offsetView(subMap,offset): sizes are not sane.\noffset == " << offset << "\nsubMap: " << subMap->description() << "\nthis->rowMap: " << this->getMap()->description());
01067     const size_t numVecs = this->getNumVectors(),
01068                 myStride = MVT::getStride(lclMV_),
01069                   newLen = subMap->getNodeNumElements();
01070     using Teuchos::ArrayRCP;
01071     // this is const, so the lclMV_ is const, so that we can only get const buffers
01072     // we will cast away the const; this is okay, because 
01073     //   a) the constructor doesn't modify the data, and 
01074     //   b) we are encapsulating in a const MV before returning
01075     ArrayRCP<Scalar> buf = MVT::getValuesNonConst(lclMV_);
01076     Teuchos::RCP<MV> subViewMV;
01077     if (isConstantStride()) {
01078       // view goes from first entry of first vector to last entry of last vector
01079       ArrayRCP<Scalar> subdata = buf.persistingView( offset, myStride * (numVecs-1) + newLen );
01080       subViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,numVecs) );
01081     }
01082     else {
01083       // use same which index, but with an offset start pointer
01084       size_t maxSubVecIndex = *std::max_element(whichVectors_.begin(), whichVectors_.end());
01085       ArrayRCP<Scalar> subdata = buf.persistingView( offset, myStride * maxSubVecIndex + newLen );
01086       subViewMV = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(subMap,subdata,myStride,whichVectors_) );
01087     }
01088     return subViewMV;      
01089   }
01090 
01091 
01092   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01093   Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01094   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subView(const Teuchos::ArrayView<const size_t> &cols) const {
01095     using Teuchos::ArrayRCP;
01096     typedef const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> CMV;
01097     TEST_FOR_EXCEPTION(cols.size() == 0, std::runtime_error,
01098         "Tpetra::MultiVector::subView(ArrayView): range must include at least one vector.");
01099     // this is const, so the lclMV_ is const, so that we can only get const buffers
01100     // we will cast away the const; this is okay, because 
01101     //   a) the constructor doesn't modify the data, and 
01102     //   b) we are encapsulating in a const MV before returning
01103     const size_t myStride = MVT::getStride(lclMV_),
01104                  myLen    = MVT::getNumRows(lclMV_),
01105               numViewCols = cols.size();
01106     // use the smallest view possible of the buffer: from the first element of the minInd vector to the last element of the maxInd vector
01107     // this minimizes overlap between views, and keeps view of the minimum amount necessary, in order to allow node to achieve maximum efficiency.
01108     // adjust the indices appropriately; shift so that smallest index is 0
01109     ArrayRCP<const Scalar> cbuf = MVT::getValues(lclMV_);
01110     ArrayRCP<Scalar>      ncbuf = Teuchos::arcp_const_cast<Scalar>(cbuf);
01111     Teuchos::Array<size_t> newCols(numViewCols);
01112     size_t minInd = Teuchos::OrdinalTraits<size_t>::max(),
01113            maxInd = Teuchos::OrdinalTraits<size_t>::zero();
01114     if (isConstantStride()) {
01115       for (size_t j=0; j < numViewCols; ++j) {
01116         newCols[j] = cols[j];
01117         if (newCols[j] < minInd) minInd = newCols[j];
01118         if (maxInd < newCols[j]) maxInd = newCols[j];
01119       }
01120     }
01121     else {
01122       for (size_t j=0; j < numViewCols; ++j) {
01123         newCols[j] = whichVectors_[cols[j]];
01124         if (newCols[j] < minInd) minInd = newCols[j];
01125         if (maxInd < newCols[j]) maxInd = newCols[j];
01126       }
01127     }
01128     ArrayRCP<Scalar> minbuf = ncbuf.persistingView(minInd * myStride, myStride * (maxInd - minInd) + myLen);
01129     for (size_t j=0; j < numViewCols; ++j) {
01130       newCols[j] -= minInd;
01131     }
01132     Teuchos::RCP<CMV> constViewMV;
01133     constViewMV = Teuchos::rcp<CMV>(
01134                     new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(), minbuf, myStride, newCols())
01135                   );
01136     return constViewMV;      
01137   }
01138 
01139 
01140   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01141   Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01142   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subView(const Teuchos::Range1D &colRng) const {
01143     TEST_FOR_EXCEPTION(colRng.size() == 0, std::runtime_error,
01144         "Tpetra::MultiVector::subView(Range1D): range must include at least one vector.");
01145     size_t numViewVecs = colRng.size();
01146     using Teuchos::ArrayRCP;
01147     // this is const, so the lclMV_ is const, so that we can only get const buffers
01148     // we will cast away the const; this is okay, because 
01149     //   a) the constructor doesn't modify the data, and 
01150     //   b) we are encapsulating in a const MV before returning
01151     ArrayRCP<const Scalar> cbuf = MVT::getValues(lclMV_);
01152     ArrayRCP<Scalar>      ncbuf = Teuchos::arcp_const_cast<Scalar>(cbuf);
01153     // resulting MultiVector is constant stride only if *this is 
01154     if (isConstantStride()) {
01155       // view goes from first entry of first vector to last entry of last vector
01156       ArrayRCP<Scalar> subdata = ncbuf.persistingView( MVT::getStride(lclMV_) * colRng.lbound(),
01157                                                        MVT::getStride(lclMV_) * (numViewVecs-1) + getLocalLength() );
01158       return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),
01159                                                                                   subdata,MVT::getStride(lclMV_),numViewVecs) );
01160     }
01161     // otherwise, use a subset of this whichVectors_ to construct new multivector
01162     Teuchos::Array<size_t> whchvecs( whichVectors_.begin()+colRng.lbound(), whichVectors_.begin()+colRng.ubound()+1 );
01163     return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),  
01164                                                                                 ncbuf,MVT::getStride(lclMV_),whchvecs) );
01165   }
01166 
01167 
01168   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01169   Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01170   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subViewNonConst(const Teuchos::ArrayView<const size_t> &cols) {
01171     TEST_FOR_EXCEPTION(cols.size() == 0, std::runtime_error,
01172         "Tpetra::MultiVector::subViewNonConst(ArrayView): range must include at least one vector.");
01173     if (isConstantStride()) {
01174       return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),
01175                                                                                   MVT::getValuesNonConst(lclMV_),MVT::getStride(lclMV_),
01176                                                                                   cols) );
01177     }
01178     // else, lookup current whichVectors_ using cols
01179     Teuchos::Array<size_t> newcols(cols.size());
01180     for (size_t j=0; j < Teuchos::as<size_t>(cols.size()); ++j) {
01181       newcols[j] = whichVectors_[cols[j]];
01182     }
01183     return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),
01184                                                                                 MVT::getValuesNonConst(lclMV_),MVT::getStride(lclMV_),
01185                                                                                 newcols()) );
01186   }
01187 
01188 
01189   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01190   Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01191   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::subViewNonConst(const Teuchos::Range1D &colRng) {
01192     TEST_FOR_EXCEPTION(colRng.size() == 0, std::runtime_error,
01193         "Tpetra::MultiVector::subViewNonConst(Range1D): range must include at least one vector.");
01194     size_t numViewVecs = colRng.size();
01195     using Teuchos::ArrayRCP;
01196     // resulting MultiVector is constant stride only if *this is 
01197     if (isConstantStride()) {
01198       // view goes from first entry of first vector to last entry of last vector
01199       const size_t stride = MVT::getStride(lclMV_);
01200       ArrayRCP<Scalar> data = MVT::getValuesNonConst(lclMV_);
01201       ArrayRCP<Scalar> subdata = data.persistingView( stride * colRng.lbound(),
01202                                                       stride * (numViewVecs-1) + getLocalLength() );
01203       return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),
01204                                                                                   subdata,stride,numViewVecs) );
01205     }
01206     // otherwise, use a subset of this whichVectors_ to construct new multivector
01207     Teuchos::Array<size_t> whchvecs( whichVectors_.begin()+colRng.lbound(), whichVectors_.begin()+colRng.ubound()+1 );
01208     const size_t stride = MVT::getStride(lclMV_);
01209     ArrayRCP<Scalar> data = MVT::getValuesNonConst(lclMV_);
01210     return Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),  
01211                                                                                 data,stride,whchvecs) );
01212   }
01213 
01214 
01215   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01216   Teuchos::RCP<const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01217   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getVector(size_t j) const {
01218 #ifdef HAVE_TPETRA_DEBUG
01219     TEST_FOR_EXCEPTION(j < 0 || j >= this->getNumVectors(), std::runtime_error,
01220         "Tpetra::MultiVector::getVector(j): index j (== " << j << ") exceeds valid column range for this multivector.");
01221 #endif
01222     // this is const, so lclMV_ is const, so we get const buff
01223     // it is safe to cast away the const because we will wrap it in a const Vector below
01224     Teuchos::ArrayRCP<Scalar> ncbuff;
01225     if (getLocalLength() > 0) {
01226       Teuchos::ArrayRCP<const Scalar> cbuff = getSubArrayRCP(MVT::getValues(lclMV_),j);
01227       ncbuff = Teuchos::arcp_const_cast<Scalar>(cbuff);
01228     }
01229     return Teuchos::rcp<const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(
01230               new Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),ncbuff)
01231            );
01232   }
01233 
01234 
01235   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01236   Teuchos::RCP<Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 
01237   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getVectorNonConst(size_t j) {
01238 #ifdef HAVE_TPETRA_DEBUG
01239     TEST_FOR_EXCEPTION(j < 0 || j >= this->getNumVectors(), std::runtime_error,
01240         "Tpetra::MultiVector::getVectorNonConst(j): index j (== " << j << ") exceeds valid column range for this multivector.");
01241 #endif
01242     Teuchos::ArrayRCP<Scalar> ncbuff;
01243     if (getLocalLength() > 0) {
01244       ncbuff = getSubArrayRCP(MVT::getValuesNonConst(lclMV_),j);
01245     }
01246     return Teuchos::rcp(new Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getMap(),ncbuff));
01247   }
01248 
01249 
01250   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01251   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dCopy(Teuchos::ArrayView<Scalar> A, size_t LDA) const
01252   {
01253     using Teuchos::ArrayRCP;
01254     TEST_FOR_EXCEPTION(LDA < getLocalLength(), std::runtime_error,
01255       "Tpetra::MultiVector::get1dCopy(A,LDA): specified stride is not large enough for local vector length.");
01256     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(A.size()) < LDA*(getNumVectors()-1)+getLocalLength(), std::runtime_error,
01257       "Tpetra::MultiVector::get1dCopy(A,LDA): specified stride/storage is not large enough for the number of vectors.");
01258     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01259     const size_t myStride = MVT::getStride(lclMV_),
01260                   numCols = getNumVectors(),
01261                   myLen   = getLocalLength();
01262     if (myLen > 0) {
01263       ArrayRCP<const Scalar> mydata = MVT::getValues(lclMV_);
01264       KOKKOS_NODE_TRACE("MultiVector::getVectorNonConst()")
01265       ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(myStride*(numCols-1)+myLen,mydata);
01266       typename Teuchos::ArrayView<Scalar>::iterator Aptr = A.begin();
01267       for (size_t j=0; j<numCols; j++) {
01268         ArrayRCP<const Scalar> myviewj = getSubArrayRCP(myview,j);
01269         std::copy(myviewj,myviewj+myLen,Aptr);
01270         Aptr += LDA;
01271       }
01272       myview = Teuchos::null;
01273       mydata = Teuchos::null;
01274     }
01275   }
01276 
01277 
01278   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01279   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get2dCopy(Teuchos::ArrayView<const Teuchos::ArrayView<Scalar> > ArrayOfPtrs) const
01280   {
01281     TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ArrayOfPtrs.size()) != getNumVectors(), std::runtime_error,
01282         "Tpetra::MultiVector::get2dCopy(ArrayOfPtrs): Array of pointers must contain as many pointers as the MultiVector has rows.");
01283     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01284     const size_t numCols = getNumVectors(),
01285                    myLen = getLocalLength();
01286     if (myLen > 0) {
01287       Teuchos::ArrayRCP<const Scalar> mybuff = MVT::getValues(lclMV_);
01288       KOKKOS_NODE_TRACE("MultiVector::get2dCopy()")
01289       Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(mybuff.size(), mybuff);
01290       for (size_t j=0; j<numCols; ++j) {
01291 #ifdef HAVE_TPETRA_DEBUG
01292         TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ArrayOfPtrs[j].size()) != getLocalLength(), std::runtime_error,
01293             "Tpetra::MultiVector::get2dCopy(ArrayOfPtrs): The ArrayView provided in ArrayOfPtrs[" << j << "] was not large enough to contain the local entries.");
01294 #endif
01295         Teuchos::ArrayRCP<const Scalar> myviewj = getSubArrayRCP(myview,j);
01296         std::copy(myviewj,myviewj+myLen,ArrayOfPtrs[j].begin());
01297       }
01298       myview = Teuchos::null;
01299       mybuff = Teuchos::null;
01300     }
01301   }
01302 
01303 
01304   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01305   Teuchos::ArrayRCP<const Scalar> MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dView() const
01306   {
01307     TEST_FOR_EXCEPTION(!isConstantStride(), std::runtime_error,
01308       "Tpetra::MultiVector::get1dView() requires that this MultiVector have constant stride.");
01309     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01310       KOKKOS_NODE_TRACE("MultiVector::get1dView()")
01311     return node->template viewBuffer<Scalar>( getStride()*(getNumVectors()-1)+getLocalLength(), MVT::getValues(lclMV_) );
01312   }
01313 
01314 
01315   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01316   Teuchos::ArrayRCP<Scalar> MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get1dViewNonConst()
01317   {
01318     TEST_FOR_EXCEPTION(!isConstantStride(), std::runtime_error,
01319       "Tpetra::MultiVector::get1dViewNonConst(): requires that this MultiVector have constant stride.");
01320     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01321       KOKKOS_NODE_TRACE("MultiVector::get1dViewNonConst()")
01322     return node->template viewBufferNonConst<Scalar>( Kokkos::ReadWrite, getStride()*(getNumVectors()-1)+getLocalLength(), MVT::getValuesNonConst(lclMV_) );
01323   }
01324 
01325 
01326   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01327   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >
01328   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get2dViewNonConst()
01329   {
01330     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01331     Teuchos::ArrayRCP< Teuchos::ArrayRCP<Scalar> > views = Teuchos::arcp<Teuchos::ArrayRCP<Scalar> >(getNumVectors());
01332     if (isConstantStride()) {
01333       const size_t myStride = MVT::getStride(lclMV_),
01334                     numCols = getNumVectors(),
01335                     myLen   = getLocalLength();
01336       if (myLen > 0) {
01337         KOKKOS_NODE_TRACE("MultiVector::get2dViewNonConst()")
01338         Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::ReadWrite,myStride*(numCols-1)+myLen,MVT::getValuesNonConst(lclMV_));
01339         for (size_t j=0; j<numCols; ++j) {
01340           views[j] = myview.persistingView(0,myLen);
01341           myview += myStride;
01342         }
01343       }
01344     }
01345     else {
01346       const size_t myStride = MVT::getStride(lclMV_),
01347                     numCols = MVT::getNumCols(lclMV_),
01348                      myCols = getNumVectors(),
01349                      myLen  = MVT::getNumRows(lclMV_);
01350       if (myLen > 0) {
01351         KOKKOS_NODE_TRACE("MultiVector::get2dViewNonConst()")
01352         Teuchos::ArrayRCP<Scalar> myview = node->template viewBufferNonConst<Scalar>(Kokkos::ReadWrite,myStride*(numCols-1)+myLen,MVT::getValuesNonConst(lclMV_));
01353         for (size_t j=0; j<myCols; ++j) {
01354           views[j] = myview.persistingView(whichVectors_[j]*myStride,myLen);
01355         }
01356       }
01357     }
01358     return views;
01359   }
01360 
01361 
01362   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01363   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > 
01364   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::get2dView() const
01365   {
01366     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01367     Teuchos::ArrayRCP< Teuchos::ArrayRCP<const Scalar> > views = Teuchos::arcp<Teuchos::ArrayRCP<const Scalar> >(getNumVectors());
01368     if (isConstantStride()) {
01369       const size_t myStride = MVT::getStride(lclMV_),
01370                     numCols = getNumVectors(),
01371                     myLen   = getLocalLength();
01372       if (myLen > 0) {
01373         KOKKOS_NODE_TRACE("MultiVector::get2dView()")
01374         Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(myStride*(numCols-1)+myLen,MVT::getValues(lclMV_));
01375         for (size_t j=0; j<numCols; ++j) {
01376           views[j] = myview.persistingView(0,myLen);
01377           myview += myStride;
01378         }
01379       }
01380     }
01381     else {
01382       const size_t myStride = MVT::getStride(lclMV_),
01383                     numCols = MVT::getNumCols(lclMV_),
01384                      myCols = getNumVectors(),
01385                      myLen  = MVT::getNumRows(lclMV_);
01386       if (myLen > 0) {
01387         KOKKOS_NODE_TRACE("MultiVector::get2dView()")
01388         Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(myStride*(numCols-1)+myLen,MVT::getValues(lclMV_));
01389         for (size_t j=0; j<myCols; ++j) {
01390           views[j] = myview.persistingView(whichVectors_[j]*myStride,myLen);
01391         }
01392       }
01393     }
01394     return views;
01395   }
01396 
01397 
01398   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01399   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::multiply(
01400       Teuchos::ETransp transA, Teuchos::ETransp transB, 
01401       const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B, 
01402       const Scalar &beta) {
01403     using Teuchos::NO_TRANS;      // enums
01404     using Teuchos::TRANS;
01405     using Teuchos::CONJ_TRANS;
01406     using Teuchos::null;
01407     using Teuchos::ScalarTraits;  // traits
01408     using Teuchos::as;
01409     using Teuchos::RCP;           // data structures
01410     using Teuchos::Array;
01411     using Teuchos::ArrayView;
01412     using Teuchos::rcp;           // initializers for data structures
01413 
01414     // This routine performs a variety of matrix-matrix multiply operations, interpreting
01415     // the MultiVector (this-aka C , A and B) as 2D matrices.  Variations are due to
01416     // the fact that A, B and C can be local replicated or global distributed
01417     // MultiVectors and that we may or may not operate with the transpose of 
01418     // A and B.  Possible cases are:
01419     //                                       Num
01420     //      OPERATIONS                        cases  Notes
01421     //  1) C(local) = A^X(local) * B^X(local)  4    (X=Trans or Not, No comm needed) 
01422     //  2) C(local) = A^T(distr) * B  (distr)  1    (2D dot product, replicate C)
01423     //  3) C(distr) = A  (distr) * B^X(local)  2    (2D vector update, no comm needed)
01424     //
01425     // The following operations are not meaningful for 1D distributions:
01426     //
01427     // u1) C(local) = A^T(distr) * B^T(distr)  1
01428     // u2) C(local) = A  (distr) * B^X(distr)  2
01429     // u3) C(distr) = A^X(local) * B^X(local)  4
01430     // u4) C(distr) = A^X(local) * B^X(distr)  4
01431     // u5) C(distr) = A^T(distr) * B^X(local)  2
01432     // u6) C(local) = A^X(distr) * B^X(local)  4
01433     // u7) C(distr) = A^X(distr) * B^X(local)  4
01434     // u8) C(local) = A^X(local) * B^X(distr)  4
01435     //
01436     // Total of 32 case (2^5).
01437 
01438     std::string errPrefix("Tpetra::MultiVector::multiply(transOpA,transOpB,alpha,A,B,beta): ");
01439 
01440     TEST_FOR_EXCEPTION( ScalarTraits<Scalar>::isComplex && (transA == TRANS || transB == TRANS), std::invalid_argument,
01441         errPrefix << "non-conjugate transpose not supported for complex types.");
01442     transA = (transA == NO_TRANS ? NO_TRANS : CONJ_TRANS);
01443     transB = (transB == NO_TRANS ? NO_TRANS : CONJ_TRANS);
01444 
01445     // Compute effective dimensions, w.r.t. transpose operations on 
01446     size_t A_nrows = (transA==CONJ_TRANS) ? A.getNumVectors() : A.getLocalLength();
01447     size_t A_ncols = (transA==CONJ_TRANS) ? A.getLocalLength() : A.getNumVectors();
01448     size_t B_nrows = (transB==CONJ_TRANS) ? B.getNumVectors() : B.getLocalLength();
01449     size_t B_ncols = (transB==CONJ_TRANS) ? B.getLocalLength() : B.getNumVectors();
01450 
01451     Scalar beta_local = beta; // local copy of beta; might be reassigned below
01452 
01453     TEST_FOR_EXCEPTION( getLocalLength() != A_nrows || getNumVectors() != B_ncols || A_ncols != B_nrows, std::runtime_error,
01454         errPrefix << "dimension of *this, op(A) and op(B) must be consistent.");
01455 
01456     bool A_is_local = !A.isDistributed();
01457     bool B_is_local = !B.isDistributed();
01458     bool C_is_local = !this->isDistributed();
01459     bool Case1 = ( C_is_local &&  A_is_local &&  B_is_local);                                           // Case 1: C(local) = A^X(local) * B^X(local)
01460     bool Case2 = ( C_is_local && !A_is_local && !B_is_local && transA==CONJ_TRANS && transB==NO_TRANS); // Case 2: C(local) = A^T(distr) * B  (distr)
01461     bool Case3 = (!C_is_local && !A_is_local &&  B_is_local && transA==NO_TRANS  );                     // Case 3: C(distr) = A  (distr) * B^X(local)
01462 
01463     // Test that we are considering a meaningful cases
01464     TEST_FOR_EXCEPTION( !Case1 && !Case2 && !Case3, std::runtime_error,
01465         errPrefix << "multiplication of op(A) and op(B) into *this is not a supported use case.");
01466 
01467     if (beta != ScalarTraits<Scalar>::zero() && Case2)
01468     {
01469       // if Case2, then C is local and contributions must be summed across all nodes
01470       // however, if beta != 0, then accumulate beta*C into the sum
01471       // when summing across all nodes, we only want to accumulate this once, so 
01472       // set beta == 0 on all nodes except node 0
01473       int MyPID = this->getMap()->getComm()->getRank();
01474       if (MyPID!=0) beta_local = ScalarTraits<Scalar>::zero();
01475     }
01476 
01477     // Check if A, B, C have constant stride, if not then make temp copy (strided)
01478     RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Atmp, Btmp; 
01479     RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >       Ctmp;
01480     if (isConstantStride() == false) Ctmp = rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(*this));
01481     else Ctmp = rcp(this,false);
01482 
01483     if (A.isConstantStride() == false) Atmp = rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A));
01484     else Atmp = rcp(&A,false);
01485 
01486     if (B.isConstantStride() == false) Btmp = rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(B));
01487     else Btmp = rcp(&B,false);
01488 
01489 #ifdef HAVE_TEUCHOS_DEBUG
01490     TEST_FOR_EXCEPTION(!Ctmp->isConstantStride() || !Btmp->isConstantStride() || !Atmp->isConstantStride(), std::logic_error,
01491         errPrefix << "failed making temporary strided copies of input multivectors.");
01492 #endif
01493 
01494     KMV &C_mv = Ctmp->lclMV_;
01495     {
01496       // get local multivectors
01497       const KMV &A_mv = Atmp->lclMV_;
01498       const KMV &B_mv = Btmp->lclMV_;
01499       // do the multiply (GEMM)
01500       MVT::GEMM(C_mv,transA,transB,alpha,A_mv,B_mv,beta_local);
01501     }
01502 
01503     // Dispose of (possibly) extra copies of A, B
01504     Atmp = null;
01505     Btmp = null;
01506 
01507     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01508     // If *this was not strided, copy the data from the strided version and then delete it
01509     if (isConstantStride() == false) {
01510       // *this is not strided, we must put data from Ctmp into *this
01511       TEST_FOR_EXCEPT(&C_mv != &lclMV_);
01512       const size_t numVecs = MVT::getNumCols(lclMV_);
01513       for (size_t j=0; j < numVecs; ++j) {
01514         KOKKOS_NODE_TRACE("MultiVector::multiply()")
01515         node->template copyBuffers<Scalar>(getLocalLength(),MVT::getValues(C_mv,j),MVT::getValuesNonConst(lclMV_,whichVectors_[j]));
01516       }
01517     }
01518 
01519     // If Case 2 then sum up *this and distribute it to all processors.
01520     if (Case2) {
01521       this->reduce();
01522     }
01523   }
01524 
01525   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01526   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::elementWiseMultiply(Scalar scalarAB, const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B, Scalar scalarThis)
01527   {
01528     using Teuchos::ArrayRCP;
01529     using Teuchos::arcp_const_cast;
01530 #ifdef HAVE_TPETRA_DEBUG
01531     TEST_FOR_EXCEPTION( !this->getMap()->isCompatible(*A.getMap()) || !this->getMap()->isCompatible(*B.getMap()), std::runtime_error,
01532         "Tpetra::MultiVector::elementWiseMultiply(): MultiVectors do not have compatible Maps:" << std::endl
01533         << "this->getMap(): " << std::endl << *this->getMap() 
01534         << "A.getMap(): " << std::endl << *A.getMap() << std::endl
01535         << "B.getMap(): " << std::endl << *B.getMap() << std::endl);
01536 #else
01537     TEST_FOR_EXCEPTION( getLocalLength() != A.getLocalLength() || getLocalLength() != B.getLocalLength(), std::runtime_error,
01538         "Tpetra::MultiVector::reciprocal(): MultiVectors do not have the same local length.");
01539 #endif
01540     TEST_FOR_EXCEPTION(B.getNumVectors() != this->getNumVectors(), std::runtime_error,
01541         "Tpetra::MultiVector::reciprocal(): MultiVectors 'this' and B must have the same number of vectors.");
01542     TEST_FOR_EXCEPTION(A.getNumVectors() != 1, std::runtime_error,
01543         "Tpetra::MultiVector::reciprocal(): MultiVectors A must have just 1 vector.");
01544     try {
01545       MVT::ElemMult(lclMV_,scalarThis,scalarAB,(const KMV &)A.lclMV_,(const KMV &)B.lclMV_);
01546     }
01547     catch (std::runtime_error &e) {
01548       TEST_FOR_EXCEPTION(true,std::runtime_error,
01549           "Tpetra::MultiVector::elementWiseMultiply(): caught exception from Kokkos:" << std::endl
01550           << e.what() << std::endl);
01551     }
01552   }
01553 
01554   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01555   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::reduce() {
01556     using Teuchos::ArrayRCP;
01557     using Teuchos::ArrayView;
01558     // this should be called only for "local" MultiVectors (!isDistributed())
01559     TEST_FOR_EXCEPTION(this->isDistributed() == true, std::runtime_error,
01560         "Tpetra::MultiVector::reduce() should only be called for non-distributed MultiVectors.");
01561     Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getMap()->getComm();
01562     if (comm->getSize() == 1) return;
01563     Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01564     // sum the data across all multivectors
01565     // need to have separate packed buffers for send and receive
01566     // if we're packed, we'll just set the receive buffer as our data, the send as a copy
01567     // if we're not packed, we'll use allocated buffers for both. 
01568     ArrayView<Scalar> target;
01569     const size_t myStride = MVT::getStride(lclMV_),
01570                  numCols  = MVT::getNumCols(lclMV_),
01571                  myLen    = MVT::getNumRows(lclMV_);
01572     Teuchos::Array<Scalar> sourceBuffer(numCols*myLen), tmparr(0);
01573     bool packed = isConstantStride() && (myStride == myLen);
01574     KOKKOS_NODE_TRACE("MultiVector::reduce()")
01575     ArrayRCP<Scalar> bufView = node->template viewBufferNonConst<Scalar>(
01576                                           Kokkos::ReadWrite,myStride*(numCols-1)+myLen,
01577                                           MVT::getValuesNonConst(lclMV_) );
01578     if (packed) {
01579       // copy data from view to sourceBuffer, reduce into view below
01580       target = bufView(0,myLen*numCols);
01581       std::copy(target.begin(),target.end(),sourceBuffer.begin());
01582     }
01583     else {
01584       // copy data into sourceBuffer, reduce from sourceBuffer into tmparr below, copy back to view after that
01585       tmparr.resize(myLen*numCols);
01586       Scalar *sptr = sourceBuffer.getRawPtr();
01587       ArrayRCP<const Scalar> vptr = bufView;
01588       for (size_t j=0; j<numCols; ++j) 
01589       {
01590         std::copy(vptr,vptr+myLen,sptr);
01591         sptr += myLen;
01592         vptr += myStride;
01593       }
01594       target = tmparr();
01595     }
01596     // reduce 
01597     Teuchos::reduceAll<int,Scalar>(*comm,Teuchos::REDUCE_SUM,Teuchos::as<int>(numCols*myLen),sourceBuffer.getRawPtr(),target.getRawPtr());
01598     if (!packed) {
01599       // copy tmparr back into view
01600       const Scalar *sptr = tmparr.getRawPtr();
01601       ArrayRCP<Scalar> vptr = bufView;
01602       for (size_t j=0; j<numCols; ++j) 
01603       {
01604         std::copy(sptr,sptr+myLen,vptr);
01605         sptr += myLen;
01606         vptr += myStride;
01607       }
01608     }
01609     bufView = Teuchos::null;
01610   }
01611 
01612 
01613   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01614   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceLocalValue(LocalOrdinal MyRow, size_t VectorIndex, const Scalar &ScalarValue)
01615   {
01616 #ifdef HAVE_TPETRA_DEBUG
01617     TEST_FOR_EXCEPTION(MyRow < this->getMap()->getMinLocalIndex() || MyRow > this->getMap()->getMaxLocalIndex(), std::runtime_error,
01618         "Tpetra::MultiVector::replaceLocalValue(): row index is invalid.");
01619     TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error,
01620         "Tpetra::MultiVector::replaceLocalValue(): vector index is invalid.");
01621 #endif
01622     getDataNonConst(VectorIndex)[MyRow] = ScalarValue;
01623   }
01624 
01625 
01626   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01627   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoLocalValue(LocalOrdinal MyRow, size_t VectorIndex, const Scalar &ScalarValue)
01628   {
01629 #ifdef HAVE_TPETRA_DEBUG
01630     TEST_FOR_EXCEPTION(MyRow < this->getMap()->getMinLocalIndex() || MyRow > this->getMap()->getMaxLocalIndex(), std::runtime_error,
01631         "Tpetra::MultiVector::sumIntoLocalValue(): row index is invalid.");
01632     TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error,
01633         "Tpetra::MultiVector::sumIntoLocalValue(): vector index is invalid.");
01634 #endif
01635     getDataNonConst(VectorIndex)[MyRow] += ScalarValue;
01636   }
01637 
01638 
01639   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01640   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceGlobalValue(GlobalOrdinal GlobalRow, size_t VectorIndex, const Scalar &ScalarValue)
01641   {
01642     LocalOrdinal MyRow = this->getMap()->getLocalElement(GlobalRow);
01643 #ifdef HAVE_TPETRA_DEBUG
01644     TEST_FOR_EXCEPTION(MyRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error,
01645         "Tpetra::MultiVector::replaceGlobalValue(): row index is not present on this processor.");
01646     TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error,
01647         "Tpetra::MultiVector::replaceGlobalValue(): vector index is invalid.");
01648 #endif
01649     getDataNonConst(VectorIndex)[MyRow] = ScalarValue;
01650   }
01651 
01652 
01653   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01654   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::sumIntoGlobalValue(GlobalOrdinal GlobalRow, size_t VectorIndex, const Scalar &ScalarValue)
01655   {
01656     LocalOrdinal MyRow = this->getMap()->getLocalElement(GlobalRow);
01657 #ifdef HAVE_TEUCHOS_DEBUG
01658     TEST_FOR_EXCEPTION(MyRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error,
01659         "Tpetra::MultiVector::sumIntoGlobalValue(): row index is not present on this processor.");
01660     TEST_FOR_EXCEPTION(VectorIndex < 0 || VectorIndex >= getNumVectors(), std::runtime_error,
01661         "Tpetra::MultiVector::sumIntoGlobalValue(): vector index is invalid.");
01662 #endif
01663     getDataNonConst(VectorIndex)[MyRow] += ScalarValue;
01664   }
01665 
01666 
01667   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01668   template <class T>
01669   Teuchos::ArrayRCP<T> MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getSubArrayRCP(Teuchos::ArrayRCP<T> arr, size_t j) const {
01670     const size_t stride = MVT::getStride(lclMV_),
01671                   myLen = getLocalLength();
01672     Teuchos::ArrayRCP<T> ret;
01673     if (isConstantStride()) {
01674       ret = arr.persistingView(j*stride,myLen);
01675     }
01676     else {
01677       ret = arr.persistingView(whichVectors_[j]*stride,myLen);
01678     }
01679     return ret;
01680   }
01681 
01682   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01683   const Kokkos::MultiVector<Scalar,Node> & 
01684   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getLocalMV() const {
01685     return lclMV_;
01686   }
01687 
01688   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01689   Kokkos::MultiVector<Scalar,Node> & 
01690   MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getLocalMVNonConst() {
01691     return lclMV_;
01692   }
01693 
01694   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01695   std::string MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const
01696   {
01697     std::ostringstream oss;
01698     oss << Teuchos::Describable::description();
01699     oss << "{length = "<<getGlobalLength()
01700         << ", getNumVectors = "<<getNumVectors()
01701         << ", isConstantStride() = "<<isConstantStride()
01702         << "}";
01703     return oss.str();
01704   }
01705 
01706   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01707   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
01708     using std::endl;
01709     using std::setw;
01710     using Teuchos::VERB_DEFAULT;
01711     using Teuchos::VERB_NONE;
01712     using Teuchos::VERB_LOW;
01713     using Teuchos::VERB_MEDIUM;
01714     using Teuchos::VERB_HIGH;
01715     using Teuchos::VERB_EXTREME;
01716     Teuchos::EVerbosityLevel vl = verbLevel;
01717     if (vl == VERB_DEFAULT) vl = VERB_LOW;
01718     Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getMap()->getComm();
01719     const int myImageID = comm->getRank(),
01720               numImages = comm->getSize();
01721     size_t width = 1;
01722     for (size_t dec=10; dec<getGlobalLength(); dec *= 10) {
01723       ++width;
01724     }
01725     Teuchos::OSTab tab(out);
01726     if (vl != VERB_NONE) {
01727       // VERB_LOW and higher prints description()
01728       if (myImageID == 0) out << this->description() << std::endl; 
01729       for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
01730         if (myImageID == imageCtr) {
01731           if (vl != VERB_LOW) {
01732             // VERB_MEDIUM and higher prints getLocalLength()
01733             out << "node " << setw(width) << myImageID << ": local length=" << getLocalLength();
01734             if (vl != VERB_MEDIUM) {
01735               // VERB_HIGH and higher prints isConstantStride() and getStride()
01736               if (isConstantStride()) out << ", constant stride=" << getStride() << endl;
01737               else out << ", non-constant stride" << endl;
01738               if (vl == VERB_EXTREME && getLocalLength() > 0) {
01739                 Teuchos::RCP<Node> node = MVT::getNode(lclMV_);
01740                 if (isConstantStride()) {
01741                   KOKKOS_NODE_TRACE("MultiVector::describe()")
01742                   Teuchos::ArrayRCP<const Scalar> myview = node->template viewBuffer<Scalar>(
01743                                           getLocalLength()+getStride()*(getNumVectors()-1), 
01744                                           MVT::getValues(lclMV_) );
01745                   // VERB_EXTREME prints values
01746                   for (size_t i=0; i<getLocalLength(); ++i) {
01747                     out << setw(width) << this->getMap()->getGlobalElement(i) << ": ";
01748                     for (size_t j=0; j<getNumVectors(); ++j) {
01749                       out << myview[j*getStride()] << "  ";
01750                     }
01751                     ++myview;
01752                     out << endl;
01753                   }
01754                   myview = Teuchos::null;
01755                 }
01756                 else {
01757                   const size_t stride = MVT::getStride(lclMV_),
01758                                rows   = MVT::getNumRows(lclMV_),
01759                                cols   = MVT::getNumCols(lclMV_);
01760                   KOKKOS_NODE_TRACE("MultiVector::describe()")
01761                   Teuchos::ArrayRCP<const Scalar> myview = 
01762                     node->template viewBuffer<Scalar>( rows + stride * (cols - 1), MVT::getValues(lclMV_) );
01763                   // VERB_EXTREME prints values
01764                   for (size_t i=0; i<getLocalLength(); ++i) {
01765                     out << setw(width) << this->getMap()->getGlobalElement(i) << ": ";
01766                     for (size_t j=0; j<getNumVectors(); ++j) {
01767                       out << myview[whichVectors_[j]*stride + i] << "  ";
01768                     }
01769                     out << endl;
01770                   }
01771                   myview = Teuchos::null;
01772                 }
01773               }
01774             }
01775             else {
01776               out << endl;
01777             }
01778           }
01779         }
01780       }
01781     }
01782   }
01783 
01784   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01785   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createViews() const {
01786     Teuchos::RCP<Node> node = this->getMap()->getNode();
01787     if (cview_ == Teuchos::null && getLocalLength() > 0) {
01788       Teuchos::ArrayRCP<const Scalar> buff = MVT::getValues(lclMV_);
01789       KOKKOS_NODE_TRACE("MultiVector::createViews()")
01790       cview_ = node->template viewBuffer<Scalar>(buff.size(),buff);
01791     }
01792   }
01793 
01794   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01795   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createViewsNonConst(Kokkos::ReadWriteOption rwo) {
01796     Teuchos::RCP<Node> node = this->getMap()->getNode();
01797     if (ncview_ == Teuchos::null && getLocalLength() > 0) {
01798       Teuchos::ArrayRCP<Scalar> buff = MVT::getValuesNonConst(lclMV_);
01799       KOKKOS_NODE_TRACE("MultiVector::createViewsNonConst()")
01800       ncview_ = node->template viewBufferNonConst<Scalar>(rwo,buff.size(),buff);
01801     }
01802   }
01803 
01804   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01805   void MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::releaseViews() const {
01806     cview_ = Teuchos::null;
01807     ncview_ = Teuchos::null;
01808   }
01809 
01810 
01811 } // namespace Tpetra
01812 
01813 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01814 Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
01815 Tpetra::createMultiVector(const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > &map, size_t numVectors) {
01816   const bool DO_INIT_TO_ZERO = true;
01817   return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map,numVectors,DO_INIT_TO_ZERO) );
01818 }
01819 
01820 //
01821 // Explicit instantiation macro
01822 //
01823 // Must be expanded from within the Tpetra namespace!
01824 //
01825 
01826 #define TPETRA_MULTIVECTOR_INSTANT(SCALAR,LO,GO,NODE) \
01827   \
01828   template class MultiVector< SCALAR , LO , GO , NODE >; \
01829   \
01830   template Teuchos::RCP< MultiVector<SCALAR,LO,GO,NODE> > \
01831   createMultiVector<SCALAR,LO,GO,NODE>(const Teuchos::RCP<const Map<LO,GO,NODE> > &map, size_t numVectors); 
01832 
01833 
01834 #endif // TPETRA_MULTIVECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines