Tpetra_DirectoryImpl_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_DIRECTORYIMPL_DEF_HPP
#define TPETRA_DIRECTORYIMPL_DEF_HPP
 
 
#include "Teuchos_TestForException.hpp"
#include "Tpetra_Distributor.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_TieBreak.hpp"
#include "Tpetra_Util.hpp"
#include "Tpetra_Details_FixedHashTable.hpp"
#include "Teuchos_Comm.hpp"
#include <memory>
#include <sstream>
 
// FIXME (mfh 16 Apr 2013) GIANT HACK BELOW
#ifdef HAVE_TPETRACORE_MPI
#include <mpi.h>
#endif  // HAVE_TPETRACORE_MPI
// FIXME (mfh 16 Apr 2013) GIANT HACK ABOVE
 
namespace Tpetra {
namespace Details {
template <class LO, class GO, class NT>
LookupStatus
Directory<LO, GO, NT>::
    getEntries(const map_type& map,
               const Teuchos::ArrayView<const GO>& globalIDs,
               const Teuchos::ArrayView<int>& nodeIDs,
               const Teuchos::ArrayView<LO>& localIDs,
               const bool computeLIDs) const {
  // Ensure that globalIDs, nodeIDs, and localIDs (if applicable)
  // all have the same size, before modifying any output arguments.
  TEUCHOS_TEST_FOR_EXCEPTION(nodeIDs.size() != globalIDs.size(),
                             std::invalid_argument, Teuchos::typeName(*this) << "::getEntries(): "
                                                                                "Output arrays do not have the right sizes.  nodeIDs.size() = "
                                                                             << nodeIDs.size() << " != globalIDs.size() = " << globalIDs.size() << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(
      computeLIDs && localIDs.size() != globalIDs.size(),
      std::invalid_argument, Teuchos::typeName(*this) << "::getEntries(): "
                                                         "Output array do not have the right sizes.  localIDs.size() = "
                                                      << localIDs.size() << " != globalIDs.size() = " << globalIDs.size() << ".");
 
  // Initially, fill nodeIDs and localIDs (if applicable) with
  // invalid values.  The "invalid" process ID is -1 (this means
  // the same thing as MPI_ANY_SOURCE to Teuchos, so it's an
  // "invalid" process ID); the invalid local ID comes from
  // OrdinalTraits.
  std::fill(nodeIDs.begin(), nodeIDs.end(), -1);
  if (computeLIDs) {
    std::fill(localIDs.begin(), localIDs.end(),
              Teuchos::OrdinalTraits<LO>::invalid());
  }
  // Actually do the work.
  return this->getEntriesImpl(map, globalIDs, nodeIDs, localIDs, computeLIDs);
}
 
template <class LO, class GO, class NT>
ReplicatedDirectory<LO, GO, NT>::
    ReplicatedDirectory(const map_type& map)
  : numProcs_(map.getComm()->getSize()) {}
 
template <class LO, class GO, class NT>
bool ReplicatedDirectory<LO, GO, NT>::
    isOneToOne(const Teuchos::Comm<int>& /* comm */) const {
  // A locally replicated Map is one-to-one only if there is no
  // replication, that is, only if the Map's communicator only has
  // one process.
  return (numProcs_ == 1);
}
 
template <class LO, class GO, class NT>
std::string
ReplicatedDirectory<LO, GO, NT>::description() const {
  std::ostringstream os;
  os << "ReplicatedDirectory"
     << "<" << Teuchos::TypeNameTraits<LO>::name()
     << ", " << Teuchos::TypeNameTraits<GO>::name()
     << ", " << Teuchos::TypeNameTraits<NT>::name() << ">";
  return os.str();
}
 
template <class LO, class GO, class NT>
ContiguousUniformDirectory<LO, GO, NT>::
    ContiguousUniformDirectory(const map_type& map) {
  TEUCHOS_TEST_FOR_EXCEPTION(!map.isContiguous(), std::invalid_argument,
                             Teuchos::typeName(*this) << " constructor: Map is not contiguous.");
  TEUCHOS_TEST_FOR_EXCEPTION(!map.isUniform(), std::invalid_argument,
                             Teuchos::typeName(*this) << " constructor: Map is not uniform.");
}
 
template <class LO, class GO, class NT>
std::string
ContiguousUniformDirectory<LO, GO, NT>::description() const {
  std::ostringstream os;
  os << "ContiguousUniformDirectory"
     << "<" << Teuchos::TypeNameTraits<LO>::name()
     << ", " << Teuchos::TypeNameTraits<GO>::name()
     << ", " << Teuchos::TypeNameTraits<NT>::name() << ">";
  return os.str();
}
 
template <class LO, class GO, class NT>
LookupStatus
ContiguousUniformDirectory<LO, GO, NT>::
    getEntriesImpl(const map_type& map,
                   const Teuchos::ArrayView<const GO>& globalIDs,
                   const Teuchos::ArrayView<int>& nodeIDs,
                   const Teuchos::ArrayView<LO>& localIDs,
                   const bool computeLIDs) const {
  using Teuchos::Comm;
  using Teuchos::RCP;
  typedef typename Teuchos::ArrayView<const GO>::size_type size_type;
  const LO invalidLid = Teuchos::OrdinalTraits<LO>::invalid();
  LookupStatus res    = AllIDsPresent;
 
  RCP<const Comm<int> > comm = map.getComm();
  const GO g_min             = map.getMinAllGlobalIndex();
 
  // Let N_G be the global number of elements in the Map,
  // and P be the number of processes in its communicator.
  // Then, N_G = P * N_L + R = R*(N_L + 1) + (P - R)*N_L.
  //
  // The first R processes own N_L+1 elements.
  // The remaining P-R processes own N_L elements.
  //
  // Let g be the current GID, g_min be the global minimum GID,
  // and g_0 = g - g_min.  If g is a valid GID in this Map, then
  // g_0 is in [0, N_G - 1].
  //
  // If g is a valid GID in this Map and g_0 < R*(N_L + 1), then
  // the rank of the process that owns g is floor(g_0 / (N_L +
  // 1)), and its corresponding local index on that process is g_0
  // mod (N_L + 1).
  //
  // Let g_R = g_0 - R*(N_L + 1).  If g is a valid GID in this Map
  // and g_0 >= R*(N_L + 1), then the rank of the process that
  // owns g is then R + floor(g_R / N_L), and its corresponding
  // local index on that process is g_R mod N_L.
 
  const size_type N_G =
      static_cast<size_type>(map.getGlobalNumElements());
  const size_type P   = static_cast<size_type>(comm->getSize());
  const size_type N_L = N_G / P;
  const size_type R   = N_G - N_L * P;  // N_G mod P
  const size_type N_R = R * (N_L + static_cast<size_type>(1));
 
#ifdef HAVE_TPETRA_DEBUG
  TEUCHOS_TEST_FOR_EXCEPTION(
      N_G != P * N_L + R, std::logic_error,
      "Tpetra::ContiguousUniformDirectory::getEntriesImpl: "
      "N_G = "
          << N_G << " != P*N_L + R = " << P << "*" << N_L << " + " << R
          << " = " << P * N_L + R << ".  "
                                     "Please report this bug to the Tpetra developers.");
#endif  // HAVE_TPETRA_DEBUG
 
  const size_type numGids = globalIDs.size();  // for const loop bound
  // Avoid signed/unsigned comparisons below, in case GO is
  // unsigned.  (Integer literals are generally signed.)
  const GO ONE = static_cast<GO>(1);
 
  if (computeLIDs) {
    for (size_type k = 0; k < numGids; ++k) {
      const GO g_0 = globalIDs[k] - g_min;
 
      // The first test is a little strange just in case GO is
      // unsigned.  Compilers raise a warning on tests like "x <
      // 0" if x is unsigned, but don't usually raise a warning if
      // the expression is a bit more complicated than that.
      if (g_0 + ONE < ONE || g_0 >= static_cast<GO>(N_G)) {
        nodeIDs[k]  = -1;
        localIDs[k] = invalidLid;
        res         = IDNotPresent;
      } else if (g_0 < static_cast<GO>(N_R)) {
        // The GID comes from the initial sequence of R processes.
        nodeIDs[k]  = static_cast<int>(g_0 / static_cast<GO>(N_L + 1));
        localIDs[k] = static_cast<LO>(g_0 % static_cast<GO>(N_L + 1));
      } else if (g_0 >= static_cast<GO>(N_R)) {
        // The GID comes from the remaining P-R processes.
        const GO g_R = g_0 - static_cast<GO>(N_R);
        nodeIDs[k]   = static_cast<int>(R + g_R / N_L);
        localIDs[k]  = static_cast<int>(g_R % N_L);
      }
#ifdef HAVE_TPETRA_DEBUG
      else {
        TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                                   "Tpetra::ContiguousUniformDirectory::getEntriesImpl: "
                                   "should never get here.  "
                                   "Please report this bug to the Tpetra developers.");
      }
#endif  // HAVE_TPETRA_DEBUG
    }
  } else {  // don't compute local indices
    for (size_type k = 0; k < numGids; ++k) {
      const GO g_0 = globalIDs[k] - g_min;
      // The first test is a little strange just in case GO is
      // unsigned.  Compilers raise a warning on tests like "x <
      // 0" if x is unsigned, but don't usually raise a warning if
      // the expression is a bit more complicated than that.
      if (g_0 + ONE < ONE || g_0 >= static_cast<GO>(N_G)) {
        nodeIDs[k] = -1;
        res        = IDNotPresent;
      } else if (g_0 < static_cast<GO>(N_R)) {
        // The GID comes from the initial sequence of R processes.
        nodeIDs[k] = static_cast<int>(g_0 / static_cast<GO>(N_L + 1));
      } else if (g_0 >= static_cast<GO>(N_R)) {
        // The GID comes from the remaining P-R processes.
        const GO g_R = g_0 - static_cast<GO>(N_R);
        nodeIDs[k]   = static_cast<int>(R + g_R / N_L);
      }
#ifdef HAVE_TPETRA_DEBUG
      else {
        TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                                   "Tpetra::ContiguousUniformDirectory::getEntriesImpl: "
                                   "should never get here.  "
                                   "Please report this bug to the Tpetra developers.");
      }
#endif  // HAVE_TPETRA_DEBUG
    }
  }
  return res;
}
 
template <class LO, class GO, class NT>
DistributedContiguousDirectory<LO, GO, NT>::
    DistributedContiguousDirectory(const map_type& map) {
  using Teuchos::arcp;
  using Teuchos::gatherAll;
  using Teuchos::RCP;
 
  RCP<const Teuchos::Comm<int> > comm = map.getComm();
 
  TEUCHOS_TEST_FOR_EXCEPTION(!map.isDistributed(), std::invalid_argument,
                             Teuchos::typeName(*this) << " constructor: Map is not distributed.");
  TEUCHOS_TEST_FOR_EXCEPTION(!map.isContiguous(), std::invalid_argument,
                             Teuchos::typeName(*this) << " constructor: Map is not contiguous.");
 
  const int numProcs = comm->getSize();
 
  // Make room for the min global ID on each process, plus one
  // entry at the end for the "max cap."
  allMinGIDs_ = arcp<GO>(numProcs + 1);
  // Get my process' min global ID.
  GO minMyGID = map.getMinGlobalIndex();
  // Gather all of the min global IDs into the first numProcs
  // entries of allMinGIDs_.
 
  // FIXME (mfh 16 Apr 2013) GIANT HACK BELOW
  //
  // The purpose of this giant hack is that gatherAll appears to
  // interpret the "receive count" argument differently than
  // MPI_Allgather does.  Matt Bettencourt reports Valgrind issues
  // (memcpy with overlapping data) with MpiComm<int>::gatherAll,
  // which could relate either to this, or to OpenMPI.
#ifdef HAVE_TPETRACORE_MPI
  MPI_Datatype rawMpiType = MPI_INT;
  bool useRawMpi          = true;
  if (typeid(GO) == typeid(int)) {
    rawMpiType = MPI_INT;
  } else if (typeid(GO) == typeid(long)) {
    rawMpiType = MPI_LONG;
  } else {
    useRawMpi = false;
  }
  if (useRawMpi) {
    using Teuchos::MpiComm;
    using Teuchos::rcp_dynamic_cast;
    RCP<const MpiComm<int> > mpiComm =
        rcp_dynamic_cast<const MpiComm<int> >(comm);
    // It could be a SerialComm instead, even in an MPI build, so
    // be sure to check.
    if (!comm.is_null()) {
      MPI_Comm rawMpiComm = *(mpiComm->getRawMpiComm());
      const int err =
          MPI_Allgather(&minMyGID, 1, rawMpiType,
                        allMinGIDs_.getRawPtr(), 1, rawMpiType,
                        rawMpiComm);
      TEUCHOS_TEST_FOR_EXCEPTION(err != MPI_SUCCESS, std::runtime_error,
                                 "Tpetra::DistributedContiguousDirectory: MPI_Allgather failed");
    } else {
      gatherAll<int, GO>(*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr());
    }
  } else {
    gatherAll<int, GO>(*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr());
  }
#else   // NOT HAVE_TPETRACORE_MPI
  gatherAll<int, GO>(*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr());
#endif  // HAVE_TPETRACORE_MPI
  // FIXME (mfh 16 Apr 2013) GIANT HACK ABOVE
 
  // gatherAll<int, GO> (*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr ());
 
  // Put the max cap at the end.  Adding one lets us write loops
  // over the global IDs with the usual strict less-than bound.
  allMinGIDs_[numProcs] = map.getMaxAllGlobalIndex() + Teuchos::OrdinalTraits<GO>::one();
}
 
template <class LO, class GO, class NT>
std::string
DistributedContiguousDirectory<LO, GO, NT>::description() const {
  std::ostringstream os;
  os << "DistributedContiguousDirectory"
     << "<" << Teuchos::TypeNameTraits<LO>::name()
     << ", " << Teuchos::TypeNameTraits<GO>::name()
     << ", " << Teuchos::TypeNameTraits<NT>::name() << ">";
  return os.str();
}
 
template <class LO, class GO, class NT>
LookupStatus
ReplicatedDirectory<LO, GO, NT>::
    getEntriesImpl(const map_type& map,
                   const Teuchos::ArrayView<const GO>& globalIDs,
                   const Teuchos::ArrayView<int>& nodeIDs,
                   const Teuchos::ArrayView<LO>& localIDs,
                   const bool computeLIDs) const {
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::Comm;
  using Teuchos::RCP;
 
  LookupStatus res                    = AllIDsPresent;
  RCP<const Teuchos::Comm<int> > comm = map.getComm();
  const int myRank                    = comm->getRank();
 
  // Map is on one process or is locally replicated.
  typename ArrayView<int>::iterator procIter = nodeIDs.begin();
  typename ArrayView<LO>::iterator lidIter   = localIDs.begin();
  typename ArrayView<const GO>::iterator gidIter;
  for (gidIter = globalIDs.begin(); gidIter != globalIDs.end(); ++gidIter) {
    if (map.isNodeGlobalElement(*gidIter)) {
      *procIter++ = myRank;
      if (computeLIDs) {
        *lidIter++ = map.getLocalElement(*gidIter);
      }
    } else {
      // Advance the pointers, leaving these values set to invalid
      procIter++;
      if (computeLIDs) {
        lidIter++;
      }
      res = IDNotPresent;
    }
  }
  return res;
}
 
template <class LO, class GO, class NT>
LookupStatus
DistributedContiguousDirectory<LO, GO, NT>::
    getEntriesImpl(const map_type& map,
                   const Teuchos::ArrayView<const GO>& globalIDs,
                   const Teuchos::ArrayView<int>& nodeIDs,
                   const Teuchos::ArrayView<LO>& localIDs,
                   const bool computeLIDs) const {
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::Comm;
  using Teuchos::RCP;
 
  RCP<const Teuchos::Comm<int> > comm = map.getComm();
  const int numProcs                  = comm->getSize();
  const LO LINVALID                   = Teuchos::OrdinalTraits<LO>::invalid();
  const GO noGIDsOnProc               = std::numeric_limits<GO>::max();
  LookupStatus res                    = AllIDsPresent;
 
  // Find the first initialized GID for search below
  int firstProcWithGIDs;
  for (firstProcWithGIDs = 0; firstProcWithGIDs < numProcs;
       firstProcWithGIDs++) {
    if (allMinGIDs_[firstProcWithGIDs] != noGIDsOnProc) break;
  }
 
  // If Map is empty, return invalid values for all requested lookups
  // This case should not happen, as an empty Map is not considered
  // Distributed.
  if (firstProcWithGIDs == numProcs) {
    // No entries in Map
    res = (globalIDs.size() > 0) ? IDNotPresent : AllIDsPresent;
    std::fill(nodeIDs.begin(), nodeIDs.end(), -1);
    if (computeLIDs)
      std::fill(localIDs.begin(), localIDs.end(), LINVALID);
    return res;
  }
 
  const GO one    = as<GO>(1);
  const GO nOverP = as<GO>(map.getGlobalNumElements() / as<global_size_t>(numProcs - firstProcWithGIDs));
 
  // Map is distributed but contiguous.
  typename ArrayView<int>::iterator procIter = nodeIDs.begin();
  typename ArrayView<LO>::iterator lidIter   = localIDs.begin();
  typename ArrayView<const GO>::iterator gidIter;
  for (gidIter = globalIDs.begin(); gidIter != globalIDs.end(); ++gidIter) {
    LO LID    = LINVALID;  // Assume not found until proven otherwise
    int image = -1;
    GO GID    = *gidIter;
    // Guess uniform distribution (TODO: maybe replace by a binary search)
    // We go through all this trouble to avoid overflow and
    // signed / unsigned casting mistakes (that were made in
    // previous versions of this code).
    int curRank;
    const GO firstGuess = firstProcWithGIDs + GID / std::max(nOverP, one);
    curRank             = as<int>(std::min(firstGuess, as<GO>(numProcs - 1)));
 
    // This while loop will stop because
    // allMinGIDs_[np] == global num elements
    while (allMinGIDs_[curRank] == noGIDsOnProc) curRank++;
 
    bool badGID = false;
    while (curRank >= firstProcWithGIDs && GID < allMinGIDs_[curRank]) {
      curRank--;
      while (curRank >= firstProcWithGIDs &&
             allMinGIDs_[curRank] == noGIDsOnProc) curRank--;
    }
    if (curRank < firstProcWithGIDs) {
      // GID is lower than all GIDs in map
      badGID = true;
    } else if (curRank == numProcs) {
      // GID is higher than all GIDs in map
      badGID = true;
    } else {
      // we have the lower bound; now limit from above
      int above = curRank + 1;
      while (allMinGIDs_[above] == noGIDsOnProc) above++;
 
      while (GID >= allMinGIDs_[above]) {
        curRank = above;
        if (curRank == numProcs) {
          // GID is higher than all GIDs in map
          badGID = true;
          break;
        }
        above++;
        while (allMinGIDs_[above] == noGIDsOnProc) above++;
      }
    }
 
    if (!badGID) {
      image = curRank;
      LID   = as<LO>(GID - allMinGIDs_[image]);
    } else {
      res = IDNotPresent;
    }
    *procIter++ = image;
    if (computeLIDs) {
      *lidIter++ = LID;
    }
  }
  return res;
}
 
template <class LO, class GO, class NT>
DistributedNoncontiguousDirectory<LO, GO, NT>::
    DistributedNoncontiguousDirectory(const map_type& map)
  : oneToOneResult_(ONE_TO_ONE_NOT_CALLED_YET)
  ,  // to be revised below
  locallyOneToOne_(true)
  ,                      // to be revised below
  useHashTables_(false)  // to be revised below
{
  initialize(map, Teuchos::null);
}
 
template <class LO, class GO, class NT>
DistributedNoncontiguousDirectory<LO, GO, NT>::
    DistributedNoncontiguousDirectory(const map_type& map,
                                      const tie_break_type& tie_break)
  : oneToOneResult_(ONE_TO_ONE_NOT_CALLED_YET)
  ,  // to be revised below
  locallyOneToOne_(true)
  ,                      // to be revised below
  useHashTables_(false)  // to be revised below
{
  initialize(map, Teuchos::ptrFromRef(tie_break));
}
 
template <class LO, class GO, class NT>
void DistributedNoncontiguousDirectory<LO, GO, NT>::
    initialize(const map_type& map,
               Teuchos::Ptr<const tie_break_type> tie_break) {
  using std::cerr;
  using std::endl;
  using Teuchos::arcp;
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::typeName;
  using Teuchos::TypeNameTraits;
  typedef Array<int>::size_type size_type;
 
  // This class' implementation of getEntriesImpl() currently
  // encodes the following assumptions:
  //
  // 1. global_size_t >= GO
  // 2. global_size_t >= int
  // 3. global_size_t >= LO
  //
  // We check these assumptions here.
  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(global_size_t) < sizeof(GO),
                             std::logic_error, typeName(*this) << ": sizeof(Tpetra::"
                                                                  "global_size_t) = "
                                                               << sizeof(global_size_t) << " < sizeof(Global"
                                                                                           "Ordinal = "
                                                               << TypeNameTraits<LO>::name() << ") = " << sizeof(GO) << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(global_size_t) < sizeof(int),
                             std::logic_error, typeName(*this) << ": sizeof(Tpetra::"
                                                                  "global_size_t) = "
                                                               << sizeof(global_size_t) << " < sizeof(int) = " << sizeof(int) << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(sizeof(global_size_t) < sizeof(LO),
                             std::logic_error, typeName(*this) << ": sizeof(Tpetra::"
                                                                  "global_size_t) = "
                                                               << sizeof(global_size_t) << " < sizeof(Local"
                                                                                           "Ordinal = "
                                                               << TypeNameTraits<LO>::name() << ") = " << sizeof(LO) << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(!map.haveGlobalConstants(), std::logic_error, "Map needs to have global constants.");
 
  RCP<const Teuchos::Comm<int> > comm = map.getComm();
  const LO LINVALID                   = Teuchos::OrdinalTraits<LO>::invalid();
  const GO minAllGID                  = map.getMinAllGlobalIndex();
  const GO maxAllGID                  = map.getMaxAllGlobalIndex();
 
  // The "Directory Map" (see below) will have a range of elements
  // from the minimum to the maximum GID of the user Map, and a
  // minimum GID of minAllGID from the user Map.  It doesn't
  // actually have to store all those entries, though do beware of
  // calling getLocalElementList on it (see Bug 5822).
  const global_size_t numGlobalEntries = maxAllGID - minAllGID + 1;
 
  // We can't afford to replicate the whole directory on each
  // process, so create the "Directory Map", a uniform contiguous
  // Map that describes how we will distribute the directory over
  // processes.
  //
  // FIXME (mfh 08 May 2012) Here we're setting minAllGID to be
  // the index base.  The index base should be separate from the
  // minimum GID.
  directoryMap_ = rcp(new map_type(numGlobalEntries, minAllGID, comm,
                                   GloballyDistributed));
  // The number of Directory elements that my process owns.
  const size_t dir_numMyEntries = directoryMap_->getLocalNumElements();
 
  // Fix for Bug 5822: If the input Map is "sparse," that is if
  // the difference between the global min and global max GID is
  // much larger than the global number of elements in the input
  // Map, then it's possible that the Directory Map might have
  // many more entries than the input Map on this process.  This
  // can cause memory scalability issues.  In that case, we switch
  // from the array-based implementation of Directory storage to
  // the hash table - based implementation.  We don't use hash
  // tables all the time, because they are slower in the common
  // case of a nonsparse Map.
  //
  // NOTE: This is a per-process decision.  Some processes may use
  // array-based storage, whereas others may use hash table -
  // based storage.
 
  // A hash table takes a constant factor more space, more or
  // less, than an array.  Thus, it's not worthwhile, even in
  // terms of memory usage, always to use a hash table.
  // Furthermore, array lookups are faster than hash table
  // lookups, so it may be worthwhile to use an array even if it
  // takes more space.  The "sparsity threshold" governs when to
  // switch to a hash table - based implementation.
  const size_t inverseSparsityThreshold = 10;
  useHashTables_ =
      (dir_numMyEntries >= inverseSparsityThreshold * map.getLocalNumElements());
 
  // Get list of process IDs that own the directory entries for the
  // Map GIDs.  These will be the targets of the sends that the
  // Distributor will do.
  const int myRank          = comm->getRank();
  const size_t numMyEntries = map.getLocalNumElements();
  Array<int> sendImageIDs(numMyEntries);
  ArrayView<const GO> myGlobalEntries = map.getLocalElementList();
  // An ID not present in this lookup indicates that it lies outside
  // of the range [minAllGID,maxAllGID] (from map_).  this means
  // something is wrong with map_, our fault.
  const LookupStatus lookupStatus =
      directoryMap_->getRemoteIndexList(myGlobalEntries, sendImageIDs);
  TEUCHOS_TEST_FOR_EXCEPTION(
      lookupStatus == IDNotPresent, std::logic_error, Teuchos::typeName(*this) << " constructor: the Directory Map could not find out where one or "
                                                                                  "more of my Map's indices should go.  The input to getRemoteIndexList "
                                                                                  "is "
                                                                               << Teuchos::toString(myGlobalEntries) << ", and the output is " << Teuchos::toString(sendImageIDs()) << ".  The input Map itself has "
                                                                                                                                                                                       "the following entries on the calling process "
                                                                               << map.getComm()->getRank() << ": " << Teuchos::toString(map.getLocalElementList()) << ", and has " << map.getGlobalNumElements() << " total global indices in [" << map.getMinAllGlobalIndex() << "," << map.getMaxAllGlobalIndex() << "].  The Directory Map has " << directoryMap_->getGlobalNumElements() << " total global indices in "
                                                                                                                                                                                                                                                                                                                                                                                                "["
                                                                               << directoryMap_->getMinAllGlobalIndex() << "," << directoryMap_->getMaxAllGlobalIndex() << "], and the calling process "
                                                                                                                                                                           "has GIDs ["
                                                                               << directoryMap_->getMinGlobalIndex() << "," << directoryMap_->getMaxGlobalIndex() << "].  "
                                                                                                                                                                     "This probably means there is a bug in Map or Directory.  "
                                                                                                                                                                     "Please report this bug to the Tpetra developers.");
 
  // Initialize the distributor using the list of process IDs to
  // which to send.  We'll use the distributor to send out triples
  // of (GID, process ID, LID).  We're sending the entries to the
  // processes that the Directory Map says should own them, which is
  // why we called directoryMap_->getRemoteIndexList() above.
  Distributor distor(comm);
  const size_t numReceives = distor.createFromSends(sendImageIDs);
 
  // NOTE (mfh 21 Mar 2012) The following code assumes that
  // sizeof(GO) >= sizeof(int) and sizeof(GO) >= sizeof(LO).
  //
  // Create and fill buffer of (GID, PID, LID) triples to send
  // out.  We pack the (GID, PID, LID) triples into a single Array
  // of GO, casting the PID from int to GO and the LID from LO to
  // GO as we do so.
  //
  // FIXME (mfh 23 Mar 2014) This assumes that sizeof(LO) <=
  // sizeof(GO) and sizeof(int) <= sizeof(GO).  The former is
  // required, and the latter is generally the case, but we should
  // still check for this.
  const int packetSize = 3;  // We're sending triples, so packet size is 3.
  Kokkos::View<GO*, Kokkos::HostSpace> exportEntries("exportEntries", packetSize * numMyEntries);
  {
    size_t exportIndex = 0;
    for (size_t i = 0; i < numMyEntries; ++i) {
      exportEntries[exportIndex++] = myGlobalEntries[i];
      exportEntries[exportIndex++] = as<GO>(myRank);
      exportEntries[exportIndex++] = as<GO>(i);
    }
  }
  // Buffer of data to receive.  The Distributor figured out for
  // us how many packets we're receiving, when we called its
  // createFromSends() method to set up the distribution plan.
  Kokkos::View<GO*, Kokkos::HostSpace> importElements("importElements", packetSize * distor.getTotalReceiveLength());
 
  // Distribute the triples of (GID, process ID, LID).
  distor.doPostsAndWaits(exportEntries, packetSize, importElements);
 
  // Unpack the redistributed data.  Both implementations of
  // Directory storage map from an LID in the Directory Map (which
  // is the LID of the GID to store) to either a PID or an LID in
  // the input Map.  Each "packet" (contiguous chunk of
  // importElements) contains a triple: (GID, PID, LID).
  if (useHashTables_) {
    // Create the hash tables.  We know exactly how many elements
    // to expect in each hash table.  FixedHashTable's constructor
    // currently requires all the keys and values at once, so we
    // have to extract them in temporary arrays.  It may be
    // possible to rewrite FixedHashTable to use a "start fill" /
    // "end fill" approach that avoids the temporary arrays, but
    // we won't try that for now.
 
    // The constructors of Array and ArrayRCP that take a number
    // of elements all initialize the arrays.  Instead, allocate
    // raw arrays, then hand them off to ArrayRCP, to avoid the
    // initial unnecessary initialization without losing the
    // benefit of exception safety (and bounds checking, in a
    // debug build).
    LO* tableKeysRaw  = NULL;
    LO* tableLidsRaw  = NULL;
    int* tablePidsRaw = NULL;
    try {
      tableKeysRaw = new LO[numReceives];
      tableLidsRaw = new LO[numReceives];
      tablePidsRaw = new int[numReceives];
    } catch (...) {
      if (tableKeysRaw != NULL) {
        delete[] tableKeysRaw;
      }
      if (tableLidsRaw != NULL) {
        delete[] tableLidsRaw;
      }
      if (tablePidsRaw != NULL) {
        delete[] tablePidsRaw;
      }
      throw;
    }
    ArrayRCP<LO> tableKeys(tableKeysRaw, 0, numReceives, true);
    ArrayRCP<LO> tableLids(tableLidsRaw, 0, numReceives, true);
    ArrayRCP<int> tablePids(tablePidsRaw, 0, numReceives, true);
 
    if (tie_break.is_null()) {
      // Fill the temporary arrays of keys and values.
      size_type importIndex = 0;
      for (size_type i = 0; i < static_cast<size_type>(numReceives); ++i) {
        const GO curGID = importElements[importIndex++];
        const LO curLID = directoryMap_->getLocalElement(curGID);
        TEUCHOS_TEST_FOR_EXCEPTION(
            curLID == LINVALID, std::logic_error,
            Teuchos::typeName(*this) << " constructor: Incoming global index "
                                     << curGID << " does not have a corresponding local index in the "
                                                  "Directory Map.  Please report this bug to the Tpetra developers.");
        tableKeys[i] = curLID;
        tablePids[i] = importElements[importIndex++];
        tableLids[i] = importElements[importIndex++];
      }
      // Set up the hash tables.  The hash tables' constructor
      // detects whether there are duplicates, so that we can set
      // locallyOneToOne_.
      lidToPidTable_ =
          rcp(new lidToPidTable_type(tableKeys(), tablePids()));
      locallyOneToOne_ = !(lidToPidTable_->hasDuplicateKeys());
      lidToLidTable_ =
          rcp(new lidToLidTable_type(tableKeys(), tableLids()));
    } else {  // tie_break is NOT null
 
      // For each directory Map LID received, collect all the
      // corresponding (PID,LID) pairs.  If the input Map is not
      // one-to-one, corresponding directory Map LIDs will have
      // more than one pair.  In that case, we will use the
      // TieBreak object to pick exactly one pair.
      typedef std::map<LO, std::vector<std::pair<int, LO> > > pair_table_type;
      pair_table_type ownedPidLidPairs;
 
      // For each directory Map LID received, collect the zero or
      // more input Map (PID,LID) pairs into ownedPidLidPairs.
      size_type importIndex = 0;
      for (size_type i = 0; i < static_cast<size_type>(numReceives); ++i) {
        const GO curGID    = importElements[importIndex++];
        const LO dirMapLid = directoryMap_->getLocalElement(curGID);
        TEUCHOS_TEST_FOR_EXCEPTION(
            dirMapLid == LINVALID, std::logic_error,
            Teuchos::typeName(*this) << " constructor: Incoming global index "
                                     << curGID << " does not have a corresponding local index in the "
                                                  "Directory Map.  Please report this bug to the Tpetra developers.");
        tableKeys[i]  = dirMapLid;
        const int PID = importElements[importIndex++];
        const int LID = importElements[importIndex++];
 
        // These may change below.  We fill them in just to ensure
        // that they won't have invalid values.
        tablePids[i] = PID;
        tableLids[i] = LID;
 
        // For every directory Map LID, we have to remember all
        // (PID, LID) pairs.  The TieBreak object will arbitrate
        // between them in the loop below.
        ownedPidLidPairs[dirMapLid].push_back(std::make_pair(PID, LID));
      }
 
      // Use TieBreak to arbitrate between (PID,LID) pairs
      // corresponding to each directory Map LID.
      //
      // FIXME (mfh 23 Mar 2014) How do I know that i is the same
      // as the directory Map LID?
      // KDD 21 Mar 2018:  It isn't, especially if the user's IDs are not
      // contiguous, but the directory map is.  Need to set tablePids[i]
      // and tableLids[i], so need to loop over numReceives (as that is
      // how those arrays are allocated).  FIXED
 
      for (size_type i = 0; i < static_cast<size_type>(numReceives); ++i) {
        const LO dirMapLid = tableKeys[i];
        const std::vector<std::pair<int, LO> >& pidLidList =
            ownedPidLidPairs[dirMapLid];
        const size_t listLen = pidLidList.size();
        if (listLen == 0) continue;  // KDD This will never happen
        const GO dirMapGid = directoryMap_->getGlobalElement(dirMapLid);
        if (listLen > 1) {
          locallyOneToOne_ = false;
        }
        // If there is some (PID,LID) pair for the current input
        // Map LID, then it makes sense to invoke the TieBreak
        // object to arbitrate between the options.  Even if
        // there is only one (PID,LID) pair, we still want to
        // give the TieBreak object a chance to do whatever it
        // likes to do, in terms of side effects (e.g., track
        // (PID,LID) pairs).
        const size_type index =
            static_cast<size_type>(tie_break->selectedIndex(dirMapGid,
                                                            pidLidList));
        tablePids[i] = pidLidList[index].first;
        tableLids[i] = pidLidList[index].second;
      }
 
      // Set up the hash tables.
      lidToPidTable_ =
          rcp(new lidToPidTable_type(tableKeys(), tablePids()));
      lidToLidTable_ =
          rcp(new lidToLidTable_type(tableKeys(), tableLids()));
    }
  } else {
    if (tie_break.is_null()) {
      // Use array-based implementation of Directory storage.
      // Allocate these arrays and fill them with invalid values,
      // in case the input Map's GID list is sparse (i.e., does
      // not populate all GIDs from minAllGID to maxAllGID).
      PIDs_ = arcp<int>(dir_numMyEntries);
      std::fill(PIDs_.begin(), PIDs_.end(), -1);
      LIDs_ = arcp<LO>(dir_numMyEntries);
      std::fill(LIDs_.begin(), LIDs_.end(), LINVALID);
      // Fill in the arrays with PIDs resp. LIDs.
      size_type importIndex = 0;
      for (size_type i = 0; i < static_cast<size_type>(numReceives); ++i) {
        const GO curGID = importElements[importIndex++];
        const LO curLID = directoryMap_->getLocalElement(curGID);
        TEUCHOS_TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error,
                                   Teuchos::typeName(*this) << " constructor: Incoming global index "
                                                            << curGID << " does not have a corresponding local index in the "
                                                                         "Directory Map.  Please report this bug to the Tpetra developers.");
 
        // If PIDs_[curLID] is not -1, then curGID is a duplicate
        // on the calling process, so the Directory is not locally
        // one-to-one.
        if (PIDs_[curLID] != -1) {
          locallyOneToOne_ = false;
        }
        PIDs_[curLID] = importElements[importIndex++];
        LIDs_[curLID] = importElements[importIndex++];
      }
    } else {
      PIDs_ = arcp<int>(dir_numMyEntries);
      LIDs_ = arcp<LO>(dir_numMyEntries);
      std::fill(PIDs_.begin(), PIDs_.end(), -1);
 
      // All received (PID, LID) pairs go into ownedPidLidPairs.
      // This is a map from the directory Map's LID to the (PID,
      // LID) pair (where the latter LID comes from the input Map,
      // not the directory Map).  If the input Map is not
      // one-to-one, corresponding LIDs will have
      // ownedPidLidPairs[curLID].size() > 1.  In that case, we
      // will use the TieBreak object to pick exactly one pair.
      Array<std::vector<std::pair<int, LO> > > ownedPidLidPairs(dir_numMyEntries);
      size_t importIndex = 0;
      for (size_t i = 0; i < numReceives; ++i) {
        const GO GID  = importElements[importIndex++];
        const int PID = importElements[importIndex++];
        const LO LID  = importElements[importIndex++];
 
        const LO dirMapLid = directoryMap_->getLocalElement(GID);
        TEUCHOS_TEST_FOR_EXCEPTION(
            dirMapLid == LINVALID, std::logic_error,
            Teuchos::typeName(*this) << " constructor: Incoming global index "
                                     << GID << " does not have a corresponding local index in the "
                                               "Directory Map.  Please report this bug to the Tpetra developers.");
        ownedPidLidPairs[dirMapLid].push_back(std::make_pair(PID, LID));
      }
 
      // Use TieBreak to arbitrate between (PID,LID) pairs
      // corresponding to each directory Map LID.
      //
      // FIXME (mfh 23 Mar 2014) How do I know that i is the same
      // as the directory Map LID?
      // KDD 21 Mar 2018:  It isn't, especially if the user's IDs are not
      // contiguous.  Loop over all ownedPidLidPairs; skip those that have
      // empty lists.  FIXED
 
      for (size_t i = 0; i < dir_numMyEntries; ++i) {
        const std::vector<std::pair<int, LO> >& pidLidList =
            ownedPidLidPairs[i];
        const size_t listLen = pidLidList.size();
        if (listLen == 0) continue;  // KDD will happen for GIDs not in
                                     // KDD the user's source map
        const LO dirMapLid = static_cast<LO>(i);
        const GO dirMapGid = directoryMap_->getGlobalElement(dirMapLid);
        if (listLen > 1) {
          locallyOneToOne_ = false;
        }
        // If there is some (PID,LID) pair for the current input
        // Map LID, then it makes sense to invoke the TieBreak
        // object to arbitrate between the options.  Even if
        // there is only one (PID,LID) pair, we still want to
        // give the TieBreak object a chance to do whatever it
        // likes to do, in terms of side effects (e.g., track
        // (PID,LID) pairs).
        const size_type index =
            static_cast<size_type>(tie_break->selectedIndex(dirMapGid,
                                                            pidLidList));
        PIDs_[i] = pidLidList[index].first;
        LIDs_[i] = pidLidList[index].second;
      }
    }
  }
}
 
template <class LO, class GO, class NT>
std::string
DistributedNoncontiguousDirectory<LO, GO, NT>::description() const {
  std::ostringstream os;
  os << "DistributedNoncontiguousDirectory"
     << "<" << Teuchos::TypeNameTraits<LO>::name()
     << ", " << Teuchos::TypeNameTraits<GO>::name()
     << ", " << Teuchos::TypeNameTraits<NT>::name() << ">";
  return os.str();
}
 
template <class LO, class GO, class NT>
LookupStatus
DistributedNoncontiguousDirectory<LO, GO, NT>::
    getEntriesImpl(const map_type& map,
                   const Teuchos::ArrayView<const GO>& globalIDs,
                   const Teuchos::ArrayView<int>& nodeIDs,
                   const Teuchos::ArrayView<LO>& localIDs,
                   const bool computeLIDs) const {
  using Details::Behavior;
  using Details::verbosePrintArray;
  using std::cerr;
  using std::endl;
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::RCP;
  using Teuchos::toString;
  using size_type = typename Array<GO>::size_type;
  const char funcPrefix[] =
      "Tpetra::"
      "DistributedNoncontiguousDirectory::getEntriesImpl: ";
  const char errSuffix[] =
      "  Please report this bug to the Tpetra developers.";
 
  RCP<const Teuchos::Comm<int> > comm = map.getComm();
  const bool verbose                  = Behavior::verbose("Directory") ||
                       Behavior::verbose("Tpetra::Directory");
  const size_t maxNumToPrint = verbose ? Behavior::verbosePrintCountThreshold() : size_t(0);
 
  std::unique_ptr<std::string> procPrefix;
  if (verbose) {
    std::ostringstream os;
    os << "Proc ";
    if (map.getComm().is_null()) {
      os << "?";
    } else {
      os << map.getComm()->getRank();
    }
    os << ": ";
    procPrefix = std::unique_ptr<std::string>(
        new std::string(os.str()));
    os << funcPrefix << "{";
    verbosePrintArray(os, globalIDs, "GIDs", maxNumToPrint);
    os << ", ";
    verbosePrintArray(os, nodeIDs, "PIDs", maxNumToPrint);
    os << ", ";
    verbosePrintArray(os, localIDs, "LIDs", maxNumToPrint);
    os << ", computeLIDs: "
       << (computeLIDs ? "true" : "false") << "}" << endl;
    cerr << os.str();
  }
 
  const size_t numEntries = globalIDs.size();
  const LO LINVALID       = Teuchos::OrdinalTraits<LO>::invalid();
  LookupStatus res        = AllIDsPresent;
 
  //
  // Set up directory structure.
  //
 
  // If we're computing LIDs, we also have to include them in each
  // packet, along with the GID and process ID.
  const int packetSize = computeLIDs ? 3 : 2;
 
  // For data distribution, we use: Surprise!  A Distributor!
  Distributor distor(comm);
 
  // Get directory locations for the requested list of entries.
  Array<int> dirImages(numEntries);
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix << "Call directoryMap_->getRemoteIndexList"
       << endl;
    cerr << os.str();
  }
  res = directoryMap_->getRemoteIndexList(globalIDs, dirImages());
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix << "Director Map getRemoteIndexList out ";
    verbosePrintArray(os, dirImages, "PIDs", maxNumToPrint);
    os << endl;
    cerr << os.str();
  }
 
  // Check for unfound globalIDs and set corresponding nodeIDs to -1
  size_t numMissing = 0;
  if (res == IDNotPresent) {
    for (size_t i = 0; i < numEntries; ++i) {
      if (dirImages[i] == -1) {
        nodeIDs[i] = -1;
        if (computeLIDs) {
          localIDs[i] = LINVALID;
        }
        numMissing++;
      }
    }
  }
 
  Array<GO> sendGIDs;
  Array<int> sendImages;
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix << "Call Distributor::createFromRecvs"
       << endl;
    cerr << os.str();
  }
  distor.createFromRecvs(globalIDs, dirImages(), sendGIDs, sendImages);
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix << "Distributor::createFromRecvs result: "
       << "{";
    verbosePrintArray(os, sendGIDs, "sendGIDs", maxNumToPrint);
    os << ", ";
    verbosePrintArray(os, sendImages, "sendPIDs", maxNumToPrint);
    os << "}" << endl;
    cerr << os.str();
  }
  const size_type numSends = sendGIDs.size();
 
  //
  // mfh 13 Nov 2012:
  //
  // The code below temporarily stores LO, GO, and int values in
  // an array of global_size_t.  If one of the signed types (LO
  // and GO should both be signed) happened to be -1 (or some
  // negative number, but -1 is the one that came up today), then
  // conversion to global_size_t will result in a huge
  // global_size_t value, and thus conversion back may overflow.
  // (Teuchos::as doesn't know that we meant it to be an LO or GO
  // all along.)
  //
  // The overflow normally would not be a problem, since it would
  // just go back to -1 again.  However, Teuchos::as does range
  // checking on conversions in a debug build, so it throws an
  // exception (std::range_error) in this case.  Range checking is
  // generally useful in debug mode, so we don't want to disable
  // this behavior globally.
  //
  // We solve this problem by forgoing use of Teuchos::as for the
  // conversions below from LO, GO, or int to global_size_t, and
  // the later conversions back from global_size_t to LO, GO, or
  // int.
  //
  // I've recorded this discussion as Bug 5760.
  //
 
  //    global_size_t >= GO
  //    global_size_t >= size_t >= int
  //    global_size_t >= size_t >= LO
  // Therefore, we can safely store all of these in a global_size_t
  Kokkos::View<global_size_t*, Kokkos::HostSpace> exports("exports", packetSize * numSends);
  {
    // Packet format:
    // - If computing LIDs: (GID, PID, LID)
    // - Otherwise:         (GID, PID)
    //
    // "PID" means "process ID" (a.k.a. "node ID," a.k.a. "rank").
 
    // Current position to which to write in exports array.  If
    // sending pairs, we pack the (GID, PID) pair for gid =
    // sendGIDs[k] in exports[2*k], exports[2*k+1].  If sending
    // triples, we pack the (GID, PID, LID) pair for gid =
    // sendGIDs[k] in exports[3*k, 3*k+1, 3*k+2].
    size_t exportsIndex = 0;
 
    if (useHashTables_) {
      if (verbose) {
        std::ostringstream os;
        os << *procPrefix << "Pack exports (useHashTables_ true)"
           << endl;
        cerr << os.str();
      }
      for (size_type gidIndex = 0; gidIndex < numSends; ++gidIndex) {
        const GO curGID = sendGIDs[gidIndex];
        // Don't use as() here (see above note).
        exports[exportsIndex++] = static_cast<global_size_t>(curGID);
        const LO curLID         = directoryMap_->getLocalElement(curGID);
        TEUCHOS_TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error, funcPrefix << "Directory Map's global index " << curGID << " lacks "
                                                                                                                                    "a corresponding local index."
                                                                                    << errSuffix);
        // Don't use as() here (see above note).
        exports[exportsIndex++] =
            static_cast<global_size_t>(lidToPidTable_->get(curLID));
        if (computeLIDs) {
          // Don't use as() here (see above note).
          exports[exportsIndex++] =
              static_cast<global_size_t>(lidToLidTable_->get(curLID));
        }
      }
    } else {
      if (verbose) {
        std::ostringstream os;
        os << *procPrefix << "Pack exports (useHashTables_ false)"
           << endl;
        cerr << os.str();
      }
      for (size_type gidIndex = 0; gidIndex < numSends; ++gidIndex) {
        const GO curGID = sendGIDs[gidIndex];
        // Don't use as() here (see above note).
        exports[exportsIndex++] = static_cast<global_size_t>(curGID);
        const LO curLID         = directoryMap_->getLocalElement(curGID);
        TEUCHOS_TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error, funcPrefix << "Directory Map's global index " << curGID << " lacks "
                                                                                                                                    "a corresponding local index."
                                                                                    << errSuffix);
        // Don't use as() here (see above note).
        exports[exportsIndex++] =
            static_cast<global_size_t>(PIDs_[curLID]);
        if (computeLIDs) {
          // Don't use as() here (see above note).
          exports[exportsIndex++] =
              static_cast<global_size_t>(LIDs_[curLID]);
        }
      }
    }
 
    TEUCHOS_TEST_FOR_EXCEPTION(exportsIndex > exports.size(), std::logic_error,
                               funcPrefix << "On Process " << comm->getRank() << ", "
                                                                                 "exportsIndex = "
                                          << exportsIndex << " > exports.size() = "
                                          << exports.size() << "." << errSuffix);
  }
 
  TEUCHOS_TEST_FOR_EXCEPTION(numEntries < numMissing, std::logic_error, funcPrefix << "On Process " << comm->getRank() << ", numEntries = " << numEntries << " < numMissing = " << numMissing << "." << errSuffix);
 
  //
  // mfh 13 Nov 2012: See note above on conversions between
  // global_size_t and LO, GO, or int.
  //
  const size_t numRecv = numEntries - numMissing;
 
  {
    const size_t importLen         = packetSize * distor.getTotalReceiveLength();
    const size_t requiredImportLen = numRecv * packetSize;
    const int myRank               = comm->getRank();
    TEUCHOS_TEST_FOR_EXCEPTION(
        importLen < requiredImportLen, std::logic_error,
        "Tpetra::Details::DistributedNoncontiguousDirectory::getEntriesImpl: "
        "On Process "
            << myRank << ": The 'imports' array must have length "
                         "at least "
            << requiredImportLen << ", but its actual length is " << importLen << ".  numRecv: " << numRecv << ", packetSize: " << packetSize << ", numEntries (# GIDs): " << numEntries << ", numMissing: " << numMissing << ": distor.getTotalReceiveLength(): "
            << distor.getTotalReceiveLength() << ".  " << std::endl
            << "Distributor description: " << distor.description() << ".  "
            << std::endl
            << "Please report this bug to the Tpetra developers.");
  }
 
  Kokkos::View<global_size_t*, Kokkos::HostSpace> imports("imports", packetSize * distor.getTotalReceiveLength());
  // FIXME (mfh 20 Mar 2014) One could overlap the sort2() below
  // with communication, by splitting this call into doPosts and
  // doWaits.  The code is still correct in this form, however.
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix << "Call doPostsAndWaits: {"
       << "packetSize: " << packetSize << ", ";
    verbosePrintArray(os, exports, "exports", maxNumToPrint);
    os << "}" << endl;
    cerr << os.str();
  }
  distor.doPostsAndWaits(exports, packetSize, imports);
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix << "doPostsAndWaits result: ";
    verbosePrintArray(os, imports, "imports", maxNumToPrint);
    os << endl;
    cerr << os.str();
  }
 
  Array<GO> sortedIDs(globalIDs);  // deep copy (for later sorting)
  Array<GO> offset(numEntries);    // permutation array (sort2 output)
  for (GO ii = 0; ii < static_cast<GO>(numEntries); ++ii) {
    offset[ii] = ii;
  }
  sort2(sortedIDs.begin(), sortedIDs.begin() + numEntries, offset.begin());
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix;
    verbosePrintArray(os, sortedIDs, "sortedIDs", maxNumToPrint);
    os << ", ";
    verbosePrintArray(os, offset, "offset", maxNumToPrint);
    os << endl;
    cerr << os.str();
  }
 
  size_t importsIndex = 0;
 
  // we know these conversions are in range, because we loaded this data
  //
  // Don't use as() for conversions here; we know they are in range.
  for (size_t i = 0; i < numRecv; ++i) {
    const GO curGID = static_cast<GO>(imports[importsIndex++]);
    auto p1         = std::equal_range(sortedIDs.begin(),
                                       sortedIDs.end(), curGID);
    if (p1.first != p1.second) {
      const size_t j = p1.first - sortedIDs.begin();
      nodeIDs[offset[j]] =
          static_cast<int>(imports[importsIndex++]);
      if (computeLIDs) {
        localIDs[offset[j]] =
            static_cast<LO>(imports[importsIndex++]);
      }
      if (nodeIDs[offset[j]] == -1) {
        res = IDNotPresent;
      }
    }
  }
 
  TEUCHOS_TEST_FOR_EXCEPTION(size_t(importsIndex) > size_t(imports.size()),
                             std::logic_error, funcPrefix << "On Process " << comm->getRank() << ": importsIndex = " << importsIndex << " > imports.size() = " << imports.size() << ".  "
                                                                                                                                                                                    "numRecv: "
                                                          << numRecv << ", packetSize: " << packetSize << ", "
                                                                                                          "numEntries (# GIDs): "
                                                          << numEntries << ", numMissing: " << numMissing << ": distor.getTotalReceiveLength(): " << distor.getTotalReceiveLength() << "." << errSuffix);
  if (verbose) {
    std::ostringstream os;
    os << *procPrefix << funcPrefix << "Done!" << endl;
    cerr << os.str();
  }
  return res;
}
 
template <class LO, class GO, class NT>
bool DistributedNoncontiguousDirectory<LO, GO, NT>::
    isOneToOne(const Teuchos::Comm<int>& comm) const {
  if (oneToOneResult_ == ONE_TO_ONE_NOT_CALLED_YET) {
    const int lcl121 = isLocallyOneToOne() ? 1 : 0;
    int gbl121       = 0;
    Teuchos::reduceAll<int, int>(comm, Teuchos::REDUCE_MIN, lcl121,
                                 Teuchos::outArg(gbl121));
    oneToOneResult_ = (gbl121 == 1) ? ONE_TO_ONE_TRUE : ONE_TO_ONE_FALSE;
  }
  return (oneToOneResult_ == ONE_TO_ONE_TRUE);
}
}  // namespace Details
}  // namespace Tpetra
 
//
// Explicit instantiation macro
//
// Must be expanded from within the Tpetra namespace!
//
#define TPETRA_DIRECTORYIMPL_INSTANT(LO, GO, NODE)                \
  namespace Details {                                             \
  template class Directory<LO, GO, NODE>;                         \
  template class ReplicatedDirectory<LO, GO, NODE>;               \
  template class ContiguousUniformDirectory<LO, GO, NODE>;        \
  template class DistributedContiguousDirectory<LO, GO, NODE>;    \
  template class DistributedNoncontiguousDirectory<LO, GO, NODE>; \
  }
 
#endif  // TPETRA_DIRECTORYIMPL_DEF_HPP