Tpetra_Details_packCrsMatrix_def.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_DETAILS_PACKCRSMATRIX_DEF_HPP
#define TPETRA_DETAILS_PACKCRSMATRIX_DEF_HPP

#include "TpetraCore_config.h"
#include "Teuchos_Array.hpp"
#include "Teuchos_ArrayView.hpp"
#include "Tpetra_Details_Behavior.hpp"
#include "Tpetra_Details_castAwayConstDualView.hpp"
#include "Tpetra_Details_createMirrorView.hpp"
#include "Tpetra_Details_getEntryOnHost.hpp"
#include "Tpetra_Details_OrdinalTraits.hpp"
#include "Tpetra_Details_PackTraits.hpp"
#include "Tpetra_Details_Profiling.hpp"
#include "Tpetra_CrsMatrix_decl.hpp"
#include <memory>
#include <sstream>
#include <stdexcept>
#include <string>


namespace Tpetra {

//
// Users must never rely on anything in the Details namespace.
//
namespace Details {

namespace PackCrsMatrixImpl {
template<class OutputOffsetsViewType,
         class CountsViewType,
         class InputOffsetsViewType,
         class InputLocalRowIndicesViewType,
         class InputLocalRowPidsViewType,
         const bool debug =
#ifdef HAVE_TPETRA_DEBUG
         true
#else
         false
#endif // HAVE_TPETRA_DEBUG
         >
class NumPacketsAndOffsetsFunctor {
public:
  typedef typename OutputOffsetsViewType::non_const_value_type output_offset_type;
  typedef typename CountsViewType::non_const_value_type count_type;
  typedef typename InputOffsetsViewType::non_const_value_type input_offset_type;
  typedef typename InputLocalRowIndicesViewType::non_const_value_type local_row_index_type;
  typedef typename InputLocalRowPidsViewType::non_const_value_type local_row_pid_type;
  // output Views drive where execution happens.
  typedef typename OutputOffsetsViewType::device_type device_type;
  static_assert (std::is_same<typename CountsViewType::device_type::execution_space,
                   typename device_type::execution_space>::value,
                 "OutputOffsetsViewType and CountsViewType must have the same execution space.");
  static_assert (Kokkos::is_view<OutputOffsetsViewType>::value,
                 "OutputOffsetsViewType must be a Kokkos::View.");
  static_assert (std::is_same<typename OutputOffsetsViewType::value_type, output_offset_type>::value,
                 "OutputOffsetsViewType must be a nonconst Kokkos::View.");
  static_assert (std::is_integral<output_offset_type>::value,
                 "The type of each entry of OutputOffsetsViewType must be a built-in integer type.");
  static_assert (Kokkos::is_view<CountsViewType>::value,
                 "CountsViewType must be a Kokkos::View.");
  static_assert (std::is_same<typename CountsViewType::value_type, output_offset_type>::value,
                 "CountsViewType must be a nonconst Kokkos::View.");
  static_assert (std::is_integral<count_type>::value,
                 "The type of each entry of CountsViewType must be a built-in integer type.");
  static_assert (Kokkos::is_view<InputOffsetsViewType>::value,
                 "InputOffsetsViewType must be a Kokkos::View.");
  static_assert (std::is_integral<input_offset_type>::value,
                 "The type of each entry of InputOffsetsViewType must be a built-in integer type.");
  static_assert (Kokkos::is_view<InputLocalRowIndicesViewType>::value,
                 "InputLocalRowIndicesViewType must be a Kokkos::View.");
  static_assert (std::is_integral<local_row_index_type>::value,
                 "The type of each entry of InputLocalRowIndicesViewType must be a built-in integer type.");

  NumPacketsAndOffsetsFunctor (const OutputOffsetsViewType& outputOffsets,
                               const CountsViewType& counts,
                               const InputOffsetsViewType& rowOffsets,
                               const InputLocalRowIndicesViewType& lclRowInds,
                               const InputLocalRowPidsViewType& lclRowPids,
                               const count_type sizeOfLclCount,
                               const count_type sizeOfGblColInd,
                               const count_type sizeOfPid,
                               const count_type sizeOfValue) :
    outputOffsets_ (outputOffsets),
    counts_ (counts),
    rowOffsets_ (rowOffsets),
    lclRowInds_ (lclRowInds),
    lclRowPids_ (lclRowPids),
    sizeOfLclCount_ (sizeOfLclCount),
    sizeOfGblColInd_ (sizeOfGblColInd),
    sizeOfPid_ (sizeOfPid),
    sizeOfValue_ (sizeOfValue),
    error_ ("error") // don't forget this, or you'll get segfaults!
  {
    if (debug) {
      const size_t numRowsToPack = static_cast<size_t> (lclRowInds_.extent (0));

      if (numRowsToPack != static_cast<size_t> (counts_.extent (0))) {
        std::ostringstream os;
        os << "lclRowInds.extent(0) = " << numRowsToPack
           << " != counts.extent(0) = " << counts_.extent (0)
           << ".";
        TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str ());
      }
      if (static_cast<size_t> (numRowsToPack + 1) !=
          static_cast<size_t> (outputOffsets_.extent (0))) {
        std::ostringstream os;
        os << "lclRowInds.extent(0) + 1 = " << (numRowsToPack + 1)
           << " != outputOffsets.extent(0) = " << outputOffsets_.extent (0)
           << ".";
        TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str ());
      }
    }
  }

  KOKKOS_INLINE_FUNCTION void
  operator() (const local_row_index_type& curInd,
              output_offset_type& update,
              const bool final) const
  {
    if (debug) {
      if (curInd < static_cast<local_row_index_type> (0)) {
        error_ () = 1;
        return;
      }
    }

    if (final) {
      if (debug) {
        if (curInd >= static_cast<local_row_index_type> (outputOffsets_.extent (0))) {
          error_ () = 2;
          return;
        }
      }
      outputOffsets_(curInd) = update;
    }

    if (curInd < static_cast<local_row_index_type> (counts_.extent (0))) {
      const auto lclRow = lclRowInds_(curInd);
      if (static_cast<size_t> (lclRow + 1) >= static_cast<size_t> (rowOffsets_.extent (0)) ||
          static_cast<local_row_index_type> (lclRow) < static_cast<local_row_index_type> (0)) {
        error_ () = 3;
        return;
      }
      // count_type could differ from the type of each row offset.
      // For example, row offsets might each be 64 bits, but if their
      // difference always fits in 32 bits, we may then safely use a
      // 32-bit count_type.
      const count_type count =
        static_cast<count_type> (rowOffsets_(lclRow+1) - rowOffsets_(lclRow));

      // We pack first the number of entries in the row, then that
      // many global column indices, then that many pids (if any),
      // then that many values.  However, if the number of entries in
      // the row is zero, we pack nothing.
      const count_type numBytes = (count == 0) ?
        static_cast<count_type> (0) :
        sizeOfLclCount_ + count * (sizeOfGblColInd_ +
                                   (lclRowPids_.size() > 0 ? sizeOfPid_ : 0) +
                                   sizeOfValue_);

      if (final) {
        counts_(curInd) = numBytes;
      }
      update += numBytes;
    }
  }

  // mfh 31 May 2017: Don't need init or join.  If you have join, MUST
  // have join both with and without volatile!  Otherwise intrawarp
  // joins are really slow on GPUs.

  int getError () const {
    typedef typename device_type::execution_space execution_space;
    auto error_h = Kokkos::create_mirror_view (error_);
    // DEEP_COPY REVIEW - DEVICE-TO-HOSTMIRROR
    Kokkos::deep_copy (execution_space(), error_h, error_);
    return error_h ();
  }

private:
  OutputOffsetsViewType outputOffsets_;
  CountsViewType counts_;
  typename InputOffsetsViewType::const_type rowOffsets_;
  typename InputLocalRowIndicesViewType::const_type lclRowInds_;
  typename InputLocalRowPidsViewType::const_type lclRowPids_;
  count_type sizeOfLclCount_;
  count_type sizeOfGblColInd_;
  count_type sizeOfPid_;
  count_type sizeOfValue_;
  Kokkos::View<int, device_type> error_;
};

template<class OutputOffsetsViewType,
         class CountsViewType,
         class InputOffsetsViewType,
         class InputLocalRowIndicesViewType,
         class InputLocalRowPidsViewType>
typename CountsViewType::non_const_value_type
computeNumPacketsAndOffsets (const OutputOffsetsViewType& outputOffsets,
                             const CountsViewType& counts,
                             const InputOffsetsViewType& rowOffsets,
                             const InputLocalRowIndicesViewType& lclRowInds,
                             const InputLocalRowPidsViewType& lclRowPids,
                             const typename CountsViewType::non_const_value_type sizeOfLclCount,
                             const typename CountsViewType::non_const_value_type sizeOfGblColInd,
                             const typename CountsViewType::non_const_value_type sizeOfPid,
                             const typename CountsViewType::non_const_value_type sizeOfValue)
{
  typedef NumPacketsAndOffsetsFunctor<OutputOffsetsViewType,
    CountsViewType, typename InputOffsetsViewType::const_type,
    typename InputLocalRowIndicesViewType::const_type,
    typename InputLocalRowPidsViewType::const_type> functor_type;
  typedef typename CountsViewType::non_const_value_type count_type;
  typedef typename OutputOffsetsViewType::size_type size_type;
  typedef typename OutputOffsetsViewType::execution_space execution_space;
  typedef typename functor_type::local_row_index_type LO;
  typedef Kokkos::RangePolicy<execution_space, LO> range_type;
  const char prefix[] = "computeNumPacketsAndOffsets: ";

  count_type count = 0;
  const count_type numRowsToPack = lclRowInds.extent (0);

  if (numRowsToPack == 0) {
    return count;
  }
  else {
    TEUCHOS_TEST_FOR_EXCEPTION
      (rowOffsets.extent (0) <= static_cast<size_type> (1),
       std::invalid_argument, prefix << "There is at least one row to pack, "
       "but the matrix has no rows.  lclRowInds.extent(0) = " <<
       numRowsToPack << ", but rowOffsets.extent(0) = " <<
       rowOffsets.extent (0) << " <= 1.");
    TEUCHOS_TEST_FOR_EXCEPTION
      (outputOffsets.extent (0) !=
       static_cast<size_type> (numRowsToPack + 1), std::invalid_argument,
       prefix << "Output dimension does not match number of rows to pack.  "
       << "outputOffsets.extent(0) = " << outputOffsets.extent (0)
       << " != lclRowInds.extent(0) + 1 = "
       << static_cast<size_type> (numRowsToPack + 1) << ".");
    TEUCHOS_TEST_FOR_EXCEPTION
      (counts.extent (0) != numRowsToPack, std::invalid_argument,
       prefix << "counts.extent(0) = " << counts.extent (0)
       << " != numRowsToPack = " << numRowsToPack << ".");

    functor_type f (outputOffsets, counts, rowOffsets,
                    lclRowInds, lclRowPids, sizeOfLclCount,
                    sizeOfGblColInd, sizeOfPid, sizeOfValue);
    Kokkos::parallel_scan (range_type (0, numRowsToPack + 1), f);

    // At least in debug mode, this functor checks for errors.
    const int errCode = f.getError ();
    TEUCHOS_TEST_FOR_EXCEPTION
      (errCode != 0, std::runtime_error, prefix << "parallel_scan error code "
       << errCode << " != 0.");

#if 0
    size_t total = 0;
    for (LO k = 0; k < numRowsToPack; ++k) {
      total += counts[k];
    }
    if (outputOffsets(numRowsToPack) != total) {
      if (errStr.get () == NULL) {
        errStr = std::unique_ptr<std::ostringstream> (new std::ostringstream ());
      }
      std::ostringstream& os = *errStr;
      os << prefix
         << "outputOffsets(numRowsToPack=" << numRowsToPack << ") "
         << outputOffsets(numRowsToPack) << " != sum of counts = "
         << total << "." << std::endl;
      if (numRowsToPack != 0) {
        // Only print the array if it's not too long.
        if (numRowsToPack < static_cast<LO> (10)) {
          os << "outputOffsets: [";
          for (LO i = 0; i <= numRowsToPack; ++i) {
            os << outputOffsets(i);
            if (static_cast<LO> (i + 1) <= numRowsToPack) {
              os << ",";
            }
          }
          os << "]" << std::endl;
          os << "counts: [";
          for (LO i = 0; i < numRowsToPack; ++i) {
            os << counts(i);
            if (static_cast<LO> (i + 1) < numRowsToPack) {
              os << ",";
            }
          }
          os << "]" << std::endl;
        }
        else {
          os << "outputOffsets(" << (numRowsToPack-1) << ") = "
             << outputOffsets(numRowsToPack-1) << "." << std::endl;
        }
      }
      count = outputOffsets(numRowsToPack);
      return {false, errStr};
    }
#endif // HAVE_TPETRA_DEBUG

    // Get last entry of outputOffsets, which is the sum of the entries
    // of counts.  Don't assume UVM.
    using Tpetra::Details::getEntryOnHost;
    return static_cast<count_type> (getEntryOnHost (outputOffsets,
                                                    numRowsToPack));
  }
}

template<class ST, class ColumnMap, class BufferDeviceType>
KOKKOS_FUNCTION
Kokkos::pair<int, size_t>
packCrsMatrixRow (const ColumnMap& col_map,
                  const Kokkos::View<char*, BufferDeviceType>& exports,
                  const typename PackTraits<typename ColumnMap::local_ordinal_type>::input_array_type& lids_in,
                  const typename PackTraits<int>::input_array_type& pids_in,
                  const typename PackTraits<ST>::input_array_type& vals_in,
                  const size_t offset,
                  const size_t num_ent,
                  const size_t num_bytes_per_value,
                  const bool pack_pids)
{
  using Kokkos::subview;
  using LO = typename ColumnMap::local_ordinal_type;
  using GO = typename ColumnMap::global_ordinal_type;
  using return_type = Kokkos::pair<int, size_t>;

  if (num_ent == 0) {
    // Empty rows always take zero bytes, to ensure sparsity.
    return return_type (0, 0);
  }

  const LO num_ent_LO = static_cast<LO> (num_ent); // packValueCount wants this
  const size_t num_ent_beg = offset;
  const size_t num_ent_len = PackTraits<LO>::packValueCount (num_ent_LO);

  const size_t gids_beg = num_ent_beg + num_ent_len;
  const size_t gids_len = num_ent * PackTraits<GO>::packValueCount (GO (0));

  const size_t pids_beg = gids_beg + gids_len;
  const size_t pids_len = pack_pids ?
    num_ent * PackTraits<int>::packValueCount (int (0)) :
    static_cast<size_t> (0);

  const size_t vals_beg = gids_beg + gids_len + pids_len;
  const size_t vals_len = num_ent * num_bytes_per_value;

  char* const num_ent_out = exports.data () + num_ent_beg;
  char* const gids_out = exports.data () + gids_beg;
  char* const pids_out = pack_pids ? exports.data () + pids_beg : NULL;
  char* const vals_out = exports.data () + vals_beg;

  size_t num_bytes_out = 0;
  int error_code = 0;
  num_bytes_out += PackTraits<LO>::packValue (num_ent_out, num_ent_LO);

  {
    // Copy column indices one at a time, so that we don't need
    // temporary storage.
    for (size_t k = 0; k < num_ent; ++k) {
      const LO lid = lids_in[k];
      const GO gid = col_map.getGlobalElement (lid);
      num_bytes_out += PackTraits<GO>::packValue (gids_out, k, gid);
    }
    // Copy PIDs one at a time, so that we don't need temporary storage.
    if (pack_pids) {
      for (size_t k = 0; k < num_ent; ++k) {
        const LO lid = lids_in[k];
        const int pid = pids_in[lid];
        num_bytes_out += PackTraits<int>::packValue (pids_out, k, pid);
      }
    }
    const auto p =
      PackTraits<ST>::packArray (vals_out, vals_in.data (), num_ent);
    error_code += p.first;
    num_bytes_out += p.second;
  }

  if (error_code != 0) {
    return return_type (10, num_bytes_out);
  }

  const size_t expected_num_bytes =
    num_ent_len + gids_len + pids_len + vals_len;
  if (num_bytes_out != expected_num_bytes) {
    return return_type (11, num_bytes_out);
  }
  return return_type (0, num_bytes_out);
}

template<class LocalMatrix, class LocalMap, class BufferDeviceType>
struct PackCrsMatrixFunctor {
  typedef LocalMatrix local_matrix_device_type;
  typedef LocalMap local_map_type;
  typedef typename local_matrix_device_type::value_type ST;
  typedef typename local_map_type::local_ordinal_type LO;
  typedef typename local_map_type::global_ordinal_type GO;
  typedef typename local_matrix_device_type::device_type DT;

  typedef Kokkos::View<const size_t*, BufferDeviceType>
    num_packets_per_lid_view_type;
  typedef Kokkos::View<const size_t*, BufferDeviceType> offsets_view_type;
  typedef Kokkos::View<char*, BufferDeviceType> exports_view_type;
  using export_lids_view_type = typename PackTraits<LO>::input_array_type;
  using source_pids_view_type = typename PackTraits<int>::input_array_type;

  typedef typename num_packets_per_lid_view_type::non_const_value_type
    count_type;
  typedef typename offsets_view_type::non_const_value_type
    offset_type;
  typedef Kokkos::pair<int, LO> value_type;

  static_assert (std::is_same<LO, typename local_matrix_device_type::ordinal_type>::value,
                 "local_map_type::local_ordinal_type and "
                 "local_matrix_device_type::ordinal_type must be the same.");

  local_matrix_device_type local_matrix;
  local_map_type local_col_map;
  exports_view_type exports;
  num_packets_per_lid_view_type num_packets_per_lid;
  export_lids_view_type export_lids;
  source_pids_view_type source_pids;
  offsets_view_type offsets;
  size_t num_bytes_per_value;
  bool pack_pids;

  PackCrsMatrixFunctor (const local_matrix_device_type& local_matrix_in,
                        const local_map_type& local_col_map_in,
                        const exports_view_type& exports_in,
                        const num_packets_per_lid_view_type& num_packets_per_lid_in,
                        const export_lids_view_type& export_lids_in,
                        const source_pids_view_type& source_pids_in,
                        const offsets_view_type& offsets_in,
                        const size_t num_bytes_per_value_in,
                        const bool pack_pids_in) :
    local_matrix (local_matrix_in),
    local_col_map (local_col_map_in),
    exports (exports_in),
    num_packets_per_lid (num_packets_per_lid_in),
    export_lids (export_lids_in),
    source_pids (source_pids_in),
    offsets (offsets_in),
    num_bytes_per_value (num_bytes_per_value_in),
    pack_pids (pack_pids_in)
  {
    const LO numRows = local_matrix_in.numRows ();
    const LO rowMapDim =
      static_cast<LO> (local_matrix.graph.row_map.extent (0));
    TEUCHOS_TEST_FOR_EXCEPTION
      (numRows != 0 && rowMapDim != numRows + static_cast<LO> (1),
       std::logic_error, "local_matrix.graph.row_map.extent(0) = "
       << rowMapDim << " != numRows (= " << numRows << " ) + 1.");
  }

  KOKKOS_INLINE_FUNCTION void init (value_type& dst) const
  {
    using ::Tpetra::Details::OrdinalTraits;
    dst = Kokkos::make_pair (0, OrdinalTraits<LO>::invalid ());
  }

  KOKKOS_INLINE_FUNCTION void
  join (value_type& dst, const value_type& src) const
  {
    // `dst` should reflect the first (least) bad index and all other
    // associated error codes and data, so prefer keeping it.
    if (src.first != 0 && dst.first == 0) {
      dst = src;
    }
  }

  KOKKOS_INLINE_FUNCTION
  void operator() (const LO i, value_type& dst) const
  {
    const size_t offset = offsets[i];
    const LO export_lid = export_lids[i];
    const size_t buf_size = exports.size();
    const size_t num_bytes = num_packets_per_lid(i);
    const size_t num_ent =
      static_cast<size_t> (local_matrix.graph.row_map[export_lid+1]
                         - local_matrix.graph.row_map[export_lid]);

    // Only pack this row's data if it has a nonzero number of
    // entries.  We can do this because receiving processes get the
    // number of packets, and will know that zero packets means zero
    // entries.
    if (num_ent == 0) {
      return;
    }

    if (export_lid >= local_matrix.numRows ()) {
      if (dst.first != 0) { // keep only the first error
        dst = Kokkos::make_pair (1, i); // invalid row
      }
      return;
    }
    else if ((offset > buf_size || offset + num_bytes > buf_size)) {
      if (dst.first != 0) { // keep only the first error
        dst = Kokkos::make_pair (2, i); // out of bounds
      }
      return;
    }

    // We can now pack this row

    // Since the matrix is locally indexed on the calling process, we
    // have to use its column Map (which it _must_ have in this case)
    // to convert to global indices.
    const auto row_beg = local_matrix.graph.row_map[export_lid];
    const auto row_end = local_matrix.graph.row_map[export_lid + 1];
    auto vals_in = subview (local_matrix.values,
                            Kokkos::make_pair (row_beg, row_end));
    auto lids_in = subview (local_matrix.graph.entries,
                            Kokkos::make_pair (row_beg, row_end));
    typedef local_map_type LMT;
    typedef BufferDeviceType BDT;
    auto p = packCrsMatrixRow<ST, LMT, BDT> (local_col_map, exports, lids_in,
                                             source_pids, vals_in, offset,
                                             num_ent, num_bytes_per_value,
                                             pack_pids);
    int error_code_this_row = p.first;
    size_t num_bytes_packed_this_row = p.second;
    if (error_code_this_row != 0) {
      if (dst.first != 0) { // keep only the first error
        dst = Kokkos::make_pair (error_code_this_row, i); // bad pack
      }
    }
    else if (num_bytes_packed_this_row != num_bytes) {
      if (dst.first != 0) { // keep only the first error
        dst = Kokkos::make_pair (3, i);
      }
    }
  }
};

template<class LocalMatrix, class LocalMap, class BufferDeviceType>
void
do_pack (const LocalMatrix& local_matrix,
         const LocalMap& local_map,
         const Kokkos::View<char*, BufferDeviceType>& exports,
         const typename PackTraits<size_t>::input_array_type& num_packets_per_lid,
         const typename PackTraits<typename LocalMap::local_ordinal_type>::input_array_type& export_lids,
         const typename PackTraits<int>::input_array_type& source_pids,
         const Kokkos::View<const size_t*, BufferDeviceType>& offsets,
         const size_t num_bytes_per_value,
         const bool pack_pids)
{
  using LO = typename LocalMap::local_ordinal_type;
  using DT = typename LocalMatrix::device_type;
  using range_type = Kokkos::RangePolicy<typename DT::execution_space, LO>;
  const char prefix[] = "Tpetra::Details::do_pack: ";

  if (export_lids.extent (0) != 0) {
    TEUCHOS_TEST_FOR_EXCEPTION
      (static_cast<size_t> (offsets.extent (0)) !=
       static_cast<size_t> (export_lids.extent (0) + 1),
       std::invalid_argument, prefix << "offsets.extent(0) = "
       << offsets.extent (0) << " != export_lids.extent(0) (= "
       << export_lids.extent (0) << ") + 1.");
    TEUCHOS_TEST_FOR_EXCEPTION
      (export_lids.extent (0) != num_packets_per_lid.extent (0),
       std::invalid_argument, prefix << "export_lids.extent(0) = " <<
       export_lids.extent (0) << " != num_packets_per_lid.extent(0) = "
       << num_packets_per_lid.extent (0) << ".");
    // If exports has nonzero length at this point, then the matrix
    // has at least one entry to pack.  Thus, if packing process
    // ranks, we had better have at least one process rank to pack.
    TEUCHOS_TEST_FOR_EXCEPTION
      (pack_pids && exports.extent (0) != 0 &&
       source_pids.extent (0) == 0, std::invalid_argument, prefix <<
       "pack_pids is true, and exports.extent(0) = " <<
       exports.extent (0)  << " != 0, meaning that we need to pack at "
       "least one matrix entry, but source_pids.extent(0) = 0.");
  }

  using pack_functor_type =
    PackCrsMatrixFunctor<LocalMatrix, LocalMap, BufferDeviceType>;
  pack_functor_type f (local_matrix, local_map, exports,
                       num_packets_per_lid, export_lids,
                       source_pids, offsets, num_bytes_per_value,
                       pack_pids);

  typename pack_functor_type::value_type result;
  range_type range (0, num_packets_per_lid.extent (0));
  Kokkos::parallel_reduce (range, f, result);

  if (result.first != 0) {
    // We can't deep_copy from AnonymousSpace Views, so we can't print
    // out any information from them in case of error.
    TEUCHOS_TEST_FOR_EXCEPTION
      (true, std::runtime_error, prefix << "PackCrsMatrixFunctor "
       "reported error code " << result.first << " for the first "
       "bad row " << result.second << ".");
  }
}

template<typename ST, typename LO, typename GO, typename NT, typename BufferDeviceType>
void
packCrsMatrix (const CrsMatrix<ST, LO, GO, NT>& sourceMatrix,
               Kokkos::DualView<char*, BufferDeviceType>& exports,
               const Kokkos::View<size_t*, BufferDeviceType>& num_packets_per_lid,
               const Kokkos::View<const LO*, BufferDeviceType>& export_lids,
               const Kokkos::View<const int*, typename NT::device_type>& export_pids,
               size_t& constant_num_packets,
               const bool pack_pids)
{
  ::Tpetra::Details::ProfilingRegion region_pack_crs_matrix(
    "Tpetra::Details::PackCrsMatrixImpl::packCrsMatrix",
    "Import/Export"
  );
  using Kokkos::View;
  typedef BufferDeviceType DT;
  typedef Kokkos::DualView<char*, BufferDeviceType> exports_view_type;
  const char prefix[] = "Tpetra::Details::PackCrsMatrixImpl::packCrsMatrix: ";
  constexpr bool debug = false;

  auto local_matrix = sourceMatrix.getLocalMatrixDevice ();
  auto local_col_map = sourceMatrix.getColMap ()->getLocalMap ();

  // Setting this to zero tells the caller to expect a possibly
  // different ("nonconstant") number of packets per local index
  // (i.e., a possibly different number of entries per row).
  constant_num_packets = 0;

  const size_t num_export_lids =
    static_cast<size_t> (export_lids.extent (0));
  TEUCHOS_TEST_FOR_EXCEPTION
    (num_export_lids !=
     static_cast<size_t> (num_packets_per_lid.extent (0)),
     std::invalid_argument, prefix << "num_export_lids.extent(0) = "
     << num_export_lids << " != num_packets_per_lid.extent(0) = "
     << num_packets_per_lid.extent (0) << ".");
  if (num_export_lids != 0) {
    TEUCHOS_TEST_FOR_EXCEPTION
      (num_packets_per_lid.data () == NULL, std::invalid_argument,
       prefix << "num_export_lids = "<< num_export_lids << " != 0, but "
       "num_packets_per_lid.data() = "
       << num_packets_per_lid.data () << " == NULL.");
  }

  const size_t num_bytes_per_lid = PackTraits<LO>::packValueCount (LO (0));
  const size_t num_bytes_per_gid = PackTraits<GO>::packValueCount (GO (0));
  const size_t num_bytes_per_pid = PackTraits<int>::packValueCount (int (0));

  size_t num_bytes_per_value = 0;
  if (PackTraits<ST>::compileTimeSize) {
    // Assume ST is default constructible; packValueCount wants an instance.
    num_bytes_per_value = PackTraits<ST>::packValueCount (ST ());
  }
  else {
    // Since the packed data come from the source matrix, we can use
    // the source matrix to get the number of bytes per Scalar value
    // stored in the matrix.  This assumes that all Scalar values in
    // the source matrix require the same number of bytes.  If the
    // source matrix has no entries on the calling process, then we
    // hope that some process does have some idea how big a Scalar
    // value is.  Of course, if no processes have any entries, then no
    // values should be packed (though this does assume that in our
    // packing scheme, rows with zero entries take zero bytes).
    size_t num_bytes_per_value_l = 0;
    if (local_matrix.values.extent(0) > 0) {
      const ST& val = local_matrix.values(0);
      num_bytes_per_value_l = PackTraits<ST>::packValueCount (val);
    }
    using Teuchos::reduceAll;
    reduceAll<int, size_t> (* (sourceMatrix.getComm ()),
                            Teuchos::REDUCE_MAX,
                            num_bytes_per_value_l,
                            Teuchos::outArg (num_bytes_per_value));
  }

  if (num_export_lids == 0) {
    exports = exports_view_type ("exports", 0);
    return;
  }

  // Array of offsets into the pack buffer.
  Kokkos::View<size_t*, DT> offsets ("offsets", num_export_lids + 1);

  // Compute number of packets per LID (row to send), as well as
  // corresponding offsets (the prefix sum of the packet counts).
  const size_t count =
    computeNumPacketsAndOffsets (offsets, num_packets_per_lid,
                                 local_matrix.graph.row_map, export_lids,
                                 export_pids,
                                 num_bytes_per_lid, num_bytes_per_gid,
                                 num_bytes_per_pid, num_bytes_per_value);

  // Resize the output pack buffer if needed.
  if (count > static_cast<size_t> (exports.extent (0))) {
    exports = exports_view_type ("exports", count);
    if (debug) {
      std::ostringstream os;
      os << "*** exports resized to " << count << std::endl;
      std::cerr << os.str ();
    }
  }
  if (debug) {
    std::ostringstream os;
    os << "*** count: " << count << ", exports.extent(0): "
       << exports.extent (0) << std::endl;
    std::cerr << os.str ();
  }

  // If exports has nonzero length at this point, then the matrix has
  // at least one entry to pack.  Thus, if packing process ranks, we
  // had better have at least one process rank to pack.
  TEUCHOS_TEST_FOR_EXCEPTION
    (pack_pids && exports.extent (0) != 0 &&
     export_pids.extent (0) == 0, std::invalid_argument, prefix <<
     "pack_pids is true, and exports.extent(0) = " <<
     exports.extent (0)  << " != 0, meaning that we need to pack at least "
     "one matrix entry, but export_pids.extent(0) = 0.");

  typedef typename std::decay<decltype (local_matrix)>::type
    local_matrix_device_type;
  typedef typename std::decay<decltype (local_col_map)>::type
    local_map_type;

  exports.modify_device ();
  auto exports_d = exports.view_device ();
  do_pack<local_matrix_device_type, local_map_type, DT>
    (local_matrix, local_col_map, exports_d, num_packets_per_lid,
     export_lids, export_pids, offsets, num_bytes_per_value,
     pack_pids);
  // If we got this far, we succeeded.
}

} // namespace PackCrsMatrixImpl

template<typename ST, typename LO, typename GO, typename NT>
void
packCrsMatrix (const CrsMatrix<ST, LO, GO, NT>& sourceMatrix,
               Teuchos::Array<char>& exports,
               const Teuchos::ArrayView<size_t>& numPacketsPerLID,
               const Teuchos::ArrayView<const LO>& exportLIDs,
               size_t& constantNumPackets)
{

  using device_type = typename CrsMatrix<ST,LO,GO,NT>::local_matrix_device_type::device_type;
  using buffer_device_type = typename DistObject<char, LO, GO, NT>::buffer_device_type;
  using host_exec_space = typename Kokkos::View<size_t*, device_type>::HostMirror::execution_space;
  using device_exec_space = typename device_type::execution_space;
  using host_dev_type = Kokkos::Device<host_exec_space, Kokkos::HostSpace>;

  // Convert all Teuchos::Array to Kokkos::View

  // This is an output array, so we don't have to copy to device here.
  // However, we'll have to remember to copy back to host when done.
  Kokkos::View<size_t*, buffer_device_type> num_packets_per_lid_d =
    create_mirror_view_from_raw_host_array (buffer_device_type (),
                                            numPacketsPerLID.getRawPtr (),
                                            numPacketsPerLID.size (), false,
                                            "num_packets_per_lid");
  // FIXME (mfh 05 Feb 2019) We should just pass the exportLIDs
  // DualView through here, instead of recreating a device View from a
  // host ArrayView that itself came from a DualView.
  //
  // This is an input array, so we have to copy to device here.
  // However, we never need to copy it back to host.
  Kokkos::View<const LO*, buffer_device_type> export_lids_d =
    create_mirror_view_from_raw_host_array (buffer_device_type (),
                                            exportLIDs.getRawPtr (),
                                            exportLIDs.size (), true,
                                            "export_lids");

  Kokkos::View<int*, device_type> export_pids_d; // output arg
  Kokkos::DualView<char*, buffer_device_type> exports_dv; // output arg
  constexpr bool pack_pids = false;
  PackCrsMatrixImpl::packCrsMatrix<ST, LO, GO, NT, buffer_device_type> (
      sourceMatrix, exports_dv, num_packets_per_lid_d, export_lids_d,
      export_pids_d, constantNumPackets, pack_pids);

  // The counts are an output of PackCrsMatrixImpl::packCrsMatrix, so we have to
  // copy them back to host.
  Kokkos::View<size_t*, host_dev_type> num_packets_per_lid_h
    (numPacketsPerLID.getRawPtr (),
     numPacketsPerLID.size ());
  // DEEP_COPY REVIEW - DEVICE-TO-HOST
  Kokkos::deep_copy (device_exec_space(), num_packets_per_lid_h, num_packets_per_lid_d);

  // FIXME (mfh 23 Aug 2017) If we're forced to use a DualView for
  // exports_dv above, then we have two host copies for exports_h.

  // The exports are an output of PackCrsMatrixImpl::packCrsMatrix, so we have
  // to copy them back to host.
  if (static_cast<size_t> (exports.size ()) !=
      static_cast<size_t> (exports_dv.extent (0))) {
    exports.resize (exports_dv.extent (0));
  }
  Kokkos::View<char*, host_dev_type> exports_h (exports.getRawPtr (),
                                                exports.size ());
  // DEEP_COPY REVIEW - DEVICE-TO-HOST
  Kokkos::deep_copy (device_exec_space(), exports_h, exports_dv.d_view);
}

template<typename ST, typename LO, typename GO, typename NT>
void
packCrsMatrixNew(
  const CrsMatrix<ST, LO, GO, NT>& sourceMatrix,
  Kokkos::DualView<char*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& exports,
  const Kokkos::DualView<size_t*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& numPacketsPerLID,
  const Kokkos::DualView<const LO*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& exportLIDs,
  size_t& constantNumPackets)
{
  using device_type = typename CrsMatrix<ST, LO, GO, NT>::device_type;
  using buffer_device_type = typename DistObject<char, LO, GO, NT>::buffer_device_type;

  // Create an empty array of PIDs, since the interface needs it.
  Kokkos::View<int*, device_type> exportPIDs_d ("exportPIDs", 0);
  constexpr bool pack_pids = false;

  // Write-only device access
  auto numPacketsPerLID_nc = numPacketsPerLID; // const DV& -> DV
  numPacketsPerLID_nc.clear_sync_state ();
  numPacketsPerLID_nc.modify_device ();
  auto numPacketsPerLID_d = numPacketsPerLID.view_device ();

  // Read-only device access
  TEUCHOS_ASSERT( ! exportLIDs.need_sync_device () );
  auto exportLIDs_d = exportLIDs.view_device ();

  ::Tpetra::Details::ProfilingRegion region_pack_crs_matrix_new(
    "Tpetra::Details::packCrsMatrixNew",
    "Import/Export"
  );
  PackCrsMatrixImpl::packCrsMatrix<ST,LO,GO,NT,buffer_device_type> (
      sourceMatrix, exports, numPacketsPerLID_d, exportLIDs_d,
      exportPIDs_d, constantNumPackets, pack_pids);
}

template<typename ST, typename LO, typename GO, typename NT>
void
packCrsMatrixWithOwningPIDs (const CrsMatrix<ST, LO, GO, NT>& sourceMatrix,
                             Kokkos::DualView<char*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& exports_dv,
                             const Teuchos::ArrayView<size_t>& numPacketsPerLID,
                             const Teuchos::ArrayView<const LO>& exportLIDs,
                             const Teuchos::ArrayView<const int>& sourcePIDs,
                             size_t& constantNumPackets)
{
  typedef typename CrsMatrix<ST,LO,GO,NT>::local_matrix_device_type local_matrix_device_type;
  typedef typename DistObject<char, LO, GO, NT>::buffer_device_type buffer_device_type;
  typedef typename Kokkos::DualView<char*, buffer_device_type>::t_host::execution_space host_exec_space;
  typedef Kokkos::Device<host_exec_space, Kokkos::HostSpace> host_dev_type;

  typename local_matrix_device_type::device_type outputDevice;
  typedef typename NT::execution_space execution_space;


  const bool verbose = ::Tpetra::Details::Behavior::verbose ();
  std::unique_ptr<std::string> prefix;
  if (verbose) {
    const int myRank = [&] () {
      auto map = sourceMatrix.getMap ();
      if (map.get () == nullptr) {
        return -1;
      }
      auto comm = map->getComm ();
      if (comm.get () == nullptr) {
        return -2;
      }
      return comm->getRank ();
    } ();
    std::ostringstream os;
    os << "Proc " << myRank << ": packCrsMatrixWithOwningPIDs: ";
    prefix = std::unique_ptr<std::string> (new std::string (os.str ()));

    std::ostringstream os2;
    os2 << *prefix << "start" << std::endl;
    std::cerr << os2.str ();
  }

  // Convert all Teuchos::Array to Kokkos::View

  // This is an output array, so we don't have to copy to device here.
  // However, we'll have to remember to copy back to host when done.
  auto num_packets_per_lid_d =
    create_mirror_view_from_raw_host_array (buffer_device_type (),
                                            numPacketsPerLID.getRawPtr (),
                                            numPacketsPerLID.size (), false,
                                            "num_packets_per_lid");

  // This is an input array, so we have to copy to device here.
  // However, we never need to copy it back to host.
  auto export_lids_d =
    create_mirror_view_from_raw_host_array (buffer_device_type (),
                                            exportLIDs.getRawPtr (),
                                            exportLIDs.size (), true,
                                            "export_lids");
  // This is an input array, so we have to copy to device here.
  // However, we never need to copy it back to host.
  auto export_pids_d =
    create_mirror_view_from_raw_host_array (outputDevice,
                                            sourcePIDs.getRawPtr (),
                                            sourcePIDs.size (), true,
                                            "export_pids");
  constexpr bool pack_pids = true;
  try {
    PackCrsMatrixImpl::packCrsMatrix
      (sourceMatrix, exports_dv, num_packets_per_lid_d, export_lids_d,
       export_pids_d, constantNumPackets, pack_pids);
  }
  catch (std::exception& e) {
    if (verbose) {
      std::ostringstream os;
      os << *prefix << "PackCrsMatrixImpl::packCrsMatrix threw: "
         << e.what () << std::endl;
      std::cerr << os.str ();
    }
    throw;
  }
  catch (...) {
    if (verbose) {
      std::ostringstream os;
      os << *prefix << "PackCrsMatrixImpl::packCrsMatrix threw an exception "
        "not a subclass of std::exception" << std::endl;
      std::cerr << os.str ();
    }
    throw;
  }

  if (numPacketsPerLID.size () != 0) {
    try {
      // The counts are an output of PackCrsMatrixImpl::packCrsMatrix,
      // so we have to copy them back to host.
      Kokkos::View<size_t*, host_dev_type> num_packets_per_lid_h
        (numPacketsPerLID.getRawPtr (), numPacketsPerLID.size ());
      // DEEP_COPY REVIEW - DEVICE-TO-HOST
      Kokkos::deep_copy (execution_space(), num_packets_per_lid_h, num_packets_per_lid_d);
    }
    catch (std::exception& e) {
      if (verbose) {
        std::ostringstream os;
        os << *prefix << "Kokkos::deep_copy threw: " << e.what () << std::endl;
        std::cerr << os.str ();
      }
      throw;
    }
    catch (...) {
      if (verbose) {
        std::ostringstream os;
        os << *prefix << "Kokkos::deep_copy threw an exception not a subclass "
          "of std::exception" << std::endl;
        std::cerr << os.str ();
      }
      throw;
    }
  }

  if (verbose) {
    std::ostringstream os;
    os << *prefix << "done" << std::endl;
    std::cerr << os.str ();
  }
}

} // namespace Details
} // namespace Tpetra

#define TPETRA_DETAILS_PACKCRSMATRIX_INSTANT( ST, LO, GO, NT ) \
  template void \
  Details::packCrsMatrix<ST, LO, GO, NT> (const CrsMatrix<ST, LO, GO, NT>&, \
    Teuchos::Array<char>&, \
    const Teuchos::ArrayView<size_t>&, \
    const Teuchos::ArrayView<const LO>&, \
    size_t&); \
  template void \
  Details::packCrsMatrixNew<ST, LO, GO, NT> (const CrsMatrix<ST, LO, GO, NT>&, \
    Kokkos::DualView<char*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
    const Kokkos::DualView<size_t*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
    const Kokkos::DualView<const LO*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
    size_t&); \
  template void \
  Details::packCrsMatrixWithOwningPIDs<ST, LO, GO, NT> (const CrsMatrix<ST, LO, GO, NT>&, \
    Kokkos::DualView<char*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
    const Teuchos::ArrayView<size_t>&, \
    const Teuchos::ArrayView<const LO>&, \
    const Teuchos::ArrayView<const int>&, \
    size_t&);

#endif // TPETRA_DETAILS_PACKCRSMATRIX_DEF_HPP