Tpetra_Import_Util2.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
// @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 <Teuchos_Array.hpp>
#include <utility>
#include <set>
 
#include "Tpetra_CrsMatrix_decl.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 <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void
lowCommunicationMakeColMapAndReindex (const Teuchos::ArrayView<const size_t> &rowPointers,
                                      const Teuchos::ArrayView<LocalOrdinal> &columnIndices_LID,
                                      const Teuchos::ArrayView<GlobalOrdinal> &columnIndices_GID,
                                      const Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > & domainMap,
                                      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(CRS_rowptr, CRS_colind, 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)
{
  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: ";
 
  // 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];
      }
    }
  }
}
 
 
 
 
// 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