contig.c

Go to the documentation of this file.

#include "metislib.h"


/*************************************************************************/
/*************************************************************************/
idx_t FindPartitionInducedComponents(graph_t *graph, idx_t *where, 
          idx_t *cptr, idx_t *cind)
{
  idx_t i, ii, j, jj, k, me=0, nvtxs, first, last, nleft, ncmps;
  idx_t *xadj, *adjncy;
  idx_t *touched, *perm, *todo;
  idx_t mustfree_ccsr=0, mustfree_where=0;

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  adjncy = graph->adjncy;

  /* Deal with NULL supplied cptr/cind vectors */
  if (cptr == NULL) {
    cptr = imalloc(nvtxs+1, "FindPartitionInducedComponents: cptr");
    cind = imalloc(nvtxs, "FindPartitionInducedComponents: cind");
    mustfree_ccsr = 1;
  }

  /* Deal with NULL supplied where vector */
  if (where == NULL) {
    where = ismalloc(nvtxs, 0, "FindPartitionInducedComponents: where");
    mustfree_where = 1;
  }

  /* Allocate memory required for the BFS traversal */
  perm    = iincset(nvtxs, 0, imalloc(nvtxs, "FindPartitionInducedComponents: perm"));
  todo    = iincset(nvtxs, 0, imalloc(nvtxs, "FindPartitionInducedComponents: todo"));
  touched = ismalloc(nvtxs, 0, "FindPartitionInducedComponents: touched");


  /* Find the connected componends induced by the partition */
  ncmps = -1;
  first = last = 0;
  nleft = nvtxs;
  while (nleft > 0) {
    if (first == last) { /* Find another starting vertex */
      cptr[++ncmps] = first;
      ASSERT(touched[todo[0]] == 0);
      i = todo[0];
      cind[last++] = i;
      touched[i] = 1;
      me = where[i];
    }

    i = cind[first++];
    k = perm[i];
    j = todo[k] = todo[--nleft];
    perm[j] = k;

    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      if (where[k] == me && !touched[k]) {
        cind[last++] = k;
        touched[k] = 1;
      }
    }
  }
  cptr[++ncmps] = first;

  if (mustfree_ccsr)
    gk_free((void **)&cptr, &cind, LTERM);
  if (mustfree_where)
    gk_free((void **)&where, LTERM);

  gk_free((void **)&perm, &todo, &touched, LTERM);

  return ncmps;
}


/*************************************************************************/
/*************************************************************************/
void ComputeBFSOrdering(ctrl_t *ctrl, graph_t *graph, idx_t *bfsperm)
{
  idx_t i, j, k, nvtxs, first, last;
  idx_t *xadj, *adjncy, *perm;

  WCOREPUSH;

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  adjncy = graph->adjncy;

  /* Allocate memory required for the BFS traversal */
  perm = iincset(nvtxs, 0, iwspacemalloc(ctrl, nvtxs));

  iincset(nvtxs, 0, bfsperm);  /* this array will also store the vertices
                                  still to be processed */

  /* Find the connected componends induced by the partition */
  first = last = 0;
  while (first < nvtxs) {
    if (first == last) { /* Find another starting vertex */
      k = bfsperm[last];
      ASSERT(perm[k] != -1);
      perm[k] = -1; /* mark node as being visited */
      last++;
    }

    i = bfsperm[first++];
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      /* if a node has been already been visited, its perm[] will be -1 */
      if (perm[k] != -1) {
        /* perm[k] is the location within bfsperm of where k resides; 
           put in that location bfsperm[last] that we are about to
           overwrite and update perm[bfsperm[last]] to reflect that. */
        bfsperm[perm[k]]    = bfsperm[last];
        perm[bfsperm[last]] = perm[k];

        bfsperm[last++] = k;  /* put node at the end of the "queue" */
        perm[k]         = -1; /* mark node as being visited */
      }
    }
  }

  WCOREPOP;
}


/*************************************************************************/
/**************************************************************************/
idx_t IsConnected(graph_t *graph, idx_t report)
{
  idx_t ncmps;

  ncmps = FindPartitionInducedComponents(graph, NULL, NULL, NULL);

  if (ncmps != 1 && report)
    printf("The graph is not connected. It has %"PRIDX" connected components.\n", ncmps);

  return (ncmps == 1);
}


/*************************************************************************/
/*************************************************************************/
idx_t IsConnectedSubdomain(ctrl_t *ctrl, graph_t *graph, idx_t pid, idx_t report)
{
  idx_t i, j, k, nvtxs, first, last, nleft, ncmps, wgt;
  idx_t *xadj, *adjncy, *where, *touched, *queue;
  idx_t *cptr;

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  adjncy = graph->adjncy;
  where  = graph->where;

  touched = ismalloc(nvtxs, 0, "IsConnected: touched");
  queue   = imalloc(nvtxs, "IsConnected: queue");
  cptr    = imalloc(nvtxs+1, "IsConnected: cptr");

  nleft = 0;
  for (i=0; i<nvtxs; i++) {
    if (where[i] == pid) 
      nleft++;
  }

  for (i=0; i<nvtxs; i++) {
    if (where[i] == pid) 
      break;
  }

  touched[i] = 1;
  queue[0] = i;
  first = 0; last = 1;

  cptr[0] = 0;  /* This actually points to queue */
  ncmps = 0;
  while (first != nleft) {
    if (first == last) { /* Find another starting vertex */
      cptr[++ncmps] = first;
      for (i=0; i<nvtxs; i++) {
        if (where[i] == pid && !touched[i])
          break;
      }
      queue[last++] = i;
      touched[i] = 1;
    }

    i = queue[first++];
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      if (where[k] == pid && !touched[k]) {
        queue[last++] = k;
        touched[k] = 1;
      }
    }
  }
  cptr[++ncmps] = first;

  if (ncmps > 1 && report) {
    printf("The graph has %"PRIDX" connected components in partition %"PRIDX":\t", ncmps, pid);
    for (i=0; i<ncmps; i++) {
      wgt = 0;
      for (j=cptr[i]; j<cptr[i+1]; j++)
        wgt += graph->vwgt[queue[j]];
      printf("[%5"PRIDX" %5"PRIDX"] ", cptr[i+1]-cptr[i], wgt);
      /*
      if (cptr[i+1]-cptr[i] == 1)
        printf("[%"PRIDX" %"PRIDX"] ", queue[cptr[i]], xadj[queue[cptr[i]]+1]-xadj[queue[cptr[i]]]);
      */
    }
    printf("\n");
  }

  gk_free((void **)&touched, &queue, &cptr, LTERM);

  return (ncmps == 1 ? 1 : 0);
}


/*************************************************************************/
/**************************************************************************/
idx_t FindSepInducedComponents(ctrl_t *ctrl, graph_t *graph, idx_t *cptr, 
          idx_t *cind)
{
  idx_t i, j, k, nvtxs, first, last, nleft, ncmps, wgt;
  idx_t *xadj, *adjncy, *where, *touched, *queue;

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  adjncy = graph->adjncy;
  where  = graph->where;

  touched = ismalloc(nvtxs, 0, "IsConnected: queue");

  for (i=0; i<graph->nbnd; i++)
    touched[graph->bndind[i]] = 1;

  queue = cind;

  nleft = 0;
  for (i=0; i<nvtxs; i++) {
    if (where[i] != 2) 
      nleft++;
  }

  for (i=0; i<nvtxs; i++) {
    if (where[i] != 2)
      break;
  }

  touched[i] = 1;
  queue[0]   = i;
  first      = 0; 
  last       = 1;
  cptr[0]    = 0;  /* This actually points to queue */
  ncmps      = 0;

  while (first != nleft) {
    if (first == last) { /* Find another starting vertex */
      cptr[++ncmps] = first;
      for (i=0; i<nvtxs; i++) {
        if (!touched[i])
          break;
      }
      queue[last++] = i;
      touched[i] = 1;
    }

    i = queue[first++];
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      if (!touched[k]) {
        queue[last++] = k;
        touched[k] = 1;
      }
    }
  }
  cptr[++ncmps] = first;

  gk_free((void **)&touched, LTERM);

  return ncmps;
}


/*************************************************************************/
/*************************************************************************/
void EliminateComponents(ctrl_t *ctrl, graph_t *graph)
{
  idx_t i, ii, j, jj, k, me, nparts, nvtxs, ncon, ncmps, other, 
        ncand, target;
  idx_t *xadj, *adjncy, *vwgt, *adjwgt, *where, *pwgts;
  idx_t *cptr, *cind, *cpvec, *pcptr, *pcind, *cwhere;
  idx_t cid, bestcid, *cwgt, *bestcwgt;
  idx_t ntodo, oldntodo, *todo;
  rkv_t *cand;
  real_t *tpwgts;
  idx_t *vmarker=NULL, *pmarker=NULL, *modind=NULL;  /* volume specific work arrays */

  WCOREPUSH;

  nvtxs  = graph->nvtxs;
  ncon   = graph->ncon;
  xadj   = graph->xadj;
  adjncy = graph->adjncy;
  vwgt   = graph->vwgt;
  adjwgt = (ctrl->objtype == METIS_OBJTYPE_VOL ? NULL : graph->adjwgt);

  where = graph->where;
  pwgts = graph->pwgts;

  nparts = ctrl->nparts;
  tpwgts = ctrl->tpwgts;

  cptr = iwspacemalloc(ctrl, nvtxs+1);
  cind = iwspacemalloc(ctrl, nvtxs);

  ncmps = FindPartitionInducedComponents(graph, where, cptr, cind);

  IFSET(ctrl->dbglvl, METIS_DBG_CONTIGINFO, 
      printf("I found %"PRIDX" components, for this %"PRIDX"-way partition\n", 
          ncmps, nparts)); 

  /* There are more components than partitions */
  if (ncmps > nparts) {
    cwgt     = iwspacemalloc(ctrl, ncon);
    bestcwgt = iwspacemalloc(ctrl, ncon);
    cpvec    = iwspacemalloc(ctrl, nparts);
    pcptr    = iset(nparts+1, 0, iwspacemalloc(ctrl, nparts+1));
    pcind    = iwspacemalloc(ctrl, ncmps);
    cwhere   = iset(nvtxs, -1, iwspacemalloc(ctrl, nvtxs));
    todo     = iwspacemalloc(ctrl, ncmps);
    cand     = (rkv_t *)wspacemalloc(ctrl, nparts*sizeof(rkv_t));

    if (ctrl->objtype == METIS_OBJTYPE_VOL) {
      /* Vol-refinement specific working arrays */
      modind  = iwspacemalloc(ctrl, nvtxs);
      vmarker = iset(nvtxs, 0, iwspacemalloc(ctrl, nvtxs));
      pmarker = iset(nparts, -1, iwspacemalloc(ctrl, nparts));
    }


    /* Get a CSR representation of the components-2-partitions mapping */
    for (i=0; i<ncmps; i++) 
      pcptr[where[cind[cptr[i]]]]++;
    MAKECSR(i, nparts, pcptr);
    for (i=0; i<ncmps; i++) 
      pcind[pcptr[where[cind[cptr[i]]]]++] = i;
    SHIFTCSR(i, nparts, pcptr);

    /* Assign the heaviest component of each partition to its original partition */
    for (ntodo=0, i=0; i<nparts; i++) {
      if (pcptr[i+1]-pcptr[i] == 1)
        bestcid = pcind[pcptr[i]];
      else {
        for (bestcid=-1, j=pcptr[i]; j<pcptr[i+1]; j++) {
          cid = pcind[j];
          iset(ncon, 0, cwgt);
          for (ii=cptr[cid]; ii<cptr[cid+1]; ii++)
            iaxpy(ncon, 1, vwgt+cind[ii]*ncon, 1, cwgt, 1);
          if (bestcid == -1 || isum(ncon, bestcwgt, 1) < isum(ncon, cwgt, 1)) {
            bestcid  = cid;
            icopy(ncon, cwgt, bestcwgt);
          }
        }
        /* Keep track of those that need to be dealt with */
        for (j=pcptr[i]; j<pcptr[i+1]; j++) {
          if (pcind[j] != bestcid)
            todo[ntodo++] = pcind[j];
        }
      }

      for (j=cptr[bestcid]; j<cptr[bestcid+1]; j++) {
        ASSERT(where[cind[j]] == i);
        cwhere[cind[j]] = i;
      }
    }


    while (ntodo > 0) {
      oldntodo = ntodo;
      for (i=0; i<ntodo; i++) {
        cid = todo[i];
        me = where[cind[cptr[cid]]];  /* Get the domain of this component */

        /* Determine the weight of the block to be moved */
        iset(ncon, 0, cwgt);
        for (j=cptr[cid]; j<cptr[cid+1]; j++) 
          iaxpy(ncon, 1, vwgt+cind[j]*ncon, 1, cwgt, 1);

        IFSET(ctrl->dbglvl, METIS_DBG_CONTIGINFO, 
            printf("Trying to move %"PRIDX" [%"PRIDX"] from %"PRIDX"\n", 
                cid, isum(ncon, cwgt, 1), me)); 

        /* Determine the connectivity */
        iset(nparts, 0, cpvec);
        for (j=cptr[cid]; j<cptr[cid+1]; j++) {
          ii = cind[j];
          for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) 
            if (cwhere[adjncy[jj]] != -1)
              cpvec[cwhere[adjncy[jj]]] += (adjwgt ? adjwgt[jj] : 1);
        }

        /* Put the neighbors into a cand[] array for sorting */
        for (ncand=0, j=0; j<nparts; j++) {
          if (cpvec[j] > 0) {
            cand[ncand].key   = cpvec[j];
            cand[ncand++].val = j;
          }
        }
        if (ncand == 0)
          continue;

        rkvsortd(ncand, cand);

        /* Limit the moves to only the top candidates, which are defined as 
           those with connectivity at least 50% of the best.
           This applies only when ncon=1, as for multi-constraint, balancing
           will be hard. */
        if (ncon == 1) {
          for (j=1; j<ncand; j++) {
            if (cand[j].key < .5*cand[0].key)
              break;
          }
          ncand = j;
        }
      
        /* Now among those, select the one with the best balance */
        target = cand[0].val;
        for (j=1; j<ncand; j++) {
          if (BetterBalanceKWay(ncon, cwgt, ctrl->ubfactors,
                1, pwgts+target*ncon, ctrl->pijbm+target*ncon,
                1, pwgts+cand[j].val*ncon, ctrl->pijbm+cand[j].val*ncon))
            target = cand[j].val;
        }

        IFSET(ctrl->dbglvl, METIS_DBG_CONTIGINFO, 
            printf("\tMoving it to %"PRIDX" [%"PRIDX"] [%"PRIDX"]\n", target, cpvec[target], ncand));

        /* Note that as a result of a previous movement, a connected component may
           now will like to stay to its original partition */
        if (target != me) {
          switch (ctrl->objtype) {
            case METIS_OBJTYPE_CUT:
              MoveGroupContigForCut(ctrl, graph, target, cid, cptr, cind);
              break;

            case METIS_OBJTYPE_VOL:
              MoveGroupContigForVol(ctrl, graph, target, cid, cptr, cind, 
                  vmarker, pmarker, modind);
              break;

            default:
              gk_errexit(SIGERR, "Unknown objtype %d\n", ctrl->objtype);
          }
        }

        /* Update the cwhere vector */
        for (j=cptr[cid]; j<cptr[cid+1]; j++) 
          cwhere[cind[j]] = target;

        todo[i] = todo[--ntodo];
      }
      if (oldntodo == ntodo) {
        IFSET(ctrl->dbglvl, METIS_DBG_CONTIGINFO, printf("Stopped at ntodo: %"PRIDX"\n", ntodo));
        break;
      }
    }

    for (i=0; i<nvtxs; i++)
      ASSERT(where[i] == cwhere[i]);

  }

  WCOREPOP;
}


/*************************************************************************/
/*************************************************************************/
void MoveGroupContigForCut(ctrl_t *ctrl, graph_t *graph, idx_t to, idx_t gid, 
         idx_t *ptr, idx_t *ind)
{
  idx_t i, ii, iii, j, jj, k, l, nvtxs, nbnd, from, me;
  idx_t *xadj, *adjncy, *adjwgt, *where, *bndptr, *bndind;
  ckrinfo_t *myrinfo;
  cnbr_t *mynbrs;

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;

  where  = graph->where;
  bndptr = graph->bndptr;
  bndind = graph->bndind;

  nbnd = graph->nbnd;

  for (iii=ptr[gid]; iii<ptr[gid+1]; iii++) {
    i    = ind[iii];
    from = where[i];

    myrinfo = graph->ckrinfo+i;
    if (myrinfo->inbr == -1) {
      myrinfo->inbr = cnbrpoolGetNext(ctrl, xadj[i+1]-xadj[i]+1);
      myrinfo->nnbrs = 0;
    }
    mynbrs = ctrl->cnbrpool + myrinfo->inbr; 

    /* find the location of 'to' in myrinfo or create it if it is not there */
    for (k=0; k<myrinfo->nnbrs; k++) {
      if (mynbrs[k].pid == to)
        break;
    }
    if (k == myrinfo->nnbrs) {
      mynbrs[k].pid = to;
      mynbrs[k].ed = 0;
      myrinfo->nnbrs++;
    }

    graph->mincut -= mynbrs[k].ed-myrinfo->id;

    /* Update ID/ED and BND related information for the moved vertex */
    iaxpy(graph->ncon,  1, graph->vwgt+i*graph->ncon, 1, graph->pwgts+to*graph->ncon,   1);
    iaxpy(graph->ncon, -1, graph->vwgt+i*graph->ncon, 1, graph->pwgts+from*graph->ncon, 1);
    UpdateMovedVertexInfoAndBND(i, from, k, to, myrinfo, mynbrs, where, nbnd, 
        bndptr, bndind, BNDTYPE_REFINE);

    /* Update the degrees of adjacent vertices */
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      ii = adjncy[j];
      me = where[ii];
      myrinfo = graph->ckrinfo+ii;

      UpdateAdjacentVertexInfoAndBND(ctrl, ii, xadj[ii+1]-xadj[ii], me,
          from, to, myrinfo, adjwgt[j], nbnd, bndptr, bndind, BNDTYPE_REFINE);
    }

    ASSERT(CheckRInfo(ctrl, graph->ckrinfo+i));
  }

  graph->nbnd = nbnd;
}


/*************************************************************************/
/*************************************************************************/
void MoveGroupContigForVol(ctrl_t *ctrl, graph_t *graph, idx_t to, idx_t gid, 
         idx_t *ptr, idx_t *ind, idx_t *vmarker, idx_t *pmarker, 
         idx_t *modind)
{
  idx_t i, ii, iii, j, jj, k, l, nvtxs, from, me, other, xgain;
  idx_t *xadj, *vsize, *adjncy, *where;
  vkrinfo_t *myrinfo, *orinfo;
  vnbr_t *mynbrs, *onbrs;

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  vsize  = graph->vsize;
  adjncy = graph->adjncy;
  where  = graph->where;

  for (iii=ptr[gid]; iii<ptr[gid+1]; iii++) {
    i    = ind[iii];
    from = where[i];

    myrinfo = graph->vkrinfo+i;
    if (myrinfo->inbr == -1) {
      myrinfo->inbr = vnbrpoolGetNext(ctrl, xadj[i+1]-xadj[i]+1);
      myrinfo->nnbrs = 0;
    }
    mynbrs = ctrl->vnbrpool + myrinfo->inbr; 

    xgain = (myrinfo->nid == 0 && myrinfo->ned > 0 ? vsize[i] : 0);

    /* find the location of 'to' in myrinfo or create it if it is not there */
    for (k=0; k<myrinfo->nnbrs; k++) {
      if (mynbrs[k].pid == to)
        break;
    }
    if (k == myrinfo->nnbrs) {
      if (myrinfo->nid > 0)
        xgain -= vsize[i];

      /* determine the volume gain resulting from that move */
      for (j=xadj[i]; j<xadj[i+1]; j++) {
        ii     = adjncy[j];
        other  = where[ii];
        orinfo = graph->vkrinfo+ii;
        onbrs  = ctrl->vnbrpool + orinfo->inbr;
        ASSERT(other != to)

        if (from == other) {
          /* Same subdomain vertex: Decrease the gain if 'to' is a new neighbor. */
          for (l=0; l<orinfo->nnbrs; l++) {
            if (onbrs[l].pid == to)
              break;
          }
          if (l == orinfo->nnbrs) 
            xgain -= vsize[ii];
        }
        else {
          /* Remote vertex: increase if 'to' is a new subdomain */
          for (l=0; l<orinfo->nnbrs; l++) {
            if (onbrs[l].pid == to)
              break;
          }
          if (l == orinfo->nnbrs) 
            xgain -= vsize[ii];

          /* Remote vertex: decrease if i is the only connection to 'from' */
          for (l=0; l<orinfo->nnbrs; l++) {
            if (onbrs[l].pid == from && onbrs[l].ned == 1) {
              xgain += vsize[ii];
              break;
            }
          }
        }
      }
      graph->minvol -= xgain;
      graph->mincut -= -myrinfo->nid;
    }
    else {
      graph->minvol -= (xgain + mynbrs[k].gv);
      graph->mincut -= mynbrs[k].ned-myrinfo->nid;
    }


    /* Update where and pwgts */
    where[i] = to;
    iaxpy(graph->ncon,  1, graph->vwgt+i*graph->ncon, 1, graph->pwgts+to*graph->ncon,   1);
    iaxpy(graph->ncon, -1, graph->vwgt+i*graph->ncon, 1, graph->pwgts+from*graph->ncon, 1);

    /* Update the id/ed/gains/bnd of potentially affected nodes */
    KWayVolUpdate(ctrl, graph, i, from, to, NULL, NULL, NULL, NULL,
        NULL, BNDTYPE_REFINE, vmarker, pmarker, modind);

    /*CheckKWayVolPartitionParams(ctrl, graph);*/
  }

  ASSERT(ComputeCut(graph, where) == graph->mincut);
  ASSERTP(ComputeVolume(graph, where) == graph->minvol,
      ("%"PRIDX" %"PRIDX"\n", ComputeVolume(graph, where), graph->minvol));

}