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