C++ idxset函数代码示例

/*************************************************************************
* This function resets the buckets
**************************************************************************/
void PQueueReset(PQueueType *queue)
{
  int i, j;
  queue->nnodes = 0;

  if (queue->type == 1) {
    queue->maxgain = -queue->ngainspan;

    j = queue->ngainspan+queue->pgainspan+1;
    queue->buckets -= queue->ngainspan;  
    for (i=0; i<j; i++)
      queue->buckets[i] = NULL;
    queue->buckets += queue->ngainspan;  
  }
  else {
    idxset(queue->maxnodes, -1, queue->locator);
  }

}

/*************************************************************************
* This function takes a graph and produces a bisection by using a region
* growing algorithm. The resulting partition is returned in
* graph->where
**************************************************************************/
void MocGrowBisection(CtrlType *ctrl, GraphType *graph, float *tpwgts, float ubfactor)
{
  int i, j, k, nvtxs, ncon, from, bestcut, mincut, nbfs;
  idxtype *bestwhere, *where;

  nvtxs = graph->nvtxs;

  MocAllocate2WayPartitionMemory(ctrl, graph);
  where = graph->where;

  bestwhere = idxmalloc(nvtxs, "BisectGraph: bestwhere");
  nbfs = 2*(nvtxs <= ctrl->CoarsenTo ? SMALLNIPARTS : LARGENIPARTS);
  bestcut = idxsum(graph->nedges, graph->adjwgt);  

  for (; nbfs>0; nbfs--) {
    idxset(nvtxs, 1, where);
    where[RandomInRange(nvtxs)] = 0;

    MocCompute2WayPartitionParams(ctrl, graph);

    MocInit2WayBalance(ctrl, graph, tpwgts);

    MocFM_2WayEdgeRefine(ctrl, graph, tpwgts, 4); 

    MocBalance2Way(ctrl, graph, tpwgts, 1.02);
    MocFM_2WayEdgeRefine(ctrl, graph, tpwgts, 4); 

    if (bestcut >= graph->mincut) {
      bestcut = graph->mincut;
      idxcopy(nvtxs, where, bestwhere);
      if (bestcut == 0)
        break;
    }
  }

  graph->mincut = bestcut;
  idxcopy(nvtxs, bestwhere, where);

  /*GKfree(&bestwhere, LTERM);*/
  GKfree1((void**)&bestwhere);
}

/*************************************************************************
* This function computes the boundary definition for balancing
**************************************************************************/
void ComputeKWayBalanceBoundary(CtrlType *ctrl, GraphType *graph, int nparts)
{
  int i, nvtxs, nbnd;
  idxtype *bndind, *bndptr;

  nvtxs = graph->nvtxs;
  bndind = graph->bndind;
  bndptr = idxset(nvtxs, -1, graph->bndptr);


  /*------------------------------------------------------------
  / Compute the new boundary
  /------------------------------------------------------------*/
  nbnd = 0;
  for (i=0; i<nvtxs; i++) {
    if (graph->rinfo[i].ed > 0) 
      BNDInsert(nbnd, bndind, bndptr, i);
  }

  graph->nbnd = nbnd;
}

/*************************************************************************
* This function takes a graph and produces a bisection by using a region
* growing algorithm. The resulting partition is returned in
* graph->where
**************************************************************************/
void MocGrowBisection2(CtrlType *ctrl, GraphType *graph, float *tpwgts, float *ubvec)
{
  int /*i, j, k,*/ nvtxs, /*ncon, from,*/ bestcut, /*mincut,*/ nbfs;
  idxtype *bestwhere, *where;

  nvtxs = graph->nvtxs;

  MocAllocate2WayPartitionMemory(ctrl, graph);
  where = graph->where;

  bestwhere = idxmalloc(nvtxs, "BisectGraph: bestwhere");
  nbfs = 2*(nvtxs <= ctrl->CoarsenTo ? SMALLNIPARTS : LARGENIPARTS);
  bestcut = idxsum(graph->nedges, graph->adjwgt);  

  for (; nbfs>0; nbfs--) {
    idxset(nvtxs, 1, where);
    where[RandomInRange(nvtxs)] = 0;

    MocCompute2WayPartitionParams(ctrl, graph);

    MocBalance2Way2(ctrl, graph, tpwgts, ubvec);

    MocFM_2WayEdgeRefine2(ctrl, graph, tpwgts, ubvec, 4); 

    MocBalance2Way2(ctrl, graph, tpwgts, ubvec);
    MocFM_2WayEdgeRefine2(ctrl, graph, tpwgts, ubvec, 4); 

    if (bestcut > graph->mincut) {
      bestcut = graph->mincut;
      idxcopy(nvtxs, where, bestwhere);
      if (bestcut == 0)
        break;
    }
  }

  graph->mincut = bestcut;
  idxcopy(nvtxs, bestwhere, where);

  GKfree((void**)&bestwhere, LTERM);
}

/*************************************************************************
* This function computes the subdomain graph
**************************************************************************/
void ComputeSubDomainGraph(GraphType *graph, int nparts, idxtype *pmat, idxtype *ndoms)
{
  int i, j, k, me, nvtxs, ndegrees;
  idxtype *xadj, *adjncy, *adjwgt, *where;
  RInfoType *rinfo;
  EDegreeType *edegrees;

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

  idxset(nparts*nparts, 0, pmat);

  for (i=0; i<nvtxs; i++) {
    if (rinfo[i].ed > 0) {
      me = where[i];
      ndegrees = rinfo[i].ndegrees;
      edegrees = rinfo[i].edegrees;

      k = me*nparts;
      for (j=0; j<ndegrees; j++) 
        pmat[k+edegrees[j].pid] += edegrees[j].ed;
    }
  }

  for (i=0; i<nparts; i++) {
    ndoms[i] = 0;
    for (j=0; j<nparts; j++) {
      if (pmat[i*nparts+j] > 0)
        ndoms[i]++;
    }
  }

}

/*************************************************************************
* This function takes a graph and produces a bisection by using a region
* growing algorithm. The resulting partition is returned in
* graph->where
**************************************************************************/
void MocGrowBisectionNew2(CtrlType *ctrl, GraphType *graph, float *tpwgts, float *ubvec)
{
  idxtype i, j, k, nvtxs, ncon, from, bestcut, mincut, nbfs, inbfs;
  idxtype *bestwhere, *where;

  nvtxs = graph->nvtxs;

  MocAllocate2WayPartitionMemory(ctrl, graph);
  where = graph->where;

  bestwhere = idxmalloc(nvtxs, "BisectGraph: bestwhere");
  nbfs = 2*(nvtxs <= ctrl->CoarsenTo ? SMALLNIPARTS : LARGENIPARTS);

  for (inbfs=0; inbfs<nbfs; inbfs++) {
    idxset(nvtxs, 1, where);
    where[RandomInRange(nvtxs)] = 0;

    MocCompute2WayPartitionParams(ctrl, graph);

    MocInit2WayBalance2(ctrl, graph, tpwgts, ubvec);

    MocFM_2WayEdgeRefine2(ctrl, graph, tpwgts, ubvec, 4); 

    if (inbfs == 0 || bestcut > graph->mincut) {
      bestcut = graph->mincut;
      idxcopy(nvtxs, where, bestwhere);
      if (bestcut == 0)
        break;
    }
  }

  graph->mincut = bestcut;
  idxcopy(nvtxs, bestwhere, where);

  gk_free((void **)&bestwhere, LTERM);
}

/*************************************************************************
* This function finds all the connected components induced by the 
* partitioning vector in wgraph->where and tries to push them around to 
* remove some of them
**************************************************************************/
void EliminateComponents(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts, float ubfactor)
{
  int i, ii, j, jj, k, me, nvtxs, tvwgt, first, last, nleft, ncmps, cwgt, other, target, deltawgt;
  idxtype *xadj, *adjncy, *vwgt, *adjwgt, *where, *pwgts, *maxpwgt;
  idxtype *cpvec, *touched, *perm, *todo, *cind, *cptr, *npcmps;

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

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

  touched = idxset(nvtxs, 0, idxwspacemalloc(ctrl, nvtxs));
  cptr = idxwspacemalloc(ctrl, nvtxs);
  cind = idxwspacemalloc(ctrl, nvtxs);
  perm = idxwspacemalloc(ctrl, nvtxs);
  todo = idxwspacemalloc(ctrl, nvtxs);
  maxpwgt = idxwspacemalloc(ctrl, nparts);
  cpvec = idxwspacemalloc(ctrl, nparts);
  npcmps = idxset(nparts, 0, idxwspacemalloc(ctrl, nparts));

  for (i=0; i<nvtxs; i++) 
    perm[i] = todo[i] = i;

  /* 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];
      npcmps[me]++;
    }

    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;

  /* printf("I found %d components, for this %d-way partition\n", ncmps, nparts); */

  if (ncmps > nparts) { /* There are more components than processors */
    /* First determine the max allowed load imbalance */
    tvwgt = idxsum(nparts, pwgts);
    for (i=0; i<nparts; i++)
      maxpwgt[i] = ubfactor*tpwgts[i]*tvwgt;

    deltawgt = 5;

    for (i=0; i<ncmps; i++) {
      me = where[cind[cptr[i]]];  /* Get the domain of this component */
      if (npcmps[me] == 1)
        continue;  /* Skip it because it is contigous */

      /*printf("Trying to move %d from %d\n", i, me); */

      /* Determine the weight of the block to be moved and abort if too high */
      for (cwgt=0, j=cptr[i]; j<cptr[i+1]; j++) 
        cwgt += vwgt[cind[j]];

      if (cwgt > .30*pwgts[me])
        continue;  /* Skip the component if it is over 30% of the weight */

      /* Determine the connectivity */
      idxset(nparts, 0, cpvec);
      for (j=cptr[i]; j<cptr[i+1]; j++) {
        ii = cind[j];
        for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) 
          cpvec[where[adjncy[jj]]] += adjwgt[jj];
      }
      cpvec[me] = 0;

      target = -1;
      for (j=0; j<nparts; j++) {
        if (cpvec[j] > 0 && (cwgt < deltawgt || pwgts[j] + cwgt < maxpwgt[j])) {
          if (target == -1 || cpvec[target] < cpvec[j])
            target = j;
//.........这里部分代码省略.........

/*************************************************************************
* This function creates the coarser graph
**************************************************************************/
void CreateCoarseGraphNoMask(CtrlType *ctrl, GraphType *graph, int cnvtxs, idxtype *match, idxtype *perm)
{
  int i, j, k, m, istart, iend, nvtxs, nedges, ncon, cnedges, v, u, dovsize;
  idxtype *xadj, *vwgt, *vsize, *adjncy, *adjwgt, *adjwgtsum, *auxadj;
  idxtype *cmap, *htable;
  idxtype *cxadj, *cvwgt, *cvsize, *cadjncy, *cadjwgt, *cadjwgtsum;
  float *nvwgt, *cnvwgt;
  GraphType *cgraph;

  dovsize = (ctrl->optype == OP_KVMETIS ? 1 : 0);

  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->ContractTmr));

  nvtxs = graph->nvtxs;
  ncon = graph->ncon;
  xadj = graph->xadj;
  vwgt = graph->vwgt;
  vsize = graph->vsize;
  nvwgt = graph->nvwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  adjwgtsum = graph->adjwgtsum;
  cmap = graph->cmap;


  /* Initialize the coarser graph */
  cgraph = SetUpCoarseGraph(graph, cnvtxs, dovsize);
  cxadj = cgraph->xadj;
  cvwgt = cgraph->vwgt;
  cvsize = cgraph->vsize;
  cnvwgt = cgraph->nvwgt;
  cadjwgtsum = cgraph->adjwgtsum;
  cadjncy = cgraph->adjncy;
  cadjwgt = cgraph->adjwgt;


  htable = idxset(cnvtxs, -1, idxwspacemalloc(ctrl, cnvtxs));

  iend = xadj[nvtxs];
  auxadj = ctrl->wspace.auxcore; 
  memcpy(auxadj, adjncy, iend*sizeof(idxtype)); 
  for (i=0; i<iend; i++)
    auxadj[i] = cmap[auxadj[i]];

  cxadj[0] = cnvtxs = cnedges = 0;
  for (i=0; i<nvtxs; i++) {
    v = perm[i];
    if (cmap[v] != cnvtxs) 
      continue;

    u = match[v];
    if (ncon == 1)
      cvwgt[cnvtxs] = vwgt[v];
    else
      scopy(ncon, nvwgt+v*ncon, cnvwgt+cnvtxs*ncon);

    if (dovsize)
      cvsize[cnvtxs] = vsize[v];

    cadjwgtsum[cnvtxs] = adjwgtsum[v];
    nedges = 0;

    istart = xadj[v];
    iend = xadj[v+1];
    for (j=istart; j<iend; j++) {
      k = auxadj[j];
      if ((m = htable[k]) == -1) {
        cadjncy[nedges] = k;
        cadjwgt[nedges] = adjwgt[j];
        htable[k] = nedges++;
      }
      else {
        cadjwgt[m] += adjwgt[j];
      }
    }

    if (v != u) { 
      if (ncon == 1)
        cvwgt[cnvtxs] += vwgt[u];
      else
        saxpy(ncon, 1.0, nvwgt+u*ncon, 1, cnvwgt+cnvtxs*ncon, 1);

      if (dovsize)
        cvsize[cnvtxs] += vsize[u];

      cadjwgtsum[cnvtxs] += adjwgtsum[u];

      istart = xadj[u];
      iend = xadj[u+1];
      for (j=istart; j<iend; j++) {
        k = auxadj[j];
        if ((m = htable[k]) == -1) {
          cadjncy[nedges] = k;
          cadjwgt[nedges] = adjwgt[j];
          htable[k] = nedges++;
        }
        else {
//.........这里部分代码省略.........

/*************************************************************************
* This function compresses a graph by merging identical vertices
* The compression should lead to at least 10% reduction.
**************************************************************************/
void CompressGraph(CtrlType *ctrl, GraphType *graph, int nvtxs, idxtype *xadj, idxtype *adjncy, idxtype *cptr, idxtype *cind)
{
  int i, ii, iii, j, jj, k, l, cnvtxs, cnedges;
  idxtype *cxadj, *cadjncy, *cvwgt, *mark, *map;
  KeyValueType *keys;

  mark = idxsmalloc(nvtxs, -1, "CompressGraph: mark");
  map = idxsmalloc(nvtxs, -1, "CompressGraph: map");
  keys = (KeyValueType *)GKmalloc(nvtxs*sizeof(KeyValueType), "CompressGraph: keys");

  /* Compute a key for each adjacency list */
  for (i=0; i<nvtxs; i++) {
    k = 0;
    for (j=xadj[i]; j<xadj[i+1]; j++)
      k += adjncy[j];
    keys[i].key = k+i; /* Add the diagonal entry as well */
    keys[i].val = i;
  }

  ikeysort(nvtxs, keys);

  l = cptr[0] = 0;
  for (cnvtxs=i=0; i<nvtxs; i++) {
    ii = keys[i].val;
    if (map[ii] == -1) { 
      mark[ii] = i;  /* Add the diagonal entry */
      for (j=xadj[ii]; j<xadj[ii+1]; j++) 
        mark[adjncy[j]] = i;

      cind[l++] = ii;
      map[ii] = cnvtxs;

      for (j=i+1; j<nvtxs; j++) {
        iii = keys[j].val;

        if (keys[i].key != keys[j].key || xadj[ii+1]-xadj[ii] != xadj[iii+1]-xadj[iii])
          break; /* Break if keys or degrees are different */

        if (map[iii] == -1) { /* Do a comparison if iii has not been mapped */ 
          for (jj=xadj[iii]; jj<xadj[iii+1]; jj++) {
            if (mark[adjncy[jj]] != i)
              break;
          }

          if (jj == xadj[iii+1]) { /* Identical adjacency structure */
            map[iii] = cnvtxs;
            cind[l++] = iii;
          }
        }
      }

      cptr[++cnvtxs] = l;
    }
  }

  /* printf("Original: %6d, Compressed: %6d\n", nvtxs, cnvtxs); */


  InitGraph(graph);

  if (cnvtxs >= COMPRESSION_FRACTION*nvtxs) {
    graph->nvtxs = nvtxs;
    graph->nedges = xadj[nvtxs];
    graph->ncon = 1;
    graph->xadj = xadj;
    graph->adjncy = adjncy;

    graph->gdata = idxmalloc(3*nvtxs+graph->nedges, "CompressGraph: gdata");
    graph->vwgt    	= graph->gdata;
    graph->adjwgtsum    = graph->gdata+nvtxs;
    graph->cmap		= graph->gdata+2*nvtxs;
    graph->adjwgt	= graph->gdata+3*nvtxs;

    idxset(nvtxs, 1, graph->vwgt);
    idxset(graph->nedges, 1, graph->adjwgt);
    for (i=0; i<nvtxs; i++)
      graph->adjwgtsum[i] = xadj[i+1]-xadj[i];

    graph->label = idxmalloc(nvtxs, "CompressGraph: label");
    for (i=0; i<nvtxs; i++)
      graph->label[i] = i;
  }
  else { /* Ok, form the compressed graph  */
    cnedges = 0;
    for (i=0; i<cnvtxs; i++) {
      ii = cind[cptr[i]];
      cnedges += xadj[ii+1]-xadj[ii];
    }

    /* Allocate memory for the compressed graph*/
    graph->gdata = idxmalloc(4*cnvtxs+1 + 2*cnedges, "CompressGraph: gdata");
    cxadj = graph->xadj		= graph->gdata;
    cvwgt = graph->vwgt         = graph->gdata + cnvtxs+1;
    graph->adjwgtsum        	= graph->gdata + 2*cnvtxs+1;
    graph->cmap                 = graph->gdata + 3*cnvtxs+1;
    cadjncy = graph->adjncy     = graph->gdata + 4*cnvtxs+1;
//.........这里部分代码省略.........

/*************************************************************************
* This function performs k-way refinement
**************************************************************************/
void MCGreedy_KWayEdgeBalanceHorizontal(CtrlType *ctrl, GraphType *graph, int nparts, 
       float *ubvec, int npasses)
{
  int i, ii, /*iii,*/ j, /*jj,*/ k, /*l,*/ pass, nvtxs, ncon, nbnd, myndegrees, oldgain, gain, nmoves; 
  int from, me, to, oldcut;
  idxtype *xadj, *adjncy, *adjwgt;
  idxtype *where, *perm, *bndptr, *bndind, *moved;
  EDegreeType *myedegrees;
  RInfoType *myrinfo;
  PQueueType queue;
  float *npwgts, *nvwgt, *minwgt, *maxwgt, tvec[MAXNCON];

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

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

  where = graph->where;
  npwgts = graph->npwgts;
  
  /* Setup the weight intervals of the various subdomains */
  minwgt =  fwspacemalloc(ctrl, ncon*nparts);
  maxwgt = fwspacemalloc(ctrl, ncon*nparts);

  for (i=0; i<nparts; i++) {
    for (j=0; j<ncon; j++) {
      maxwgt[i*ncon+j] = ubvec[j]/nparts;
      minwgt[i*ncon+j] = 1.0/(ubvec[j]*nparts);
    }
  }

  perm = idxwspacemalloc(ctrl, nvtxs);
  moved = idxwspacemalloc(ctrl, nvtxs);

  PQueueInit(ctrl, &queue, nvtxs, graph->adjwgtsum[idxamax(nvtxs, graph->adjwgtsum)]);

  if (ctrl->dbglvl&DBG_REFINE) {
    printf("Partitions: [%5.4f %5.4f], Nv-Nb[%6d %6d]. Cut: %6d, LB: ",
            npwgts[samin(ncon*nparts, npwgts)], npwgts[samax(ncon*nparts, npwgts)], 
            graph->nvtxs, graph->nbnd, graph->mincut);
    ComputeHKWayLoadImbalance(ncon, nparts, npwgts, tvec);
    for (i=0; i<ncon; i++)
      printf("%.3f ", tvec[i]);
    printf("[B]\n");
  }


  for (pass=0; pass<npasses; pass++) {
    ASSERT(ComputeCut(graph, where) == graph->mincut);

    /* Check to see if things are out of balance, given the tolerance */
    if (MocIsHBalanced(ncon, nparts, npwgts, ubvec))
      break;

    PQueueReset(&queue);
    idxset(nvtxs, -1, moved);

    oldcut = graph->mincut;
    nbnd = graph->nbnd;

    RandomPermute(nbnd, perm, 1);
    for (ii=0; ii<nbnd; ii++) {
      i = bndind[perm[ii]];
      PQueueInsert(&queue, i, graph->rinfo[i].ed - graph->rinfo[i].id);
      moved[i] = 2;
    }

    nmoves = 0;
    for (;;) {
      if ((i = PQueueGetMax(&queue)) == -1) 
        break;
      moved[i] = 1;

      myrinfo = graph->rinfo+i;
      from = where[i];
      nvwgt = graph->nvwgt+i*ncon;

      if (AreAllHVwgtsBelow(ncon, 1.0, npwgts+from*ncon, -1.0, nvwgt, minwgt+from*ncon))
        continue;   /* This cannot be moved! */

      myedegrees = myrinfo->edegrees;
      myndegrees = myrinfo->ndegrees;

      for (k=0; k<myndegrees; k++) {
        to = myedegrees[k].pid;
        if (IsHBalanceBetterFT(ncon, nparts, npwgts+from*ncon, npwgts+to*ncon, nvwgt, ubvec))
          break;
      }
      if (k == myndegrees) 
        continue;  /* break out if you did not find a candidate */

      for (j=k+1; j<myndegrees; j++) {
        to = myedegrees[j].pid;
//.........这里部分代码省略.........

/*************************************************************************
* This function projects a partition, and at the same time computes the
* parameters for refinement.
**************************************************************************/
void MocProject2WayPartition(CtrlType *ctrl, GraphType *graph)
{
    int i, j, k, nvtxs, nbnd, me;
    idxtype *xadj, *adjncy, *adjwgt, *adjwgtsum;
    idxtype *cmap, *where, *id, *ed, *bndptr, *bndind;
    idxtype *cwhere, *cid, *ced, *cbndptr;
    GraphType *cgraph;

    cgraph = graph->coarser;
    cwhere = cgraph->where;
    cid = cgraph->id;
    ced = cgraph->ed;
    cbndptr = cgraph->bndptr;

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

    MocAllocate2WayPartitionMemory(ctrl, graph);

    where = graph->where;
    id = idxset(nvtxs, 0, graph->id);
    ed = idxset(nvtxs, 0, graph->ed);
    bndptr = idxset(nvtxs, -1, graph->bndptr);
    bndind = graph->bndind;


    /* Go through and project partition and compute id/ed for the nodes */
    for (i=0; i<nvtxs; i++) {
        k = cmap[i];
        where[i] = cwhere[k];
        cmap[i] = cbndptr[k];
    }

    for (nbnd=0, i=0; i<nvtxs; i++) {
        me = where[i];

        id[i] = adjwgtsum[i];

        if (xadj[i] == xadj[i+1]) {
            bndptr[i] = nbnd;
            bndind[nbnd++] = i;
        }
        else {
            if (cmap[i] != -1) { /* If it is an interface node. Note that cmap[i] = cbndptr[cmap[i]] */
                for (j=xadj[i]; j<xadj[i+1]; j++) {
                    if (me != where[adjncy[j]])
                        ed[i] += adjwgt[j];
                }
                id[i] -= ed[i];

                if (ed[i] > 0 || xadj[i] == xadj[i+1]) {
                    bndptr[i] = nbnd;
                    bndind[nbnd++] = i;
                }
            }
        }
    }

    graph->mincut = cgraph->mincut;
    graph->nbnd = nbnd;
    scopy(2*graph->ncon, cgraph->npwgts, graph->npwgts);

    FreeGraph(graph->coarser);
    graph->coarser = NULL;

}

/*************************************************************************
* This function performs k-way refinement
**************************************************************************/
void Greedy_KWayEdgeRefine(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts, float ubfactor, int npasses)
{
    int i, ii, iii, j, jj, k, l, pass, nvtxs, nbnd, tvwgt, myndegrees, oldgain, gain;
    int from, me, to, oldcut, vwgt;
    idxtype *xadj, *adjncy, *adjwgt;
    idxtype *where, *pwgts, *perm, *bndptr, *bndind, *minwgt, *maxwgt, *moved, *itpwgts;
    EDegreeType *myedegrees;
    RInfoType *myrinfo;
    PQueueType queue;

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

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

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

    /* Setup the weight intervals of the various subdomains */
    minwgt =  idxwspacemalloc(ctrl, nparts);
    maxwgt = idxwspacemalloc(ctrl, nparts);
    itpwgts = idxwspacemalloc(ctrl, nparts);
    tvwgt = idxsum(nparts, pwgts);
    ASSERT(tvwgt == idxsum(nvtxs, graph->vwgt));

    for (i=0; i<nparts; i++) {
        itpwgts[i] = tpwgts[i]*tvwgt;
        maxwgt[i] = tpwgts[i]*tvwgt*ubfactor;
        minwgt[i] = tpwgts[i]*tvwgt*(1.0/ubfactor);
    }

    perm = idxwspacemalloc(ctrl, nvtxs);
    moved = idxwspacemalloc(ctrl, nvtxs);

    PQueueInit(ctrl, &queue, nvtxs, graph->adjwgtsum[idxamax(nvtxs, graph->adjwgtsum)]);

    IFSET(ctrl->dbglvl, DBG_REFINE,
          printf("Partitions: [%6d %6d]-[%6d %6d], Balance: %5.3f, Nv-Nb[%6d %6d]. Cut: %6d\n",
                 pwgts[idxamin(nparts, pwgts)], pwgts[idxamax(nparts, pwgts)], minwgt[0], maxwgt[0],
                 1.0*nparts*pwgts[idxamax(nparts, pwgts)]/tvwgt, graph->nvtxs, graph->nbnd,
                 graph->mincut));

    for (pass=0; pass<npasses; pass++) {
        ASSERT(ComputeCut(graph, where) == graph->mincut);

        PQueueReset(&queue);
        idxset(nvtxs, -1, moved);

        oldcut = graph->mincut;
        nbnd = graph->nbnd;

        RandomPermute(nbnd, perm, 1);
        for (ii=0; ii<nbnd; ii++) {
            i = bndind[perm[ii]];
            PQueueInsert(&queue, i, graph->rinfo[i].ed - graph->rinfo[i].id);
            moved[i] = 2;
        }

        for (iii=0;; iii++) {
            if ((i = PQueueGetMax(&queue)) == -1)
                break;
            moved[i] = 1;

            myrinfo = graph->rinfo+i;
            from = where[i];
            vwgt = graph->vwgt[i];

            if (pwgts[from]-vwgt < minwgt[from])
                continue;   /* This cannot be moved! */

            myedegrees = myrinfo->edegrees;
            myndegrees = myrinfo->ndegrees;

            j = myrinfo->id;
            for (k=0; k<myndegrees; k++) {
                to = myedegrees[k].pid;
                gain = myedegrees[k].ed-j; /* j = myrinfo->id. Allow good nodes to move */
                if (pwgts[to]+vwgt <= maxwgt[to]+gain && gain >= 0)
                    break;
            }
            if (k == myndegrees)
                continue;  /* break out if you did not find a candidate */

            for (j=k+1; j<myndegrees; j++) {
                to = myedegrees[j].pid;
                if ((myedegrees[j].ed > myedegrees[k].ed && pwgts[to]+vwgt <= maxwgt[to]) ||
                        (myedegrees[j].ed == myedegrees[k].ed &&
                         itpwgts[myedegrees[k].pid]*pwgts[to] < itpwgts[to]*pwgts[myedegrees[k].pid]))
                    k = j;
            }

            to = myedegrees[k].pid;

            j = 0;
//.........这里部分代码省略.........

/*************************************************************************
* This function projects a partition, and at the same time computes the
* parameters for refinement.
**************************************************************************/
void ProjectKWayPartition(CtrlType *ctrl, GraphType *graph, int nparts)
{
  int i, j, k, nvtxs, nbnd, me, other, istart, iend, ndegrees;
  idxtype *xadj, *adjncy, *adjwgt, *adjwgtsum;
  idxtype *cmap, *where, *bndptr, *bndind;
  idxtype *cwhere;
  GraphType *cgraph;
  RInfoType *crinfo, *rinfo, *myrinfo;
  EDegreeType *myedegrees;
  idxtype *htable;

  cgraph = graph->coarser;
  cwhere = cgraph->where;
  crinfo = cgraph->rinfo;

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

  AllocateKWayPartitionMemory(ctrl, graph, nparts);
  where = graph->where;
  rinfo = graph->rinfo;
  bndind = graph->bndind;
  bndptr = idxset(nvtxs, -1, graph->bndptr);

  /* Go through and project partition and compute id/ed for the nodes */
  for (i=0; i<nvtxs; i++) {
    k = cmap[i];
    where[i] = cwhere[k];
    cmap[i] = crinfo[k].ed;  /* For optimization */
  }

  htable = idxset(nparts, -1, idxwspacemalloc(ctrl, nparts));

  ctrl->wspace.cdegree = 0;
  for (nbnd=0, i=0; i<nvtxs; i++) {
    me = where[i];

    myrinfo = rinfo+i;
    myrinfo->id = myrinfo->ed = myrinfo->ndegrees = 0;
    myrinfo->edegrees = NULL;

    myrinfo->id = adjwgtsum[i];

    if (cmap[i] > 0) { /* If it is an interface node. Note cmap[i] = crinfo[cmap[i]].ed */
      istart = xadj[i];
      iend = xadj[i+1];

      myedegrees = myrinfo->edegrees = ctrl->wspace.edegrees+ctrl->wspace.cdegree;
      ctrl->wspace.cdegree += iend-istart;

      ndegrees = 0;
      for (j=istart; j<iend; j++) {
        other = where[adjncy[j]];
        if (me != other) {
          myrinfo->ed += adjwgt[j];
          if ((k = htable[other]) == -1) {
            htable[other] = ndegrees;
            myedegrees[ndegrees].pid = other;
            myedegrees[ndegrees++].ed = adjwgt[j];
          }
          else {
            myedegrees[k].ed += adjwgt[j];
          }
        }
      }
      myrinfo->id -= myrinfo->ed;

      /* Remove space for edegrees if it was interior */
      if (myrinfo->ed == 0) { 
        myrinfo->edegrees = NULL;
        ctrl->wspace.cdegree -= iend-istart;
      }
      else {
        if (myrinfo->ed-myrinfo->id >= 0) 
          BNDInsert(nbnd, bndind, bndptr, i); 

        myrinfo->ndegrees = ndegrees;

        for (j=0; j<ndegrees; j++)
          htable[myedegrees[j].pid] = -1;
      }
    }
  }

  idxcopy(nparts, cgraph->pwgts, graph->pwgts);
  graph->mincut = cgraph->mincut;
  graph->nbnd = nbnd;

  FreeGraph(graph->coarser);
  graph->coarser = NULL;

  idxwspacefree(ctrl, nparts);
//.........这里部分代码省略.........

/*************************************************************************
* This function computes the initial id/ed 
**************************************************************************/
void ComputeKWayPartitionParams(CtrlType *ctrl, GraphType *graph, int nparts)
{
  int i, j, k, nvtxs, nbnd, mincut, me, other;
  idxtype *xadj, *vwgt, *adjncy, *adjwgt, *pwgts, *where, *bndind, *bndptr;
  RInfoType *rinfo, *myrinfo;
  EDegreeType *myedegrees;

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

  where = graph->where;
  pwgts = idxset(nparts, 0, graph->pwgts);
  bndind = graph->bndind;
  bndptr = idxset(nvtxs, -1, graph->bndptr);
  rinfo = graph->rinfo;


  /*------------------------------------------------------------
  / Compute now the id/ed degrees
  /------------------------------------------------------------*/
  ctrl->wspace.cdegree = 0;
  nbnd = mincut = 0;
  for (i=0; i<nvtxs; i++) {
    me = where[i];
    pwgts[me] += vwgt[i];

    myrinfo = rinfo+i;
    myrinfo->id = myrinfo->ed = myrinfo->ndegrees = 0;
    myrinfo->edegrees = NULL;

    for (j=xadj[i]; j<xadj[i+1]; j++) {
      if (me != where[adjncy[j]])
        myrinfo->ed += adjwgt[j];
    }
    myrinfo->id = graph->adjwgtsum[i] - myrinfo->ed;

    if (myrinfo->ed > 0) 
      mincut += myrinfo->ed;

    if (myrinfo->ed-myrinfo->id >= 0)
      BNDInsert(nbnd, bndind, bndptr, i);

    /* Time to compute the particular external degrees */
    if (myrinfo->ed > 0) { 
      myedegrees = myrinfo->edegrees = ctrl->wspace.edegrees+ctrl->wspace.cdegree;
      ctrl->wspace.cdegree += xadj[i+1]-xadj[i];

      for (j=xadj[i]; j<xadj[i+1]; j++) {
        other = where[adjncy[j]];
        if (me != other) {
          for (k=0; k<myrinfo->ndegrees; k++) {
            if (myedegrees[k].pid == other) {
              myedegrees[k].ed += adjwgt[j];
              break;
            }
          }
          if (k == myrinfo->ndegrees) {
            myedegrees[myrinfo->ndegrees].pid = other;
            myedegrees[myrinfo->ndegrees++].ed = adjwgt[j];
          }
        }
      }

      ASSERT(myrinfo->ndegrees <= xadj[i+1]-xadj[i]);
    }
  }

  graph->mincut = mincut/2;
  graph->nbnd = nbnd;

}

void METIS_PartMeshDual_WV(int *ne, int *nn, idxtype *elmnts, int *etype, int *numflag, 
                        int *nparts, int *edgecut, idxtype *epart, idxtype *npart, idxtype *vwgts)
{
  int i, j, k, me;
  idxtype *xadj, *adjncy, *pwgts, *nptr, *nind;
  int options[10], pnumflag=0, wgtflag=2;
  int nnbrs, nbrind[200], nbrwgt[200], maxpwgt;
  int esize, esizes[] = {-1, 3, 4, 8, 4};

  esize = esizes[*etype];

  if (*numflag == 1)
    ChangeMesh2CNumbering((*ne)*esize, elmnts);

  xadj = idxmalloc(*ne+1, "METIS_MESHPARTNODAL: xadj");
  adjncy = idxmalloc(esize*(*ne), "METIS_MESHPARTNODAL: adjncy");

  METIS_MeshToDual(ne, nn, elmnts, etype, &pnumflag, xadj, adjncy);

  options[0] = 0;
  METIS_PartGraphKway(ne, xadj, adjncy, vwgts, NULL, &wgtflag, &pnumflag, nparts, options, edgecut, epart);

  /* Construct the node-element list */
  nptr = idxsmalloc(*nn+1, 0, "METIS_MESHPARTDUAL: nptr");
  for (j=esize*(*ne), i=0; i<j; i++) 
    nptr[elmnts[i]]++;
  MAKECSR(i, *nn, nptr);

  nind = idxmalloc(nptr[*nn], "METIS_MESHPARTDUAL: nind");
  for (k=i=0; i<(*ne); i++) {
    for (j=0; j<esize; j++, k++) 
      nind[nptr[elmnts[k]]++] = i;
  }
  for (i=(*nn); i>0; i--)
    nptr[i] = nptr[i-1];
  nptr[0] = 0;


  /* OK, now compute a nodal partition based on the element partition npart */
  idxset(*nn, -1, npart);
  pwgts = idxsmalloc(*nparts, 0, "METIS_MESHPARTDUAL: pwgts");
  for (i=0; i<*nn; i++) {
    me = epart[nind[nptr[i]]];
    for (j=nptr[i]+1; j<nptr[i+1]; j++) {
      if (epart[nind[j]] != me)
        break;
    }
    if (j == nptr[i+1]) {
      npart[i] = me;
      pwgts[me]++;
    }
  }

  maxpwgt = 1.03*(*nn)/(*nparts);
  for (i=0; i<*nn; i++) {
    if (npart[i] == -1) { /* Assign the boundary element */
      nnbrs = 0;
      for (j=nptr[i]; j<nptr[i+1]; j++) {
        me = epart[nind[j]];
        for (k=0; k<nnbrs; k++) {
          if (nbrind[k] == me) {
            nbrwgt[k]++;
            break;
          }
        }
        if (k == nnbrs) {
          nbrind[nnbrs] = me;
          nbrwgt[nnbrs++] = 1;
        }
      }
      /* Try to assign it first to the domain with most things in common */
      j = iamax(nnbrs, nbrwgt);
      if (pwgts[nbrind[j]] < maxpwgt) {
        npart[i] = nbrind[j];
      }
      else {
        /* If that fails, assign it to a light domain */
        npart[i] = nbrind[0];
        for (j=0; j<nnbrs; j++) {
          if (pwgts[nbrind[j]] < maxpwgt) {
            npart[i] = nbrind[j];
            break;
          }
        }
      }
      pwgts[npart[i]]++;
    }
  }

  if (*numflag == 1)
    ChangeMesh2FNumbering2((*ne)*esize, elmnts, *ne, *nn, epart, npart);

  GKfree(&xadj, &adjncy, &pwgts, &nptr, &nind, LTERM);

}