Tpetra_Import_Util2.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_IMPORT_UTIL2_HPP
#define TPETRA_IMPORT_UTIL2_HPP


#include "Tpetra_ConfigDefs.hpp"
#include "Tpetra_Import.hpp"
#include "Tpetra_HashTable.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_Util.hpp"
#include "Tpetra_Distributor.hpp"
#include "Tpetra_Details_reallocDualViewIfNeeded.hpp"
#include "Tpetra_Details_MpiTypeTraits.hpp"
#include "Tpetra_Vector.hpp"
#include "Kokkos_DualView.hpp"
#include "KokkosSparse_SortCrs.hpp"
#include <Teuchos_Array.hpp>
#include "Tpetra_Details_createMirrorView.hpp"
#include <Kokkos_UnorderedMap.hpp>
#include <unordered_map>
#include <utility>
#include <set>

#include "Tpetra_CrsMatrix_decl.hpp"

#include <Kokkos_Core.hpp>
#include <Kokkos_Sort.hpp>

namespace Tpetra {
namespace Import_Util {

template <typename Scalar, typename Ordinal>
void sortCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                    const Teuchos::ArrayView<Ordinal>& CRS_colind,
                    const Teuchos::ArrayView<Scalar>& CRS_vals);

template <typename Ordinal>
void sortCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                    const Teuchos::ArrayView<Ordinal>& CRS_colind);

template <typename rowptr_array_type, typename colind_array_type, typename vals_array_type>
void sortCrsEntries(const rowptr_array_type& CRS_rowptr,
                    const colind_array_type& CRS_colind,
                    const vals_array_type& CRS_vals);

template <typename rowptr_array_type, typename colind_array_type>
void sortCrsEntries(const rowptr_array_type& CRS_rowptr,
                    const colind_array_type& CRS_colind);

template <typename Scalar, typename Ordinal>
void sortAndMergeCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                            const Teuchos::ArrayView<Ordinal>& CRS_colind,
                            const Teuchos::ArrayView<Scalar>& CRS_vals);

template <typename Ordinal>
void sortAndMergeCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                            const Teuchos::ArrayView<Ordinal>& CRS_colind);

template <class rowptr_view_type, class colind_view_type, class vals_view_type>
void sortAndMergeCrsEntries(const rowptr_view_type& CRS_rowptr,
                            const colind_view_type& CRS_colind,
                            const vals_view_type& CRS_vals);

template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void lowCommunicationMakeColMapAndReindex(
    const Teuchos::ArrayView<const size_t>& rowptr,
    const Teuchos::ArrayView<LocalOrdinal>& colind_LID,
    const Teuchos::ArrayView<GlobalOrdinal>& colind_GID,
    const Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& domainMapRCP,
    const Teuchos::ArrayView<const int>& owningPIDs,
    Teuchos::Array<int>& remotePIDs,
    Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& colMap);

template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void lowCommunicationMakeColMapAndReindex(
    const Kokkos::View<size_t*, typename Node::device_type> rowptr_view,
    const Kokkos::View<LocalOrdinal*, typename Node::device_type> colind_LID_view,
    const Kokkos::View<GlobalOrdinal*, typename Node::device_type> colind_GID_view,
    const Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& domainMapRCP,
    const Teuchos::ArrayView<const int>& owningPIDs,
    Teuchos::Array<int>& remotePIDs,
    Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& colMap);

template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void getTwoTransferOwnershipVector(const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& transferThatDefinesOwnership,
                                   bool useReverseModeForOwnership,
                                   const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& transferForMigratingData,
                                   bool useReverseModeForMigration,
                                   Tpetra::Vector<int, LocalOrdinal, GlobalOrdinal, Node>& owningPIDs);

}  // namespace Import_Util
}  // namespace Tpetra

//
// Implementations
//

namespace Tpetra {
namespace Import_Util {

template <typename PID, typename GlobalOrdinal>
bool sort_PID_then_GID(const std::pair<PID, GlobalOrdinal>& a,
                       const std::pair<PID, GlobalOrdinal>& b) {
  if (a.first != b.first)
    return (a.first < b.first);
  return (a.second < b.second);
}

template <typename PID,
          typename GlobalOrdinal,
          typename LocalOrdinal>
bool sort_PID_then_pair_GID_LID(const std::pair<PID, std::pair<GlobalOrdinal, LocalOrdinal>>& a,
                                const std::pair<PID, std::pair<GlobalOrdinal, LocalOrdinal>>& b) {
  if (a.first != b.first)
    return a.first < b.first;
  else
    return (a.second.first < b.second.first);
}

template <typename Scalar,
          typename LocalOrdinal,
          typename GlobalOrdinal,
          typename Node>
void reverseNeighborDiscovery(const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& SourceMatrix,
                              const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::row_ptrs_host_view_type& rowptr,
                              const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_inds_host_view_type& colind,
                              const Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& RowTransfer,
                              Teuchos::RCP<const Tpetra::Import<LocalOrdinal, GlobalOrdinal, Node>> MyImporter,
                              Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> MyDomainMap,
                              Teuchos::ArrayRCP<int>& type3PIDs,
                              Teuchos::ArrayRCP<LocalOrdinal>& type3LIDs,
                              Teuchos::RCP<const Teuchos::Comm<int>>& rcomm) {
#ifdef HAVE_TPETRACORE_MPI
  using Teuchos::TimeMonitor;
  using ::Tpetra::Details::Behavior;
  typedef LocalOrdinal LO;
  typedef GlobalOrdinal GO;
  typedef std::pair<GO, GO> pidgidpair_t;
  using Teuchos::RCP;
  const std::string prefix{" Import_Util2::ReverseND:: "};
  const std::string label("IU2::Neighbor");

  // There can be no neighbor discovery if you don't have an importer
  if (MyImporter.is_null()) return;

  std::ostringstream errstr;
  bool error      = false;
  auto const comm = MyDomainMap->getComm();

  MPI_Comm rawComm = getRawMpiComm(*comm);
  const int MyPID  = rcomm->getRank();

  // Things related to messages I am sending in forward mode (RowTransfer)
  // *** Note: this will be incorrect for transferAndFillComplete if it is in reverse mode. FIXME cbl.
  auto ExportPIDs    = RowTransfer.getExportPIDs();
  auto ExportLIDs    = RowTransfer.getExportLIDs();
  auto NumExportLIDs = RowTransfer.getNumExportIDs();

  Distributor& Distor   = MyImporter->getDistributor();
  const size_t NumRecvs = Distor.getNumReceives();
  const size_t NumSends = Distor.getNumSends();
  auto RemoteLIDs       = MyImporter->getRemoteLIDs();
  auto const ProcsFrom  = Distor.getProcsFrom();
  auto const ProcsTo    = Distor.getProcsTo();

  auto LengthsFrom                                                 = Distor.getLengthsFrom();
  auto MyColMap                                                    = SourceMatrix.getColMap();
  const size_t numCols                                             = MyColMap->getLocalNumElements();
  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> target = MyImporter->getTargetMap();

  // Get the owning pids in a special way,
  // s.t. ProcsFrom[RemotePIDs[i]] is the proc that owns RemoteLIDs[j]....
  Teuchos::Array<int> RemotePIDOrder(numCols, -1);

  // For each remote ID, record index into ProcsFrom, who owns it.
  for (size_t i = 0, j = 0; i < NumRecvs; ++i) {
    for (size_t k = 0; k < LengthsFrom[i]; ++k) {
      const int pid = ProcsFrom[i];
      if (pid != MyPID) {
        RemotePIDOrder[RemoteLIDs[j]] = i;
      }
      j++;
    }
  }

  // Step One: Start tacking the (GID,PID) pairs on the std sets
  //
  // For each index in ProcsFrom, we will insert into a set of (PID,
  // GID) pairs, in order to build a list of such pairs for each of
  // those processes.  Since this is building a reverse, we will send
  // to these processes.
  Teuchos::Array<int> ReverseSendSizes(NumRecvs, 0);
  // do this as C array to avoid Teuchos::Array value initialization of all reserved memory
  Teuchos::Array<Teuchos::ArrayRCP<pidgidpair_t>> RSB(NumRecvs);

  {
#ifdef HAVE_TPETRA_MMM_TIMINGS
    TimeMonitor set_all(*TimeMonitor::getNewTimer(prefix + std::string("isMMallSetRSB")));
#endif

    // 25 Jul 2018: CBL
    // todo:std::unordered_set (hash table),
    // with an adequate prereservation ("bucket count").
    // An onordered_set has to have a custom hasher for pid/gid pair
    // However, when pidsets is copied to RSB, it will be in key
    // order _not_ in pid,gid order. (unlike std::set).
    // Impliment this with a reserve, and time BOTH building pidsets
    // _and_ the sort after the receive. Even if unordered_set saves
    // time, if it causes the sort to be longer, it's not a win.

    Teuchos::Array<std::set<pidgidpair_t>> pidsets(NumRecvs);
    {
#ifdef HAVE_TPETRA_MMM_TIMINGS
      TimeMonitor set_insert(*TimeMonitor::getNewTimer(prefix + std::string("isMMallSetRSBinsert")));
#endif
      for (size_t i = 0; i < NumExportLIDs; i++) {
        LO lid     = ExportLIDs[i];
        GO exp_pid = ExportPIDs[i];
        for (auto j = rowptr[lid]; j < rowptr[lid + 1]; j++) {
          int pid_order = RemotePIDOrder[colind[j]];
          if (pid_order != -1) {
            GO gid     = MyColMap->getGlobalElement(colind[j]);  // Epetra SM.GCID46 =>sm->graph-> {colmap(colind)}
            auto tpair = pidgidpair_t(exp_pid, gid);
            pidsets[pid_order].insert(pidsets[pid_order].end(), tpair);
          }
        }
      }
    }

    {
#ifdef HAVE_TPETRA_MMM_TIMINGS
      TimeMonitor set_cpy(*TimeMonitor::getNewTimer(prefix + std::string("isMMallSetRSBcpy")));
#endif
      int jj = 0;
      for (auto&& ps : pidsets) {
        auto s  = ps.size();
        RSB[jj] = Teuchos::arcp(new pidgidpair_t[s], 0, s, true);
        std::copy(ps.begin(), ps.end(), RSB[jj]);
        ReverseSendSizes[jj] = s;
        ++jj;
      }
    }
  }  // end of set based packing.

  Teuchos::Array<int> ReverseRecvSizes(NumSends, -1);
  Teuchos::Array<MPI_Request> rawBreq(ProcsFrom.size() + ProcsTo.size(), MPI_REQUEST_NULL);
  // 25 Jul 2018: MPI_TAG_UB is the largest tag value; could be < 32768.
  const int mpi_tag_base_ = 3;

  int mpireq_idx = 0;
  for (int i = 0; i < ProcsTo.size(); ++i) {
    int Rec_Tag            = mpi_tag_base_ + ProcsTo[i];
    int* thisrecv          = (int*)(&ReverseRecvSizes[i]);
    MPI_Request rawRequest = MPI_REQUEST_NULL;
    MPI_Irecv(const_cast<int*>(thisrecv),
              1,
              MPI_INT,
              ProcsTo[i],
              Rec_Tag,
              rawComm,
              &rawRequest);
    rawBreq[mpireq_idx++] = rawRequest;
  }
  for (int i = 0; i < ProcsFrom.size(); ++i) {
    int Send_Tag           = mpi_tag_base_ + MyPID;
    int* mysend            = (int*)(&ReverseSendSizes[i]);
    MPI_Request rawRequest = MPI_REQUEST_NULL;
    MPI_Isend(mysend,
              1,
              MPI_INT,
              ProcsFrom[i],
              Send_Tag,
              rawComm,
              &rawRequest);
    rawBreq[mpireq_idx++] = rawRequest;
  }
  Teuchos::Array<MPI_Status> rawBstatus(rawBreq.size());
#ifdef HAVE_TPETRA_DEBUG
  const int err1 =
#endif
      MPI_Waitall(rawBreq.size(), rawBreq.getRawPtr(),
                  rawBstatus.getRawPtr());

#ifdef HAVE_TPETRA_DEBUG
  if (err1) {
    errstr << MyPID << "sE1 reverseNeighborDiscovery Mpi_Waitall error on send ";
    error = true;
    std::cerr << errstr.str() << std::flush;
  }
#endif

  int totalexportpairrecsize = 0;
  for (size_t i = 0; i < NumSends; ++i) {
    totalexportpairrecsize += ReverseRecvSizes[i];
#ifdef HAVE_TPETRA_DEBUG
    if (ReverseRecvSizes[i] < 0) {
      errstr << MyPID << "E4 reverseNeighborDiscovery: got < 0 for receive size " << ReverseRecvSizes[i] << std::endl;
      error = true;
    }
#endif
  }
  Teuchos::ArrayRCP<pidgidpair_t> AllReverseRecv = Teuchos::arcp(new pidgidpair_t[totalexportpairrecsize], 0, totalexportpairrecsize, true);
  int offset                                     = 0;
  mpireq_idx                                     = 0;
  for (int i = 0; i < ProcsTo.size(); ++i) {
    int recv_data_size     = ReverseRecvSizes[i] * 2;
    int recvData_MPI_Tag   = mpi_tag_base_ * 2 + ProcsTo[i];
    MPI_Request rawRequest = MPI_REQUEST_NULL;
    GO* rec_bptr           = (GO*)(&AllReverseRecv[offset]);
    offset += ReverseRecvSizes[i];
    MPI_Irecv(rec_bptr,
              recv_data_size,
              ::Tpetra::Details::MpiTypeTraits<GO>::getType(rec_bptr[0]),
              ProcsTo[i],
              recvData_MPI_Tag,
              rawComm,
              &rawRequest);
    rawBreq[mpireq_idx++] = rawRequest;
  }
  for (int ii = 0; ii < ProcsFrom.size(); ++ii) {
    GO* send_bptr          = (GO*)(RSB[ii].getRawPtr());
    MPI_Request rawSequest = MPI_REQUEST_NULL;
    int send_data_size     = ReverseSendSizes[ii] * 2;  // 2 == count of pair
    int sendData_MPI_Tag   = mpi_tag_base_ * 2 + MyPID;
    MPI_Isend(send_bptr,
              send_data_size,
              ::Tpetra::Details::MpiTypeTraits<GO>::getType(send_bptr[0]),
              ProcsFrom[ii],
              sendData_MPI_Tag,
              rawComm,
              &rawSequest);

    rawBreq[mpireq_idx++] = rawSequest;
  }
#ifdef HAVE_TPETRA_DEBUG
  const int err =
#endif
      MPI_Waitall(rawBreq.size(),
                  rawBreq.getRawPtr(),
                  rawBstatus.getRawPtr());
#ifdef HAVE_TPETRA_DEBUG
  if (err) {
    errstr << MyPID << "E3.r reverseNeighborDiscovery Mpi_Waitall error on receive ";
    error = true;
    std::cerr << errstr.str() << std::flush;
  }
#endif
  std::sort(AllReverseRecv.begin(), AllReverseRecv.end(), Tpetra::Import_Util::sort_PID_then_GID<GlobalOrdinal, GlobalOrdinal>);

  auto newEndOfPairs = std::unique(AllReverseRecv.begin(), AllReverseRecv.end());
  // don't resize to remove non-unique, just use the end-of-unique iterator
  if (AllReverseRecv.begin() == newEndOfPairs) return;
  int ARRsize = std::distance(AllReverseRecv.begin(), newEndOfPairs);
  auto rPIDs  = Teuchos::arcp(new int[ARRsize], 0, ARRsize, true);
  auto rLIDs  = Teuchos::arcp(new LocalOrdinal[ARRsize], 0, ARRsize, true);

  int tsize = 0;
  for (auto itr = AllReverseRecv.begin(); itr != newEndOfPairs; ++itr) {
    if ((int)(itr->first) != MyPID) {
      rPIDs[tsize]     = (int)itr->first;
      LocalOrdinal lid = MyDomainMap->getLocalElement(itr->second);
      rLIDs[tsize]     = lid;
      tsize++;
    }
  }

  type3PIDs = rPIDs.persistingView(0, tsize);
  type3LIDs = rLIDs.persistingView(0, tsize);

  if (error) {
    std::cerr << errstr.str() << std::flush;
    comm->barrier();
    comm->barrier();
    comm->barrier();
    MPI_Abort(MPI_COMM_WORLD, -1);
  }
#endif
}

// Note: This should get merged with the other Tpetra sort routines eventually.
template <typename Scalar, typename Ordinal>
void sortCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                    const Teuchos::ArrayView<Ordinal>& CRS_colind,
                    const Teuchos::ArrayView<Scalar>& CRS_vals) {
  // For each row, sort column entries from smallest to largest.
  // Use shell sort. Stable sort so it is fast if indices are already sorted.
  // Code copied from  Epetra_CrsMatrix::SortEntries()
  size_t NumRows = CRS_rowptr.size() - 1;
  size_t nnz     = CRS_colind.size();

  const bool permute_values_array = CRS_vals.size() > 0;

  for (size_t i = 0; i < NumRows; i++) {
    size_t start = CRS_rowptr[i];
    if (start >= nnz) continue;

    size_t NumEntries = CRS_rowptr[i + 1] - start;
    Teuchos::ArrayRCP<Scalar> locValues;
    if (permute_values_array)
      locValues = Teuchos::arcp<Scalar>(&CRS_vals[start], 0, NumEntries, false);
    Teuchos::ArrayRCP<Ordinal> locIndices(&CRS_colind[start], 0, NumEntries, false);

    Ordinal n = NumEntries;
    Ordinal m = 1;
    while (m < n) m = m * 3 + 1;
    m /= 3;

    while (m > 0) {
      Ordinal max = n - m;
      for (Ordinal j = 0; j < max; j++) {
        for (Ordinal k = j; k >= 0; k -= m) {
          if (locIndices[k + m] >= locIndices[k])
            break;
          if (permute_values_array) {
            Scalar dtemp     = locValues[k + m];
            locValues[k + m] = locValues[k];
            locValues[k]     = dtemp;
          }
          Ordinal itemp     = locIndices[k + m];
          locIndices[k + m] = locIndices[k];
          locIndices[k]     = itemp;
        }
      }
      m = m / 3;
    }
  }
}

template <typename Ordinal>
void sortCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                    const Teuchos::ArrayView<Ordinal>& CRS_colind) {
  // Generate dummy values array
  Teuchos::ArrayView<Tpetra::Details::DefaultTypes::scalar_type> CRS_vals;
  sortCrsEntries(CRS_rowptr, CRS_colind, CRS_vals);
}

namespace Impl {

template <class RowOffsetsType, class ColumnIndicesType, class ValuesType>
class SortCrsEntries {
 private:
  typedef typename ColumnIndicesType::non_const_value_type ordinal_type;
  typedef typename ValuesType::non_const_value_type scalar_type;

 public:
  SortCrsEntries(const RowOffsetsType& ptr,
                 const ColumnIndicesType& ind,
                 const ValuesType& val)
    : ptr_(ptr)
    , ind_(ind)
    , val_(val) {
    static_assert(std::is_signed<ordinal_type>::value,
                  "The type of each "
                  "column index -- that is, the type of each entry of ind "
                  "-- must be signed in order for this functor to work.");
  }

  KOKKOS_FUNCTION void operator()(const size_t i) const {
    const size_t nnz                = ind_.extent(0);
    const size_t start              = ptr_(i);
    const bool permute_values_array = val_.extent(0) > 0;

    if (start < nnz) {
      const size_t NumEntries = ptr_(i + 1) - start;

      const ordinal_type n = static_cast<ordinal_type>(NumEntries);
      ordinal_type m       = 1;
      while (m < n) m = m * 3 + 1;
      m /= 3;

      while (m > 0) {
        ordinal_type max = n - m;
        for (ordinal_type j = 0; j < max; j++) {
          for (ordinal_type k = j; k >= 0; k -= m) {
            const size_t sk = start + k;
            if (ind_(sk + m) >= ind_(sk)) {
              break;
            }
            if (permute_values_array) {
              const scalar_type dtemp = val_(sk + m);
              val_(sk + m)            = val_(sk);
              val_(sk)                = dtemp;
            }
            const ordinal_type itemp = ind_(sk + m);
            ind_(sk + m)             = ind_(sk);
            ind_(sk)                 = itemp;
          }
        }
        m = m / 3;
      }
    }
  }

  static void
  sortCrsEntries(const RowOffsetsType& ptr,
                 const ColumnIndicesType& ind,
                 const ValuesType& val) {
    // For each row, sort column entries from smallest to largest.
    // Use shell sort. Stable sort so it is fast if indices are already sorted.
    // Code copied from  Epetra_CrsMatrix::SortEntries()
    // NOTE: This should not be taken as a particularly efficient way to sort
    // rows of matrices in parallel.  But it is correct, so that's something.
    if (ptr.extent(0) == 0) {
      return;  // no rows, so nothing to sort
    }
    const size_t NumRows = ptr.extent(0) - 1;

    typedef size_t index_type;  // what this function was using; not my choice
    typedef typename ValuesType::execution_space execution_space;
    // Specify RangePolicy explicitly, in order to use correct execution
    // space.  See #1345.
    typedef Kokkos::RangePolicy<execution_space, index_type> range_type;
    Kokkos::parallel_for("sortCrsEntries", range_type(0, NumRows),
                         SortCrsEntries(ptr, ind, val));
  }

 private:
  RowOffsetsType ptr_;
  ColumnIndicesType ind_;
  ValuesType val_;
};

}  // namespace Impl

template <typename rowptr_array_type, typename colind_array_type, typename vals_array_type>
void sortCrsEntries(const rowptr_array_type& CRS_rowptr,
                    const colind_array_type& CRS_colind,
                    const vals_array_type& CRS_vals) {
  Impl::SortCrsEntries<rowptr_array_type, colind_array_type,
                       vals_array_type>::sortCrsEntries(CRS_rowptr, CRS_colind, CRS_vals);
}

template <typename rowptr_array_type, typename colind_array_type>
void sortCrsEntries(const rowptr_array_type& CRS_rowptr,
                    const colind_array_type& CRS_colind) {
  // Generate dummy values array
  typedef typename colind_array_type::execution_space execution_space;
  typedef Tpetra::Details::DefaultTypes::scalar_type scalar_type;
  typedef typename Kokkos::View<scalar_type*, execution_space> scalar_view_type;
  scalar_view_type CRS_vals;
  sortCrsEntries<rowptr_array_type, colind_array_type,
                 scalar_view_type>(CRS_rowptr, CRS_colind, CRS_vals);
}

// Note: This should get merged with the other Tpetra sort routines eventually.
template <typename Scalar, typename Ordinal>
void sortAndMergeCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                            const Teuchos::ArrayView<Ordinal>& CRS_colind,
                            const Teuchos::ArrayView<Scalar>& CRS_vals) {
  // For each row, sort column entries from smallest to largest,
  // merging column ids that are identify by adding values.  Use shell
  // sort. Stable sort so it is fast if indices are already sorted.
  // Code copied from Epetra_CrsMatrix::SortEntries()

  if (CRS_rowptr.size() == 0) {
    return;  // no rows, so nothing to sort
  }
  const size_t NumRows = CRS_rowptr.size() - 1;
  const size_t nnz     = CRS_colind.size();
  size_t new_curr      = CRS_rowptr[0];
  size_t old_curr      = CRS_rowptr[0];

  const bool permute_values_array = CRS_vals.size() > 0;

  for (size_t i = 0; i < NumRows; i++) {
    const size_t old_rowptr_i = CRS_rowptr[i];
    CRS_rowptr[i]             = old_curr;
    if (old_rowptr_i >= nnz) continue;

    size_t NumEntries = CRS_rowptr[i + 1] - old_rowptr_i;
    Teuchos::ArrayRCP<Scalar> locValues;
    if (permute_values_array)
      locValues = Teuchos::arcp<Scalar>(&CRS_vals[old_rowptr_i], 0, NumEntries, false);
    Teuchos::ArrayRCP<Ordinal> locIndices(&CRS_colind[old_rowptr_i], 0, NumEntries, false);

    // Sort phase
    Ordinal n = NumEntries;
    Ordinal m = n / 2;

    while (m > 0) {
      Ordinal max = n - m;
      for (Ordinal j = 0; j < max; j++) {
        for (Ordinal k = j; k >= 0; k -= m) {
          if (locIndices[k + m] >= locIndices[k])
            break;
          if (permute_values_array) {
            Scalar dtemp     = locValues[k + m];
            locValues[k + m] = locValues[k];
            locValues[k]     = dtemp;
          }
          Ordinal itemp     = locIndices[k + m];
          locIndices[k + m] = locIndices[k];
          locIndices[k]     = itemp;
        }
      }
      m = m / 2;
    }

    // Merge & shrink
    for (size_t j = old_rowptr_i; j < CRS_rowptr[i + 1]; j++) {
      if (j > old_rowptr_i && CRS_colind[j] == CRS_colind[new_curr - 1]) {
        if (permute_values_array) CRS_vals[new_curr - 1] += CRS_vals[j];
      } else if (new_curr == j) {
        new_curr++;
      } else {
        CRS_colind[new_curr] = CRS_colind[j];
        if (permute_values_array) CRS_vals[new_curr] = CRS_vals[j];
        new_curr++;
      }
    }
    old_curr = new_curr;
  }

  CRS_rowptr[NumRows] = new_curr;
}

template <typename Ordinal>
void sortAndMergeCrsEntries(const Teuchos::ArrayView<size_t>& CRS_rowptr,
                            const Teuchos::ArrayView<Ordinal>& CRS_colind) {
  Teuchos::ArrayView<Tpetra::Details::DefaultTypes::scalar_type> CRS_vals;
  return sortAndMergeCrsEntries<Tpetra::Details::DefaultTypes::scalar_type, Ordinal>(CRS_rowptr, CRS_colind, CRS_vals);
}

template <class rowptr_view_type, class colind_view_type, class vals_view_type>
void sortAndMergeCrsEntries(rowptr_view_type& CRS_rowptr,
                            colind_view_type& CRS_colind,
                            vals_view_type& CRS_vals) {
  using execution_space = typename vals_view_type::execution_space;

  auto CRS_rowptr_in = CRS_rowptr;
  auto CRS_colind_in = CRS_colind;
  auto CRS_vals_in   = CRS_vals;

  KokkosSparse::sort_and_merge_matrix<execution_space, rowptr_view_type,
                                      colind_view_type, vals_view_type>(CRS_rowptr_in, CRS_colind_in, CRS_vals_in,
                                                                        CRS_rowptr, CRS_colind, CRS_vals);
}

template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void lowCommunicationMakeColMapAndReindexSerial(const Teuchos::ArrayView<const size_t>& rowptr,
                                                const Teuchos::ArrayView<LocalOrdinal>& colind_LID,
                                                const Teuchos::ArrayView<GlobalOrdinal>& colind_GID,
                                                const Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& domainMapRCP,
                                                const Teuchos::ArrayView<const int>& owningPIDs,
                                                Teuchos::Array<int>& remotePIDs,
                                                Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& colMap) {
  using Teuchos::rcp;
  typedef LocalOrdinal LO;
  typedef GlobalOrdinal GO;
  typedef Tpetra::global_size_t GST;
  typedef Tpetra::Map<LO, GO, Node> map_type;
  const char prefix[] = "lowCommunicationMakeColMapAndReindexSerial: ";

  // The domainMap is an RCP because there is a shortcut for a
  // (common) special case to return the columnMap = domainMap.
  const map_type& domainMap = *domainMapRCP;

  // Scan all column indices and sort into two groups:
  // Local:  those whose GID matches a GID of the domain map on this processor and
  // Remote: All others.
  const size_t numDomainElements = domainMap.getLocalNumElements();
  Teuchos::Array<bool> LocalGIDs;
  if (numDomainElements > 0) {
    LocalGIDs.resize(numDomainElements, false);  // Assume domain GIDs are not local
  }

  // In principle it is good to have RemoteGIDs and RemotGIDList be as
  // long as the number of remote GIDs on this processor, but this
  // would require two passes through the column IDs, so we make it
  // the max of 100 and the number of block rows.
  //
  // FIXME (mfh 11 Feb 2015) Tpetra::Details::HashTable can hold at
  // most INT_MAX entries, but it's possible to have more rows than
  // that (if size_t is 64 bits and int is 32 bits).
  const size_t numMyRows = rowptr.size() - 1;
  const int hashsize     = std::max(static_cast<int>(numMyRows), 100);

  Tpetra::Details::HashTable<GO, LO> RemoteGIDs(hashsize);
  Teuchos::Array<GO> RemoteGIDList;
  RemoteGIDList.reserve(hashsize);
  Teuchos::Array<int> PIDList;
  PIDList.reserve(hashsize);

  // Here we start using the *LocalOrdinal* colind_LID array.  This is
  // safe even if both columnIndices arrays are actually the same
  // (because LocalOrdinal==GO).  For *local* GID's set
  // colind_LID with with their LID in the domainMap.  For *remote*
  // GIDs, we set colind_LID with (numDomainElements+NumRemoteColGIDs)
  // before the increment of the remote count.  These numberings will
  // be separate because no local LID is greater than
  // numDomainElements.

  size_t NumLocalColGIDs = 0;
  LO NumRemoteColGIDs    = 0;
  for (size_t i = 0; i < numMyRows; ++i) {
    for (size_t j = rowptr[i]; j < rowptr[i + 1]; ++j) {
      const GO GID = colind_GID[j];
      // Check if GID matches a row GID
      const LO LID = domainMap.getLocalElement(GID);
      if (LID != -1) {
        const bool alreadyFound = LocalGIDs[LID];
        if (!alreadyFound) {
          LocalGIDs[LID] = true;  // There is a column in the graph associated with this domain map GID
          NumLocalColGIDs++;
        }
        colind_LID[j] = LID;
      } else {
        const LO hash_value = RemoteGIDs.get(GID);
        if (hash_value == -1) {  // This means its a new remote GID
          const int PID = owningPIDs[j];
          TEUCHOS_TEST_FOR_EXCEPTION(
              PID == -1, std::invalid_argument, prefix << "Cannot figure out if "
                                                          "PID is owned.");
          colind_LID[j] = static_cast<LO>(numDomainElements + NumRemoteColGIDs);
          RemoteGIDs.add(GID, NumRemoteColGIDs);
          RemoteGIDList.push_back(GID);
          PIDList.push_back(PID);
          NumRemoteColGIDs++;
        } else {
          colind_LID[j] = static_cast<LO>(numDomainElements + hash_value);
        }
      }
    }
  }

  // Possible short-circuit: If all domain map GIDs are present as
  // column indices, then set ColMap=domainMap and quit.
  if (domainMap.getComm()->getSize() == 1) {
    // Sanity check: When there is only one process, there can be no
    // remoteGIDs.
    TEUCHOS_TEST_FOR_EXCEPTION(
        NumRemoteColGIDs != 0, std::runtime_error, prefix << "There is only one "
                                                             "process in the domain Map's communicator, which means that there are no "
                                                             "\"remote\" indices.  Nevertheless, some column indices are not in the "
                                                             "domain Map.");
    if (static_cast<size_t>(NumLocalColGIDs) == numDomainElements) {
      // In this case, we just use the domainMap's indices, which is,
      // not coincidently, what we clobbered colind with up above
      // anyway.  No further reindexing is needed.
      colMap = domainMapRCP;
      return;
    }
  }

  // Now build the array containing column GIDs
  // Build back end, containing remote GIDs, first
  const LO numMyCols = NumLocalColGIDs + NumRemoteColGIDs;
  Teuchos::Array<GO> ColIndices;
  GO* RemoteColIndices = NULL;
  if (numMyCols > 0) {
    ColIndices.resize(numMyCols);
    if (NumLocalColGIDs != static_cast<size_t>(numMyCols)) {
      RemoteColIndices = &ColIndices[NumLocalColGIDs];  // Points to back half of ColIndices
    }
  }

  for (LO i = 0; i < NumRemoteColGIDs; ++i) {
    RemoteColIndices[i] = RemoteGIDList[i];
  }

  // Build permute array for *remote* reindexing.
  Teuchos::Array<LO> RemotePermuteIDs(NumRemoteColGIDs);
  for (LO i = 0; i < NumRemoteColGIDs; ++i) {
    RemotePermuteIDs[i] = i;
  }

  // Sort External column indices so that all columns coming from a
  // given remote processor are contiguous.  This is a sort with two
  // auxilary arrays: RemoteColIndices and RemotePermuteIDs.
  Tpetra::sort3(PIDList.begin(), PIDList.end(),
                ColIndices.begin() + NumLocalColGIDs,
                RemotePermuteIDs.begin());

  // Stash the RemotePIDs.
  //
  // Note: If Teuchos::Array had a shrink_to_fit like std::vector,
  // we'd call it here.
  remotePIDs = PIDList;

  // Sort external column indices so that columns from a given remote
  // processor are not only contiguous but also in ascending
  // order. NOTE: I don't know if the number of externals associated
  // with a given remote processor is known at this point ... so I
  // count them here.

  // NTS: Only sort the RemoteColIndices this time...
  LO StartCurrent = 0, StartNext = 1;
  while (StartNext < NumRemoteColGIDs) {
    if (PIDList[StartNext] == PIDList[StartNext - 1]) {
      StartNext++;
    } else {
      Tpetra::sort2(ColIndices.begin() + NumLocalColGIDs + StartCurrent,
                    ColIndices.begin() + NumLocalColGIDs + StartNext,
                    RemotePermuteIDs.begin() + StartCurrent);
      StartCurrent = StartNext;
      StartNext++;
    }
  }
  Tpetra::sort2(ColIndices.begin() + NumLocalColGIDs + StartCurrent,
                ColIndices.begin() + NumLocalColGIDs + StartNext,
                RemotePermuteIDs.begin() + StartCurrent);

  // Reverse the permutation to get the information we actually care about
  Teuchos::Array<LO> ReverseRemotePermuteIDs(NumRemoteColGIDs);
  for (LO i = 0; i < NumRemoteColGIDs; ++i) {
    ReverseRemotePermuteIDs[RemotePermuteIDs[i]] = i;
  }

  // Build permute array for *local* reindexing.
  bool use_local_permute = false;
  Teuchos::Array<LO> LocalPermuteIDs(numDomainElements);

  // Now fill front end. Two cases:
  //
  // (1) If the number of Local column GIDs is the same as the number
  //     of Local domain GIDs, we can simply read the domain GIDs into
  //     the front part of ColIndices, otherwise
  //
  // (2) We step through the GIDs of the domainMap, checking to see if
  //     each domain GID is a column GID.  we want to do this to
  //     maintain a consistent ordering of GIDs between the columns
  //     and the domain.
  Teuchos::ArrayView<const GO> domainGlobalElements = domainMap.getLocalElementList();
  if (static_cast<size_t>(NumLocalColGIDs) == numDomainElements) {
    if (NumLocalColGIDs > 0) {
      // Load Global Indices into first numMyCols elements column GID list
      std::copy(domainGlobalElements.begin(), domainGlobalElements.end(),
                ColIndices.begin());
    }
  } else {
    LO NumLocalAgain  = 0;
    use_local_permute = true;
    for (size_t i = 0; i < numDomainElements; ++i) {
      if (LocalGIDs[i]) {
        LocalPermuteIDs[i]          = NumLocalAgain;
        ColIndices[NumLocalAgain++] = domainGlobalElements[i];
      }
    }
    TEUCHOS_TEST_FOR_EXCEPTION(
        static_cast<size_t>(NumLocalAgain) != NumLocalColGIDs,
        std::runtime_error, prefix << "Local ID count test failed.");
  }

  // Make column Map
  const GST minus_one = Teuchos::OrdinalTraits<GST>::invalid();
  colMap              = rcp(new map_type(minus_one, ColIndices, domainMap.getIndexBase(),
                                         domainMap.getComm()));

  // Low-cost reindex of the matrix
  for (size_t i = 0; i < numMyRows; ++i) {
    for (size_t j = rowptr[i]; j < rowptr[i + 1]; ++j) {
      const LO ID = colind_LID[j];
      if (static_cast<size_t>(ID) < numDomainElements) {
        if (use_local_permute) {
          colind_LID[j] = LocalPermuteIDs[colind_LID[j]];
        }
        // In the case where use_local_permute==false, we just copy
        // the DomainMap's ordering, which it so happens is what we
        // put in colind_LID to begin with.
      } else {
        colind_LID[j] = NumLocalColGIDs + ReverseRemotePermuteIDs[colind_LID[j] - numDomainElements];
      }
    }
  }
}

template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void lowCommunicationMakeColMapAndReindex(
    const Teuchos::ArrayView<const size_t>& rowptr,
    const Teuchos::ArrayView<LocalOrdinal>& colind_LID,
    const Teuchos::ArrayView<GlobalOrdinal>& colind_GID,
    const Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& domainMapRCP,
    const Teuchos::ArrayView<const int>& owningPIDs,
    Teuchos::Array<int>& remotePIDs,
    Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& colMap) {
  using Teuchos::rcp;
  typedef LocalOrdinal LO;
  typedef GlobalOrdinal GO;
  typedef Tpetra::global_size_t GST;
  typedef Tpetra::Map<LO, GO, Node> map_type;
  const char prefix[] = "lowCommunicationMakeColMapAndReindex: ";

  typedef typename Node::device_type DT;
  using execution_space = typename DT::execution_space;
  execution_space exec;
  using team_policy = Kokkos::TeamPolicy<execution_space, Kokkos::Schedule<Kokkos::Dynamic>>;
  typedef typename map_type::local_map_type local_map_type;

  // Create device mirror and host mirror views from function parameters
  // When we pass in views instead of Teuchos::ArrayViews, we can avoid copying views
  auto colind_LID_view = Details::create_mirror_view_from_raw_host_array(exec, colind_LID.getRawPtr(), colind_LID.size(), true, "colind_LID");
  auto rowptr_view     = Details::create_mirror_view_from_raw_host_array(exec, rowptr.getRawPtr(), rowptr.size(), true, "rowptr");
  auto colind_GID_view = Details::create_mirror_view_from_raw_host_array(exec, colind_GID.getRawPtr(), colind_GID.size(), true, "colind_GID");
  auto owningPIDs_view = Details::create_mirror_view_from_raw_host_array(exec, owningPIDs.getRawPtr(), owningPIDs.size(), true, "owningPIDs");

  typename decltype(colind_LID_view)::host_mirror_type colind_LID_host(colind_LID.getRawPtr(), colind_LID.size());
  typename decltype(colind_GID_view)::host_mirror_type colind_GID_host(colind_GID.getRawPtr(), colind_GID.size());

  Kokkos::deep_copy(colind_LID_view, colind_LID_host);
  Kokkos::deep_copy(colind_GID_view, colind_GID_host);

  // The domainMap is an RCP because there is a shortcut for a
  // (common) special case to return the columnMap = domainMap.
  const map_type& domainMap = *domainMapRCP;

  Kokkos::UnorderedMap<LO, bool, DT> LocalGIDs_view_map(colind_LID.size());
  Kokkos::UnorderedMap<GO, LO, DT> RemoteGIDs_view_map(colind_LID.size());

  const size_t numMyRows         = rowptr.size() - 1;
  local_map_type domainMap_local = domainMap.getLocalMap();

  const size_t numDomainElements = domainMap.getLocalNumElements();
  Kokkos::View<bool*, DT> LocalGIDs_view("LocalGIDs", numDomainElements);
  auto LocalGIDs_host = Kokkos::create_mirror_view(LocalGIDs_view);

  size_t NumLocalColGIDs = 0;

  // Scan all column indices and sort into two groups:
  // Local:  those whose GID matches a GID of the domain map on this processor and
  // Remote: All others.
  // Kokkos::Parallel_reduce sums up NumLocalColGIDs, while we use the size of the Remote GIDs map to find NumRemoteColGIDs
  Kokkos::parallel_reduce(
      team_policy(numMyRows, Kokkos::AUTO), KOKKOS_LAMBDA(const typename team_policy::member_type& member, size_t& update) {
        const int i                 = member.league_rank();
        size_t NumLocalColGIDs_temp = 0;
        size_t rowptr_start         = rowptr_view[i];
        size_t rowptr_end           = rowptr_view[i + 1];
        Kokkos::parallel_reduce(
            Kokkos::TeamThreadRange(member, rowptr_start, rowptr_end), [&](const size_t j, size_t& innerUpdate) {
              const GO GID = colind_GID_view[j];
              // Check if GID matches a row GID in local domain map
              const LO LID = domainMap_local.getLocalElement(GID);
              if (LID != -1) {
                auto outcome = LocalGIDs_view_map.insert(LID);
                // Fresh insert
                if (outcome.success()) {
                  LocalGIDs_view[LID] = true;
                  innerUpdate++;
                }
              } else {
                const int PID = owningPIDs_view[j];
                auto outcome  = RemoteGIDs_view_map.insert(GID, PID);
                if (outcome.success() && PID == -1) {
                  Kokkos::abort("Cannot figure out if ID is owned.\n");
                }
              }
            },
            NumLocalColGIDs_temp);
        if (member.team_rank() == 0) update += NumLocalColGIDs_temp;
      },
      NumLocalColGIDs);

  LO NumRemoteColGIDs = RemoteGIDs_view_map.size();

  Kokkos::View<int*, DT> PIDList_view("PIDList", NumRemoteColGIDs);
  auto PIDList_host = Kokkos::create_mirror_view(PIDList_view);

  Kokkos::View<GO*, DT> RemoteGIDList_view("RemoteGIDList", NumRemoteColGIDs);
  auto RemoteGIDList_host = Kokkos::create_mirror_view(RemoteGIDList_view);

  // For each index in RemoteGIDs_map that contains a GID, use "update" to indicate the number of GIDs "before" this GID
  // This maps each element in the RemoteGIDs hash table to an index in RemoteGIDList / PIDList without any overwriting or empty spaces between indices
  Kokkos::parallel_scan(
      Kokkos::RangePolicy<execution_space>(0, RemoteGIDs_view_map.capacity()), KOKKOS_LAMBDA(const int i, GO& update, const bool final) {
        if (final && RemoteGIDs_view_map.valid_at(i)) {
          RemoteGIDList_view[update] = RemoteGIDs_view_map.key_at(i);
          PIDList_view[update]       = RemoteGIDs_view_map.value_at(i);
        }
        if (RemoteGIDs_view_map.valid_at(i)) {
          update += 1;
        }
      });

  // Possible short-circuit: If all domain map GIDs are present as
  // column indices, then set ColMap=domainMap and quit.
  if (domainMap.getComm()->getSize() == 1) {
    // Sanity check: When there is only one process, there can be no
    // remoteGIDs.
    TEUCHOS_TEST_FOR_EXCEPTION(
        NumRemoteColGIDs != 0, std::runtime_error, prefix << "There is only one "
                                                             "process in the domain Map's communicator, which means that there are no "
                                                             "\"remote\" indices.  Nevertheless, some column indices are not in the "
                                                             "domain Map.");
    if (static_cast<size_t>(NumLocalColGIDs) == numDomainElements) {
      // In this case, we just use the domainMap's indices, which is,
      // not coincidently, what we clobbered colind with up above
      // anyway.  No further reindexing is needed.
      colMap = domainMapRCP;

      // Fill out local colMap (which should only contain local GIDs)
      auto localColMap = colMap->getLocalMap();
      Kokkos::parallel_for(
          Kokkos::RangePolicy<execution_space>(0, colind_GID.size()), KOKKOS_LAMBDA(const int i) {
            colind_LID_view[i] = localColMap.getLocalElement(colind_GID_view[i]);
          });
      Kokkos::deep_copy(execution_space(), colind_LID_host, colind_LID_view);
      return;
    }
  }

  // Now build the array containing column GIDs
  // Build back end, containing remote GIDs, first
  const LO numMyCols = NumLocalColGIDs + NumRemoteColGIDs;
  Kokkos::View<GO*, DT> ColIndices_view("ColIndices", numMyCols);

  // We don't need to load the backend of ColIndices or sort if there are no remote GIDs
  if (NumRemoteColGIDs > 0) {
    if (NumLocalColGIDs != static_cast<size_t>(numMyCols)) {
      Kokkos::parallel_for(
          Kokkos::RangePolicy<execution_space>(0, NumRemoteColGIDs), KOKKOS_LAMBDA(const int i) {
            ColIndices_view[NumLocalColGIDs + i] = RemoteGIDList_view[i];
          });
    }

    // Find the largest PID for bin sorting purposes
    int PID_max = 0;
    Kokkos::parallel_reduce(
        Kokkos::RangePolicy<execution_space>(0, PIDList_host.size()), KOKKOS_LAMBDA(const int i, int& max) {
          if (max < PIDList_view[i]) max = PIDList_view[i];
        },
        Kokkos::Max<int>(PID_max));

    using KeyViewTypePID = decltype(PIDList_view);
    using BinSortOpPID   = Kokkos::BinOp1D<KeyViewTypePID>;

    // Make a subview of ColIndices for remote GID sorting
    auto ColIndices_subview = Kokkos::subview(ColIndices_view, Kokkos::make_pair(NumLocalColGIDs, ColIndices_view.size()));

    // Make binOp with bins = PID_max + 1, min = 0, max = PID_max
    BinSortOpPID binOp2(PID_max + 1, 0, PID_max);

    // Sort External column indices so that all columns coming from a
    // given remote processor are contiguous.  This is a sort with one
    // auxilary array: RemoteColIndices
    Kokkos::BinSort<KeyViewTypePID, BinSortOpPID> bin_sort2(PIDList_view, 0, PIDList_view.size(), binOp2, false);
    bin_sort2.create_permute_vector(exec);
    bin_sort2.sort(exec, PIDList_view);
    bin_sort2.sort(exec, ColIndices_subview);

    // Deep copy back from device to host
    Kokkos::deep_copy(exec, PIDList_host, PIDList_view);

    // Stash the RemotePIDs. Once remotePIDs is changed to become a Kokkos view, we can remove this and copy directly.
    // Note: If Teuchos::Array had a shrink_to_fit like std::vector,
    // we'd call it here.

    exec.fence("fence before setting PIDList");
    Teuchos::Array<int> PIDList(NumRemoteColGIDs);
    for (LO i = 0; i < NumRemoteColGIDs; ++i) {
      PIDList[i] = PIDList_host[i];
    }

    remotePIDs = PIDList;

    // Sort external column indices so that columns from a given remote
    // processor are not only contiguous but also in ascending
    // order. NOTE: I don't know if the number of externals associated
    // with a given remote processor is known at this point ... so I
    // count them here.
    LO StartCurrent = 0, StartNext = 1;
    while (StartNext < NumRemoteColGIDs) {
      if (PIDList_host[StartNext] == PIDList_host[StartNext - 1]) {
        StartNext++;
      } else {
        Kokkos::sort(ColIndices_view, NumLocalColGIDs + StartCurrent, NumLocalColGIDs + StartNext);
        StartCurrent = StartNext;
        StartNext++;
      }
    }

    Kokkos::sort(ColIndices_view, NumLocalColGIDs + StartCurrent, NumLocalColGIDs + StartNext);
  }

  // Build permute array for *local* reindexing.

  // Now fill front end. Two cases:
  //
  // (1) If the number of Local column GIDs is the same as the number
  //     of Local domain GIDs, we can simply read the domain GIDs into
  //     the front part of ColIndices, otherwise
  //
  // (2) We step through the GIDs of the domainMap, checking to see if
  //     each domain GID is a column GID.  we want to do this to
  //     maintain a consistent ordering of GIDs between the columns
  //     and the domain.
  if (static_cast<size_t>(NumLocalColGIDs) == numDomainElements) {
    if (NumLocalColGIDs > 0) {
      // Load Global Indices into first numMyCols elements column GID list
      Kokkos::parallel_for(
          Kokkos::RangePolicy<execution_space>(0, numDomainElements), KOKKOS_LAMBDA(const int i) {
            ColIndices_view[i] = domainMap_local.getGlobalElement(i);
          });
    }
  } else {
    // This part isn't actually tested in the unit tests
    LO NumLocalAgain = 0;
    Kokkos::parallel_scan(
        Kokkos::RangePolicy<execution_space>(0, numDomainElements), KOKKOS_LAMBDA(const int i, LO& update, const bool final) {
          if (final && LocalGIDs_view[i]) {
            ColIndices_view[update] = domainMap_local.getGlobalElement(i);
          }
          if (LocalGIDs_view[i]) {
            update++;
          }
        },
        NumLocalAgain);

    TEUCHOS_TEST_FOR_EXCEPTION(
        static_cast<size_t>(NumLocalAgain) != NumLocalColGIDs,
        std::runtime_error, prefix << "Local ID count test failed.");
  }

  // Make column Map
  const GST minus_one = Teuchos::OrdinalTraits<GST>::invalid();

  colMap = rcp(new map_type(minus_one, ColIndices_view, domainMap.getIndexBase(),
                            domainMap.getComm()));

  // Fill out colind_LID using local map
  auto localColMap = colMap->getLocalMap();
  Kokkos::parallel_for(
      Kokkos::RangePolicy<execution_space>(0, colind_GID.size()), KOKKOS_LAMBDA(const int i) {
        colind_LID_view[i] = localColMap.getLocalElement(colind_GID_view[i]);
      });

  // For now, we copy back into colind_LID_host (which also overwrites the colind_LID Tuechos array)
  // When colind_LID becomes a Kokkos View we can delete this
  Kokkos::deep_copy(exec, colind_LID_host, colind_LID_view);
}

template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void lowCommunicationMakeColMapAndReindex(
    const Kokkos::View<size_t*, typename Node::device_type> rowptr_view,
    const Kokkos::View<LocalOrdinal*, typename Node::device_type> colind_LID_view,
    const Kokkos::View<GlobalOrdinal*, typename Node::device_type> colind_GID_view,
    const Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& domainMapRCP,
    const Kokkos::View<int*, typename Node::device_type> owningPIDs_view,
    Teuchos::Array<int>& remotePIDs,
    Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>& colMap) {
  using Teuchos::rcp;
  typedef LocalOrdinal LO;
  typedef GlobalOrdinal GO;
  typedef Tpetra::global_size_t GST;
  typedef Tpetra::Map<LO, GO, Node> map_type;
  const char prefix[] = "lowCommunicationMakeColMapAndReindex: ";

  typedef typename Node::device_type DT;
  using execution_space = typename DT::execution_space;
  execution_space exec;
  using team_policy = Kokkos::TeamPolicy<execution_space, Kokkos::Schedule<Kokkos::Dynamic>>;
  typedef typename map_type::local_map_type local_map_type;

  // The domainMap is an RCP because there is a shortcut for a
  // (common) special case to return the columnMap = domainMap.
  const map_type& domainMap = *domainMapRCP;

  Kokkos::UnorderedMap<LO, bool, DT> LocalGIDs_view_map(colind_LID_view.size());
  Kokkos::UnorderedMap<GO, LO, DT> RemoteGIDs_view_map(colind_LID_view.size());

  const size_t numMyRows         = rowptr_view.size() - 1;
  local_map_type domainMap_local = domainMap.getLocalMap();

  const size_t numDomainElements = domainMap.getLocalNumElements();
  Kokkos::View<bool*, DT> LocalGIDs_view("LocalGIDs", numDomainElements);
  auto LocalGIDs_host = Kokkos::create_mirror_view(LocalGIDs_view);

  size_t NumLocalColGIDs = 0;

  // Scan all column indices and sort into two groups:
  // Local:  those whose GID matches a GID of the domain map on this processor and
  // Remote: All others.
  // Kokkos::Parallel_reduce sums up NumLocalColGIDs, while we use the size of the Remote GIDs map to find NumRemoteColGIDs
  Kokkos::parallel_reduce(
      team_policy(numMyRows, Kokkos::AUTO), KOKKOS_LAMBDA(const typename team_policy::member_type& member, size_t& update) {
        const int i                 = member.league_rank();
        size_t NumLocalColGIDs_temp = 0;
        size_t rowptr_start         = rowptr_view[i];
        size_t rowptr_end           = rowptr_view[i + 1];
        Kokkos::parallel_reduce(
            Kokkos::TeamThreadRange(member, rowptr_start, rowptr_end), [&](const size_t j, size_t& innerUpdate) {
              const GO GID = colind_GID_view[j];
              // Check if GID matches a row GID in local domain map
              const LO LID = domainMap_local.getLocalElement(GID);
              if (LID != -1) {
                auto outcome = LocalGIDs_view_map.insert(LID);
                // Fresh insert
                if (outcome.success()) {
                  LocalGIDs_view[LID] = true;
                  innerUpdate++;
                }
              } else {
                const int PID = owningPIDs_view[j];
                auto outcome  = RemoteGIDs_view_map.insert(GID, PID);
                if (outcome.success() && PID == -1) {
                  Kokkos::abort("Cannot figure out if ID is owned.\n");
                }
              }
            },
            NumLocalColGIDs_temp);
        if (member.team_rank() == 0) update += NumLocalColGIDs_temp;
      },
      NumLocalColGIDs);

  LO NumRemoteColGIDs = RemoteGIDs_view_map.size();

  Kokkos::View<int*, DT> PIDList_view("PIDList_d", NumRemoteColGIDs);

  Kokkos::View<GO*, DT> RemoteGIDList_view("RemoteGIDList", NumRemoteColGIDs);
  auto RemoteGIDList_host = Kokkos::create_mirror_view(RemoteGIDList_view);

  // For each index in RemoteGIDs_map that contains a GID, use "update" to indicate the number of GIDs "before" this GID
  // This maps each element in the RemoteGIDs hash table to an index in RemoteGIDList / PIDList without any overwriting or empty spaces between indices
  Kokkos::parallel_scan(
      Kokkos::RangePolicy<execution_space>(0, RemoteGIDs_view_map.capacity()), KOKKOS_LAMBDA(const int i, GO& update, const bool final) {
        if (final && RemoteGIDs_view_map.valid_at(i)) {
          RemoteGIDList_view[update] = RemoteGIDs_view_map.key_at(i);
          PIDList_view[update]       = RemoteGIDs_view_map.value_at(i);
        }
        if (RemoteGIDs_view_map.valid_at(i)) {
          update += 1;
        }
      });

  // Possible short-circuit: If all domain map GIDs are present as
  // column indices, then set ColMap=domainMap and quit.
  if (domainMap.getComm()->getSize() == 1) {
    // Sanity check: When there is only one process, there can be no
    // remoteGIDs.
    TEUCHOS_TEST_FOR_EXCEPTION(
        NumRemoteColGIDs != 0, std::runtime_error, prefix << "There is only one "
                                                             "process in the domain Map's communicator, which means that there are no "
                                                             "\"remote\" indices.  Nevertheless, some column indices are not in the "
                                                             "domain Map.");
    if (static_cast<size_t>(NumLocalColGIDs) == numDomainElements) {
      // In this case, we just use the domainMap's indices, which is,
      // not coincidently, what we clobbered colind with up above
      // anyway.  No further reindexing is needed.
      colMap = domainMapRCP;

      // Fill out local colMap (which should only contain local GIDs)
      auto localColMap = colMap->getLocalMap();
      Kokkos::parallel_for(
          Kokkos::RangePolicy<execution_space>(0, colind_GID_view.size()), KOKKOS_LAMBDA(const int i) {
            colind_LID_view[i] = localColMap.getLocalElement(colind_GID_view[i]);
          });
      return;
    }
  }

  // Now build the array containing column GIDs
  // Build back end, containing remote GIDs, first
  const LO numMyCols = NumLocalColGIDs + NumRemoteColGIDs;
  Kokkos::View<GO*, DT> ColIndices_view("ColIndices", numMyCols);

  // We don't need to load the backend of ColIndices or sort if there are no remote GIDs
  if (NumRemoteColGIDs > 0) {
    if (NumLocalColGIDs != static_cast<size_t>(numMyCols)) {
      Kokkos::parallel_for(
          Kokkos::RangePolicy<execution_space>(0, NumRemoteColGIDs), KOKKOS_LAMBDA(const int i) {
            ColIndices_view[NumLocalColGIDs + i] = RemoteGIDList_view[i];
          });
    }

    // Find the largest PID for bin sorting purposes
    int PID_max = 0;
    Kokkos::parallel_reduce(
        Kokkos::RangePolicy<execution_space>(0, PIDList_view.size()), KOKKOS_LAMBDA(const int i, int& max) {
          if (max < PIDList_view[i]) max = PIDList_view[i];
        },
        Kokkos::Max<int>(PID_max));

    using KeyViewTypePID = decltype(PIDList_view);
    using BinSortOpPID   = Kokkos::BinOp1D<KeyViewTypePID>;

    // Make a subview of ColIndices for remote GID sorting
    auto ColIndices_subview = Kokkos::subview(ColIndices_view, Kokkos::make_pair(NumLocalColGIDs, ColIndices_view.size()));

    // Make binOp with bins = PID_max + 1, min = 0, max = PID_max
    BinSortOpPID binOp2(PID_max + 1, 0, PID_max);

    // Sort External column indices so that all columns coming from a
    // given remote processor are contiguous.  This is a sort with one
    // auxilary array: RemoteColIndices
    Kokkos::BinSort<KeyViewTypePID, BinSortOpPID> bin_sort2(PIDList_view, 0, PIDList_view.size(), binOp2, false);
    bin_sort2.create_permute_vector(exec);
    bin_sort2.sort(exec, PIDList_view);
    bin_sort2.sort(exec, ColIndices_subview);

    // Deep copy back from device to host
    // Stash the RemotePIDs. Once remotePIDs is changed to become a Kokkos view, we can remove this and copy directly.
    Teuchos::Array<int> PIDList(NumRemoteColGIDs);
    Kokkos::View<int*, Kokkos::HostSpace> PIDList_host(PIDList.data(), PIDList.size());
    Kokkos::deep_copy(exec, PIDList_host, PIDList_view);
    exec.fence();

    remotePIDs = PIDList;

    // Sort external column indices so that columns from a given remote
    // processor are not only contiguous but also in ascending
    // order. NOTE: I don't know if the number of externals associated
    // with a given remote processor is known at this point ... so I
    // count them here.
    LO StartCurrent = 0, StartNext = 1;
    while (StartNext < NumRemoteColGIDs) {
      if (PIDList_host[StartNext] == PIDList_host[StartNext - 1]) {
        StartNext++;
      } else {
        Kokkos::sort(ColIndices_view, NumLocalColGIDs + StartCurrent, NumLocalColGIDs + StartNext);
        StartCurrent = StartNext;
        StartNext++;
      }
    }

    Kokkos::sort(ColIndices_view, NumLocalColGIDs + StartCurrent, NumLocalColGIDs + StartNext);
  }

  // Build permute array for *local* reindexing.

  // Now fill front end. Two cases:
  //
  // (1) If the number of Local column GIDs is the same as the number
  //     of Local domain GIDs, we can simply read the domain GIDs into
  //     the front part of ColIndices, otherwise
  //
  // (2) We step through the GIDs of the domainMap, checking to see if
  //     each domain GID is a column GID.  we want to do this to
  //     maintain a consistent ordering of GIDs between the columns
  //     and the domain.
  if (static_cast<size_t>(NumLocalColGIDs) == numDomainElements) {
    if (NumLocalColGIDs > 0) {
      // Load Global Indices into first numMyCols elements column GID list
      Kokkos::parallel_for(
          Kokkos::RangePolicy<execution_space>(0, numDomainElements), KOKKOS_LAMBDA(const int i) {
            ColIndices_view[i] = domainMap_local.getGlobalElement(i);
          });
    }
  } else {
    // This part isn't actually tested in the unit tests
    LO NumLocalAgain = 0;
    Kokkos::parallel_scan(
        Kokkos::RangePolicy<execution_space>(0, numDomainElements), KOKKOS_LAMBDA(const int i, LO& update, const bool final) {
          if (final && LocalGIDs_view[i]) {
            ColIndices_view[update] = domainMap_local.getGlobalElement(i);
          }
          if (LocalGIDs_view[i]) {
            update++;
          }
        },
        NumLocalAgain);

    TEUCHOS_TEST_FOR_EXCEPTION(
        static_cast<size_t>(NumLocalAgain) != NumLocalColGIDs,
        std::runtime_error, prefix << "Local ID count test failed.");
  }

  // Make column Map
  const GST minus_one = Teuchos::OrdinalTraits<GST>::invalid();

  colMap = rcp(new map_type(minus_one, ColIndices_view, domainMap.getIndexBase(),
                            domainMap.getComm()));

  // Fill out colind_LID using local map
  auto localColMap = colMap->getLocalMap();
  Kokkos::parallel_for(
      Kokkos::RangePolicy<execution_space>(0, colind_GID_view.size()), KOKKOS_LAMBDA(const int i) {
        colind_LID_view[i] = localColMap.getLocalElement(colind_GID_view[i]);
      });
}

// Generates an list of owning PIDs based on two transfer (aka import/export objects)
// Let:
//   OwningMap = useReverseModeForOwnership ? transferThatDefinesOwnership.getTargetMap() : transferThatDefinesOwnership.getSourceMap();
//   MapAo     = useReverseModeForOwnership ? transferThatDefinesOwnership.getSourceMap() : transferThatDefinesOwnership.getTargetMap();
//   MapAm     = useReverseModeForMigration ? transferThatDefinesMigration.getTargetMap() : transferThatDefinesMigration.getSourceMap();
//   VectorMap = useReverseModeForMigration ? transferThatDefinesMigration.getSourceMap() : transferThatDefinesMigration.getTargetMap();
// Precondition:
//  1) MapAo.isSameAs(*MapAm)                      - map compatibility between transfers
//  2) VectorMap->isSameAs(*owningPIDs->getMap())  - map compabibility between transfer & vector
//  3) OwningMap->isOneToOne()                     - owning map is 1-to-1
//  --- Precondition 3 is only checked in DEBUG mode ---
// Postcondition:
//   owningPIDs[VectorMap->getLocalElement(GID i)] =   j iff  (OwningMap->isLocalElement(GID i) on rank j)
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void getTwoTransferOwnershipVector(const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& transferThatDefinesOwnership,
                                   bool useReverseModeForOwnership,
                                   const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& transferThatDefinesMigration,
                                   bool useReverseModeForMigration,
                                   Tpetra::Vector<int, LocalOrdinal, GlobalOrdinal, Node>& owningPIDs) {
  typedef Tpetra::Import<LocalOrdinal, GlobalOrdinal, Node> import_type;
  typedef Tpetra::Export<LocalOrdinal, GlobalOrdinal, Node> export_type;

  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> OwningMap = useReverseModeForOwnership ? transferThatDefinesOwnership.getTargetMap() : transferThatDefinesOwnership.getSourceMap();
  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> MapAo     = useReverseModeForOwnership ? transferThatDefinesOwnership.getSourceMap() : transferThatDefinesOwnership.getTargetMap();
  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> MapAm     = useReverseModeForMigration ? transferThatDefinesMigration.getTargetMap() : transferThatDefinesMigration.getSourceMap();
  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> VectorMap = useReverseModeForMigration ? transferThatDefinesMigration.getSourceMap() : transferThatDefinesMigration.getTargetMap();

  TEUCHOS_TEST_FOR_EXCEPTION(!MapAo->isSameAs(*MapAm), std::runtime_error, "Tpetra::Import_Util::getTwoTransferOwnershipVector map mismatch between transfers");
  TEUCHOS_TEST_FOR_EXCEPTION(!VectorMap->isSameAs(*owningPIDs.getMap()), std::runtime_error, "Tpetra::Import_Util::getTwoTransferOwnershipVector map mismatch transfer and vector");
#ifdef HAVE_TPETRA_DEBUG
  TEUCHOS_TEST_FOR_EXCEPTION(!OwningMap->isOneToOne(), std::runtime_error, "Tpetra::Import_Util::getTwoTransferOwnershipVector owner must be 1-to-1");
#endif

  int rank = OwningMap->getComm()->getRank();
  // Generate "A" vector and fill it with owning information.  We can read this from transferThatDefinesOwnership w/o communication
  // Note:  Due to the 1-to-1 requirement, several of these options throw
  Tpetra::Vector<int, LocalOrdinal, GlobalOrdinal, Node> temp(MapAo);
  const import_type* ownAsImport = dynamic_cast<const import_type*>(&transferThatDefinesOwnership);
  const export_type* ownAsExport = dynamic_cast<const export_type*>(&transferThatDefinesOwnership);

  Teuchos::ArrayRCP<int> pids    = temp.getDataNonConst();
  Teuchos::ArrayView<int> v_pids = pids();
  if (ownAsImport && useReverseModeForOwnership) {
    TEUCHOS_TEST_FOR_EXCEPTION(1, std::runtime_error, "Tpetra::Import_Util::getTwoTransferOwnershipVector owner must be 1-to-1");
  } else if (ownAsImport && !useReverseModeForOwnership)
    getPids(*ownAsImport, v_pids, false);
  else if (ownAsExport && useReverseModeForMigration) {
    TEUCHOS_TEST_FOR_EXCEPTION(1, std::runtime_error, "Tpetra::Import_Util::getTwoTransferOwnershipVector this option not yet implemented");
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(1, std::runtime_error, "Tpetra::Import_Util::getTwoTransferOwnershipVector owner must be 1-to-1");
  }

  const import_type* xferAsImport = dynamic_cast<const import_type*>(&transferThatDefinesMigration);
  const export_type* xferAsExport = dynamic_cast<const export_type*>(&transferThatDefinesMigration);
  TEUCHOS_TEST_FOR_EXCEPTION(!xferAsImport && !xferAsExport, std::runtime_error, "Tpetra::Import_Util::getTwoTransferOwnershipVector transfer undefined");

  // Migrate from "A" vector to output vector
  owningPIDs.putScalar(rank);
  if (xferAsImport && useReverseModeForMigration)
    owningPIDs.doExport(temp, *xferAsImport, Tpetra::REPLACE);
  else if (xferAsImport && !useReverseModeForMigration)
    owningPIDs.doImport(temp, *xferAsImport, Tpetra::REPLACE);
  else if (xferAsExport && useReverseModeForMigration)
    owningPIDs.doImport(temp, *xferAsExport, Tpetra::REPLACE);
  else
    owningPIDs.doExport(temp, *xferAsExport, Tpetra::REPLACE);
}

}  // namespace Import_Util
}  // namespace Tpetra

#endif  // TPETRA_IMPORT_UTIL_HPP