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