Tpetra_MultiVector_def.hpp

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

Generated on Tue Jul 13 09:39:07 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7