/*************************************************************************
* This function computes the balance of the partitioning
**************************************************************************/
void ComputePartitionBalance(GraphType *graph, int nparts, idxtype *where, float *ubvec)
{
  int i, j, nvtxs, ncon;
  idxtype *kpwgts, *vwgt;
  /*float balance;*/

  nvtxs = graph->nvtxs;
  ncon = graph->ncon;
  vwgt = graph->vwgt;

  kpwgts = idxsmalloc(nparts, 0, "ComputePartitionInfo: kpwgts");

  if (vwgt == NULL) {
    for (i=0; i<nvtxs; i++)
      kpwgts[where[i]]++;
    ubvec[0] = 1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*nvtxs);
  }
  else {
    for (j=0; j<ncon; j++) {
      idxset(nparts, 0, kpwgts);
      for (i=0; i<graph->nvtxs; i++)
        kpwgts[where[i]] += vwgt[i*ncon+j];

      ubvec[j] = 1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts));
    }
  }

  free(kpwgts);

}

/*************************************************************************
* This function is the driver for the partition refinement mode of ParMETIS
**************************************************************************/
void Order_Partition(CtrlType *ctrl, GraphType *graph, WorkSpaceType *wspace)
{

  SetUp(ctrl, graph, wspace);
  graph->ncon = 1;

  IFSET(ctrl->dbglvl, DBG_PROGRESS, rprintf(ctrl, "[%6d %8d %5d %5d][%d][%d]\n",
        graph->gnvtxs, GlobalSESum(ctrl, graph->nedges), GlobalSEMin(ctrl, graph->nvtxs),
        GlobalSEMax(ctrl, graph->nvtxs), ctrl->CoarsenTo,
        GlobalSEMax(ctrl, graph->vwgt[idxamax(graph->nvtxs, graph->vwgt)])));

  if (graph->gnvtxs < 1.3*ctrl->CoarsenTo || (graph->finer != NULL && graph->gnvtxs > graph->finer->gnvtxs*COARSEN_FRACTION)) {
    /* Compute the initial npart-way multisection */
    InitMultisection(ctrl, graph, wspace);

    if (graph->finer == NULL) { /* Do that only of no-coarsening took place */
      ComputeNodePartitionParams(ctrl, graph, wspace);
      KWayNodeRefine(ctrl, graph, wspace, 2*NGR_PASSES, ORDER_UNBALANCE_FRACTION);
    }
  }
  else { /* Coarsen it and the partition it */
    Mc_LocalMatch_HEM(ctrl, graph, wspace);

    Order_Partition(ctrl, graph->coarser, wspace);

    Moc_ProjectPartition(ctrl, graph, wspace);
    ComputeNodePartitionParams(ctrl, graph, wspace);
    KWayNodeRefine(ctrl, graph, wspace, 2*NGR_PASSES, ORDER_UNBALANCE_FRACTION);
  }
}

/*************************************************************************
* This function computes movement statistics for adaptive refinement
* schemes
**************************************************************************/
void Mc_ComputeMoveStatistics(CtrlType *ctrl, GraphType *graph, int *nmoved, int *maxin, int *maxout)
{
  int i, nvtxs, nparts, myhome;
  idxtype *vwgt, *where;
  idxtype *lend, *gend, *lleft, *gleft, *lstart, *gstart;

  nvtxs = graph->nvtxs;
  vwgt = graph->vwgt;
  where = graph->where;
  nparts = ctrl->nparts;

  lstart = idxsmalloc(nparts, 0, "ComputeMoveStatistics: lstart");
  gstart = idxsmalloc(nparts, 0, "ComputeMoveStatistics: gstart");
  lleft = idxsmalloc(nparts, 0, "ComputeMoveStatistics: lleft");
  gleft = idxsmalloc(nparts, 0, "ComputeMoveStatistics: gleft");
  lend = idxsmalloc(nparts, 0, "ComputeMoveStatistics: lend");
  gend = idxsmalloc(nparts, 0, "ComputeMoveStatistics: gend");

  for (i=0; i<nvtxs; i++) {
    myhome = (ctrl->ps_relation == COUPLED) ? ctrl->mype : graph->home[i];
    lstart[myhome] += (graph->vsize == NULL) ? 1 : graph->vsize[i];
    lend[where[i]] += (graph->vsize == NULL) ? 1 : graph->vsize[i];
    if (where[i] != myhome)
      lleft[myhome] += (graph->vsize == NULL) ? 1 : graph->vsize[i];
  }

  /* PrintVector(ctrl, ctrl->npes, 0, lend, "Lend: "); */

  MPI_Allreduce((void *)lstart, (void *)gstart, nparts, IDX_DATATYPE, MPI_SUM, ctrl->comm);
  MPI_Allreduce((void *)lleft, (void *)gleft, nparts, IDX_DATATYPE, MPI_SUM, ctrl->comm);
  MPI_Allreduce((void *)lend, (void *)gend, nparts, IDX_DATATYPE, MPI_SUM, ctrl->comm);

  *nmoved = idxsum(nparts, gleft);
  *maxout = gleft[idxamax(nparts, gleft)];
  for (i=0; i<nparts; i++)
    lstart[i] = gend[i]+gleft[i]-gstart[i];
  *maxin = lstart[idxamax(nparts, lstart)];

  GKfree((void **)&lstart, (void **)&gstart, (void **)&lleft, (void **)&gleft, (void **)&lend, (void **)&gend, LTERM);
}

/*************************************************************************
* This function reads the element node array of a mesh
**************************************************************************/
idxtype *ReadMesh(char *filename, int *ne, int *nn, int *etype)
{
  int i, j, k, esize;
  idxtype *elmnts;
  FILE *fpin;

  if ((fpin = fopen(filename, "r")) == NULL) {
    printf("Failed to open file %s\n", filename);
    exit(0);
  }

  if (fscanf(fpin, "%d %d", ne, etype) != 2) {
    printf("Header line of input file does not contain two numbers.\n");
    exit(0);
  }

  switch (*etype) {
    case 1:
      esize = 3;
      break;
    case 2:
      esize = 4;
      break;
    case 3:
      esize = 8;
      break;
    case 4:
      esize = 4;
      break;
    default:
      errexit("Unknown mesh-element type: %d\n", *etype);
  }

  elmnts = idxmalloc(esize*(*ne), "ReadMesh: elmnts");

  for (j=esize*(*ne), i=0; i<j; i++) {
    if (fscanf(fpin, "%d", elmnts+i) != 1) {
      printf("Missing node number %d for element %d\n", i%esize+1, i/esize);
      exit(0);
    }
    elmnts[i]--;
  }

  fclose(fpin);

  *nn = elmnts[idxamax(j, elmnts)]+1;

  return elmnts;
}

/*************************************************************************
* This function computes the balance of the element partitioning
**************************************************************************/
float ComputeElementBalance(int ne, int nparts, idxtype *where)
{
  int i;
  idxtype *kpwgts;
  float balance;

  kpwgts = idxsmalloc(nparts, 0, "ComputeElementBalance: kpwgts");

  for (i=0; i<ne; i++)
    kpwgts[where[i]]++;

  balance = 1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts));

  free(kpwgts);

  return balance;

}

/*************************************************************************
* This function computes the subdomain graph
**************************************************************************/
void EliminateSubDomainEdges(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts)
{
  int i, ii, j, k, me, other, nvtxs, total, max, avg, totalout, nind, ncand, ncand2, target, target2, nadd;
  int min, move, cpwgt, tvwgt;
  idxtype *xadj, *adjncy, *vwgt, *adjwgt, *pwgts, *where, *maxpwgt, *pmat, *ndoms, *mypmat, *otherpmat, *ind;
  KeyValueType *cand, *cand2;

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

  where = graph->where;
  pwgts = graph->pwgts;  /* We assume that this is properly initialized */

  maxpwgt = idxwspacemalloc(ctrl, nparts);
  ndoms = idxwspacemalloc(ctrl, nparts);
  otherpmat = idxwspacemalloc(ctrl, nparts);
  ind = idxwspacemalloc(ctrl, nvtxs);
  pmat = ctrl->wspace.pmat;

  cand = (KeyValueType *)GKmalloc(nparts*sizeof(KeyValueType), "EliminateSubDomainEdges: cand");
  cand2 = (KeyValueType *)GKmalloc(nparts*sizeof(KeyValueType), "EliminateSubDomainEdges: cand");

  /* Compute the pmat matrix and ndoms */
  ComputeSubDomainGraph(graph, nparts, pmat, ndoms);


  /* Compute the maximum allowed weight for each domain */
  tvwgt = idxsum(nparts, pwgts);
  for (i=0; i<nparts; i++)
    maxpwgt[i] = 1.25*tpwgts[i]*tvwgt;


  /* Get into the loop eliminating subdomain connections */
  for (;;) {
    total = idxsum(nparts, ndoms);
    avg = total/nparts;
    max = ndoms[idxamax(nparts, ndoms)];

    /* printf("Adjacent Subdomain Stats: Total: %3d, Max: %3d, Avg: %3d [%5d]\n", total, max, avg, idxsum(nparts*nparts, pmat)); */

    if (max < 1.4*avg)
      break;

    me = idxamax(nparts, ndoms);
    mypmat = pmat + me*nparts;
    totalout = idxsum(nparts, mypmat);

    /*printf("Me: %d, TotalOut: %d,\n", me, totalout);*/

    /* Sort the connections according to their cut */
    for (ncand2=0, i=0; i<nparts; i++) {
      if (mypmat[i] > 0) {
        cand2[ncand2].key = mypmat[i];
        cand2[ncand2++].val = i;
      }
    }
    ikeysort(ncand2, cand2);

    move = 0;
    for (min=0; min<ncand2; min++) {
      if (cand2[min].key > totalout/(2*ndoms[me])) 
        break;

      other = cand2[min].val;

      /*printf("\tMinOut: %d to %d\n", mypmat[other], other);*/

      idxset(nparts, 0, otherpmat);

      /* Go and find the vertices in 'other' that are connected in 'me' */
      for (nind=0, i=0; i<nvtxs; i++) {
        if (where[i] == other) {
          for (j=xadj[i]; j<xadj[i+1]; j++) {
            if (where[adjncy[j]] == me) {
              ind[nind++] = i;
              break;
            }
          }
        }
      }

      /* Go and construct the otherpmat to see where these nind vertices are connected to */
      for (cpwgt=0, ii=0; ii<nind; ii++) {
        i = ind[ii];
        cpwgt += vwgt[i];

        for (j=xadj[i]; j<xadj[i+1]; j++) 
          otherpmat[where[adjncy[j]]] += adjwgt[j];
      }
      otherpmat[other] = 0;

      for (ncand=0, i=0; i<nparts; i++) {
        if (otherpmat[i] > 0) {
          cand[ncand].key = -otherpmat[i];
//.........这里部分代码省略.........

/*************************************************************************
* This function performs k-way refinement
**************************************************************************/
void Greedy_KWayEdgeBalanceMConn(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, nmoves; 
  int from, me, to, oldcut, vwgt, maxndoms, nadd;
  idxtype *xadj, *adjncy, *adjwgt;
  idxtype *where, *pwgts, *perm, *bndptr, *bndind, *minwgt, *maxwgt, *moved, *itpwgts;
  idxtype *phtable, *pmat, *pmatptr, *ndoms;
  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;
  
  pmat = ctrl->wspace.pmat;
  phtable = idxwspacemalloc(ctrl, nparts);
  ndoms = idxwspacemalloc(ctrl, nparts);

  ComputeSubDomainGraph(graph, nparts, pmat, ndoms);


  /* 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 [B]\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);

    /* Check to see if things are out of balance, given the tolerance */
    for (i=0; i<nparts; i++) {
      if (pwgts[i] > maxwgt[i])
        break;
    }
    if (i == nparts) /* Things are balanced. Return right away */
      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;
    }

    maxndoms = ndoms[idxamax(nparts, ndoms)];

    for (nmoves=0;;) {
      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;

      /* Determine the valid domains */
      for (j=0; j<myndegrees; j++) {
        to = myedegrees[j].pid;
//.........这里部分代码省略.........

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

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

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

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

  pmat = ctrl->wspace.pmat;
  phtable = idxwspacemalloc(ctrl, nparts);
  ndoms = idxwspacemalloc(ctrl, nparts);

  ComputeSubDomainGraph(graph, nparts, pmat, ndoms);

  /* 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);

  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);

    maxndoms = ndoms[idxamax(nparts, ndoms)];

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

    RandomPermute(nbnd, perm, 1);
    for (nmoves=iii=0; iii<graph->nbnd; iii++) {
      ii = perm[iii];
      if (ii >= nbnd)
        continue;
      i = bndind[ii];

      myrinfo = graph->rinfo+i;

      if (myrinfo->ed >= myrinfo->id) { /* Total ED is too high */
        from = where[i];
        vwgt = graph->vwgt[i];

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

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

        /* Determine the valid domains */
        for (j=0; j<myndegrees; j++) {
          to = myedegrees[j].pid;
          phtable[to] = 1;
          pmatptr = pmat + to*nparts;
          for (nadd=0, k=0; k<myndegrees; k++) {
            if (k == j)
              continue;

            l = myedegrees[k].pid;
            if (pmatptr[l] == 0) {
              if (ndoms[l] > maxndoms-1) {
                phtable[to] = 0;
                nadd = maxndoms;
                break;
              }
              nadd++;
            }
          }
          if (ndoms[to]+nadd > maxndoms)
            phtable[to] = 0;
//.........这里部分代码省略.........

/*************************************************************************
* 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 converts a mesh into a dual graph
**************************************************************************/
void ParMETIS_V3_Mesh2Dual(idxtype *elmdist, idxtype *eptr, idxtype *eind, 
                 int *numflag, int *ncommonnodes, idxtype **xadj, 
		 idxtype **adjncy, MPI_Comm *comm)
{
  int i, j, jj, k, kk, m;
  int npes, mype, pe, count, mask, pass;
  int nelms, lnns, my_nns, node;
  int firstelm, firstnode, lnode, nrecv, nsend;
  int *scounts, *rcounts, *sdispl, *rdispl;
  idxtype *nodedist, *nmap, *auxarray;
  idxtype *gnptr, *gnind, *nptr, *nind, *myxadj, *myadjncy = NULL;
  idxtype *sbuffer, *rbuffer, *htable;
  KeyValueType *nodelist, *recvbuffer;
  idxtype ind[200], wgt[200];
  int gmaxnode, gminnode;
  CtrlType ctrl;


  SetUpCtrl(&ctrl, -1, 0, *comm);

  npes = ctrl.npes;
  mype = ctrl.mype;

  nelms = elmdist[mype+1]-elmdist[mype];

  if (*numflag == 1) 
    ChangeNumberingMesh2(elmdist, eptr, eind, NULL, NULL, NULL, npes, mype, 1);

  mask = (1<<11)-1;

  /*****************************/
  /* Determine number of nodes */
  /*****************************/
  gminnode = GlobalSEMin(&ctrl, eind[idxamin(eptr[nelms], eind)]);
  for (i=0; i<eptr[nelms]; i++)
    eind[i] -= gminnode;

  gmaxnode = GlobalSEMax(&ctrl, eind[idxamax(eptr[nelms], eind)]);


  /**************************/
  /* Check for input errors */
  /**************************/
  ASSERTS(nelms > 0);

  /* construct node distribution array */
  nodedist = idxsmalloc(npes+1, 0, "nodedist");
  for (nodedist[0]=0, i=0,j=gmaxnode+1; i<npes; i++) {
    k = j/(npes-i);
    nodedist[i+1] = nodedist[i]+k;
    j -= k;
  }
  my_nns = nodedist[mype+1]-nodedist[mype];
  firstnode = nodedist[mype];

  nodelist = (KeyValueType *)GKmalloc(eptr[nelms]*sizeof(KeyValueType), "nodelist");
  auxarray = idxmalloc(eptr[nelms], "auxarray");
  htable   = idxsmalloc(amax(my_nns, mask+1), -1, "htable");
  scounts  = imalloc(4*npes+2, "scounts");
  rcounts  = scounts+npes;
  sdispl   = scounts+2*npes;
  rdispl   = scounts+3*npes+1;


  /*********************************************/
  /* first find a local numbering of the nodes */
  /*********************************************/
  for (i=0; i<nelms; i++) {
    for (j=eptr[i]; j<eptr[i+1]; j++) {
      nodelist[j].key = eind[j];
      nodelist[j].val = j;
      auxarray[j]     = i; /* remember the local element ID that uses this node */
    }
  }
  ikeysort(eptr[nelms], nodelist);

  for (count=1, i=1; i<eptr[nelms]; i++) {
    if (nodelist[i].key > nodelist[i-1].key)
      count++;
  }

  lnns = count;
  nmap = idxmalloc(lnns, "nmap");

  /* renumber the nodes of the elements array */
  count = 1;
  nmap[0] = nodelist[0].key;
  eind[nodelist[0].val] = 0;
  nodelist[0].val = auxarray[nodelist[0].val];  /* Store the local element ID */
  for (i=1; i<eptr[nelms]; i++) {
    if (nodelist[i].key > nodelist[i-1].key) {
      nmap[count] = nodelist[i].key;
      count++;
    }
    eind[nodelist[i].val] = count-1;
    nodelist[i].val = auxarray[nodelist[i].val];  /* Store the local element ID */
  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	part = idxmalloc(graph.nvtxs, "main: part");

	printf("------------------------------------------------\n");
	printf("Clustering....\n");
	fflush(stdout);
	starttimer(METISTmr);         //YK: main algorithm starts here!

	if ( noOfSingletons > 0 )
	{
		
		mlmcl(&(noSingletonGraph->nvtxs), noSingletonGraph->xadj, noSingletonGraph->adjncy,
			noSingletonGraph->vwgt,noSingletonGraph->adjwgt, &wgtflag, part, opt ); 
	}
	else	
	{
		mlmcl(&graph.nvtxs, graph.xadj, graph.adjncy,graph.vwgt,
			graph.adjwgt, &wgtflag, part, opt ); 
	}

	stoptimer(METISTmr); 

	printf("------------------------------------------------\n");
	if ( noOfSingletons > 0 )
	{
		npart=mapPartition(part,noSingletonGraph->nvtxs);
		ncut=ComputeNCut(noSingletonGraph, part,npart);
//		printf("In graph that does not include singletons,");
//		printf("No. of Clusters: %d, N-Cut value: %.2f\n",npart,ncut);


		idxtype *clusterSizes = histogram(part,
					graph.nvtxs-noOfSingletons, npart);

		int maxSize = clusterSizes[idxamax(npart, clusterSizes)];
		float avgClusterSize =
						(graph.nvtxs-noOfSingletons)*1.0/(npart);
		float balance =	(maxSize*1.0) /
				((graph.nvtxs-noOfSingletons)*1.0/npart);
		float stdDevn = stdDeviation(clusterSizes, npart);
		float avgNcut = ncut * 1.0/npart;
		float normStdDevn = stdDevn/avgClusterSize;

	// Warning: This computation only works if the singletons
	// have been placed in their own clusters. This works for
	// MLR-MCL, in other words, because it is guaranteed to
	// place singletons in their own clusters.
		printf("Output statistics for graph without singletons\n");
		printf("Clusters: %d N-Cut: %.3f", 
					npart, ncut);
		printf(" AvgN-Cut: %.3f Balance in cluster sizes: %.2f ",avgNcut,
					balance); 
		printf("Std_Deviation in cluster sizes: %.2f ", stdDevn);
		printf("Coefficient_of_Variation: %.2f\n", normStdDevn);

		free( clusterSizes );

		npart += noOfSingletons;
	//	ncut += noOfSingletons;
		printf("Output statistics for original graph\n");

		mapIndices(part, nodeMap, graph.nvtxs, npart-noOfSingletons);
	}
	else
	{
		npart=mapPartition(part,graph.nvtxs);
		ncut=ComputeNCut(&graph, part,npart);

/*************************************************************************
* 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;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    }

    MPI_Bcast(&ier, 1, MPI_INT, npes-1, comm);
    if (ier > 0){
      fclose(fpin);
      MPI_Finalize();
      exit(0);
    }

    line = (char *)GKmalloc(sizeof(char)*(MAXLINE+1), "line");

    fgets(line, MAXLINE, fpin);
    sscanf(line, "%d %d", &gnelms, &etype);

    /* Construct elmdist and send it to all the processors */
    elmdist[0] = 0;
    for (i=0,j=gnelms; i<npes; i++) {
      k = j/(npes-i);
      elmdist[i+1] = elmdist[i]+k;
      j -= k;
    }

    MPI_Bcast((void *)elmdist, npes+1, IDX_DATATYPE, npes-1, comm);
  }
  else {
    MPI_Bcast(&ier, 1, MPI_INT, npes-1, comm);
    if (ier > 0){
      MPI_Finalize();
      exit(0);
    }

    MPI_Bcast((void *)elmdist, npes+1, IDX_DATATYPE, npes-1, comm);
  }

  MPI_Bcast((void *)(&etype), 1, MPI_INT, npes-1, comm);

  gnelms = mesh->gnelms = elmdist[npes];
  nelms = mesh->nelms = elmdist[mype+1]-elmdist[mype];
  mesh->etype = etype;
  esize = esizes[etype];
  mgcnum = mgcnums[etype];

  elements = mesh->elements = idxmalloc(nelms*esize, "ParallelReadMesh: elements");

  if (mype == npes-1) {
    maxnelms = 0;
    for (i=0; i<npes; i++) {
      maxnelms = (maxnelms > elmdist[i+1]-elmdist[i]) ?
      maxnelms : elmdist[i+1]-elmdist[i];
    }

    your_elements = idxmalloc(maxnelms*esize, "your_elements");

    for (pe=0; pe<npes; pe++) {
      your_nelms = elmdist[pe+1]-elmdist[pe];
      for (i=0; i<your_nelms; i++) {

        fgets(line, MAXLINE, fpin);
        oldstr = line;
        newstr = NULL;

        /*************************************/
        /* could get element weigts here too */
        /*************************************/

        for (j=0; j<esize; j++) {
          your_elements[i*esize+j] = (int)strtol(oldstr, &newstr, 10);
          oldstr = newstr;
        }
      }

      if (pe < npes-1) {
        MPI_Send((void *)your_elements, your_nelms*esize, IDX_DATATYPE, pe, 0, comm);
      }
      else {
        for (i=0; i<your_nelms*esize; i++)
          elements[i] = your_elements[i];
      }
    }
    fclose(fpin);
    free(your_elements);
  }
  else {
    MPI_Recv((void *)elements, nelms*esize, IDX_DATATYPE, npes-1, 0, comm, &stat);
  }

  /*********************************/
  /* now check for number of nodes */
  /*********************************/
  minnode = elements[idxamin(nelms*esize, elements)];
  MPI_Allreduce((void *)&minnode, (void *)&gminnode, 1, MPI_INT, MPI_MIN, comm);
  for (i=0; i<nelms*esize; i++)
    elements[i] -= gminnode;

  maxnode = elements[idxamax(nelms*esize, elements)];
  MPI_Allreduce((void *)&maxnode, (void *)&gmaxnode, 1, MPI_INT, MPI_MAX, comm);
  mesh->gnns = gmaxnode+1;

  if (mype==0) printf("Nelements: %d, Nnodes: %d, EType: %d\n", gnelms, mesh->gnns, etype);
}

/*************************************************************************
* This function computes cuts and balance information
**************************************************************************/
void ComputePartitionInfo(GraphType *graph, int nparts, idxtype *where)
{
  int i, j, /*k,*/ nvtxs, ncon, mustfree=0;
  idxtype *xadj, *adjncy, *vwgt, *adjwgt, *kpwgts, *tmpptr;
  idxtype *padjncy, *padjwgt, *padjcut;

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

  if (vwgt == NULL) {
    vwgt = graph->vwgt = idxsmalloc(nvtxs, 1, "vwgt");
    mustfree = 1;
  }
  if (adjwgt == NULL) {
    adjwgt = graph->adjwgt = idxsmalloc(xadj[nvtxs], 1, "adjwgt");
    mustfree += 2;
  }

  printf("%d-way Cut: %5d, Vol: %5d, ", nparts, ComputeCut(graph, where), ComputeVolume(graph, where));

  /* Compute balance information */
  kpwgts = idxsmalloc(ncon*nparts, 0, "ComputePartitionInfo: kpwgts");

  for (i=0; i<nvtxs; i++) {
    for (j=0; j<ncon; j++) 
      kpwgts[where[i]*ncon+j] += vwgt[i*ncon+j];
  }

  if (ncon == 1) {
    printf("\tBalance: %5.3f out of %5.3f\n", 
            1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts)),
            1.0*nparts*vwgt[idxamax(nvtxs, vwgt)]/(1.0*idxsum(nparts, kpwgts)));
  }
  else {
    printf("\tBalance:");
    for (j=0; j<ncon; j++) 
      printf(" (%5.3f out of %5.3f)", 
            1.0*nparts*kpwgts[ncon*idxamax_strd(nparts, kpwgts+j, ncon)+j]/(1.0*idxsum_strd(nparts, kpwgts+j, ncon)),
            1.0*nparts*vwgt[ncon*idxamax_strd(nvtxs, vwgt+j, ncon)+j]/(1.0*idxsum_strd(nparts, kpwgts+j, ncon)));
    printf("\n");
  }


  /* Compute p-adjncy information */
  padjncy = idxsmalloc(nparts*nparts, 0, "ComputePartitionInfo: padjncy");
  padjwgt = idxsmalloc(nparts*nparts, 0, "ComputePartitionInfo: padjwgt");
  padjcut = idxsmalloc(nparts*nparts, 0, "ComputePartitionInfo: padjwgt");

  idxset(nparts, 0, kpwgts);
  for (i=0; i<nvtxs; i++) {
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      if (where[i] != where[adjncy[j]]) {
        padjncy[where[i]*nparts+where[adjncy[j]]] = 1;
        padjcut[where[i]*nparts+where[adjncy[j]]] += adjwgt[j];
        if (kpwgts[where[adjncy[j]]] == 0) {
          padjwgt[where[i]*nparts+where[adjncy[j]]]++;
          kpwgts[where[adjncy[j]]] = 1;
        }
      }
    }
    for (j=xadj[i]; j<xadj[i+1]; j++) 
      kpwgts[where[adjncy[j]]] = 0;
  }

  for (i=0; i<nparts; i++)
    kpwgts[i] = idxsum(nparts, padjncy+i*nparts);
  printf("Min/Max/Avg/Bal # of adjacent     subdomains: %5d %5d %5.2f %7.3f\n",
    kpwgts[idxamin(nparts, kpwgts)], kpwgts[idxamax(nparts, kpwgts)], 
    1.0*idxsum(nparts, kpwgts)/(1.0*nparts), 
    1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts)));

  for (i=0; i<nparts; i++)
    kpwgts[i] = idxsum(nparts, padjcut+i*nparts);
  printf("Min/Max/Avg/Bal # of adjacent subdomain cuts: %5d %5d %5d %7.3f\n",
    kpwgts[idxamin(nparts, kpwgts)], kpwgts[idxamax(nparts, kpwgts)], idxsum(nparts, kpwgts)/nparts, 
    1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts)));

  for (i=0; i<nparts; i++)
    kpwgts[i] = idxsum(nparts, padjwgt+i*nparts);
  printf("Min/Max/Avg/Bal/Frac # of interface    nodes: %5d %5d %5d %7.3f %7.3f\n",
    kpwgts[idxamin(nparts, kpwgts)], kpwgts[idxamax(nparts, kpwgts)], idxsum(nparts, kpwgts)/nparts, 
    1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts)), 1.0*idxsum(nparts, kpwgts)/(1.0*nvtxs));

  tmpptr = graph->where;
  graph->where = where;
  for (i=0; i<nparts; i++)
    IsConnectedSubdomain(NULL, graph, i, 1);
  graph->where = tmpptr;

  if (mustfree == 1 || mustfree == 3) {
    free(vwgt);
    graph->vwgt = NULL;
  }
//.........这里部分代码省略.........

/*************************************************************************
* This function performs an edge-based FM refinement
**************************************************************************/
void FM_2WayEdgeRefine(CtrlType *ctrl, GraphType *graph, int *tpwgts, int npasses)
{
  int i, ii, j, k, kwgt, nvtxs, nbnd, nswaps, from, to, pass, me, limit, tmp;
  idxtype *xadj, *vwgt, *adjncy, *adjwgt, *where, *id, *ed, *bndptr, *bndind, *pwgts;
  idxtype *moved, *swaps, *perm;
  PQueueType parts[2];
  int higain, oldgain, mincut, mindiff, origdiff, initcut, newcut, mincutorder, avgvwgt;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  vwgt = graph->vwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  where = graph->where;
  id = graph->id;
  ed = graph->ed;
  pwgts = graph->pwgts;
  bndptr = graph->bndptr;
  bndind = graph->bndind;

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

  limit = (int) amin(amax(0.01*nvtxs, 15), 100);
  avgvwgt = amin((pwgts[0]+pwgts[1])/20, 2*(pwgts[0]+pwgts[1])/nvtxs);

  tmp = graph->adjwgtsum[idxamax(nvtxs, graph->adjwgtsum)];
  PQueueInit(ctrl, &parts[0], nvtxs, tmp);
  PQueueInit(ctrl, &parts[1], nvtxs, tmp);

  IFSET(ctrl->dbglvl, DBG_REFINE, 
     printf("Partitions: [%6d %6d] T[%6d %6d], Nv-Nb[%6d %6d]. ICut: %6d\n",
             pwgts[0], pwgts[1], tpwgts[0], tpwgts[1], graph->nvtxs, graph->nbnd, graph->mincut));

  origdiff = abs(tpwgts[0]-pwgts[0]);
  idxset(nvtxs, -1, moved);
  for (pass=0; pass<npasses; pass++) { /* Do a number of passes */
    PQueueReset(&parts[0]);
    PQueueReset(&parts[1]);

    mincutorder = -1;
    newcut = mincut = initcut = graph->mincut;
    mindiff = abs(tpwgts[0]-pwgts[0]);

    ASSERT(ComputeCut(graph, where) == graph->mincut);
    ASSERT(CheckBnd(graph));

    /* Insert boundary nodes in the priority queues */
    nbnd = graph->nbnd;
    RandomPermute(nbnd, perm, 1);
    for (ii=0; ii<nbnd; ii++) {
      i = perm[ii];
      ASSERT(ed[bndind[i]] > 0 || id[bndind[i]] == 0);
      ASSERT(bndptr[bndind[i]] != -1);
      PQueueInsert(&parts[where[bndind[i]]], bndind[i], ed[bndind[i]]-id[bndind[i]]);
    }

    for (nswaps=0; nswaps<nvtxs; nswaps++) {
      from = (tpwgts[0]-pwgts[0] < tpwgts[1]-pwgts[1] ? 0 : 1);
      to = (from+1)%2;

      if ((higain = PQueueGetMax(&parts[from])) == -1)
        break;
      ASSERT(bndptr[higain] != -1);

      newcut -= (ed[higain]-id[higain]);
      INC_DEC(pwgts[to], pwgts[from], vwgt[higain]);

      if ((newcut < mincut && abs(tpwgts[0]-pwgts[0]) <= origdiff+avgvwgt) || 
          (newcut == mincut && abs(tpwgts[0]-pwgts[0]) < mindiff)) {
        mincut = newcut;
        mindiff = abs(tpwgts[0]-pwgts[0]);
        mincutorder = nswaps;
      }
      else if (nswaps-mincutorder > limit) { /* We hit the limit, undo last move */
        newcut += (ed[higain]-id[higain]);
        INC_DEC(pwgts[from], pwgts[to], vwgt[higain]);
        break;
      }

      where[higain] = to;
      moved[higain] = nswaps;
      swaps[nswaps] = higain;

      IFSET(ctrl->dbglvl, DBG_MOVEINFO, 
        printf("Moved %6d from %d. [%3d %3d] %5d [%4d %4d]\n", higain, from, ed[higain]-id[higain], vwgt[higain], newcut, pwgts[0], pwgts[1]));

      /**************************************************************
      * Update the id[i]/ed[i] values of the affected nodes
      ***************************************************************/
      SWAP(id[higain], ed[higain], tmp);
      if (ed[higain] == 0 && xadj[higain] < xadj[higain+1]) 
        BNDDelete(nbnd, bndind,  bndptr, higain);

      for (j=xadj[higain]; j<xadj[higain+1]; j++) {
        k = adjncy[j];
//.........这里部分代码省略.........

/*************************************************************************
* This function balances two partitions by moving boundary nodes
* from the domain that is overweight to the one that is underweight.
**************************************************************************/
void Bnd2WayBalance(CtrlType *ctrl, GraphType *graph, int *tpwgts)
{
  int i, ii, j, k, kwgt, nvtxs, nbnd, nswaps, from, to, pass, me, tmp;
  idxtype *xadj, *vwgt, *adjncy, *adjwgt, *where, *id, *ed, *bndptr, *bndind, *pwgts;
  idxtype *moved, *perm;
  PQueueType parts;
  int higain, oldgain, mincut, mindiff;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  vwgt = graph->vwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  where = graph->where;
  id = graph->id;
  ed = graph->ed;
  pwgts = graph->pwgts;
  bndptr = graph->bndptr;
  bndind = graph->bndind;

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

  /* Determine from which domain you will be moving data */
  mindiff = abs(tpwgts[0]-pwgts[0]);
  from = (pwgts[0] < tpwgts[0] ? 1 : 0);
  to = (from+1)%2;

  IFSET(ctrl->dbglvl, DBG_REFINE, 
     printf("Partitions: [%6d %6d] T[%6d %6d], Nv-Nb[%6d %6d]. ICut: %6d [B]\n",
             pwgts[0], pwgts[1], tpwgts[0], tpwgts[1], graph->nvtxs, graph->nbnd, graph->mincut));

  tmp = graph->adjwgtsum[idxamax(nvtxs, graph->adjwgtsum)];
  PQueueInit(ctrl, &parts, nvtxs, tmp);

  idxset(nvtxs, -1, moved);

  ASSERT(ComputeCut(graph, where) == graph->mincut);
  ASSERT(CheckBnd(graph));

  /* Insert the boundary nodes of the proper partition whose size is OK in the priority queue */
  nbnd = graph->nbnd;
  RandomPermute(nbnd, perm, 1);
  for (ii=0; ii<nbnd; ii++) {
    i = perm[ii];
    ASSERT(ed[bndind[i]] > 0 || id[bndind[i]] == 0);
    ASSERT(bndptr[bndind[i]] != -1);
    if (where[bndind[i]] == from && vwgt[bndind[i]] <= mindiff)
      PQueueInsert(&parts, bndind[i], ed[bndind[i]]-id[bndind[i]]);
  }

  mincut = graph->mincut;
  for (nswaps=0; nswaps<nvtxs; nswaps++) {
    if ((higain = PQueueGetMax(&parts)) == -1)