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