Zoltan2_CommGraphModel.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//   Zoltan2: A package of combinatorial algorithms for scientific computing
//                  Copyright 2012 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 Karen Devine      (kddevin@sandia.gov)
//                    Erik Boman        (egboman@sandia.gov)
//                    Siva Rajamanickam (srajama@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
 
#ifndef _ZOLTAN2_COMMGRAPHMODEL_HPP_
#define _ZOLTAN2_COMMGRAPHMODEL_HPP_
 
#include <Zoltan2_Model.hpp>
#include <Zoltan2_ModelHelpers.hpp>
#include <Zoltan2_InputTraits.hpp>
#include <Zoltan2_MatrixAdapter.hpp>
#include <Zoltan2_GraphAdapter.hpp>
#include <Zoltan2_XpetraCrsGraphAdapter.hpp>
#include <Zoltan2_IdentifierAdapter.hpp>
#include <Zoltan2_VectorAdapter.hpp>
#include <Zoltan2_MeshAdapter.hpp>
#include <Zoltan2_StridedData.hpp>
#include <unordered_map>
 
namespace Zoltan2 {
 
 
template <typename Adapter>
class CommGraphModel : public Model<Adapter>
{
public:
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
  typedef typename Adapter::scalar_t    scalar_t;
  typedef typename Adapter::gno_t       gno_t;
  typedef typename Adapter::lno_t       lno_t;
  typedef typename Adapter::node_t      node_t;
  typedef typename Adapter::user_t      user_t;
  typedef typename Adapter::userCoord_t userCoord_t;
  typedef StridedData<lno_t, scalar_t>  input_t;
  typedef typename Adapter::offset_t    offset_t;
#endif
 
  ~CommGraphModel() { }
 
  CommGraphModel(const RCP<const MatrixAdapter<user_t,userCoord_t> > &ia,
     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm)
  {
    throw std::runtime_error("CommGraphModel is not implemented for MatrixAdapter yet.");
  }
 
  CommGraphModel(const RCP<const GraphAdapter<user_t,userCoord_t> > &ia,
     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm);
 
  CommGraphModel(const RCP<const MeshAdapter<user_t> > &ia,
     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm)
  {
    throw std::runtime_error("CommGraphModel is not implemented for MeshAdapter yet.");
  }
 
  CommGraphModel(const RCP<const VectorAdapter<userCoord_t> > &/* ia */,
     const RCP<const Environment> &/* env */, const RCP<const Comm<int> > &/* comm */)
  {
    throw std::runtime_error("cannot build CommGraphModel from VectorAdapter");
  }
 
  CommGraphModel(const RCP<const IdentifierAdapter<user_t> > &ia,
     const RCP<const Environment> &env, const RCP<const Comm<int> > &comm)
  {
    throw std::runtime_error("cannot build GraphModel from IdentifierAdapter");
  }
 
  const RCP<const Comm<int> > getComm() { return comm_; }
 
  size_t getLocalNumVertices() const { return nLocalVertices_; }
 
  size_t getGlobalNumVertices() const { return nGlobalVertices_; }
 
  size_t getLocalNumEdges() const { return nLocalEdges_; }
 
  size_t getGlobalNumEdges() const { return nGlobalEdges_; }
 
  int getNumWeightsPerVertex() const { return nWeightsPerVertex_; }
 
  int getNumWeightsPerEdge() const { return nWeightsPerEdge_; }
 
 
  size_t getVertexList(
    ArrayView<const gno_t> &Ids,
    ArrayView<input_t> &wgts) const
  {
    Ids = vGids_.view(0, nLocalVertices_);
    wgts = vWeights_.view(0, nWeightsPerVertex_);
    return nLocalVertices_;
  }
 
  // Implied Vertex LNOs from getVertexList are used as indices to offsets
  // array.
  // Vertex GNOs are returned as neighbors in edgeIds.
 
  size_t getEdgeList( ArrayView<const gno_t> &edgeIds,
    ArrayView<const offset_t> &offsets,
    ArrayView<input_t> &wgts) const
  {
    edgeIds = eGids_.view(0, nLocalEdges_);
    offsets = eOffsets_.view(0, nLocalVertices_+1);
    wgts = eWeights_.view(0, nWeightsPerEdge_);
    return nLocalEdges_;
  }
 
  inline void getVertexDist(ArrayView<size_t> &vtxdist) const
  {
    vtxdist = vtxDist_();
    if (vtxDist_.size() == 0) {
      throw std::runtime_error("getVertexDist is available only "
                               "when consecutiveIdsRequired");
    }
  }
 
  // The Model interface.
 
  size_t getLocalNumObjects() const { return nLocalVertices_; }
 
  size_t getGlobalNumObjects() const { return nGlobalVertices_; }
 
  // Migration-related functions.
 
  int getNumActiveRanks() const { return nActiveRanks_; }
 
  int getDestinationRank() const { return destRank_; }
 
  int getStartRank() const { return startRank_; }
 
  int getEndRank() const { return endRank_; }
 
private:
 
  void print(); // For debugging
  void migrateGraph(); // For debugging
 
  const RCP<const Environment > env_;
  const RCP<const Comm<int> > comm_;
 
  int threshold_;        // threshold on #vertices each rank stores post-migration 
  int nActiveRanks_ ;    // # ranks for the small graph to be partitioned on
  int destRank_, startRank_, endRank_;
                       
 
  size_t nLocalVertices_;                // # local vertices in built graph
  size_t nGlobalVertices_;               // # global vertices in built graph
  ArrayRCP<gno_t> vGids_;                  // vertices of graph built in model;
                                           // may be same as adapter's input
                                           // or may be renumbered 0 to (N-1).
 
  int nWeightsPerVertex_;
  ArrayRCP<input_t> vWeights_;
 
  // Note: in some cases, size of these arrays
  // may be larger than nLocalEdges_.  So do not use .size().
  // Use nLocalEdges_, nGlobalEdges_
 
  size_t nLocalEdges_;                  // # local edges in built graph
  size_t nGlobalEdges_;                 // # global edges in built graph
  ArrayRCP<gno_t> eGids_;                 // edges of graph built in model
  ArrayRCP<offset_t> eOffsets_;           // edge offsets build in model
                                          // May be same as adapter's input 
                                          // or may differ
                                          // due to renumbering, self-edge
                                          // removal, or local graph.
 
  int nWeightsPerEdge_;
  ArrayRCP<input_t> eWeights_;            // edge weights in built graph
                                          // May be same as adapter's input
                                          // or may differ due to self-edge
                                          // removal, or local graph.
 
  ArrayRCP<size_t> vtxDist_;              // If consecutiveIdsRequired,
                                          // vtxDist (as needed by ParMETIS
                                          // and Scotch) is also created.
                                          // Otherwise, it is Teuchos::null.
};
 
 
// GraphModel from GraphAdapter
template <typename Adapter>
CommGraphModel<Adapter>::CommGraphModel(
  const RCP<const GraphAdapter<user_t,userCoord_t> > &bia,
  const RCP<const Environment> &env,
  const RCP<const Comm<int> > &comm):
        env_(env),
        comm_(comm),
        nLocalVertices_(0),
        nGlobalVertices_(0),
        vGids_(),
        nWeightsPerVertex_(0),
        vWeights_(),
        nLocalEdges_(0),
        nGlobalEdges_(0),
        eGids_(),
        eOffsets_(),
        nWeightsPerEdge_(0),
        eWeights_(),
        vtxDist_()
{
  int commSize = comm_->getSize();
 
  // Get XpetraCrsGraphAdapter from GraphAdapter 
  RCP<XpetraCrsGraphAdapter<user_t, userCoord_t>> ia; 
  try{
    RCP<GraphAdapter<user_t, userCoord_t>> tmp = 
       rcp_const_cast<GraphAdapter<user_t, userCoord_t>>(bia);
    ia = rcp_dynamic_cast<XpetraCrsGraphAdapter<user_t, userCoord_t>>(tmp);
  }
  Z2_FORWARD_EXCEPTIONS;
 
  // Get the graph from the input adapter
  auto inGraph = ia->getXpetraGraph();
 
  // Get the importer of the graph
  auto imp = inGraph->getImporter();
 
  // Identify nbor PIDs and number of entries sent per PID
  std::map<int,double> exportpidmap;
  auto exportpids = imp->getExportPIDs();
  size_t nexportpids = imp->getNumExportIDs();
  for (size_t i = 0; i < nexportpids; i++) {
    int k = exportpids[i];
    if (exportpidmap.find(k) != exportpidmap.end())
      exportpidmap[k] = exportpidmap[k] + 1.;
    else
      exportpidmap[k] = 1.;
  }
 
  // Set sizes
  // There is only one vertex in each rank 
  nLocalVertices_ = 1;
  nLocalEdges_ = exportpidmap.size();
 
  // Allocate space
  vGids_ = arcp(new gno_t[nLocalVertices_],
    0, nLocalVertices_, true);
  eGids_ = arcp(new gno_t[nLocalEdges_],
    0, nLocalEdges_, true);  
  eOffsets_ = arcp(new offset_t[nLocalVertices_+1],
       0, nLocalVertices_+1, true);
  scalar_t *wgts2 = new scalar_t [nLocalEdges_];
 
  // Form the vertices
  vGids_[0] = comm->getRank();
  
  // Form the edges
  size_t ptr = 0;
  eOffsets_[0] = ptr;
  for (std::map<int,double>::iterator it = exportpidmap.begin();
       it != exportpidmap.end(); it++) {
    eGids_[ptr] = it->first;
    wgts2[ptr++] = it->second;
  }
  eOffsets_[nLocalVertices_] = ptr;
    
  // Edge weights
  nWeightsPerEdge_ = 1;
  input_t *wgts = new input_t [nWeightsPerEdge_];
  eWeights_ = arcp(wgts, 0, nWeightsPerEdge_, true);
 
  for (int w=0; w < nWeightsPerEdge_; w++){
    int stride=0;
    ArrayRCP<const scalar_t> wgtArray = arcp(wgts2, 0, nLocalEdges_, true);
    eWeights_[w] = input_t(wgtArray, stride);
  }
 
  // Vertex weights
  nWeightsPerVertex_ = 1;
  input_t *weightInfo = new input_t [nWeightsPerVertex_];
 
  for (int idx=0; idx < nWeightsPerVertex_; idx++){
    scalar_t *wgt = new scalar_t [nLocalVertices_];
    wgt[0] = inGraph->getNodeNumEntries();
    ArrayRCP<const scalar_t> wgtArray = arcp(wgt, 0, nLocalVertices_, true);
    weightInfo[idx] = input_t(wgtArray, 1);
  }
 
  vWeights_ = arcp<input_t>(weightInfo, 0, nWeightsPerVertex_, true);
 
  reduceAll<int, size_t>(*comm_, Teuchos::REDUCE_SUM, 1,
       &nLocalVertices_, &nGlobalVertices_);
  reduceAll<int, size_t>(*comm_, Teuchos::REDUCE_SUM, 1,
       &nLocalEdges_, &nGlobalEdges_);
 
  // Build vtxDist_ array starting with vGid on each rank
  vtxDist_ = arcp(new size_t[commSize+1], 0, commSize+1, true);
  vtxDist_[0] = 0;
  Teuchos::gatherAll(*comm_, 1, &nLocalVertices_, commSize, &vtxDist_[1]);
  for (int i = 0; i < commSize; i++)
    vtxDist_[i+1] += vtxDist_[i];
 
  // Migrate the quotient graph into smaller number of MPI ranks (active ranks)
  migrateGraph();
}
 
template <typename Adapter>
void CommGraphModel<Adapter>::migrateGraph()
{
 
  // Set default threshold for migration
  threshold_ = 1024;  
 
  // Check if the user set the threshold value 
  const ParameterList &pl = env_->getParameters();
  const Teuchos::ParameterEntry *pe = pl.getEntryPtr("quotient_threshold");
  if (pe)
    threshold_ = pe->getValue<int>(&threshold_);
 
  // Compute the sizes of/in the new distribution
  nActiveRanks_ = std::ceil((double) nGlobalVertices_ / threshold_);
  size_t avgVertexShare = nGlobalVertices_ / nActiveRanks_;
  size_t myVertexShare = 0;
 
  int me = comm_->getRank();
  int commSize = comm_->getSize();
 
  // Save the original pointers 
  ArrayRCP<offset_t> old_eOffsets_ = eOffsets_;
  ArrayRCP<gno_t> old_eGids_ = eGids_;
  size_t old_nLocalEdges_ = nLocalEdges_;
  ArrayRCP<input_t> old_vWeights_ = vWeights_;
  ArrayRCP<input_t> old_eWeights_ = eWeights_;
 
  // Compute whom to send to
  destRank_ = me / (int) avgVertexShare;
  if(destRank_ >= nActiveRanks_)
    destRank_ = nActiveRanks_ - 1;
  
  // Start with sending the size of the edge list 
  RCP<CommRequest<int>> *requests;
  if(me < nActiveRanks_) {
 
    // Determine the range of ranks to receive edges from
    // Needs to be updated when chunks are introduced 
    startRank_ = me * static_cast<int>(avgVertexShare);
    endRank_ = (me+1) * static_cast<int>(avgVertexShare);
    if(me == nActiveRanks_ - 1 ) // Last rank gets the surplus 
      endRank_ = static_cast<int>(nGlobalVertices_);
    myVertexShare = endRank_ - startRank_;
 
    eOffsets_ = arcp(new offset_t[myVertexShare+1], 0, myVertexShare+1, true);
    eOffsets_[0] = 0;
 
    // Receive the sizes of their edge list
    requests = new RCP<CommRequest<int>>[myVertexShare];
    for(int i = startRank_; i < endRank_; i++) {
      requests[i-startRank_] = Teuchos::ireceive<int, offset_t>(*comm_, 
                 arcp(&eOffsets_[i-startRank_+1], 0, 1, false), 
                 i);
    } 
 
    // Send adjacency size  even though this rank will remain active
    Teuchos::send<int, offset_t>(*comm_, 1, &old_eOffsets_[nLocalVertices_], destRank_);
     
    // Wait 
    Teuchos::waitAll<int>(*comm_, Teuchos::arrayView(requests, myVertexShare)); 
    
    // Prefix sum over the offsets
    for(size_t i = 1; i <= myVertexShare; i++)
      eOffsets_[i] += eOffsets_[i-1];
 
    // Recompute the number of local edges
    nLocalEdges_ = eOffsets_[myVertexShare];
 
    // Reallocate the adjacency array
    eGids_ = arcp(new gno_t[nLocalEdges_], 0, nLocalEdges_, true);
 
    
    // Receive the adjacency lists
    for(int i = startRank_; i < endRank_; i++) {
      offset_t adjStartRank_ = eOffsets_[i-startRank_];
      offset_t adjSize = eOffsets_[i-startRank_+1] - adjStartRank_;
      requests[i-startRank_] = Teuchos::ireceive<int, gno_t>(*comm_, 
              arcp(&eGids_[adjStartRank_], 0, adjSize, false), 
              i);
    } 
   
    // Send adjacency even though this rank will remain active
    Teuchos::send<int, gno_t>(*comm_, old_nLocalEdges_, &old_eGids_[0], destRank_);
    Teuchos::waitAll<int>(*comm_, Teuchos::arrayView(requests, myVertexShare)); 
 
    
    // Migrate vertex weights arrays
    scalar_t *wgts = new scalar_t [myVertexShare];
    for(int i = startRank_; i < endRank_; i++) {
      requests[i-startRank_] = Teuchos::ireceive<int, scalar_t>(*comm_, 
                 arcp(&wgts[i-startRank_], 0, 1, false), // assumes one vertex per rank  
                 i);
    } 
 
    const scalar_t *wPtr;
    size_t wLen = 0;
    int stride = 0;
    old_vWeights_[0].getStridedList(wLen, wPtr, stride);
    Teuchos::send<int, scalar_t>(*comm_, nLocalVertices_, wPtr, destRank_);
    
    Teuchos::waitAll<int>(*comm_, Teuchos::arrayView(requests, myVertexShare));     
 
    input_t *weightInfo = new input_t [nWeightsPerVertex_];
    for (int idx=0; idx < nWeightsPerVertex_; idx++){
      ArrayRCP<const scalar_t> wgtArray = arcp(wgts, 0, myVertexShare, true);
      weightInfo[idx] = input_t(wgtArray, 1);
    }
    vWeights_ = arcp<input_t>(weightInfo, 0, nWeightsPerVertex_, true);
 
    // Migrate edge weights arrays
    scalar_t *ewgts = new scalar_t [nLocalEdges_];
    for(int i = startRank_; i < endRank_; i++) {
      offset_t adjStartRank_ = eOffsets_[i-startRank_];
      offset_t adjSize = eOffsets_[i-startRank_+1] - adjStartRank_;
      requests[i-startRank_] = Teuchos::ireceive<int, scalar_t>(*comm_, 
                 arcp(&ewgts[adjStartRank_], 0, adjSize, false), // assumes one vertex per rank  
                 i);
    } 
 
    old_eWeights_[0].getStridedList(wLen, wPtr, stride);
    Teuchos::send<int, scalar_t>(*comm_, old_nLocalEdges_, wPtr, destRank_);
    
    Teuchos::waitAll<int>(*comm_, Teuchos::arrayView(requests, myVertexShare));     
 
    input_t *eweightInfo = new input_t [nWeightsPerEdge_];
    for (int idx=0; idx < nWeightsPerEdge_; idx++){
      ArrayRCP<const scalar_t> ewgtArray = arcp(ewgts, 0, nLocalEdges_, true);
      eweightInfo[idx] = input_t(ewgtArray, 1);
    }
    eWeights_ = arcp<input_t>(eweightInfo, 0, nWeightsPerEdge_, true);
 
    
    // Finalize the migration
    vGids_ = arcp(new gno_t[myVertexShare], 0, myVertexShare, true);
    for(int i = startRank_; i < endRank_; i++)
      vGids_[i-startRank_] = i;
    
    nLocalVertices_ = myVertexShare;
    
 
  }
  else {
 
    // Send adjacency size 
    Teuchos::send<int, offset_t>(*comm_, 1, &eOffsets_[nLocalVertices_], destRank_);
    
    // Send adjacency list
    Teuchos::send<int, gno_t>(*comm_, nLocalEdges_, &eGids_[0], destRank_);
    
    // Send vertex weights list
    const scalar_t *wPtr;
    size_t wLen = 0;
    int stride = 0;
    vWeights_[0].getStridedList(wLen, wPtr, stride);
    Teuchos::send<int, scalar_t>(*comm_, nLocalVertices_, wPtr, destRank_);
 
    // Send edge weights list
    eWeights_[0].getStridedList(wLen, wPtr, stride);
    Teuchos::send<int, scalar_t>(*comm_, nLocalEdges_, wPtr, destRank_);
    
    nLocalVertices_ = 0;
  }
 
  for (int i = 0; i <= commSize; i++)
    vtxDist_[i] = 0;
 
  Teuchos::gatherAll(*comm_, 1, &nLocalVertices_, commSize, &vtxDist_[1]);
  for (int i = 0; i < commSize; i++)
    vtxDist_[i+1] += vtxDist_[i];
  
}
 
template <typename Adapter>
void CommGraphModel<Adapter>::print()
{
  std::ostream *os = env_->getDebugOStream();
  
  int me = comm_->getRank();
  
  *os << me
      << " Nvtx  " << nLocalVertices_
      << " Nedge " << nLocalEdges_
      << " NVWgt " << nWeightsPerVertex_
      << " NEWgt " << nWeightsPerEdge_
      << std::endl;
 
  for (size_t i = 0; i < nLocalVertices_; i++) {
    *os << me << " " << i << " GID " << vGids_[i] << ": ";
    for (offset_t j = eOffsets_[i]; j < eOffsets_[i+1]; j++)
      *os << eGids_[j] << " " ;
    *os << std::endl;
  }
  
  if (nWeightsPerVertex_) {
    for (size_t i = 0; i < nLocalVertices_; i++) {
      *os << me << " " << i << " VWGTS " << vGids_[i] << ": ";
      for (int j = 0; j < nWeightsPerVertex_; j++)
        *os << vWeights_[j][i] << " ";
      *os << std::endl;
    }
  }
 
  if (nWeightsPerEdge_) {
    for (size_t i = 0; i < nLocalVertices_; i++) {
      *os << me << " " << i << " EWGTS " << vGids_[i] << ": ";
      for (offset_t j = eOffsets_[i]; j < eOffsets_[i+1]; j++) {
        *os << eGids_[j] << " (";
        for (int w = 0; w < nWeightsPerEdge_; w++)
          *os << eWeights_[w][j] << " ";
        *os << ") ";
      }
      *os << std::endl;
    }
  }
 
}
 
}   // namespace Zoltan2
 
 
#endif