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