Tpetra_CrsMatrix_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.
//
// ************************************************************************
// @HEADER
 
// clang-format off
#ifndef TPETRA_CRSMATRIX_DECL_HPP
#define TPETRA_CRSMATRIX_DECL_HPP

 
#include "Tpetra_CrsMatrix_fwd.hpp"
#include "Tpetra_LocalCrsMatrixOperator.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 "KokkosSparse_Utils.hpp"
#include "KokkosSparse_CrsMatrix.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
  {
  // clang-format on
private:
  using dist_object_type =
      DistObject<char, LocalOrdinal, GlobalOrdinal,
                 Node>; 
  // clang-format off
 
  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 Kokkos::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::HostMirror;
 

    using local_multiply_op_type =
      LocalCrsMatrixOperator<scalar_type,
                             scalar_type,
                             device_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);

    // 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:
    local_matrix_device_type getLocalMatrixDevice () const;
    local_matrix_host_type getLocalMatrixHost () const;

    std::shared_ptr<local_multiply_op_type> getLocalMultiplyOperator () const;

    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 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;
 
    using ordinal_rowptrs_type = typename local_multiply_op_type::ordinal_view_type;
    mutable ordinal_rowptrs_type ordinalRowptrs;
 
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:
 
  // clang-format on
  using dist_object_type::
      copyAndPermute; 
  // clang-format off
 
    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;
 
  // clang-format on
  using dist_object_type::packAndPrepare; 
                                          // clang-format off
 
  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;
 
  // clang-format on
  using dist_object_type::unpackAndCombine; 
                                            // clang-format off

    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 Kokkos::ArithTraits<impl_scalar_type> STS;
    typedef Kokkos::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);

    static size_t
    mergeRowIndicesAndValues (size_t rowLen, local_ordinal_type* cols, impl_scalar_type* vals);

    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);
 
 
  protected:
 
    // Friend the tester for CrsMatrix::swap
    friend class Tpetra::crsMatrix_Swap_Tester<Scalar, LocalOrdinal, GlobalOrdinal, Node>;

    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_;
 
  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