Tpetra_MultiVector_decl.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
// @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 "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 <type_traits>
 
namespace Tpetra {
 
  namespace Details {
    template<class DstMultiVectorType, class SrcMultiVectorType>
    struct MultiVectorCloner {
      typedef DstMultiVectorType dst_mv_type;
      typedef SrcMultiVectorType src_mv_type;
 
      static Teuchos::RCP<dst_mv_type>
      clone (const src_mv_type& X,
             const Teuchos::RCP<typename dst_mv_type::node_type>& node2);
    };
  } // namespace Details
 
  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);
 
  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);
 
namespace Classes {
  // 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 = ::Tpetra::Details::DefaultTypes::scalar_type,
            class LocalOrdinal = ::Tpetra::Details::DefaultTypes::local_ordinal_type,
            class GlobalOrdinal = ::Tpetra::Details::DefaultTypes::global_ordinal_type,
            class Node = ::Tpetra::Details::DefaultTypes::node_type>
  class MultiVector :
    public DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node>
  {
  public:
 
 
    typedef Scalar scalar_type;
    typedef typename Kokkos::Details::ArithTraits<Scalar>::val_type impl_scalar_type;
    typedef LocalOrdinal local_ordinal_type;
    typedef GlobalOrdinal global_ordinal_type;
    typedef Node node_type;
 
    typedef typename Node::device_type device_type;
 
  private:
    typedef DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> base_type;
 
  public:
    typedef typename Kokkos::Details::InnerProductSpaceTraits<impl_scalar_type>::dot_type dot_type;
 
    typedef typename Kokkos::Details::ArithTraits<impl_scalar_type>::mag_type mag_type;
 
    typedef typename Node::execution_space execution_space;
 
    typedef Kokkos::DualView<impl_scalar_type**, Kokkos::LayoutLeft,
      typename execution_space::execution_space> dual_view_type;
 
    typedef Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
 
 
 
    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);
 
    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);
 
  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 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 map_type& subMap,
                 const size_t offset = 0);
 
    template <class Node2>
    Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node2> >
    clone (const Teuchos::RCP<Node2>& node2) const;
 
    void swap(MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & mv);
 
    virtual ~MultiVector ();
 
 
 
  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) const;
 
    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) const
    {
      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) const;
 
    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) const
    {
      sumIntoGlobalValue (gblRow, col, static_cast<impl_scalar_type> (val), atomic);
    }
 
    void
    replaceLocalValue (const LocalOrdinal lclRow,
                       const size_t col,
                       const impl_scalar_type& value) const;
 
    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) const
    {
      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) const;
 
    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) const
    {
      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();
 
    MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
    operator= (const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& source);
 
 
 
 
    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 ();
 
    dual_view_type getDualView () const;
 
    template<class TargetDeviceType>
    void sync () {
      getDualView ().template sync<TargetDeviceType> ();
    }
 
    template<class TargetDeviceType>
    bool need_sync () const {
      return getDualView ().template need_sync<TargetDeviceType> ();
    }
 
    template<class TargetDeviceType>
    void modify () {
      getDualView ().template modify<TargetDeviceType> ();
    }
 
    template<class TargetDeviceType>
    typename Kokkos::Impl::if_c<
      std::is_same<
        typename device_type::memory_space,
        typename TargetDeviceType::memory_space>::value,
      typename dual_view_type::t_dev,
      typename dual_view_type::t_host>::type
    getLocalView () const {
      return getDualView ().template view<TargetDeviceType> ();
    }
 
 
 
    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);
      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 .....
      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);
 
    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 {
      typedef typename Kokkos::View<mag_type*, Kokkos::HostSpace> host_norms_view_type;
      host_norms_view_type h_norms("Tpetra::MV::h_norms",norms.extent(0));
      this->normImpl (h_norms, NORM_ONE);
      Kokkos::deep_copy(norms,h_norms);
    }
    void norm1 (const Kokkos::View<mag_type*, Kokkos::HostSpace>& norms) const;
 
    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.
      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];
      }
    }
 
    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 {
      typedef typename Kokkos::View<mag_type*, Kokkos::HostSpace> host_norms_view_type;
      host_norms_view_type h_norms("Tpetra::MV::h_norms",norms.extent(0));
      this->normImpl (h_norms, NORM_TWO);
      Kokkos::deep_copy(norms,h_norms);
    }
    void norm2 (const Kokkos::View<mag_type*, Kokkos::HostSpace>& norms) const;
 
    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.
      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];
      }
    }
 
    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 {
      typedef typename Kokkos::View<mag_type*, Kokkos::HostSpace> host_norms_view_type;
      host_norms_view_type h_norms("Tpetra::MV::h_norms",norms.extent(0));
      this->normImpl (h_norms, NORM_INF);
      Kokkos::deep_copy(norms,h_norms);
    }
    void normInf (const Kokkos::View<mag_type*, Kokkos::HostSpace>& norms) const;
 
    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.
      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 TPETRA_DEPRECATED
    normWeighted (const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& weights,
                  const Teuchos::ArrayView<mag_type>& norms) const;
 
    template <typename T>
    typename std::enable_if< ! (std::is_same<mag_type,T>::value), void >::type
    TPETRA_DEPRECATED
    normWeighted (const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& weights,
                  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->normWeighted (weights, 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;
 
 
 
 
    virtual std::string description() const;
 
    virtual void
    describe (Teuchos::FancyOStream& out,
              const Teuchos::EVerbosityLevel verbLevel =
              Teuchos::Describable::verbLevel_default) const;
 
    virtual void
    removeEmptyProcessesInPlace (const Teuchos::RCP<const map_type>& newMap);
 
    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);
 
 
    // \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;
 
  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 dual_view_type view_;
 
    mutable dual_view_type origView_;
 
    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);
 
 
 
    enum EWhichNorm {
      NORM_ONE, //<! Use the one-norm
      NORM_TWO, //<! Use the two-norm
      NORM_INF  //<! Use the infinity-norm
    };
 
    void
    normImpl (const Kokkos::View<mag_type*, Kokkos::HostSpace>& norms,
              const EWhichNorm whichNorm) const;
 
 
 
 
 
 
    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;
 
 
 
    virtual bool
    checkSizes (const SrcDistObject& sourceObj);
 
    virtual size_t constantNumberOfPackets () const;
 
    virtual bool useNewInterface () { return true; }
 
    typedef typename DistObject<Scalar, LocalOrdinal, GlobalOrdinal,
                                Node>::buffer_device_type buffer_device_type;
 
    virtual void
    copyAndPermuteNew (const SrcDistObject& sourceObj,
                       const size_t numSameIDs,
                       const Kokkos::DualView<const local_ordinal_type*, device_type>& permuteToLIDs,
                       const Kokkos::DualView<const local_ordinal_type*, device_type>& permuteFromLIDs);
 
    virtual void
    packAndPrepareNew (const SrcDistObject& sourceObj,
                       const Kokkos::DualView<const local_ordinal_type*, device_type>& exportLIDs,
                       Kokkos::DualView<impl_scalar_type*, buffer_device_type>& exports,
                       const Kokkos::DualView<size_t*, buffer_device_type>& /* numPacketsPerLID */,
                       size_t& constantNumPackets,
                       Distributor& /* distor */);
 
    virtual void
    unpackAndCombineNew (const Kokkos::DualView<const LocalOrdinal*, device_type>& importLIDs,
                         const Kokkos::DualView<const impl_scalar_type*, buffer_device_type>& imports,
                         const Kokkos::DualView<const size_t*, buffer_device_type>& /* numPacketsPerLID */,
                         const size_t constantNumPackets,
                         Distributor& /* distor */,
                         const CombineMode CM);
  }; // 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_ ();
  }
} // namespace Classes
 
  namespace Details {
 
    template<class DstMultiVectorType,
             class SrcMultiVectorType>
    Teuchos::RCP<typename MultiVectorCloner<DstMultiVectorType, SrcMultiVectorType>::dst_mv_type>
    MultiVectorCloner<DstMultiVectorType, SrcMultiVectorType>::
    clone (const src_mv_type& X,
           const Teuchos::RCP<typename dst_mv_type::node_type>& node2)
    {
      using Teuchos::RCP;
      typedef typename src_mv_type::map_type src_map_type;
      typedef typename dst_mv_type::map_type dst_map_type;
      typedef typename dst_mv_type::node_type dst_node_type;
      typedef typename dst_mv_type::dual_view_type dst_dual_view_type;
 
      // Clone X's Map to have the new Node type.
      RCP<const src_map_type> map1 = X.getMap ();
      RCP<const dst_map_type> map2 = map1.is_null () ?
        Teuchos::null : map1->template clone<dst_node_type> (node2);
 
      const size_t lclNumRows = X.getLocalLength ();
      const size_t numCols = X.getNumVectors ();
      dst_dual_view_type Y_view ("MV::dual_view", lclNumRows, numCols);
 
      RCP<dst_mv_type> Y (new dst_mv_type (map2, Y_view));
      // Let deep_copy do the work for us, to avoid code duplication.
      ::Tpetra::deep_copy (*Y, X);
      return Y ;
    }
 
  } // namespace Details
 
  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)
  {
    typedef typename DN::device_type DD;
    //typedef typename SN::device_type SD;
    using ::Tpetra::Classes::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).");
 
    if (src.isConstantStride () && dst.isConstantStride ()) {
      typedef typename MultiVector<DS, DL, DG, DN>::dual_view_type::t_host::execution_space HES;
      typedef typename MultiVector<DS, DL, DG, DN>::dual_view_type::t_dev::execution_space DES;
 
      // If we need sync to device, then host has the most recent version.
      const bool useHostVersion = src.template need_sync<typename SN::device_type> ();
 
      if (! useHostVersion) {
        // Device memory has the most recent version of src.
        dst.template modify<DES> (); // We are about to modify dst on device.
        // Copy from src to dst on device.
        Details::localDeepCopyConstStride (dst.template getLocalView<DES> (),
                                           src.template getLocalView<typename SN::device_type> ());
        dst.template sync<HES> (); // Sync dst from device to host.
      }
      else { // Host memory has the most recent version of src.
        dst.template modify<HES> (); // We are about to modify dst on host.
        // Copy from src to dst on host.
        Details::localDeepCopyConstStride (dst.template getLocalView<Kokkos::HostSpace> (),
                                           src.template getLocalView<Kokkos::HostSpace> ());
        dst.template sync<DES> (); // Sync dst from host to device.
      }
    }
    else {
      typedef Kokkos::DualView<SL*, DD> whichvecs_type;
      typedef typename whichvecs_type::t_dev::execution_space DES;
      typedef typename whichvecs_type::t_host::execution_space HES;
 
      if (dst.isConstantStride ()) {
        const SL numWhichVecs =
          static_cast<SL> (getMultiVectorWhichVectors (src).size ());
        const std::string whichVecsLabel ("MV::deep_copy::whichVecs");
 
        // We can't sync src, since it is only an input argument.
        // Thus, we have to use the most recently modified version of
        // src, which could be either the device or host version.
        //
        // If we need sync to device, then host has the most recent version.
        const bool useHostVersion = src.template need_sync<typename SN::device_type> ();
 
        if (! useHostVersion) { // Copy from the device version of src.
          // whichVecs tells the kernel which vectors (columns) of src
          // to copy.  Fill whichVecs on the host, and sync to device.
          whichvecs_type whichVecs (whichVecsLabel, numWhichVecs);
          whichVecs.template modify<HES> ();
 
          Teuchos::ArrayView<const size_t> src_whichVectors =
            getMultiVectorWhichVectors (src);
          for (SL i = 0; i < numWhichVecs; ++i) {
            whichVecs.h_view(i) = static_cast<SL> (src_whichVectors[i]);
          }
          // Sync the host version of whichVecs to the device.
          whichVecs.template sync<DES> ();
 
          // Mark the device version of dst's DualView as modified.
          dst.template modify<DES> ();
          // Copy from the selected vectors of src to dst, on the device.
          Details::localDeepCopy (dst.template getLocalView<typename DN::device_type> (),
                                  src.template getLocalView<typename SN::device_type> (),
                                  dst.isConstantStride (),
                                  src.isConstantStride (),
                                  whichVecs.d_view,
                                  whichVecs.d_view);
          // Sync dst's DualView to the host.  This is cheaper than
          // repeating the above copy from src to dst on the host.
          dst.template sync<HES> ();
        }
        else { // host version of src was the most recently modified
          // Copy from the host version of src.
          //
          // whichVecs tells the kernel which vectors (columns) of src
          // to copy.  Fill whichVecs on the host, and use it there.
          typedef Kokkos::View<SL*, HES> the_whichvecs_type;
          the_whichvecs_type whichVecs (whichVecsLabel, numWhichVecs);
          Teuchos::ArrayView<const size_t> src_whichVectors =
            getMultiVectorWhichVectors (src);
          for (SL i = 0; i < numWhichVecs; ++i) {
            whichVecs(i) = static_cast<SL> (src_whichVectors[i]);
          }
          // Copy from the selected vectors of src to dst, on the
          // host.  The function ignores the first instance of
          // 'whichVecs' in this case.
          Details::localDeepCopy (dst.template getLocalView<Kokkos::HostSpace> (),
                                  src.template getLocalView<Kokkos::HostSpace> (),
                                  dst.isConstantStride (),
                                  src.isConstantStride (),
                                  whichVecs, whichVecs);
          // Sync dst back to the device, since we only copied on the host.
          dst.template sync<DES> ();
        }
      }
      else { // dst is NOT constant stride
        if (src.isConstantStride ()) {
 
          // If we need sync to device, then host has the most recent version.
          const bool useHostVersion = src.template need_sync<typename SN::device_type> ();
 
          if (! useHostVersion) { // Copy from the device version of src.
            // whichVecs tells the kernel which vectors (columns) of dst
            // to copy.  Fill whichVecs on the host, and sync to device.
            typedef Kokkos::DualView<DL*, DES> the_whichvecs_type;
            const std::string whichVecsLabel ("MV::deep_copy::whichVecs");
            Teuchos::ArrayView<const size_t> dst_whichVectors =
              getMultiVectorWhichVectors (dst);
            const DL numWhichVecs = static_cast<DL> (dst_whichVectors.size ());
            the_whichvecs_type whichVecs (whichVecsLabel, numWhichVecs);
            whichVecs.template modify<HES> ();
            for (DL i = 0; i < numWhichVecs; ++i) {
              whichVecs.h_view(i) = dst_whichVectors[i];
            }
            // Sync the host version of whichVecs to the device.
            whichVecs.template sync<DES> ();
 
            // Copy src to the selected vectors of dst, on the device.
            Details::localDeepCopy (dst.template getLocalView<typename DN::device_type> (),
                                    src.template getLocalView<typename SN::device_type> (),
                                    dst.isConstantStride (),
                                    src.isConstantStride (),
                                    whichVecs.d_view,
                                    whichVecs.d_view);
            // We can't sync src and repeat the above copy on the
            // host, so sync dst back to the host.
            //
            // FIXME (mfh 29 Jul 2014) This may overwrite columns that
            // don't actually belong to dst's view.
            dst.template sync<HES> ();
          }
          else { // host version of src was the most recently modified
            // Copy from the host version of src.
            //
            // whichVecs tells the kernel which vectors (columns) of src
            // to copy.  Fill whichVecs on the host, and use it there.
            typedef Kokkos::View<DL*, HES> the_whichvecs_type;
            Teuchos::ArrayView<const size_t> dst_whichVectors =
              getMultiVectorWhichVectors (dst);
            const DL numWhichVecs = static_cast<DL> (dst_whichVectors.size ());
            the_whichvecs_type whichVecs ("MV::deep_copy::whichVecs",
                                          numWhichVecs);
            for (DL i = 0; i < numWhichVecs; ++i) {
              whichVecs(i) = static_cast<DL> (dst_whichVectors[i]);
            }
            // Copy from src to the selected vectors of dst, on the host.
            Details::localDeepCopy (dst.template getLocalView<Kokkos::HostSpace> (),
                                    src.template getLocalView<Kokkos::HostSpace> (),
                                    dst.isConstantStride (),
                                    src.isConstantStride (),
                                    whichVecs, whichVecs);
            // Sync dst back to the device, since we only copied on the host.
            //
            // FIXME (mfh 29 Jul 2014) This may overwrite columns that
            // don't actually belong to dst's view.
            dst.template sync<DES> ();
          }
        }
        else { // neither src nor dst have constant stride
 
          // If we need sync to device, then host has the most recent version.
          const bool useHostVersion = src.template need_sync<typename SN::device_type> ();
 
          if (! useHostVersion) { // Copy from the device version of src.
            // whichVectorsDst tells the kernel which columns of dst
            // to copy.  Fill it on the host, and sync to device.
            Teuchos::ArrayView<const size_t> dst_whichVectors =
              getMultiVectorWhichVectors (dst);
            const DL dstNumWhichVecs =
              static_cast<DL> (dst_whichVectors.size ());
            Kokkos::DualView<DL*, DES> whichVecsDst ("MV::deep_copy::whichVecsDst",
                                                     dstNumWhichVecs);
            whichVecsDst.template modify<HES> ();
            for (DL i = 0; i < dstNumWhichVecs; ++i) {
              whichVecsDst.h_view(i) = static_cast<DL> (dst_whichVectors[i]);
            }
            // Sync the host version of whichVecsDst to the device.
            whichVecsDst.template sync<DES> ();
 
            // whichVectorsSrc tells the kernel which vectors
            // (columns) of src to copy.  Fill it on the host, and
            // sync to device.  Use the destination MultiVector's
            // LocalOrdinal type here.
            Teuchos::ArrayView<const size_t> src_whichVectors =
              getMultiVectorWhichVectors (src);
            const DL srcNumWhichVecs =
              static_cast<DL> (src_whichVectors.size ());
            Kokkos::DualView<DL*, DES> whichVecsSrc ("MV::deep_copy::whichVecsSrc",
                                                     srcNumWhichVecs);
            whichVecsSrc.template modify<HES> ();
            for (DL i = 0; i < srcNumWhichVecs; ++i) {
              whichVecsSrc.h_view(i) = static_cast<DL> (src_whichVectors[i]);
            }
            // Sync the host version of whichVecsSrc to the device.
            whichVecsSrc.template sync<DES> ();
 
            // Copy from the selected vectors of src to the selected
            // vectors of dst, on the device.
            Details::localDeepCopy (dst.template getLocalView<typename DN::device_type> (),
                                    src.template getLocalView<typename SN::device_type> (),
                                    dst.isConstantStride (),
                                    src.isConstantStride (),
                                    whichVecsDst.d_view,
                                    whichVecsSrc.d_view);
          }
          else {
            Teuchos::ArrayView<const size_t> dst_whichVectors =
              getMultiVectorWhichVectors (dst);
            const DL dstNumWhichVecs =
              static_cast<DL> (dst_whichVectors.size ());
            Kokkos::View<DL*, HES> whichVectorsDst ("dstWhichVecs",
                                                    dstNumWhichVecs);
            for (DL i = 0; i < dstNumWhichVecs; ++i) {
              whichVectorsDst(i) = dst_whichVectors[i];
            }
 
            Teuchos::ArrayView<const size_t> src_whichVectors =
              getMultiVectorWhichVectors (src);
            // Use the destination MultiVector's LocalOrdinal type here.
            const DL srcNumWhichVecs =
              static_cast<DL> (src_whichVectors.size ());
            Kokkos::View<DL*, HES> whichVectorsSrc ("srcWhichVecs",
                                                    srcNumWhichVecs);
            for (DL i = 0; i < srcNumWhichVecs; ++i) {
              whichVectorsSrc(i) = src_whichVectors[i];
            }
 
            // Copy from the selected vectors of src to the selected
            // vectors of dst, on the host.
            Details::localDeepCopy (dst.template getLocalView<Kokkos::HostSpace> (),
                                    src.template getLocalView<Kokkos::HostSpace> (),
                                    dst.isConstantStride (),
                                    src.isConstantStride (),
                                    whichVectorsDst, whichVectorsSrc);
            // We can't sync src and repeat the above copy on the
            // host, so sync dst back to the host.
            //
            // FIXME (mfh 29 Jul 2014) This may overwrite columns that
            // don't actually belong to dst's view.
            dst.template sync<HES> ();
          }
        }
      }
    }
  }
} // 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