Tpetra_CrsMatrix_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_CRSMATRIX_DECL_HPP
#define TPETRA_CRSMATRIX_DECL_HPP

 
#include "Tpetra_CrsMatrix_fwd.hpp"
#include "TpetraExt_MatrixMatrix_fwd.hpp"
#include "KokkosSparse_Utils.hpp"
#include "KokkosSparse_CrsMatrix.hpp"
#include "Tpetra_Details_MatrixApplyHelper.hpp"
#include "Tpetra_RowMatrix_decl.hpp"
#include "Tpetra_Exceptions.hpp"
#include "Tpetra_DistObject.hpp"
#include "Tpetra_CrsGraph.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_Details_PackTraits.hpp"  // unused here, could delete
#include "Tpetra_Details_ExecutionSpacesUser.hpp"
#include "Teuchos_DataAccess.hpp"
 
#include <memory>  // std::shared_ptr
 
namespace Tpetra {
 
// Forward declaration for CrsMatrix::swap() test
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
class crsMatrix_Swap_Tester;

template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Import<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& importer,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap  = Teuchos::null,
                               const Teuchos::RCP<Teuchos::ParameterList>& params                           = Teuchos::null);

template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Import<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& rowImporter,
                               const Import<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& domainImporter,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap,
                               const Teuchos::RCP<Teuchos::ParameterList>& params);

template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Export<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& exporter,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap  = Teuchos::null,
                               const Teuchos::RCP<Teuchos::ParameterList>& params                           = Teuchos::null);

template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Export<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& rowExporter,
                               const Export<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& domainExporter,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap,
                               const Teuchos::RCP<Teuchos::ParameterList>& params);

namespace Details {
template <class SC, class LO, class GO, class NO>
void residual(const Operator<SC, LO, GO, NO>& A,
              const MultiVector<SC, LO, GO, NO>& X,
              const MultiVector<SC, LO, GO, NO>& B,
              MultiVector<SC, LO, GO, NO>& R);
}

template <class Scalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node>
class CrsMatrix : public RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>,
                  public DistObject<char, LocalOrdinal, GlobalOrdinal, Node>,
                  public Details::Spaces::User {
 private:
  using dist_object_type =
      DistObject<char, LocalOrdinal, GlobalOrdinal,
                 Node>;  
 
 public:


  using scalar_type = Scalar;
  using local_ordinal_type = LocalOrdinal;
  using global_ordinal_type = GlobalOrdinal;
  using device_type = typename Node::device_type;
  using execution_space = typename device_type::execution_space;
  using memory_space = typename device_type::memory_space;

  using node_type = Node;

  using map_type = Map<LocalOrdinal, GlobalOrdinal, Node>;

  using import_type = Import<LocalOrdinal, GlobalOrdinal, Node>;

  using export_type = Export<LocalOrdinal, GlobalOrdinal, Node>;

  using row_matrix_type = RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;

  using impl_scalar_type = typename row_matrix_type::impl_scalar_type;
  using mag_type = typename KokkosKernels::ArithTraits<impl_scalar_type>::mag_type;

  using crs_graph_type = CrsGraph<LocalOrdinal, GlobalOrdinal, Node>;

  using local_graph_device_type = typename crs_graph_type::local_graph_device_type;
  using local_graph_host_type   = typename crs_graph_type::local_graph_host_type;

  using local_matrix_device_type =
      KokkosSparse::CrsMatrix<impl_scalar_type,
                              local_ordinal_type,
                              device_type,
                              void,
                              typename local_graph_device_type::size_type>;
  using local_matrix_host_type =
      typename local_matrix_device_type::host_mirror_type;
 
  using row_ptrs_device_view_type =
      typename row_matrix_type::row_ptrs_device_view_type;
  using row_ptrs_host_view_type =
      typename row_matrix_type::row_ptrs_host_view_type;
 
  using local_inds_device_view_type =
      typename row_matrix_type::local_inds_device_view_type;
  using local_inds_host_view_type =
      typename row_matrix_type::local_inds_host_view_type;
  using nonconst_local_inds_host_view_type =
      typename row_matrix_type::nonconst_local_inds_host_view_type;
 
  using global_inds_device_view_type =
      typename row_matrix_type::global_inds_device_view_type;
  using global_inds_host_view_type =
      typename row_matrix_type::global_inds_host_view_type;
  using nonconst_global_inds_host_view_type =
      typename row_matrix_type::nonconst_global_inds_host_view_type;
 
  using values_device_view_type =
      typename row_matrix_type::values_device_view_type;
  using values_host_view_type =
      typename row_matrix_type::values_host_view_type;
  using nonconst_values_host_view_type =
      typename row_matrix_type::nonconst_values_host_view_type;



  CrsMatrix(const CrsMatrix<Scalar, LocalOrdinal,
                            GlobalOrdinal, Node>&) = default;

  CrsMatrix(CrsMatrix<Scalar, LocalOrdinal,
                      GlobalOrdinal, Node>&&) = default;

  CrsMatrix&
  operator=(const CrsMatrix<Scalar, LocalOrdinal,
                            GlobalOrdinal, Node>&) = default;

  CrsMatrix&
  operator=(CrsMatrix<Scalar, LocalOrdinal,
                      GlobalOrdinal, Node>&&) = default;

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const size_t maxNumEntriesPerRow,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::ArrayView<const size_t>& numEntPerRowToAlloc,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const size_t maxNumEntPerRow,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const Teuchos::ArrayView<const size_t>& numEntPerRowToAlloc,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  explicit CrsMatrix(CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& matrix,
                     const Teuchos::RCP<const crs_graph_type>& graph,
                     const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  explicit CrsMatrix(const Teuchos::RCP<const crs_graph_type>& graph,
                     const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  explicit CrsMatrix(const Teuchos::RCP<const crs_graph_type>& graph,
                     const typename local_matrix_device_type::values_type& values,
                     const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const typename local_graph_device_type::row_map_type& rowPointers,
            const typename local_graph_device_type::entries_type::non_const_type& columnIndices,
            const typename local_matrix_device_type::values_type& values,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const Teuchos::ArrayRCP<size_t>& rowPointers,
            const Teuchos::ArrayRCP<LocalOrdinal>& columnIndices,
            const Teuchos::ArrayRCP<Scalar>& values,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const local_matrix_device_type& lclMatrix,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const local_matrix_device_type& lclMatrix,
            const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const Teuchos::RCP<const map_type>& domainMap      = Teuchos::null,
            const Teuchos::RCP<const map_type>& rangeMap       = Teuchos::null,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  CrsMatrix(const local_matrix_device_type& lclMatrix,
            const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const Teuchos::RCP<const map_type>& domainMap,
            const Teuchos::RCP<const map_type>& rangeMap,
            const Teuchos::RCP<const import_type>& importer,
            const Teuchos::RCP<const export_type>& exporter,
            const Teuchos::RCP<Teuchos::ParameterList>& params =
                Teuchos::null);

  CrsMatrix(const Teuchos::RCP<const map_type>& rowMap,
            const Teuchos::RCP<const map_type>& colMap,
            const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& source,
            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  // This function in 'Copy' mode is only guaranteed to work correctly for matrices
  // which are fillComplete.
  CrsMatrix(const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& source,
            const Teuchos::DataAccess copyOrView);

  virtual ~CrsMatrix() = default;
 
  // This friend declaration makes the clone() method work.
  template <class S2, class LO2, class GO2, class N2>
  friend class CrsMatrix;
 
  // This friend declaration allows for fused residual calculation
  template <class S2, class LO2, class GO2, class N2>
  friend void Details::residual(const Operator<S2, LO2, GO2, N2>& A,
                                const MultiVector<S2, LO2, GO2, N2>& X,
                                const MultiVector<S2, LO2, GO2, N2>& B,
                                MultiVector<S2, LO2, GO2, N2>& R);
 
  // This friend declaration allows for batching of apply calls
  template <class MatrixArray, class MultiVectorArray>
  friend void batchedApply(const MatrixArray& Matrices,
                           const typename std::remove_pointer<typename MultiVectorArray::value_type>::type& X,
                           MultiVectorArray& Y,
                           typename std::remove_pointer<typename MatrixArray::value_type>::type::scalar_type alpha,
                           typename std::remove_pointer<typename MatrixArray::value_type>::type::scalar_type beta,
                           Teuchos::RCP<Teuchos::ParameterList> params);
 
 public:


  //
  void
  insertGlobalValues(const GlobalOrdinal globalRow,
                     const Teuchos::ArrayView<const GlobalOrdinal>& cols,
                     const Teuchos::ArrayView<const Scalar>& vals);

  void
  insertGlobalValues(const GlobalOrdinal globalRow,
                     const LocalOrdinal numEnt,
                     const Scalar vals[],
                     const GlobalOrdinal inds[]);

  void
  insertLocalValues(const LocalOrdinal localRow,
                    const Teuchos::ArrayView<const LocalOrdinal>& cols,
                    const Teuchos::ArrayView<const Scalar>& vals,
                    const CombineMode CM = ADD);

  void
  insertLocalValues(const LocalOrdinal localRow,
                    const LocalOrdinal numEnt,
                    const Scalar vals[],
                    const LocalOrdinal cols[],
                    const CombineMode CM = ADD);
 
 protected:
  virtual LocalOrdinal
  replaceGlobalValuesImpl(impl_scalar_type rowVals[],
                          const crs_graph_type& graph,
                          const RowInfo& rowInfo,
                          const GlobalOrdinal inds[],
                          const impl_scalar_type newVals[],
                          const LocalOrdinal numElts);
 
 public:
  local_ordinal_type
  replaceGlobalValues(
      const global_ordinal_type globalRow,
      const Kokkos::View<const global_ordinal_type*, Kokkos::AnonymousSpace>& inputInds,
      const Kokkos::View<const impl_scalar_type*, Kokkos::AnonymousSpace>& inputVals);

  LocalOrdinal
  replaceGlobalValues(const GlobalOrdinal globalRow,
                      const Teuchos::ArrayView<const GlobalOrdinal>& cols,
                      const Teuchos::ArrayView<const Scalar>& vals);

  LocalOrdinal
  replaceGlobalValues(const GlobalOrdinal globalRow,
                      const LocalOrdinal numEnt,
                      const Scalar vals[],
                      const GlobalOrdinal cols[]);
 
 protected:
  virtual LocalOrdinal
  replaceLocalValuesImpl(impl_scalar_type rowVals[],
                         const crs_graph_type& graph,
                         const RowInfo& rowInfo,
                         const LocalOrdinal inds[],
                         const impl_scalar_type newVals[],
                         const LocalOrdinal numElts);
 
 public:
  local_ordinal_type
  replaceLocalValues(
      const local_ordinal_type localRow,
      const Kokkos::View<const local_ordinal_type*, Kokkos::AnonymousSpace>& inputInds,
      const Kokkos::View<const impl_scalar_type*, Kokkos::AnonymousSpace>& inputVals);

  LocalOrdinal
  replaceLocalValues(const LocalOrdinal localRow,
                     const Teuchos::ArrayView<const LocalOrdinal>& cols,
                     const Teuchos::ArrayView<const Scalar>& vals);

  LocalOrdinal
  replaceLocalValues(const LocalOrdinal localRow,
                     const LocalOrdinal numEnt,
                     const Scalar inputVals[],
                     const LocalOrdinal inputCols[]);
 
 private:
  static const bool useAtomicUpdatesByDefault =
#ifdef KOKKOS_ENABLE_SERIAL
      !std::is_same<execution_space, Kokkos::Serial>::value;
#else
      true;
#endif  // KOKKOS_ENABLE_SERIAL

 protected:
  virtual LocalOrdinal
  sumIntoGlobalValuesImpl(impl_scalar_type rowVals[],
                          const crs_graph_type& graph,
                          const RowInfo& rowInfo,
                          const GlobalOrdinal inds[],
                          const impl_scalar_type newVals[],
                          const LocalOrdinal numElts,
                          const bool atomic = useAtomicUpdatesByDefault);
 
 public:
  LocalOrdinal
  sumIntoGlobalValues(const GlobalOrdinal globalRow,
                      const Teuchos::ArrayView<const GlobalOrdinal>& cols,
                      const Teuchos::ArrayView<const Scalar>& vals,
                      const bool atomic = useAtomicUpdatesByDefault);

  LocalOrdinal
  sumIntoGlobalValues(const GlobalOrdinal globalRow,
                      const LocalOrdinal numEnt,
                      const Scalar vals[],
                      const GlobalOrdinal cols[],
                      const bool atomic = useAtomicUpdatesByDefault);
 
 protected:
  virtual LocalOrdinal
  sumIntoLocalValuesImpl(impl_scalar_type rowVals[],
                         const crs_graph_type& graph,
                         const RowInfo& rowInfo,
                         const LocalOrdinal inds[],
                         const impl_scalar_type newVals[],
                         const LocalOrdinal numElts,
                         const bool atomic = useAtomicUpdatesByDefault);
 
 public:
  local_ordinal_type
  sumIntoLocalValues(
      const local_ordinal_type localRow,
      const Kokkos::View<const local_ordinal_type*, Kokkos::AnonymousSpace>& inputInds,
      const Kokkos::View<const impl_scalar_type*, Kokkos::AnonymousSpace>& inputVals,
      const bool atomic = useAtomicUpdatesByDefault);

  LocalOrdinal
  sumIntoLocalValues(const LocalOrdinal localRow,
                     const Teuchos::ArrayView<const LocalOrdinal>& cols,
                     const Teuchos::ArrayView<const Scalar>& vals,
                     const bool atomic = useAtomicUpdatesByDefault);

  LocalOrdinal
  sumIntoLocalValues(const LocalOrdinal localRow,
                     const LocalOrdinal numEnt,
                     const Scalar vals[],
                     const LocalOrdinal cols[],
                     const bool atomic = useAtomicUpdatesByDefault);
 
 private:
  LocalOrdinal
  transformLocalValues(impl_scalar_type rowVals[],
                       const crs_graph_type& graph,
                       const RowInfo& rowInfo,
                       const LocalOrdinal inds[],
                       const impl_scalar_type newVals[],
                       const LocalOrdinal numElts,
                       std::function<impl_scalar_type(const impl_scalar_type&, const impl_scalar_type&)> f,
                       const bool atomic = useAtomicUpdatesByDefault);

  LocalOrdinal
  transformGlobalValues(impl_scalar_type rowVals[],
                        const crs_graph_type& graph,
                        const RowInfo& rowInfo,
                        const GlobalOrdinal inds[],
                        const impl_scalar_type newVals[],
                        const LocalOrdinal numElts,
                        std::function<impl_scalar_type(const impl_scalar_type&, const impl_scalar_type&)> f,
                        const bool atomic = useAtomicUpdatesByDefault);

  LocalOrdinal
  transformLocalValues(const LocalOrdinal lclRow,
                       const LocalOrdinal numInputEnt,
                       const impl_scalar_type inputVals[],
                       const LocalOrdinal inputCols[],
                       std::function<impl_scalar_type(const impl_scalar_type&, const impl_scalar_type&)> f,
                       const bool atomic = useAtomicUpdatesByDefault);

  LocalOrdinal
  transformGlobalValues(const GlobalOrdinal gblRow,
                        const LocalOrdinal numInputEnt,
                        const impl_scalar_type inputVals[],
                        const GlobalOrdinal inputCols[],
                        std::function<impl_scalar_type(const impl_scalar_type&, const impl_scalar_type&)> f,
                        const bool atomic = useAtomicUpdatesByDefault);
 
 public:
  template <class LocalIndicesViewType,
            class ImplScalarViewType,
            class BinaryFunction>
  LocalOrdinal
  transformLocalValues(const LocalOrdinal lclRow,
                       const typename UnmanagedView<LocalIndicesViewType>::type& inputInds,
                       const typename UnmanagedView<ImplScalarViewType>::type& inputVals,
                       BinaryFunction f,
                       const bool atomic = useAtomicUpdatesByDefault) {
    // We use static_assert here to check the template parameters,
    // rather than std::enable_if (e.g., on the return value, to
    // enable compilation only if the template parameters match the
    // desired attributes).  This turns obscure link errors into
    // clear compilation errors.  It also makes the return value a
    // lot easier to see.
    static_assert(Kokkos::is_view<LocalIndicesViewType>::value,
                  "First template parameter LocalIndicesViewType must be "
                  "a Kokkos::View.");
    static_assert(Kokkos::is_view<ImplScalarViewType>::value,
                  "Second template parameter ImplScalarViewType must be a "
                  "Kokkos::View.");
    static_assert(static_cast<int>(LocalIndicesViewType::rank) == 1,
                  "First template parameter LocalIndicesViewType must "
                  "have rank 1.");
    static_assert(static_cast<int>(ImplScalarViewType::rank) == 1,
                  "Second template parameter ImplScalarViewType must have "
                  "rank 1.");
    static_assert(std::is_same<
                      typename LocalIndicesViewType::non_const_value_type,
                      local_ordinal_type>::value,
                  "First template parameter LocalIndicesViewType must "
                  "contain values of type local_ordinal_type.");
    static_assert(std::is_same<
                      typename ImplScalarViewType::non_const_value_type,
                      impl_scalar_type>::value,
                  "Second template parameter ImplScalarViewType must "
                  "contain values of type impl_scalar_type.");
    typedef LocalOrdinal LO;
    const LO numInputEnt = inputInds.extent(0);
    if (static_cast<LO>(inputVals.extent(0)) != numInputEnt) {
      return Teuchos::OrdinalTraits<LO>::invalid();
    }
    return this->transformLocalValues(lclRow,
                                      numInputEnt,
                                      inputVals.data(),
                                      inputInds.data(),
                                      f,
                                      atomic);
  }

  template <class BinaryFunction, class InputMemorySpace>
  LocalOrdinal
  transformGlobalValues(const GlobalOrdinal gblRow,
                        const Kokkos::View<const GlobalOrdinal*,
                                           InputMemorySpace,
                                           Kokkos::MemoryUnmanaged>& inputInds,
                        const Kokkos::View<const impl_scalar_type*,
                                           InputMemorySpace,
                                           Kokkos::MemoryUnmanaged>& inputVals,
                        BinaryFunction f,
                        const bool atomic = useAtomicUpdatesByDefault) {
    typedef LocalOrdinal LO;
    const LO numInputEnt = inputInds.extent(0);
    if (static_cast<LO>(inputVals.extent(0)) != numInputEnt) {
      return Teuchos::OrdinalTraits<LO>::invalid();
    }
    return this->transformGlobalValues(gblRow,
                                       numInputEnt,
                                       inputVals.data(),
                                       inputInds.data(),
                                       f,
                                       atomic);
  }

  void setAllToScalar(const Scalar& alpha);

  void scale(const Scalar& alpha);

  void
  setAllValues(const typename local_graph_device_type::row_map_type& ptr,
               const typename local_graph_device_type::entries_type::non_const_type& ind,
               const typename local_matrix_device_type::values_type& val);

  void
  setAllValues(const local_matrix_device_type& localMatrix);

  void
  setAllValues(const Teuchos::ArrayRCP<size_t>& ptr,
               const Teuchos::ArrayRCP<LocalOrdinal>& ind,
               const Teuchos::ArrayRCP<Scalar>& val);

  row_ptrs_host_view_type getLocalRowPtrsHost() const { return getCrsGraph()->getLocalRowPtrsHost(); }

  row_ptrs_device_view_type getLocalRowPtrsDevice() const { return getCrsGraph()->getLocalRowPtrsDevice(); }

  local_inds_host_view_type getLocalIndicesHost() const { return getCrsGraph()->getLocalIndicesHost(); }

  local_inds_device_view_type getLocalIndicesDevice() const { return getCrsGraph()->getLocalIndicesDevice(); }



  void globalAssemble();

  void resumeFill(const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  void
  fillComplete(const Teuchos::RCP<const map_type>& domainMap,
               const Teuchos::RCP<const map_type>& rangeMap,
               const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  void
  fillComplete(const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  void
  expertStaticFillComplete(const Teuchos::RCP<const map_type>& domainMap,
                           const Teuchos::RCP<const map_type>& rangeMap,
                           const Teuchos::RCP<const import_type>& importer    = Teuchos::null,
                           const Teuchos::RCP<const export_type>& exporter    = Teuchos::null,
                           const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

  void
  replaceColMap(const Teuchos::RCP<const map_type>& newColMap);

  void
  reindexColumns(crs_graph_type* const graph,
                 const Teuchos::RCP<const map_type>& newColMap,
                 const Teuchos::RCP<const import_type>& newImport = Teuchos::null,
                 const bool sortEachRow                           = true);

  void
  replaceDomainMap(const Teuchos::RCP<const map_type>& newDomainMap);

  void
  replaceDomainMapAndImporter(const Teuchos::RCP<const map_type>& newDomainMap,
                              Teuchos::RCP<const import_type>& newImporter);

  void
  replaceRangeMap(const Teuchos::RCP<const map_type>& newRangeMap);

  void
  replaceRangeMapAndExporter(const Teuchos::RCP<const map_type>& newRangeMap,
                             Teuchos::RCP<const export_type>& newExporter);

  virtual void
  removeEmptyProcessesInPlace(const Teuchos::RCP<const map_type>& newMap) override;



  Teuchos::RCP<const Teuchos::Comm<int> > getComm() const override;

  Teuchos::RCP<const map_type> getRowMap() const override;

  Teuchos::RCP<const map_type> getColMap() const override;

  Teuchos::RCP<const RowGraph<LocalOrdinal, GlobalOrdinal, Node> >
  getGraph() const override;

  Teuchos::RCP<const crs_graph_type> getCrsGraph() const;
 
 private:
  const crs_graph_type& getCrsGraphRef() const;
 
 public:
#if __armclang_major__ == 22 && __armclang_minor__ == 1
  // On Stria, PR 13052 caused a 25% performance regression in the
  // CGSolve performance test that is fixed by forcing
  // getLocalMatrixDevice to always be inlined. Restrict the fix
  // to the specific toolchain where the problem was observed
#define TPETRA_DETAILS_ALWAYS_INLINE __attribute__((always_inline))
#else
#define TPETRA_DETAILS_ALWAYS_INLINE
#endif
  TPETRA_DETAILS_ALWAYS_INLINE local_matrix_device_type
  getLocalMatrixDevice() const;
  local_matrix_host_type getLocalMatrixHost() const;
#undef TPETRA_DETAILS_ALWAYS_INLINE

  global_size_t getGlobalNumRows() const override;

  global_size_t getGlobalNumCols() const override;

  size_t getLocalNumRows() const override;

  size_t getLocalNumCols() const override;

  GlobalOrdinal getIndexBase() const override;

  global_size_t getGlobalNumEntries() const override;

  size_t getLocalNumEntries() const override;

  size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const override;

  size_t getNumEntriesInLocalRow(local_ordinal_type localRow) const override;

  size_t getGlobalMaxNumRowEntries() const override;

  size_t getLocalMaxNumRowEntries() const override;

  virtual LocalOrdinal getBlockSize() const override { return 1; }

  bool hasColMap() const override;

  bool isLocallyIndexed() const override;

  bool isGloballyIndexed() const override;

  bool isFillComplete() const override;

  bool isFillActive() const;


  bool isStorageOptimized() const;

  bool isStaticGraph() const;

  mag_type getNormInf() const;

  //
  mag_type getNorm1(bool assumeSymmetric = false) const;

  mag_type getFrobeniusNorm() const override;

  virtual bool supportsRowViews() const override;
 
 protected:
  using values_dualv_type =
      Kokkos::DualView<impl_scalar_type*, device_type>;
  using values_wdv_type =
      Details::WrappedDualView<values_dualv_type>;
  values_wdv_type valuesUnpacked_wdv;
  mutable values_wdv_type valuesPacked_wdv;
 
 public:
  void
  getGlobalRowCopy(GlobalOrdinal GlobalRow,
                   nonconst_global_inds_host_view_type& Indices,
                   nonconst_values_host_view_type& Values,
                   size_t& NumEntries) const override;
  void
  getLocalRowCopy(LocalOrdinal LocalRow,
                  nonconst_local_inds_host_view_type& Indices,
                  nonconst_values_host_view_type& Values,
                  size_t& NumEntries) const override;

 
  void
  getGlobalRowView(GlobalOrdinal GlobalRow,
                   global_inds_host_view_type& indices,
                   values_host_view_type& values) const override;

  void
  getLocalRowView(LocalOrdinal LocalRow,
                  local_inds_host_view_type& indices,
                  values_host_view_type& values) const override;


  void
  getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& diag) const override;

  void getLocalDiagOffsets(Teuchos::ArrayRCP<size_t>& offsets) const;

  void
  getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& diag,
                   const Kokkos::View<const size_t*, device_type,
                                      Kokkos::MemoryUnmanaged>& offsets) const;

  void
  getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& diag,
                   const Teuchos::ArrayView<const size_t>& offsets) const;

  void
  leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) override;

  void
  rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) override;



  void
  localApply(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& X,
             MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Y,
             const Teuchos::ETransp mode = Teuchos::NO_TRANS,
             const Scalar& alpha         = Teuchos::ScalarTraits<Scalar>::one(),
             const Scalar& beta          = Teuchos::ScalarTraits<Scalar>::zero()) const;

  template <class T>
  Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node> >
  convert() const;



  void
  apply(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& X,
        MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Y,
        Teuchos::ETransp mode = Teuchos::NO_TRANS,
        Scalar alpha          = Teuchos::ScalarTraits<Scalar>::one(),
        Scalar beta           = Teuchos::ScalarTraits<Scalar>::zero()) const override;

  bool hasTransposeApply() const override;

  Teuchos::RCP<const map_type> getDomainMap() const override;

  Teuchos::RCP<const map_type> getRangeMap() const override;



  virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
  add(const Scalar& alpha,
      const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
      const Scalar& beta,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& domainMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& rangeMap,
      const Teuchos::RCP<Teuchos::ParameterList>& params) const override;



  std::string description() const override;

  void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const override;



  typedef typename DistObject<Scalar, LocalOrdinal, GlobalOrdinal,
                              Node>::buffer_device_type buffer_device_type;
 
  virtual bool
  checkSizes(const SrcDistObject& source) override;
 
  void
  applyCrsPadding(
      const typename crs_graph_type::padding_type& padding,
      const bool verbose);
 
 private:
  void
  copyAndPermuteStaticGraph(
      const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& source,
      const size_t numSameIDs,
      const LocalOrdinal permuteToLIDs[],
      const LocalOrdinal permuteFromLIDs[],
      const size_t numPermutes);
 
  void
  copyAndPermuteNonStaticGraph(
      const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& source,
      const size_t numSameIDs,
      const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs_dv,
      const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs_dv,
      const size_t numPermutes);
 
 protected:
  using dist_object_type::
      copyAndPermute;  
 
  virtual void
  copyAndPermute(const SrcDistObject& source,
                 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& source,
                 const Kokkos::DualView<
                     const local_ordinal_type*,
                     buffer_device_type>& exportLIDs,
                 Kokkos::DualView<char*, buffer_device_type>& exports,
                 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
                 size_t& constantNumPackets) override;
 
  using dist_object_type::packAndPrepare;  
 
 private:
  void
  unpackAndCombineImpl(
      const Kokkos::DualView<const local_ordinal_type*,
                             buffer_device_type>& importLIDs,
      Kokkos::DualView<char*, buffer_device_type> imports,
      Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
      const size_t constantNumPackets,
      const CombineMode combineMode,
      const bool verbose);

  void
  unpackAndCombineImplNonStatic(
      const Kokkos::DualView<const local_ordinal_type*,
                             buffer_device_type>& importLIDs,
      Kokkos::DualView<char*, buffer_device_type> imports,
      Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
      const size_t constantNumPackets,
      const CombineMode combineMode);
 
 public:
  void
  unpackAndCombine(const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& importLIDs,
                   Kokkos::DualView<char*, buffer_device_type> imports,
                   Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
                   const size_t constantNumPackets,
                   const CombineMode CM) override;
 
  using dist_object_type::unpackAndCombine;  

  void
  packNew(const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
          Kokkos::DualView<char*, buffer_device_type>& exports,
          const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
          size_t& constantNumPackets) const;
 
 private:
  void
  packNonStaticNew(const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
                   Kokkos::DualView<char*, buffer_device_type>& exports,
                   const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
                   size_t& constantNumPackets) const;

  size_t
  packRow(char exports[],
          const size_t offset,
          const size_t numEnt,
          const GlobalOrdinal gidsIn[],
          const impl_scalar_type valsIn[],
          const size_t numBytesPerValue) const;

  bool
  packRowStatic(char* const numEntOut,
                char* const valOut,
                char* const indOut,
                const size_t numEnt,
                const LocalOrdinal lclRow) const;

  size_t
  unpackRow(GlobalOrdinal gidsOut[],
            impl_scalar_type valsOut[],
            const char imports[],
            const size_t offset,
            const size_t numBytes,
            const size_t numEnt,
            const size_t numBytesPerValue);

  void
  allocatePackSpaceNew(Kokkos::DualView<char*, buffer_device_type>& exports,
                       size_t& totalNumEntries,
                       const Kokkos::DualView<const local_ordinal_type*,
                                              buffer_device_type>& exportLIDs) const;
 
 public:
  typename local_matrix_host_type::values_type::const_type
  getLocalValuesHost(Access::ReadOnlyStruct s) const {
    return valuesPacked_wdv.getHostView(s);
  }

  typename local_matrix_host_type::values_type
  getLocalValuesHost(Access::ReadWriteStruct s) {
    return valuesPacked_wdv.getHostView(s);
  }

  typename local_matrix_host_type::values_type
  getLocalValuesHost(Access::OverwriteAllStruct s) {
    return valuesPacked_wdv.getHostView(s);
  }

  typename local_matrix_device_type::values_type::const_type
  getLocalValuesDevice(Access::ReadOnlyStruct s) const {
    return valuesPacked_wdv.getDeviceView(s);
  }

  typename local_matrix_device_type::values_type
  getLocalValuesDevice(Access::ReadWriteStruct s) {
    return valuesPacked_wdv.getDeviceView(s);
  }

  typename local_matrix_device_type::values_type
  getLocalValuesDevice(Access::OverwriteAllStruct s) {
    return valuesPacked_wdv.getDeviceView(s);
  }
 
 private:
  // Friend declaration for nonmember function.
  template <class CrsMatrixType>
  friend Teuchos::RCP<CrsMatrixType>
  Tpetra::importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                         const Import<typename CrsMatrixType::local_ordinal_type,
                                                      typename CrsMatrixType::global_ordinal_type,
                                                      typename CrsMatrixType::node_type>& importer,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& domainMap,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& rangeMap,
                                         const Teuchos::RCP<Teuchos::ParameterList>& params);
 
  // Friend declaration for nonmember function.
  template <class CrsMatrixType>
  friend Teuchos::RCP<CrsMatrixType>
  Tpetra::importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                         const Import<typename CrsMatrixType::local_ordinal_type,
                                                      typename CrsMatrixType::global_ordinal_type,
                                                      typename CrsMatrixType::node_type>& rowImporter,
                                         const Import<typename CrsMatrixType::local_ordinal_type,
                                                      typename CrsMatrixType::global_ordinal_type,
                                                      typename CrsMatrixType::node_type>& domainImporter,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& domainMap,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& rangeMap,
                                         const Teuchos::RCP<Teuchos::ParameterList>& params);
 
  // Friend declaration for nonmember function.
  template <class CrsMatrixType>
  friend Teuchos::RCP<CrsMatrixType>
  Tpetra::exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                         const Export<typename CrsMatrixType::local_ordinal_type,
                                                      typename CrsMatrixType::global_ordinal_type,
                                                      typename CrsMatrixType::node_type>& exporter,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& domainMap,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& rangeMap,
                                         const Teuchos::RCP<Teuchos::ParameterList>& params);
 
  // Friend declaration for nonmember function.
  template <class CrsMatrixType>
  friend Teuchos::RCP<CrsMatrixType>
  Tpetra::exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                         const Export<typename CrsMatrixType::local_ordinal_type,
                                                      typename CrsMatrixType::global_ordinal_type,
                                                      typename CrsMatrixType::node_type>& rowExporter,
                                         const Export<typename CrsMatrixType::local_ordinal_type,
                                                      typename CrsMatrixType::global_ordinal_type,
                                                      typename CrsMatrixType::node_type>& domainExporter,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& domainMap,
                                         const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                      typename CrsMatrixType::global_ordinal_type,
                                                                      typename CrsMatrixType::node_type> >& rangeMap,
                                         const Teuchos::RCP<Teuchos::ParameterList>& params);
 
 public:
  void
  importAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& destMatrix,
                        const import_type& importer,
                        const Teuchos::RCP<const map_type>& domainMap,
                        const Teuchos::RCP<const map_type>& rangeMap,
                        const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;

  void
  importAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& destMatrix,
                        const import_type& rowImporter,
                        const import_type& domainImporter,
                        const Teuchos::RCP<const map_type>& domainMap,
                        const Teuchos::RCP<const map_type>& rangeMap,
                        const Teuchos::RCP<Teuchos::ParameterList>& params) const;

  void
  exportAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& destMatrix,
                        const export_type& exporter,
                        const Teuchos::RCP<const map_type>& domainMap      = Teuchos::null,
                        const Teuchos::RCP<const map_type>& rangeMap       = Teuchos::null,
                        const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;

  void
  exportAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& destMatrix,
                        const export_type& rowExporter,
                        const export_type& domainExporter,
                        const Teuchos::RCP<const map_type>& domainMap,
                        const Teuchos::RCP<const map_type>& rangeMap,
                        const Teuchos::RCP<Teuchos::ParameterList>& params) const;
 
 private:
  void
  transferAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& destMatrix,
                          const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer,
                          const Teuchos::RCP<const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node> >& domainTransfer,
                          const Teuchos::RCP<const map_type>& domainMap      = Teuchos::null,
                          const Teuchos::RCP<const map_type>& rangeMap       = Teuchos::null,
                          const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;

 protected:
  virtual void
  insertGlobalValuesImpl(crs_graph_type& graph,
                         RowInfo& rowInfo,
                         const GlobalOrdinal gblColInds[],
                         const impl_scalar_type vals[],
                         const size_t numInputEnt);
 
 private:
  void
  insertGlobalValuesFiltered(
      const GlobalOrdinal globalRow,
      const Teuchos::ArrayView<const GlobalOrdinal>& indices,
      const Teuchos::ArrayView<const Scalar>& values,
      const bool debug);

  void
  insertGlobalValuesFilteredChecked(
      const GlobalOrdinal globalRow,
      const Teuchos::ArrayView<const GlobalOrdinal>& indices,
      const Teuchos::ArrayView<const Scalar>& values,
      const char* const prefix,
      const bool debug,
      const bool verbose);

  void
  combineGlobalValues(
      const GlobalOrdinal globalRowIndex,
      const Teuchos::ArrayView<const GlobalOrdinal>& columnIndices,
      const Teuchos::ArrayView<const Scalar>& values,
      const Tpetra::CombineMode combineMode,
      const char* const prefix,
      const bool debug,
      const bool verbose);

  LocalOrdinal
  combineGlobalValuesRaw(const LocalOrdinal lclRow,
                         const LocalOrdinal numEnt,
                         const impl_scalar_type vals[],
                         const GlobalOrdinal cols[],
                         const Tpetra::CombineMode combineMode,
                         const char* const prefix,
                         const bool debug,
                         const bool verbose);

  template <class BinaryFunction>
  LocalOrdinal
  transformGlobalValues(const GlobalOrdinal globalRow,
                        const Teuchos::ArrayView<const GlobalOrdinal>& indices,
                        const Teuchos::ArrayView<const Scalar>& values,
                        BinaryFunction f,
                        const bool atomic = useAtomicUpdatesByDefault) {
    typedef impl_scalar_type IST;
    typedef LocalOrdinal LO;
    typedef GlobalOrdinal GO;
 
    const LO numInputEnt = static_cast<LO>(indices.size());
    if (static_cast<LO>(values.size()) != numInputEnt) {
      return Teuchos::OrdinalTraits<LO>::invalid();
    }
 
    const GO* const inputCols = indices.getRawPtr();
    const IST* const inputVals =
        reinterpret_cast<const IST*>(values.getRawPtr());
    return this->transformGlobalValues(globalRow, numInputEnt, inputVals,
                                       inputCols, f, atomic);
  }

  void
  insertNonownedGlobalValues(const GlobalOrdinal globalRow,
                             const Teuchos::ArrayView<const GlobalOrdinal>& indices,
                             const Teuchos::ArrayView<const Scalar>& values);

  void
  insertIndicesAndValues(crs_graph_type& graph,
                         RowInfo& rowInfo,
                         const typename crs_graph_type::SLocalGlobalViews& newInds,
                         const Teuchos::ArrayView<impl_scalar_type>& oldRowVals,
                         const Teuchos::ArrayView<const impl_scalar_type>& newRowVals,
                         const ELocalGlobal lg,
                         const ELocalGlobal I);
 
 protected:
  // useful typedefs
  typedef Teuchos::OrdinalTraits<LocalOrdinal> OTL;
  typedef KokkosKernels::ArithTraits<impl_scalar_type> STS;
  typedef KokkosKernels::ArithTraits<mag_type> STM;
  typedef MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> MV;
  typedef Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> V;
  typedef crs_graph_type Graph;
 
  // Enums
  enum GraphAllocationStatus {
    GraphAlreadyAllocated,
    GraphNotYetAllocated
  };
 
 protected:
  void allocateValues(ELocalGlobal lg, GraphAllocationStatus gas,
                      const bool verbose);

  void
  sortAndMergeIndicesAndValues(const bool sorted,
                               const bool merged);
 
 public:
  bool haveGlobalConstants() const;
 
 protected:
  mutable Teuchos::RCP<MV> importMV_;

  mutable Teuchos::RCP<MV> exportMV_;

  Teuchos::RCP<MV>
  getColumnMapMultiVector(const MV& X_domainMap,
                          const bool force = false) const;

  Teuchos::RCP<MV>
  getRowMapMultiVector(const MV& Y_rangeMap,
                       const bool force = false) const;

  void
  applyNonTranspose(const MV& X_in,
                    MV& Y_in,
                    Scalar alpha,
                    Scalar beta) const;

  void
  applyTranspose(const MV& X_in,
                 MV& Y_in,
                 const Teuchos::ETransp mode,
                 Scalar alpha,
                 Scalar beta) const;
 
  // matrix data accessors

  typename values_dualv_type::t_host::const_type
  getValuesViewHost(const RowInfo& rowinfo) const;

  typename values_dualv_type::t_dev::const_type
  getValuesViewDevice(const RowInfo& rowinfo) const;

  typename values_dualv_type::t_host
  getValuesViewHostNonConst(const RowInfo& rowinfo);

  typename values_dualv_type::t_dev
  getValuesViewDeviceNonConst(const RowInfo& rowinfo);
 
 private:
  // TODO: When KokkosKernels 4.4 is released, local_matrix_device_type can be permanently modified to use the default_size_type
  // of KK. This is always a type that is enabled by KK's ETI (preferring int if both or neither int and size_t are enabled).
  //
  // At that point the ApplyHelper can be replaced with just a SPMVHandle.
  using local_matrix_int_rowptrs_device_type =
      KokkosSparse::CrsMatrix<impl_scalar_type,
                              local_ordinal_type,
                              device_type,
                              void,
                              int>;

  using ApplyHelper = Details::MatrixApplyHelper<
      local_matrix_device_type,
      local_matrix_int_rowptrs_device_type,
      typename MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::device_view_type>;
 
  std::shared_ptr<ApplyHelper> getApplyHelper() const {
    if (!applyHelper) {
      auto A_lcl  = getLocalMatrixDevice();
      applyHelper = std::make_shared<ApplyHelper>(A_lcl.nnz(), A_lcl.graph.row_map);
    }
    return applyHelper;
  }
 
 protected:
  // Friend the tester for CrsMatrix::swap
  friend class Tpetra::crsMatrix_Swap_Tester<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
 
  // Friend the matrix multiply kernels so they can access internally-cached integer
  // row pointers without making them part of the CrsMatrix interface
  template <typename S, typename LO, typename GO, typename NODE, typename LOV>
  friend struct Tpetra::MMdetails::KernelWrappers;
  template <typename S, typename LO, typename GO, typename NODE, typename LOV>
  friend struct Tpetra::MMdetails::KernelWrappers2;
  // Really we want to friend:
  // template <typename S, typename LO, typename GO, typename NODE, typename LOV>
  // friend void Tpetra::MMdetails::kokkos_kernels_mult_A_B_newmatrix(...),
  // but nvcc seemingly doesn't reliably honor this friendship (because it's a template class?)
  // Instead, we defined a little non-template wrapper class over by kokkos_kernels_mult_A_B_newmatrix so we
  // can friend that here instead, which seems to work.
  // This is needed so kokkos_kernels_mult_A_B_newmatrix can access the lazy integer row pointer
  // helpers
  friend struct Tpetra::MMdetails::CrsMatrixApplyHelperAccess;
 
  // friend Matrix Matrix utility function that needs to access integer-typed rowptrs
  friend void Tpetra::MMdetails::import_and_extract_views<Scalar, LocalOrdinal, GlobalOrdinal, Node>(
      const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
      Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > targetMap,
      CrsMatrixStruct<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Aview,
      Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > prototypeImporter,
      bool userAssertsThereAreNoRemotes,
      const std::string& label,
      const Teuchos::RCP<Teuchos::ParameterList>& params);

  void swap(CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& matrix);
 
 protected:
  void fillLocalMatrix(const Teuchos::RCP<Teuchos::ParameterList>& params);

  void fillLocalGraphAndMatrix(const Teuchos::RCP<Teuchos::ParameterList>& params);

  void checkInternalState() const;


  Teuchos::RCP<const Graph> staticGraph_;
  Teuchos::RCP<Graph> myGraph_;
 
 protected:
  Details::EStorageStatus storageStatus_ =
      Details::STORAGE_1D_UNPACKED;

  bool fillComplete_ = false;

  std::map<GlobalOrdinal, std::pair<Teuchos::Array<GlobalOrdinal>,
                                    Teuchos::Array<Scalar> > >
      nonlocals_;
 
 private:
  mutable std::shared_ptr<ApplyHelper> applyHelper;
 
 public:
  // FIXME (mfh 24 Feb 2014) Is it _really_ necessary to make this a
  // public inner class of CrsMatrix?  It looks like it doesn't
  // depend on any implementation details of CrsMatrix at all.  It
  // should really be declared and defined outside of CrsMatrix.
  template <class DestViewType, class SrcViewType,
            class DestOffsetViewType, class SrcOffsetViewType>
  struct pack_functor {
    typedef typename DestViewType::execution_space execution_space;
    SrcViewType src_;
    DestViewType dst_;
    SrcOffsetViewType src_offset_;
    DestOffsetViewType dst_offset_;
    typedef typename DestOffsetViewType::non_const_value_type scalar_index_type;
 
    pack_functor(DestViewType dst,
                 const SrcViewType src,
                 DestOffsetViewType dst_offset,
                 const SrcOffsetViewType src_offset)
      : src_(src)
      , dst_(dst)
      , src_offset_(src_offset)
      , dst_offset_(dst_offset) {}
 
    KOKKOS_INLINE_FUNCTION
    void operator()(const LocalOrdinal row) const {
      scalar_index_type srcPos       = src_offset_(row);
      const scalar_index_type dstEnd = dst_offset_(row + 1);
      scalar_index_type dstPos       = dst_offset_(row);
      for (; dstPos < dstEnd; ++dstPos, ++srcPos) {
        dst_(dstPos) = src_(srcPos);
      }
    }
  };
};  // class CrsMatrix

template <class Scalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node>
Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
createCrsMatrix(
    const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& map,
    const size_t maxNumEntriesPerRow                   = 0,
    const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) {
  using matrix_type =
      CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
  return Teuchos::rcp(new matrix_type(map, maxNumEntriesPerRow,
                                      params));
}
 
template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Import<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& importer,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap,
                               const Teuchos::RCP<Teuchos::ParameterList>& params) {
  Teuchos::RCP<CrsMatrixType> destMatrix;
  sourceMatrix->importAndFillComplete(destMatrix, importer, domainMap, rangeMap, params);
  return destMatrix;
}
 
template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Import<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& rowImporter,
                               const Import<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& domainImporter,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap,
                               const Teuchos::RCP<Teuchos::ParameterList>& params) {
  Teuchos::RCP<CrsMatrixType> destMatrix;
  sourceMatrix->importAndFillComplete(destMatrix, rowImporter, domainImporter, domainMap, rangeMap, params);
  return destMatrix;
}
 
template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Export<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& exporter,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap,
                               const Teuchos::RCP<Teuchos::ParameterList>& params) {
  Teuchos::RCP<CrsMatrixType> destMatrix;
  sourceMatrix->exportAndFillComplete(destMatrix, exporter, domainMap, rangeMap, params);
  return destMatrix;
}
 
template <class CrsMatrixType>
Teuchos::RCP<CrsMatrixType>
exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                               const Export<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& rowExporter,
                               const Export<typename CrsMatrixType::local_ordinal_type,
                                            typename CrsMatrixType::global_ordinal_type,
                                            typename CrsMatrixType::node_type>& domainExporter,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& domainMap,
                               const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                            typename CrsMatrixType::global_ordinal_type,
                                                            typename CrsMatrixType::node_type> >& rangeMap,
                               const Teuchos::RCP<Teuchos::ParameterList>& params) {
  Teuchos::RCP<CrsMatrixType> destMatrix;
  sourceMatrix->exportAndFillComplete(destMatrix, rowExporter, domainExporter, domainMap, rangeMap, params);
  return destMatrix;
}

template <class CrsMatrixType>
void removeCrsMatrixZeros(CrsMatrixType& matrix,
                          typename Teuchos::ScalarTraits<typename CrsMatrixType::scalar_type>::magnitudeType const& threshold =
                              Teuchos::ScalarTraits<typename CrsMatrixType::scalar_type>::magnitude(Teuchos::ScalarTraits<typename CrsMatrixType::scalar_type>::zero())) {
  auto localMatrix       = matrix.getLocalMatrixDevice();
  size_t nnzBefore       = localMatrix.nnz();
  localMatrix            = KokkosSparse::removeCrsMatrixZeros(localMatrix, threshold);
  size_t localNNZRemoved = nnzBefore - localMatrix.nnz();
  // Skip the expertStaticFillComplete if no entries were removed on any process.
  // The fill complete can perform MPI collectives, so it can only be skipped on all processes or none.
  size_t globalNNZRemoved = 0;
  Teuchos::reduceAll<int, size_t>(*(matrix.getComm()), Teuchos::REDUCE_SUM, 1, &localNNZRemoved, &globalNNZRemoved);
  if (globalNNZRemoved != size_t(0)) {
    matrix.resumeFill();
    matrix.setAllValues(localMatrix);
    matrix.expertStaticFillComplete(matrix.getDomainMap(), matrix.getRangeMap());
  }
}
 
}  // namespace Tpetra

 
#endif  // TPETRA_CRSMATRIX_DECL_HPP