Tpetra_MultiVector_decl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef TPETRA_MULTIVECTOR_DECL_HPP
#define TPETRA_MULTIVECTOR_DECL_HPP
 
 
#include "Tpetra_MultiVector_fwd.hpp"
#include "Tpetra_Vector_fwd.hpp"
#include "Tpetra_FEMultiVector_fwd.hpp"
#include "Tpetra_DistObject.hpp"
#include "Tpetra_Map_fwd.hpp"
#include "Tpetra_Details_Behavior.hpp"
#include "Kokkos_DualView.hpp"
#include "Teuchos_BLAS_types.hpp"
#include "Teuchos_DataAccess.hpp"
#include "Teuchos_Range1D.hpp"
#include "KokkosKernels_ArithTraits.hpp"
#ifndef KOKKOS_ENABLE_DEPRECATED_CODE_5
#include "KokkosKernels_InnerProductSpaceTraits.hpp"
#else
#include "Kokkos_InnerProductSpaceTraits.hpp"
#endif
#include "Tpetra_KokkosRefactor_Details_MultiVectorLocalDeepCopy.hpp"
#include "Tpetra_Access.hpp"
#include "Tpetra_Details_WrappedDualView.hpp"
#include <type_traits>
 
#ifdef HAVE_TPETRACORE_TEUCHOSNUMERICS
#ifndef DOXYGEN_SHOULD_SKIP_THIS
namespace Teuchos {
template <class OrdinalType, class ScalarType>
class SerialDenseMatrix;  // forward declaration
}
#endif  // DOXYGEN_SHOULD_SKIP_THIS
#endif  // HAVE_TPETRACORE_TEUCHOSNUMERICS
 
namespace Tpetra {
 
template <class DS, class DL, class DG, class DN,
          class SS, class SL, class SG, class SN>
void deep_copy(MultiVector<DS, DL, DG, DN>& dst,
               const MultiVector<SS, SL, SG, SN>& src);
 
#ifdef HAVE_TPETRACORE_TEUCHOSNUMERICS
template <class ST, class LO, class GO, class NT>
void deep_copy(MultiVector<ST, LO, GO, NT>& dst,
               const Teuchos::SerialDenseMatrix<int, ST>& src);
 
template <class ST, class LO, class GO, class NT>
void deep_copy(Teuchos::SerialDenseMatrix<int, ST>& dst,
               const MultiVector<ST, LO, GO, NT>& src);
#endif  // HAVE_TPETRACORE_TEUCHOSNUMERICS
 
template <class ST, class LO, class GO, class NT>
MultiVector<ST, LO, GO, NT>
createCopy(const MultiVector<ST, LO, GO, NT>& src);
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
createMultiVector(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node>>& map,
                  const size_t numVectors);
 
// WARNING NOT FOR USERS
// This means we don't need to make MultiVector a friend of
// Vector or of itself (with different template parameters).
template <class SC, class LO, class GO, class NT>
Teuchos::ArrayView<const size_t>
getMultiVectorWhichVectors(const MultiVector<SC, LO, GO, NT>& X);
 
template <class Scalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node>
class MultiVector : public DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
 public:
 
 
  using scalar_type = Scalar;
  using impl_scalar_type =
      typename KokkosKernels::ArithTraits<Scalar>::val_type;
 
  using map_type = Map<LocalOrdinal, GlobalOrdinal, Node>;
  using local_ordinal_type = typename map_type::local_ordinal_type;
  using global_ordinal_type = typename map_type::global_ordinal_type;
  using device_type = typename map_type::device_type;
  using node_type = typename map_type::node_type;
 
#ifndef KOKKOS_ENABLE_DEPRECATED_CODE_5
  using dot_type =
      typename KokkosKernels::Details::InnerProductSpaceTraits<impl_scalar_type>::dot_type;
#else
  using dot_type =
      typename Kokkos::Details::InnerProductSpaceTraits<impl_scalar_type>::dot_type;
#endif
 
  using mag_type = typename KokkosKernels::ArithTraits<impl_scalar_type>::mag_type;
 
  using execution_space = typename device_type::execution_space;
 
  using dual_view_type         = Kokkos::DualView<impl_scalar_type**,
                                          Kokkos::LayoutLeft,
                                          device_type>;
  using wrapped_dual_view_type = Details::WrappedDualView<dual_view_type>;
 
  using host_view_type   = typename dual_view_type::t_host;
  using device_view_type = typename dual_view_type::t_dev;
 
 
 
  MultiVector();
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const size_t numVecs,
              const bool zeroOut = true);
 
  MultiVector(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& source,
              const Teuchos::DataAccess copyOrView);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const Teuchos::ArrayView<const Scalar>& A,
              const size_t LDA,
              const size_t NumVectors);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar>>& ArrayOfPtrs,
              const size_t NumVectors);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const dual_view_type& view);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const typename dual_view_type::t_dev& d_view);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const dual_view_type& view,
              const dual_view_type& origView);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const wrapped_dual_view_type& d_view);
 
 protected:
  MultiVector(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& X,
              const size_t j);
 
 public:
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const dual_view_type& view,
              const Teuchos::ArrayView<const size_t>& whichVectors);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const wrapped_dual_view_type& view,
              const Teuchos::ArrayView<const size_t>& whichVectors);
 
  MultiVector(const Teuchos::RCP<const map_type>& map,
              const dual_view_type& view,
              const dual_view_type& origView,
              const Teuchos::ArrayView<const size_t>& whichVectors);
 
  MultiVector(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& X,
              const Teuchos::RCP<const map_type>& subMap,
              const local_ordinal_type rowOffset = 0);
 
  MultiVector(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& X,
              const map_type& subMap,
              const size_t offset = 0);
 
  MultiVector(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&) = default;
 
  MultiVector(MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&&) = default;
 
  MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
  operator=(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&) = default;
 
  MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
  operator=(MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&&) = default;
 
  virtual ~MultiVector() = default;
 
  void swap(MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& mv);
 
 
 
 protected:
  static const bool useAtomicUpdatesByDefault =
#ifdef KOKKOS_ENABLE_SERIAL
      !std::is_same<execution_space, Kokkos::Serial>::value;
#else
      true;
#endif  // KOKKOS_ENABLE_SERIAL
 
 public:
  void
  replaceGlobalValue(const GlobalOrdinal gblRow,
                     const size_t col,
                     const impl_scalar_type& value);
 
  template <typename T>
  typename std::enable_if<!std::is_same<T, impl_scalar_type>::value && std::is_convertible<T, impl_scalar_type>::value, void>::type
  replaceGlobalValue(GlobalOrdinal globalRow,
                     size_t col,
                     const T& value) {
    replaceGlobalValue(globalRow, col, static_cast<impl_scalar_type>(value));
  }
 
  void
  sumIntoGlobalValue(const GlobalOrdinal gblRow,
                     const size_t col,
                     const impl_scalar_type& value,
                     const bool atomic = useAtomicUpdatesByDefault);
 
  template <typename T>
  typename std::enable_if<!std::is_same<T, impl_scalar_type>::value && std::is_convertible<T, impl_scalar_type>::value, void>::type
  sumIntoGlobalValue(const GlobalOrdinal gblRow,
                     const size_t col,
                     const T& val,
                     const bool atomic = useAtomicUpdatesByDefault) {
    sumIntoGlobalValue(gblRow, col, static_cast<impl_scalar_type>(val), atomic);
  }
 
  void
  replaceLocalValue(const LocalOrdinal lclRow,
                    const size_t col,
                    const impl_scalar_type& value);
 
  template <typename T>
  typename std::enable_if<!std::is_same<T, impl_scalar_type>::value && std::is_convertible<T, impl_scalar_type>::value, void>::type
  replaceLocalValue(const LocalOrdinal lclRow,
                    const size_t col,
                    const T& val) {
    replaceLocalValue(lclRow, col, static_cast<impl_scalar_type>(val));
  }
 
  void
  sumIntoLocalValue(const LocalOrdinal lclRow,
                    const size_t col,
                    const impl_scalar_type& val,
                    const bool atomic = useAtomicUpdatesByDefault);
 
  template <typename T>
  typename std::enable_if<!std::is_same<T, impl_scalar_type>::value && std::is_convertible<T, impl_scalar_type>::value, void>::type
  sumIntoLocalValue(const LocalOrdinal lclRow,
                    const size_t col,
                    const T& val,
                    const bool atomic = useAtomicUpdatesByDefault) {
    sumIntoLocalValue(lclRow, col, static_cast<impl_scalar_type>(val), atomic);
  }
 
  void putScalar(const Scalar& value);
 
  template <typename T>
  typename std::enable_if<!std::is_same<T, impl_scalar_type>::value && std::is_convertible<T, impl_scalar_type>::value, void>::type
  putScalar(const T& value) {
    putScalar(static_cast<impl_scalar_type>(value));
  }
 
  void randomize();
 
  void randomize(const Scalar& minVal, const Scalar& maxVal);
 
  void replaceMap(const Teuchos::RCP<const map_type>& map);
 
  void reduce();
 
 
 
 
  Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  subCopy(const Teuchos::Range1D& colRng) const;
 
  Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  subCopy(const Teuchos::ArrayView<const size_t>& cols) const;
 
  Teuchos::RCP<const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  subView(const Teuchos::Range1D& colRng) const;
 
  Teuchos::RCP<const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  subView(const Teuchos::ArrayView<const size_t>& cols) const;
 
  Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  subViewNonConst(const Teuchos::Range1D& colRng);
 
  Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  subViewNonConst(const Teuchos::ArrayView<const size_t>& cols);
 
  Teuchos::RCP<const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  offsetView(const Teuchos::RCP<const map_type>& subMap,
             const size_t offset) const;
 
  Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  offsetViewNonConst(const Teuchos::RCP<const map_type>& subMap,
                     const size_t offset);
 
  Teuchos::RCP<const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  getVector(const size_t j) const;
 
  Teuchos::RCP<Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  getVectorNonConst(const size_t j);
 
  Teuchos::ArrayRCP<const Scalar> getData(size_t j) const;
 
  Teuchos::ArrayRCP<Scalar> getDataNonConst(size_t j);
 
  void
  get1dCopy(const Teuchos::ArrayView<Scalar>& A,
            const size_t LDA) const;
 
  void
  get2dCopy(const Teuchos::ArrayView<const Teuchos::ArrayView<Scalar>>& ArrayOfPtrs) const;
 
  Teuchos::ArrayRCP<const Scalar> get1dView() const;
 
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar>> get2dView() const;
 
  Teuchos::ArrayRCP<Scalar> get1dViewNonConst();
 
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar>> get2dViewNonConst();
 
  typename dual_view_type::t_host::const_type getLocalViewHost(Access::ReadOnlyStruct) const;
 
  typename dual_view_type::t_host getLocalViewHost(Access::ReadWriteStruct);
 
  typename dual_view_type::t_host getLocalViewHost(Access::OverwriteAllStruct);
 
  typename dual_view_type::t_dev::const_type getLocalViewDevice(Access::ReadOnlyStruct) const;
 
  typename dual_view_type::t_dev getLocalViewDevice(Access::ReadWriteStruct);
 
  typename dual_view_type::t_dev getLocalViewDevice(Access::OverwriteAllStruct);
 
  wrapped_dual_view_type getWrappedDualView() const;
 
  template <class TargetDeviceType>
  bool need_sync() const {
    return view_.getDualView().template need_sync<TargetDeviceType>();
  }
 
  bool need_sync_host() const;
 
  bool need_sync_device() const;
 
  template <class TargetDeviceType>
  typename std::remove_reference<decltype(std::declval<dual_view_type>().template view<TargetDeviceType>())>::type::const_type
  getLocalView(Access::ReadOnlyStruct s) const {
    return view_.template getView<TargetDeviceType>(s);
  }
 
  template <class TargetDeviceType>
  typename std::remove_reference<decltype(std::declval<dual_view_type>().template view<TargetDeviceType>())>::type
  getLocalView(Access::ReadWriteStruct s) {
    return view_.template getView<TargetDeviceType>(s);
  }
 
  template <class TargetDeviceType>
  typename std::remove_reference<decltype(std::declval<dual_view_type>().template view<TargetDeviceType>())>::type
  getLocalView(Access::OverwriteAllStruct s) {
    return view_.template getView<TargetDeviceType>(s);
  }
 
 
 
  void
  dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
      const Teuchos::ArrayView<dot_type>& dots) const;
 
  template <typename T>
  typename std::enable_if<!(std::is_same<dot_type, T>::value), void>::type
  dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
      const Teuchos::ArrayView<T>& dots) const {
    const size_t sz = static_cast<size_t>(dots.size());
    Teuchos::Array<dot_type> dts(sz);
    this->dot(A, dts);
    for (size_t i = 0; i < sz; ++i) {
      // If T and dot_type differ, this does an implicit conversion.
      dots[i] = dts[i];
    }
  }
 
  template <typename T>
  typename std::enable_if<!(std::is_same<dot_type, T>::value), void>::type
  dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
      std::vector<T>& dots) const {
    const size_t sz = dots.size();
    Teuchos::Array<dot_type> dts(sz);
    this->dot(A, dts);
    for (size_t i = 0; i < sz; ++i) {
      // If T and dot_type differ, this does an implicit conversion.
      dots[i] = dts[i];
    }
  }
 
  void
  dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
      const Kokkos::View<dot_type*, Kokkos::HostSpace>& dots) const;
 
  template <class ViewType>
  void
  dot(typename std::enable_if<std::is_same<typename ViewType::value_type, dot_type>::value &&
                                  std::is_same<typename ViewType::memory_space, typename device_type::memory_space>::value,
                              const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>::type& A,
      const ViewType& dots) const {
    const Kokkos::View<dot_type*, Kokkos::HostSpace> h_dots("Tpetra::Dots", dots.extent(0));
    this->dot(A, h_dots);
    // DEEP_COPY REVIEW - NOT TESTED
    Kokkos::deep_copy(dots, h_dots);
  }
 
  template <typename T>
  typename std::enable_if<!(std::is_same<dot_type, T>::value), void>::type
  dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
      const Kokkos::View<T*, device_type>& dots) const {
    const size_t numDots = dots.extent(0);
    Kokkos::View<dot_type*, device_type> dts("MV::dot tmp", numDots);
    // Call overload that takes a Kokkos::View<dot_type*, device_type>.
    this->dot(A, dts);
    // FIXME (mfh 14 Jul 2014) Does this actually work if dot_type
    // and T differ?  We would need a test for this, but only the
    // Sacado and Stokhos packages are likely to care about this use
    // case.  It could also come up for Kokkos::complex ->
    // std::complex conversions, but those two implementations
    // should generally be bitwise compatible.
    // CT: no this can't possible work .....
    // DEEP_COPY REVIEW - NOT TESTED
    Kokkos::deep_copy(dots, dts);
  }
 
  void abs(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A);
 
  void reciprocal(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A);
 
  void scale(const Scalar& alpha);
 
  void scale(const Teuchos::ArrayView<const Scalar>& alpha);
 
  void scale(const Kokkos::View<const impl_scalar_type*, device_type>& alpha);
 
  void
  scale(const Scalar& alpha,
        const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A);
 
  void
  update(const Scalar& alpha,
         const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
         const Scalar& beta);
 
  void
  update(const Scalar& alpha,
         const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
         const Scalar& beta,
         const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& B,
         const Scalar& gamma);
 
  void
  norm1(const Kokkos::View<mag_type*, Kokkos::HostSpace>& norms) const;
 
  template <class ViewType>
  typename std::enable_if<std::is_same<typename ViewType::value_type, mag_type>::value &&
                          std::is_same<typename ViewType::memory_space, typename device_type::memory_space>::value>::type
  norm1(const ViewType& norms) const {
    // FIXME (mfh 11 Apr 2019) The enable_ifs make it useless for
    // this method to be templated.  (It only exists in case
    // HostSpace = device_type::memory_space.)
    using host_norms_view_type = Kokkos::View<mag_type*, Kokkos::HostSpace>;
    host_norms_view_type h_norms("Tpetra::MV::h_norms", norms.extent(0));
    this->norm1(h_norms);
    // DEEP_COPY REVIEW - HOST-TO-DEVICE
    Kokkos::deep_copy(execution_space(), norms, h_norms);
  }
 
  template <typename T>
  typename std::enable_if<!(std::is_same<mag_type, T>::value), void>::type
  norm1(const Kokkos::View<T*, device_type>& norms) const {
    const size_t numNorms = norms.extent(0);
    Kokkos::View<mag_type*, device_type> tmpNorms("MV::norm1 tmp", numNorms);
    // Call overload that takes a Kokkos::View<mag_type*, device_type>.
    this->norm1(tmpNorms);
    // FIXME (mfh 15 Jul 2014) Does this actually work if mag_type
    // and T differ?  We would need a test for this, but only the
    // Sacado and Stokhos packages are likely to care about this use
    // case.  It could also come up with Kokkos::complex ->
    // std::complex conversion.
    // DEEP_COPY REVIEW - NOT TESTED
    Kokkos::deep_copy(norms, tmpNorms);
  }
 
  void norm1(const Teuchos::ArrayView<mag_type>& norms) const;
 
  template <typename T>
  typename std::enable_if<!(std::is_same<mag_type, T>::value), void>::type
  norm1(const Teuchos::ArrayView<T>& norms) const {
    typedef typename Teuchos::ArrayView<T>::size_type size_type;
    const size_type sz = norms.size();
    Teuchos::Array<mag_type> theNorms(sz);
    this->norm1(theNorms);
    for (size_type i = 0; i < sz; ++i) {
      // If T and mag_type differ, this does an implicit conversion.
      norms[i] = theNorms[i];
    }
  }
 
  void
  norm2(const Kokkos::View<mag_type*, Kokkos::HostSpace>& norms) const;
 
  template <class ViewType>
  typename std::enable_if<std::is_same<typename ViewType::value_type, mag_type>::value &&
                          std::is_same<typename ViewType::memory_space, typename device_type::memory_space>::value>::type
  norm2(const ViewType& norms) const {
    // FIXME (mfh 11 Apr 2019) The enable_ifs make it useless for
    // this method to be templated.  (It only exists in case
    // HostSpace = device_type::memory_space.)
    using host_norms_view_type = Kokkos::View<mag_type*, Kokkos::HostSpace>;
    host_norms_view_type h_norms("Tpetra::MV::h_norms", norms.extent(0));
    this->norm2(h_norms);
    // DEEP_COPY REVIEW - NOT TESTED
    Kokkos::deep_copy(norms, h_norms);
  }
 
  template <typename T>
  typename std::enable_if<!(std::is_same<mag_type, T>::value), void>::type
  norm2(const Kokkos::View<T*, device_type>& norms) const {
    const size_t numNorms = norms.extent(0);
    Kokkos::View<mag_type*, device_type> theNorms("MV::norm2 tmp", numNorms);
    // Call overload that takes a Kokkos::View<mag_type*, device_type>.
    this->norm2(theNorms);
    // FIXME (mfh 14 Jul 2014) Does this actually work if mag_type
    // and T differ?  We would need a test for this, but only the
    // Sacado and Stokhos packages are likely to care about this use
    // case.  This could also come up with Kokkos::complex ->
    // std::complex conversion.
    // DEEP_COPY REVIEW - NOT TESTED
    Kokkos::deep_copy(norms, theNorms);
  }
 
  void norm2(const Teuchos::ArrayView<mag_type>& norms) const;
 
  template <typename T>
  typename std::enable_if<!(std::is_same<mag_type, T>::value), void>::type
  norm2(const Teuchos::ArrayView<T>& norms) const {
    typedef typename Teuchos::ArrayView<T>::size_type size_type;
    const size_type sz = norms.size();
    Teuchos::Array<mag_type> theNorms(sz);
    this->norm2(theNorms);
    for (size_type i = 0; i < sz; ++i) {
      // If T and mag_type differ, this does an implicit conversion.
      norms[i] = theNorms[i];
    }
  }
 
  void normInf(const Kokkos::View<mag_type*, Kokkos::HostSpace>& norms) const;
 
  template <class ViewType>
  typename std::enable_if<std::is_same<typename ViewType::value_type, mag_type>::value &&
                          std::is_same<typename ViewType::memory_space, typename device_type::memory_space>::value>::type
  normInf(const ViewType& norms) const {
    // FIXME (mfh 11 Apr 2019) The enable_ifs make it useless for
    // this method to be templated.  (It only exists in case
    // HostSpace = device_type::memory_space.)
    using host_norms_view_type = Kokkos::View<mag_type*, Kokkos::HostSpace>;
    host_norms_view_type h_norms("Tpetra::MV::h_norms", norms.extent(0));
    this->normInf(h_norms);
    // DEEP_COPY REVIEW - HOST-TO-DEVICE
    Kokkos::deep_copy(execution_space(), norms, h_norms);
  }
 
  template <typename T>
  typename std::enable_if<!(std::is_same<mag_type, T>::value), void>::type
  normInf(const Kokkos::View<T*, device_type>& norms) const {
    const size_t numNorms = norms.extent(0);
    Kokkos::View<mag_type*, device_type> theNorms("MV::normInf tmp", numNorms);
    // Call overload that takes a Kokkos::View<mag_type*, device_type>.
    this->normInf(theNorms);
    // FIXME (mfh 15 Jul 2014) Does this actually work if mag_type
    // and T differ?  We would need a test for this, but only the
    // Sacado and Stokhos packages are likely to care about this use
    // case.  This could also come up with Kokkos::complex ->
    // std::complex conversion.
    // DEEP_COPY REVIEW - NOT TESTED
    Kokkos::deep_copy(norms, theNorms);
  }
 
  void normInf(const Teuchos::ArrayView<mag_type>& norms) const;
 
  template <typename T>
  typename std::enable_if<!(std::is_same<mag_type, T>::value), void>::type
  normInf(const Teuchos::ArrayView<T>& norms) const {
    typedef typename Teuchos::ArrayView<T>::size_type size_type;
    const size_type sz = norms.size();
    Teuchos::Array<mag_type> theNorms(sz);
    this->norm2(theNorms);
    for (size_type i = 0; i < sz; ++i) {
      // If T and mag_type differ, this does an implicit conversion.
      norms[i] = theNorms[i];
    }
  }
 
  void meanValue(const Teuchos::ArrayView<impl_scalar_type>& means) const;
 
  template <typename T>
  typename std::enable_if<!std::is_same<impl_scalar_type, T>::value, void>::type
  meanValue(const Teuchos::ArrayView<T>& means) const {
    typedef typename Teuchos::Array<T>::size_type size_type;
    const size_type numMeans = means.size();
 
    Teuchos::Array<impl_scalar_type> theMeans(numMeans);
    this->meanValue(theMeans());
    for (size_type k = 0; k < numMeans; ++k) {
      means[k] = static_cast<T>(theMeans[k]);
    }
  }
 
  void
  multiply(Teuchos::ETransp transA,
           Teuchos::ETransp transB,
           const Scalar& alpha,
           const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
           const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& B,
           const Scalar& beta);
 
  void
  elementWiseMultiply(Scalar scalarAB,
                      const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
                      const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& B,
                      Scalar scalarThis);
 
 
  size_t getNumVectors() const;
 
  size_t getLocalLength() const;
 
  global_size_t getGlobalLength() const;
 
  size_t getStride() const;
 
  bool isConstantStride() const;
 
  bool aliases(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& other) const;
 
 
 
 
  virtual std::string description() const override;
 
  virtual void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const override;
 
  virtual void
  removeEmptyProcessesInPlace(const Teuchos::RCP<const map_type>& newMap) override;
 
  void setCopyOrView(const Teuchos::DataAccess copyOrView) {
    TEUCHOS_TEST_FOR_EXCEPTION(
        copyOrView == Teuchos::Copy, std::invalid_argument,
        "Tpetra::MultiVector::setCopyOrView: The Kokkos refactor version of "
        "MultiVector _only_ implements view semantics.  You may not call this "
        "method with copyOrView = Teuchos::Copy.  The only valid argument is "
        "Teuchos::View.");
  }
 
  // This method ONLY exists for the circa 2014 "Kokkos refactor"
  // effort.  It ALWAYS returns Teuchos::View.
  Teuchos::DataAccess getCopyOrView() const {
    return Teuchos::View;
  }
 
  void
  assign(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& src);
 
  template <class T>
  Teuchos::RCP<MultiVector<T, LocalOrdinal, GlobalOrdinal, Node>>
  convert() const;
 
  // \brief Checks to see if the local length, number of vectors and size of Scalar type match
  bool isSameSize(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) const;
 
 private:
  using base_type = DistObject<scalar_type, local_ordinal_type,
                               global_ordinal_type, node_type>;
 
 protected:
  template <class DS, class DL, class DG, class DN,
            class SS, class SL, class SG, class SN>
  friend void ::Tpetra::deep_copy(MultiVector<DS, DL, DG, DN>& dst,
                                  const MultiVector<SS, SL, SG, SN>& src);
 
  mutable wrapped_dual_view_type view_;
 
  Teuchos::Array<size_t> whichVectors_;
 
  template <class SC, class LO, class GO, class NT>
  friend ::Teuchos::ArrayView<const size_t> getMultiVectorWhichVectors(const ::Tpetra::MultiVector<SC, LO, GO, NT>& X);
 
 
 
  std::string
  descriptionImpl(const std::string& className) const;
 
  std::string
  localDescribeToString(const Teuchos::EVerbosityLevel vl) const;
 
  void
  describeImpl(Teuchos::FancyOStream& out,
               const std::string& className,
               const Teuchos::EVerbosityLevel verbLevel =
                   Teuchos::Describable::verbLevel_default) const;
 
  // Return true if and only if VectorIndex is a valid column index.
  bool vectorIndexOutOfRange(const size_t VectorIndex) const;
 
  template <class T>
  Teuchos::ArrayRCP<T>
  getSubArrayRCP(Teuchos::ArrayRCP<T> arr, size_t j) const;
 
  size_t getOrigNumLocalRows() const;
 
  size_t getOrigNumLocalCols() const;
 
 
 
  using buffer_device_type =
      typename DistObject<scalar_type,
                          local_ordinal_type,
                          global_ordinal_type,
                          node_type>::buffer_device_type;
 
  virtual bool
  checkSizes(const SrcDistObject& sourceObj) override;
 
  virtual size_t constantNumberOfPackets() const override;
 
  virtual void copyAndPermute(
      const SrcDistObject& sourceObj, const size_t numSameIDs,
      const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs,
      const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs,
      const CombineMode CM, const execution_space& space) override;
 
  virtual void copyAndPermute(
      const SrcDistObject& sourceObj, const size_t numSameIDs,
      const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs,
      const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs,
      const CombineMode CM) override;
 
  virtual void
  packAndPrepare(const SrcDistObject& sourceObj,
                 const Kokkos::DualView<
                     const local_ordinal_type*,
                     buffer_device_type>& exportLIDs,
                 Kokkos::DualView<
                     impl_scalar_type*,
                     buffer_device_type>& exports,
                 Kokkos::DualView<
                     size_t*,
                     buffer_device_type> /* numPacketsPerLID */,
                 size_t& constantNumPackets,
                 const execution_space& space) override;
 
  virtual void
  packAndPrepare(const SrcDistObject& sourceObj,
                 const Kokkos::DualView<
                     const local_ordinal_type*,
                     buffer_device_type>& exportLIDs,
                 Kokkos::DualView<
                     impl_scalar_type*,
                     buffer_device_type>& exports,
                 Kokkos::DualView<
                     size_t*,
                     buffer_device_type> /* numPacketsPerLID */,
                 size_t& constantNumPackets) override;
 
  virtual void
  unpackAndCombine(const Kokkos::DualView<
                       const local_ordinal_type*,
                       buffer_device_type>& importLIDs,
                   Kokkos::DualView<
                       impl_scalar_type*,
                       buffer_device_type>
                       imports,
                   Kokkos::DualView<
                       size_t*,
                       buffer_device_type> /* numPacketsPerLID */,
                   const size_t constantNumPackets,
                   const CombineMode CM,
                   const execution_space& space) override;
 
  virtual void
  unpackAndCombine(const Kokkos::DualView<
                       const local_ordinal_type*,
                       buffer_device_type>& importLIDs,
                   Kokkos::DualView<
                       impl_scalar_type*,
                       buffer_device_type>
                       imports,
                   Kokkos::DualView<
                       size_t*,
                       buffer_device_type> /* numPacketsPerLID */,
                   const size_t constantNumPackets,
                   const CombineMode CM) override;
 
 private:
  // If comm buffers can be aliased to the data view, use this
  // implementation.
  template <class NO = Node>
  typename std::enable_if<std::is_same<typename Tpetra::Details::DefaultTypes::CommBufferMemorySpace<typename NO::execution_space>::type,
                                       typename NO::device_type::memory_space>::value,
                          bool>::type
  reallocImportsIfNeededImpl(const size_t newSize,
                             const bool verbose,
                             const std::string* prefix,
                             const bool areRemoteLIDsContiguous,
                             const CombineMode CM);
 
  // If comm buffers cannot be aliased to the data view, use this
  // implementation. (Just calls DistObject::reallocImportsIfNeeded.)
  template <class NO = Node>
  typename std::enable_if<!std::is_same<typename Tpetra::Details::DefaultTypes::CommBufferMemorySpace<typename NO::execution_space>::type,
                                        typename NO::device_type::memory_space>::value,
                          bool>::type
  reallocImportsIfNeededImpl(const size_t newSize,
                             const bool verbose,
                             const std::string* prefix,
                             const bool areRemoteLIDsContiguous,
                             const CombineMode CM);
 
 protected:
  virtual bool
  reallocImportsIfNeeded(const size_t newSize,
                         const bool verbose,
                         const std::string* prefix,
                         const bool areRemoteLIDsContiguous = false,
                         const CombineMode CM               = INSERT) override;
 
 public:
  bool importsAreAliased();
 
 protected:
  Kokkos::DualView<impl_scalar_type*, buffer_device_type> unaliased_imports_;
 
};  // class MultiVector
 
template <class SC, class LO, class GO, class NT>
Teuchos::ArrayView<const size_t>
getMultiVectorWhichVectors(const MultiVector<SC, LO, GO, NT>& X) {
  return X.whichVectors_();
}
 
template <class ST, class LO, class GO, class NT>
void deep_copy(MultiVector<ST, LO, GO, NT>& dst,
               const MultiVector<ST, LO, GO, NT>& src) {
  // NOTE (mfh 11 Sep 2014) We can't implement deep_copy with
  // shallow-copy operator=, because that would invalidate existing
  // views of dst!
  dst.assign(src);
}
 
// Implementation of the most generic version of MultiVector deep_copy.
template <class DS, class DL, class DG, class DN,
          class SS, class SL, class SG, class SN>
void deep_copy(MultiVector<DS, DL, DG, DN>& dst,
               const MultiVector<SS, SL, SG, SN>& src) {
  using ::Tpetra::getMultiVectorWhichVectors;
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      dst.getGlobalLength() != src.getGlobalLength() ||
          dst.getNumVectors() != src.getNumVectors(),
      std::invalid_argument,
      "Tpetra::deep_copy: Global dimensions of the two Tpetra::MultiVector "
      "objects do not match.  src has dimensions ["
          << src.getGlobalLength()
          << "," << src.getNumVectors() << "], and dst has dimensions ["
          << dst.getGlobalLength() << "," << dst.getNumVectors() << "].");
 
  // FIXME (mfh 28 Jul 2014) Don't throw; just set a local error flag.
  TEUCHOS_TEST_FOR_EXCEPTION(
      dst.getLocalLength() != src.getLocalLength(), std::invalid_argument,
      "Tpetra::deep_copy: The local row counts of the two Tpetra::MultiVector "
      "objects do not match.  src has "
          << src.getLocalLength() << " row(s) "
          << " and dst has " << dst.getLocalLength() << " row(s).");
 
  const bool srcMostUpToDateOnDevice = !src.need_sync_device();
 
  if (src.isConstantStride() && dst.isConstantStride()) {
    if (srcMostUpToDateOnDevice) {
      Details::localDeepCopyConstStride(
          dst.getLocalViewDevice(Access::OverwriteAll),
          src.getLocalViewDevice(Access::ReadOnly));
    } else {
      Details::localDeepCopyConstStride(
          dst.getLocalViewDevice(Access::OverwriteAll),
          src.getLocalViewHost(Access::ReadOnly));
    }
  } else {
    auto dstWhichVecs = getMultiVectorWhichVectors(dst);
    auto srcWhichVecs = getMultiVectorWhichVectors(src);
 
    if (srcMostUpToDateOnDevice) {
      Details::localDeepCopy(dst.getLocalViewDevice(Access::OverwriteAll),
                             src.getLocalViewDevice(Access::ReadOnly),
                             dst.isConstantStride(),
                             src.isConstantStride(),
                             dstWhichVecs,
                             srcWhichVecs);
    } else {
      Details::localDeepCopy(dst.getLocalViewDevice(Access::OverwriteAll),
                             src.getLocalViewHost(Access::ReadOnly),
                             dst.isConstantStride(),
                             src.isConstantStride(),
                             dstWhichVecs,
                             srcWhichVecs);
    }
  }
}
}  // namespace Tpetra
 
namespace Teuchos {
 
// Give Teuchos::TypeNameTraits<Tpetra::MultiVector<...> > a
// human-readable definition.
template <class SC, class LO, class GO, class NT>
class TypeNameTraits<Tpetra::MultiVector<SC, LO, GO, NT>> {
 public:
  static std::string name() {
    return std::string("Tpetra::MultiVector<") +
           TypeNameTraits<SC>::name() + "," +
           TypeNameTraits<LO>::name() + "," +
           TypeNameTraits<GO>::name() + "," +
           TypeNameTraits<NT>::name() + ">";
  }
 
  static std::string
  concreteName(const Tpetra::MultiVector<SC, LO, GO, NT>&) {
    return name();
  }
};
}  // namespace Teuchos
 
#endif  // TPETRA_MULTIVECTOR_DECL_HPP