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

// clang-format off
#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 "Kokkos_ArithTraits.hpp"
#include "Kokkos_InnerProductSpaceTraits.hpp"
#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 Kokkos::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;

    using dot_type =
      typename Kokkos::Details::InnerProductSpaceTraits<impl_scalar_type>::dot_type;

    using mag_type = typename Kokkos::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>& norms) 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;

  // clang-format on
  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;
  // clang-format off

    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