creating branch for adding support for CUDA 3.x and Fermi

[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 6829]

1 files changed, 1778 insertions, 0 deletions

diff --git a/benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/METISLib/kwayvolfm.c b/benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/METISLib/kwayvolfm.c
new file mode 100644
index 0000000..4e1112d
--- /dev/null
+++ b/benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/METISLib/kwayvolfm.c
@@ -0,0 +1,1778 @@
+/*
+ * kwayvolfm.c
+ *
+ * This file contains code that implements the multilevel k-way refinement
+ *
+ * Started 7/8/98
+ * George
+ *
+ * $Id: kwayvolfm.c,v 1.2 2003/07/31 06:14:01 karypis Exp $
+ *
+ */
+
+#include <metis.h>
+
+
+/*************************************************************************
+* This function performs k-way refinement
+**************************************************************************/
+void Random_KWayVolRefine(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts, 
+                          float ubfactor, int npasses, int ffactor)
+{
+  int i, ii, iii, j, jj, k, kk, l, u, pass, nvtxs, nmoves, tvwgt, myndegrees, xgain; 
+  int from, me, to, oldcut, oldvol, vwgt;
+  idxtype *xadj, *adjncy, *adjwgt;
+  idxtype *where, *pwgts, *perm, *bndptr, *bndind, *minwgt, *maxwgt, *itpwgts, *updind, *marker, *phtable;
+  VEDegreeType *myedegrees;
+  VRInfoType *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;
+  
+  /* 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));
+
+  updind = idxmalloc(nvtxs, "Random_KWayVolRefine: updind");
+  marker = idxsmalloc(nvtxs, 0, "Random_KWayVolRefine: marker");
+  phtable = idxsmalloc(nparts, -1, "Random_KWayVolRefine: phtable");
+
+  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("VolPart: [%5d %5d]-[%5d %5d], Balance: %3.2f, Nv-Nb[%5d %5d]. Cut: %5d, Vol: %5d\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, graph->minvol));
+
+  for (pass=0; pass<npasses; pass++) {
+    ASSERT(ComputeCut(graph, where) == graph->mincut);
+
+    oldcut = graph->mincut;
+    oldvol = graph->minvol;
+
+    RandomPermute(graph->nbnd, perm, 1);
+    for (nmoves=iii=0; iii<graph->nbnd; iii++) {
+      ii = perm[iii];
+      if (ii >= graph->nbnd)
+        continue;
+      i = bndind[ii];
+      myrinfo = graph->vrinfo+i;
+
+      if (myrinfo->gv >= 0) { /* Total volume gain is too high */
+        from = where[i];
+        vwgt = graph->vwgt[i];
+
+        if (myrinfo->id > 0 && pwgts[from]-vwgt < minwgt[from]) 
+          continue;   /* This cannot be moved! */
+
+        xgain = (myrinfo->id == 0 && myrinfo->ed > 0 ? graph->vsize[i] : 0);
+
+        myedegrees = myrinfo->edegrees;
+        myndegrees = myrinfo->ndegrees;
+
+        for (k=0; k<myndegrees; k++) {
+          to = myedegrees[k].pid;
+          if (pwgts[to]+vwgt <= maxwgt[to]+ffactor*myedegrees[k].gv && xgain+myedegrees[k].gv >= 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 (pwgts[to]+vwgt > maxwgt[to])
+            continue;
+          if (myedegrees[j].gv > myedegrees[k].gv ||
+              (myedegrees[j].gv == myedegrees[k].gv && myedegrees[j].ed > myedegrees[k].ed) ||
+              (myedegrees[j].gv == myedegrees[k].gv && 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;
+        if (xgain+myedegrees[k].gv > 0 || myedegrees[k].ed-myrinfo->id > 0)
+          j = 1;
+        else if (myedegrees[k].ed-myrinfo->id == 0) {
+          if ((iii&5) == 0 || pwgts[from] >= maxwgt[from] || itpwgts[from]*(pwgts[to]+vwgt) < itpwgts[to]*pwgts[from])
+            j = 1;
+        }
+        if (j == 0)
+          continue;
+          
+        /*=====================================================================
+        * If we got here, we can now move the vertex from 'from' to 'to' 
+        *======================================================================*/
+        INC_DEC(pwgts[to], pwgts[from], vwgt);
+        graph->mincut -= myedegrees[k].ed-myrinfo->id;
+        graph->minvol -= (xgain+myedegrees[k].gv);
+        where[i] = to;
+
+        IFSET(ctrl->dbglvl, DBG_MOVEINFO, printf("\t\tMoving %6d from %3d to %3d. Gain: [%4d %4d]. Cut: %6d, Vol: %6d\n", 
+              i, from, to, xgain+myedegrees[k].gv, myedegrees[k].ed-myrinfo->id, graph->mincut, graph->minvol));
+
+        KWayVolUpdate(ctrl, graph, i, from, to, marker, phtable, updind);
+
+        nmoves++;
+
+        /* CheckVolKWayPartitionParams(ctrl, graph, nparts); */
+      }
+    }
+
+    IFSET(ctrl->dbglvl, DBG_REFINE,
+       printf("\t[%6d %6d], Balance: %5.3f, Nb: %6d. Nmoves: %5d, Cut: %6d, Vol: %6d\n",
+               pwgts[idxamin(nparts, pwgts)], pwgts[idxamax(nparts, pwgts)],
+               1.0*nparts*pwgts[idxamax(nparts, pwgts)]/tvwgt, graph->nbnd, nmoves, graph->mincut, 
+               graph->minvol));
+
+    if (graph->minvol == oldvol && graph->mincut == oldcut)
+      break;
+  }
+
+  GKfree(&marker, &updind, &phtable, LTERM);
+
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nvtxs);
+}
+
+
+/*************************************************************************
+* This function performs k-way refinement
+**************************************************************************/
+void Random_KWayVolRefineMConn(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts, 
+            float ubfactor, int npasses, int ffactor)
+{
+  int i, ii, iii, j, jj, k, kk, l, u, pass, nvtxs, nmoves, tvwgt, myndegrees, xgain; 
+  int from, me, to, oldcut, oldvol, vwgt, nadd, maxndoms;
+  idxtype *xadj, *adjncy, *adjwgt;
+  idxtype *where, *pwgts, *perm, *bndptr, *bndind, *minwgt, *maxwgt, *itpwgts, *updind, *marker, *phtable;
+  idxtype *pmat, *pmatptr, *ndoms;
+  VEDegreeType *myedegrees;
+  VRInfoType *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;
+  
+  /* 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));
+
+  updind = idxmalloc(nvtxs, "Random_KWayVolRefine: updind");
+  marker = idxsmalloc(nvtxs, 0, "Random_KWayVolRefine: marker");
+  phtable = idxsmalloc(nparts, -1, "Random_KWayVolRefine: phtable");
+
+  pmat = ctrl->wspace.pmat;
+  ndoms = idxwspacemalloc(ctrl, nparts);
+
+  ComputeVolSubDomainGraph(graph, nparts, pmat, ndoms);
+
+  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("VolPart: [%5d %5d]-[%5d %5d], Balance: %3.2f, Nv-Nb[%5d %5d]. Cut: %5d, Vol: %5d\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, graph->minvol));
+
+  for (pass=0; pass<npasses; pass++) {
+    ASSERT(ComputeCut(graph, where) == graph->mincut);
+
+    maxndoms = ndoms[idxamax(nparts, ndoms)];
+
+    oldcut = graph->mincut;
+    oldvol = graph->minvol;
+
+    RandomPermute(graph->nbnd, perm, 1);
+    for (nmoves=iii=0; iii<graph->nbnd; iii++) {
+      ii = perm[iii];
+      if (ii >= graph->nbnd)
+        continue;
+      i = bndind[ii];
+      myrinfo = graph->vrinfo+i;
+
+      if (myrinfo->gv >= 0) { /* Total volume gain is too high */
+        from = where[i];
+        vwgt = graph->vwgt[i];
+
+        if (myrinfo->id > 0 && pwgts[from]-vwgt < minwgt[from]) 
+          continue;   /* This cannot be moved! */
+
+        xgain = (myrinfo->id == 0 && myrinfo->ed > 0 ? graph->vsize[i] : 0);
+
+        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;
+          if (nadd == 0)
+            phtable[to] = 2;
+        }
+
+        for (k=0; k<myndegrees; k++) {
+          to = myedegrees[k].pid;
+          if (!phtable[to])
+            continue;
+          if (pwgts[to]+vwgt <= maxwgt[to]+ffactor*myedegrees[k].gv && xgain+myedegrees[k].gv >= 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 (!phtable[to] || pwgts[to]+vwgt > maxwgt[to])
+            continue;
+          if (myedegrees[j].gv > myedegrees[k].gv ||
+              (myedegrees[j].gv == myedegrees[k].gv && myedegrees[j].ed > myedegrees[k].ed) ||
+              (myedegrees[j].gv == myedegrees[k].gv && 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;
+        if (xgain+myedegrees[k].gv > 0 || myedegrees[k].ed-myrinfo->id > 0)
+          j = 1;
+        else if (myedegrees[k].ed-myrinfo->id == 0) {
+          if ((iii&5) == 0 || phtable[myedegrees[k].pid] == 2 || pwgts[from] >= maxwgt[from] || itpwgts[from]*(pwgts[to]+vwgt) < itpwgts[to]*pwgts[from])
+            j = 1;
+        }
+
+        if (j == 0)
+          continue;
+
+        for (j=0; j<myndegrees; j++) 
+          phtable[myedegrees[j].pid] = -1;
+
+          
+        /*=====================================================================
+        * If we got here, we can now move the vertex from 'from' to 'to' 
+        *======================================================================*/
+        INC_DEC(pwgts[to], pwgts[from], vwgt);
+        graph->mincut -= myedegrees[k].ed-myrinfo->id;
+        graph->minvol -= (xgain+myedegrees[k].gv);
+        where[i] = to;
+
+        IFSET(ctrl->dbglvl, DBG_MOVEINFO, printf("\t\tMoving %6d from %3d to %3d. Gain: [%4d %4d]. Cut: %6d, Vol: %6d\n", 
+              i, from, to, xgain+myedegrees[k].gv, myedegrees[k].ed-myrinfo->id, graph->mincut, graph->minvol));
+
+        /* Update pmat to reflect the move of 'i' */
+        pmat[from*nparts+to] += (myrinfo->id-myedegrees[k].ed);
+        pmat[to*nparts+from] += (myrinfo->id-myedegrees[k].ed);
+        if (pmat[from*nparts+to] == 0) {
+          ndoms[from]--;
+          if (ndoms[from]+1 == maxndoms)
+            maxndoms = ndoms[idxamax(nparts, ndoms)];
+        }
+        if (pmat[to*nparts+from] == 0) {
+          ndoms[to]--;
+          if (ndoms[to]+1 == maxndoms)
+            maxndoms = ndoms[idxamax(nparts, ndoms)];
+        }
+
+        for (j=xadj[i]; j<xadj[i+1]; j++) {
+          ii = adjncy[j];
+          me = where[ii];
+
+          /* Update pmat to reflect the move of 'i' for domains other than 'from' and 'to' */
+          if (me != from && me != to) {
+            pmat[me*nparts+from] -= adjwgt[j];
+            pmat[from*nparts+me] -= adjwgt[j];
+            if (pmat[me*nparts+from] == 0) {
+              ndoms[me]--;
+              if (ndoms[me]+1 == maxndoms)
+                maxndoms = ndoms[idxamax(nparts, ndoms)];
+            }
+            if (pmat[from*nparts+me] == 0) {
+              ndoms[from]--;
+              if (ndoms[from]+1 == maxndoms)
+                maxndoms = ndoms[idxamax(nparts, ndoms)];
+            }
+
+            if (pmat[me*nparts+to] == 0) {
+              ndoms[me]++;
+              if (ndoms[me] > maxndoms) {
+                IFSET(ctrl->dbglvl, DBG_REFINE, printf("You just increased the maxndoms: %d %d\n", ndoms[me], maxndoms));
+                maxndoms = ndoms[me];
+              }
+            }
+            if (pmat[to*nparts+me] == 0) {
+              ndoms[to]++;
+              if (ndoms[to] > maxndoms) {
+                IFSET(ctrl->dbglvl, DBG_REFINE, printf("You just increased the maxndoms: %d %d\n", ndoms[to], maxndoms));
+                maxndoms = ndoms[to];
+              }
+            }
+            pmat[me*nparts+to] += adjwgt[j];
+            pmat[to*nparts+me] += adjwgt[j];
+          }
+        }
+
+        KWayVolUpdate(ctrl, graph, i, from, to, marker, phtable, updind);
+
+        nmoves++;
+
+        /* CheckVolKWayPartitionParams(ctrl, graph, nparts); */
+      }
+    }
+
+    IFSET(ctrl->dbglvl, DBG_REFINE,
+       printf("\t[%6d %6d], Balance: %5.3f, Nb: %6d. Nmoves: %5d, Cut: %6d, Vol: %6d\n",
+               pwgts[idxamin(nparts, pwgts)], pwgts[idxamax(nparts, pwgts)],
+               1.0*nparts*pwgts[idxamax(nparts, pwgts)]/tvwgt, graph->nbnd, nmoves, graph->mincut, 
+               graph->minvol));
+
+    if (graph->minvol == oldvol && graph->mincut == oldcut)
+      break;
+  }
+
+  GKfree(&marker, &updind, &phtable, LTERM);
+
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nvtxs);
+}
+
+
+
+
+/*************************************************************************
+* This function performs k-way refinement
+**************************************************************************/
+void Greedy_KWayVolBalance(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts, 
+                           float ubfactor, int npasses)
+{
+  int i, ii, iii, j, jj, k, kk, l, u, pass, nvtxs, nmoves, tvwgt, myndegrees, xgain; 
+  int from, me, to, vwgt, gain;
+  idxtype *xadj, *adjncy, *adjwgt;
+  idxtype *where, *pwgts, *perm, *moved, *bndptr, *bndind, *minwgt, *maxwgt, *itpwgts, *updind, *marker, *phtable;
+  VEDegreeType *myedegrees;
+  VRInfoType *myrinfo; 
+  PQueueType queue;
+
+  nvtxs = graph->nvtxs;
+  xadj = graph->xadj;
+  adjncy = graph->adjncy;
+  adjwgt = graph->adjwgt;
+
+  bndptr = graph->bndptr;
+  bndind = graph->bndind;
+
+  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));
+
+  updind = idxmalloc(nvtxs, "Random_KWayVolRefine: updind");
+  marker = idxsmalloc(nvtxs, 0, "Random_KWayVolRefine: marker");
+  phtable = idxsmalloc(nparts, -1, "Random_KWayVolRefine: phtable");
+
+  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("VolPart: [%5d %5d]-[%5d %5d], Balance: %3.2f, Nv-Nb[%5d %5d]. Cut: %5d, Vol: %5d [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, graph->minvol));
+
+
+  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);
+
+    RandomPermute(graph->nbnd, perm, 1);
+    for (ii=0; ii<graph->nbnd; ii++) {
+      i = bndind[perm[ii]];
+      PQueueInsert(&queue, i, graph->vrinfo[i].gv);
+      moved[i] = 2;
+    }
+
+    for (nmoves=0;;) {
+      if ((i = PQueueGetMax(&queue)) == -1) 
+        break;
+      moved[i] = 1;
+
+      myrinfo = graph->vrinfo+i;
+      from = where[i];
+      vwgt = graph->vwgt[i];
+
+      if (pwgts[from]-vwgt < minwgt[from]) 
+        continue;   /* This cannot be moved! */
+
+      xgain = (myrinfo->id == 0 && myrinfo->ed > 0 ? graph->vsize[i] : 0);
+
+      myedegrees = myrinfo->edegrees;
+      myndegrees = myrinfo->ndegrees;
+
+      for (k=0; k<myndegrees; k++) {
+        to = myedegrees[k].pid;
+        if (pwgts[to]+vwgt <= maxwgt[to] || 
+            itpwgts[from]*(pwgts[to]+vwgt) <= itpwgts[to]*pwgts[from]) 
+          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 (itpwgts[myedegrees[k].pid]*pwgts[to] < itpwgts[to]*pwgts[myedegrees[k].pid]) 
+          k = j;
+      }
+
+      to = myedegrees[k].pid;
+
+      if (pwgts[from] < maxwgt[from] && pwgts[to] > minwgt[to] && 
+          (xgain+myedegrees[k].gv < 0 || 
+           (xgain+myedegrees[k].gv == 0 &&  myedegrees[k].ed-myrinfo->id < 0))
+         )
+        continue;
+  
+
+      /*=====================================================================
+      * If we got here, we can now move the vertex from 'from' to 'to' 
+      *======================================================================*/
+      INC_DEC(pwgts[to], pwgts[from], vwgt);
+      graph->mincut -= myedegrees[k].ed-myrinfo->id;
+      graph->minvol -= (xgain+myedegrees[k].gv);
+      where[i] = to;
+
+      IFSET(ctrl->dbglvl, DBG_MOVEINFO, printf("\t\tMoving %6d from %3d to %3d. Gain: [%4d %4d]. Cut: %6d, Vol: %6d\n", 
+            i, from, to, xgain+myedegrees[k].gv, myedegrees[k].ed-myrinfo->id, graph->mincut, graph->minvol));
+
+      KWayVolUpdate(ctrl, graph, i, from, to, marker, phtable, updind);
+
+      nmoves++;
+
+      /*CheckVolKWayPartitionParams(ctrl, graph, nparts); */
+    }
+
+    IFSET(ctrl->dbglvl, DBG_REFINE,
+       printf("\t[%6d %6d], Balance: %5.3f, Nb: %6d. Nmoves: %5d, Cut: %6d, Vol: %6d\n",
+               pwgts[idxamin(nparts, pwgts)], pwgts[idxamax(nparts, pwgts)],
+               1.0*nparts*pwgts[idxamax(nparts, pwgts)]/tvwgt, graph->nbnd, nmoves, graph->mincut, 
+               graph->minvol));
+
+  }
+
+  GKfree(&marker, &updind, &phtable, LTERM);
+
+  PQueueFree(ctrl, &queue);
+
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nvtxs);
+  idxwspacefree(ctrl, nvtxs);
+}
+
+
+
+/*************************************************************************
+* This function performs k-way refinement
+**************************************************************************/
+void Greedy_KWayVolBalanceMConn(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts, 
+                                float ubfactor, int npasses)
+{
+  int i, ii, iii, j, jj, k, kk, l, u, pass, nvtxs, nmoves, tvwgt, myndegrees, xgain; 
+  int from, me, to, vwgt, gain, maxndoms, nadd;
+  idxtype *xadj, *adjncy, *adjwgt;
+  idxtype *where, *pwgts, *perm, *moved, *bndptr, *bndind, *minwgt, *maxwgt, *itpwgts, *updind, *marker, *phtable;
+  idxtype *pmat, *pmatptr, *ndoms;
+  VEDegreeType *myedegrees;
+  VRInfoType *myrinfo; 
+  PQueueType queue;
+
+  nvtxs = graph->nvtxs;
+  xadj = graph->xadj;
+  adjncy = graph->adjncy;
+  adjwgt = graph->adjwgt;
+
+  bndptr = graph->bndptr;
+  bndind = graph->bndind;
+
+  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));
+
+  updind = idxmalloc(nvtxs, "Random_KWayVolRefine: updind");
+  marker = idxsmalloc(nvtxs, 0, "Random_KWayVolRefine: marker");
+  phtable = idxsmalloc(nparts, -1, "Random_KWayVolRefine: phtable");
+
+  pmat = ctrl->wspace.pmat;
+  ndoms = idxwspacemalloc(ctrl, nparts);
+
+  ComputeVolSubDomainGraph(graph, nparts, pmat, ndoms);
+
+  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("VolPart: [%5d %5d]-[%5d %5d], Balance: %3.2f, Nv-Nb[%5d %5d]. Cut: %5d, Vol: %5d [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, graph->minvol));
+
+
+  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);
+
+    RandomPermute(graph->nbnd, perm, 1);
+    for (ii=0; ii<graph->nbnd; ii++) {
+      i = bndind[perm[ii]];
+      PQueueInsert(&queue, i, graph->vrinfo[i].gv);
+      moved[i] = 2;
+    }
+
+    maxndoms = ndoms[idxamax(nparts, ndoms)];
+
+    for (nmoves=0;;) {
+      if ((i = PQueueGetMax(&queue)) == -1) 
+        break;
+      moved[i] = 1;
+
+      myrinfo = graph->vrinfo+i;
+      from = where[i];
+      vwgt = graph->vwgt[i];
+
+      if (pwgts[from]-vwgt < minwgt[from]) 
+        continue;   /* This cannot be moved! */
+
+      xgain = (myrinfo->id == 0 && myrinfo->ed > 0 ? graph->vsize[i] : 0);
+
+      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;
+      }
+
+      for (k=0; k<myndegrees; k++) {
+        to = myedegrees[k].pid;
+        if (!phtable[to])
+          continue;
+        if (pwgts[to]+vwgt <= maxwgt[to] || 
+            itpwgts[from]*(pwgts[to]+vwgt) <= itpwgts[to]*pwgts[from]) 
+          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 (!phtable[to])
+          continue;
+        if (itpwgts[myedegrees[k].pid]*pwgts[to] < itpwgts[to]*pwgts[myedegrees[k].pid]) 
+          k = j;
+      }
+
+      to = myedegrees[k].pid;
+
+      for (j=0; j<myndegrees; j++) 
+        phtable[myedegrees[j].pid] = -1;
+
+      if (pwgts[from] < maxwgt[from] && pwgts[to] > minwgt[to] && 
+          (xgain+myedegrees[k].gv < 0 || 
+           (xgain+myedegrees[k].gv == 0 &&  myedegrees[k].ed-myrinfo->id < 0))
+         )
+        continue;
+  
+
+      /*=====================================================================
+      * If we got here, we can now move the vertex from 'from' to 'to' 
+      *======================================================================*/
+      INC_DEC(pwgts[to], pwgts[from], vwgt);
+      graph->mincut -= myedegrees[k].ed-myrinfo->id;
+      graph->minvol -= (xgain+myedegrees[k].gv);
+      where[i] = to;
+
+      IFSET(ctrl->dbglvl, DBG_MOVEINFO, printf("\t\tMoving %6d from %3d to %3d. Gain: [%4d %4d]. Cut: %6d, Vol: %6d\n", 
+            i, from, to, xgain+myedegrees[k].gv, myedegrees[k].ed-myrinfo->id, graph->mincut, graph->minvol));
+
+      /* Update pmat to reflect the move of 'i' */
+      pmat[from*nparts+to] += (myrinfo->id-myedegrees[k].ed);
+      pmat[to*nparts+from] += (myrinfo->id-myedegrees[k].ed);
+      if (pmat[from*nparts+to] == 0) {
+        ndoms[from]--;
+        if (ndoms[from]+1 == maxndoms)
+          maxndoms = ndoms[idxamax(nparts, ndoms)];
+      }
+      if (pmat[to*nparts+from] == 0) {
+        ndoms[to]--;
+        if (ndoms[to]+1 == maxndoms)
+          maxndoms = ndoms[idxamax(nparts, ndoms)];
+      }
+
+      for (j=xadj[i]; j<xadj[i+1]; j++) {
+        ii = adjncy[j];
+        me = where[ii];
+
+        /* Update pmat to reflect the move of 'i' for domains other than 'from' and 'to' */
+        if (me != from && me != to) {
+          pmat[me*nparts+from] -= adjwgt[j];
+          pmat[from*nparts+me] -= adjwgt[j];
+          if (pmat[me*nparts+from] == 0) {
+            ndoms[me]--;
+            if (ndoms[me]+1 == maxndoms)
+              maxndoms = ndoms[idxamax(nparts, ndoms)];
+          }
+          if (pmat[from*nparts+me] == 0) {
+            ndoms[from]--;
+            if (ndoms[from]+1 == maxndoms)
+              maxndoms = ndoms[idxamax(nparts, ndoms)];
+          }
+
+          if (pmat[me*nparts+to] == 0) {
+            ndoms[me]++;
+            if (ndoms[me] > maxndoms) {
+              IFSET(ctrl->dbglvl, DBG_REFINE, printf("You just increased the maxndoms: %d %d\n", ndoms[me], maxndoms));
+              maxndoms = ndoms[me];
+            }
+          }
+          if (pmat[to*nparts+me] == 0) {
+            ndoms[to]++;
+            if (ndoms[to] > maxndoms) {
+              IFSET(ctrl->dbglvl, DBG_REFINE, printf("You just increased the maxndoms: %d %d\n", ndoms[to], maxndoms));
+              maxndoms = ndoms[to];
+            }
+          }
+          pmat[me*nparts+to] += adjwgt[j];
+          pmat[to*nparts+me] += adjwgt[j];
+        }
+      }
+
+      KWayVolUpdate(ctrl, graph, i, from, to, marker, phtable, updind);
+
+      nmoves++;
+
+      /*CheckVolKWayPartitionParams(ctrl, graph, nparts); */
+    }
+
+    IFSET(ctrl->dbglvl, DBG_REFINE,
+       printf("\t[%6d %6d], Balance: %5.3f, Nb: %6d. Nmoves: %5d, Cut: %6d, Vol: %6d\n",
+               pwgts[idxamin(nparts, pwgts)], pwgts[idxamax(nparts, pwgts)],
+               1.0*nparts*pwgts[idxamax(nparts, pwgts)]/tvwgt, graph->nbnd, nmoves, graph->mincut, 
+               graph->minvol));
+
+  }
+
+  GKfree(&marker, &updind, &phtable, LTERM);
+
+  PQueueFree(ctrl, &queue);
+
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nvtxs);
+  idxwspacefree(ctrl, nvtxs);
+}
+
+
+
+
+/*************************************************************************
+* This function updates the edge and volume gains as a result of moving
+* v from 'from' to 'to'.
+* The working arrays marker and phtable are assumed to be initialized to
+* -1, and they left to -1 upon return
+**************************************************************************/
+void KWayVolUpdate(CtrlType *ctrl, GraphType *graph, int v, int from, int to,
+                   idxtype *marker, idxtype *phtable, idxtype *updind)
+{
+  int ii, iii, j, jj, k, kk, l, u, nupd, other, me, myidx; 
+  idxtype *xadj, *vsize, *adjncy, *adjwgt, *where;
+  VEDegreeType *myedegrees, *oedegrees;
+  VRInfoType *myrinfo, *orinfo;
+
+  xadj = graph->xadj;
+  adjncy = graph->adjncy;
+  adjwgt = graph->adjwgt;
+  vsize = graph->vsize;
+  where = graph->where;
+
+  myrinfo = graph->vrinfo+v;
+  myedegrees = myrinfo->edegrees;
+
+
+  /*======================================================================
+   * Remove the contributions on the gain made by 'v'. 
+   *=====================================================================*/
+  for (k=0; k<myrinfo->ndegrees; k++)
+    phtable[myedegrees[k].pid] = k;
+  phtable[from] = k;
+
+  myidx = phtable[to];  /* Keep track of the index in myedegrees of the 'to' domain */
+
+  for (j=xadj[v]; j<xadj[v+1]; j++) {
+    ii = adjncy[j];
+    other = where[ii];
+    orinfo = graph->vrinfo+ii;
+    oedegrees = orinfo->edegrees;
+
+    if (other == from) {
+      for (k=0; k<orinfo->ndegrees; k++) {
+        if (phtable[oedegrees[k].pid] == -1) 
+          oedegrees[k].gv += vsize[v];
+      }
+    }
+    else {
+      ASSERT(phtable[other] != -1);
+
+      if (myedegrees[phtable[other]].ned > 1) {
+        for (k=0; k<orinfo->ndegrees; k++) {
+          if (phtable[oedegrees[k].pid] == -1) 
+            oedegrees[k].gv += vsize[v];
+        }
+      }
+      else { /* There is only one connection */
+        for (k=0; k<orinfo->ndegrees; k++) {
+          if (phtable[oedegrees[k].pid] != -1) 
+            oedegrees[k].gv -= vsize[v];
+        }
+      }
+    }
+  }
+
+  for (k=0; k<myrinfo->ndegrees; k++)
+    phtable[myedegrees[k].pid] = -1;
+  phtable[from] = -1;
+
+
+  /*======================================================================
+   * Update the id/ed of vertex 'v'
+   *=====================================================================*/
+  myrinfo->ed += myrinfo->id-myedegrees[myidx].ed;
+  SWAP(myrinfo->id, myedegrees[myidx].ed, j);
+  SWAP(myrinfo->nid, myedegrees[myidx].ned, j);
+  if (myedegrees[myidx].ed == 0) 
+    myedegrees[myidx] = myedegrees[--myrinfo->ndegrees];
+  else
+    myedegrees[myidx].pid = from;
+
+  /*======================================================================
+   * Update the degrees of adjacent vertices and their volume gains
+   *=====================================================================*/
+  marker[v] = 1;
+  updind[0] = v;
+  nupd = 1;
+  for (j=xadj[v]; j<xadj[v+1]; j++) {
+    ii = adjncy[j];
+    me = where[ii];
+
+    if (!marker[ii]) {  /* The marking is done for boundary and max gv calculations */
+      marker[ii] = 2;
+      updind[nupd++] = ii;
+    }
+
+    myrinfo = graph->vrinfo+ii;
+    if (myrinfo->edegrees == NULL) {
+      myrinfo->edegrees = ctrl->wspace.vedegrees+ctrl->wspace.cdegree;
+      ctrl->wspace.cdegree += xadj[ii+1]-xadj[ii];
+    }
+    myedegrees = myrinfo->edegrees;
+
+    if (me == from) {
+      INC_DEC(myrinfo->ed, myrinfo->id, adjwgt[j]);
+      myrinfo->nid--;
+    } 
+    else if (me == to) {
+      INC_DEC(myrinfo->id, myrinfo->ed, adjwgt[j]);
+      myrinfo->nid++;
+    }
+
+    /* Remove the edgeweight from the 'pid == from' entry of the vertex */
+    if (me != from) {
+      for (k=0; k<myrinfo->ndegrees; k++) {
+        if (myedegrees[k].pid == from) {
+          if (myedegrees[k].ned == 1) {
+            myedegrees[k] = myedegrees[--myrinfo->ndegrees];
+            marker[ii] = 1;  /* You do a complete .gv calculation */
+
+            /* All vertices adjacent to 'ii' need to be updated */
+            for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) {
+              u = adjncy[jj];
+              other = where[u];
+              orinfo = graph->vrinfo+u;
+              oedegrees = orinfo->edegrees;
+
+              for (kk=0; kk<orinfo->ndegrees; kk++) {
+                if (oedegrees[kk].pid == from) {
+                  oedegrees[kk].gv -= vsize[ii];
+                  break;
+                }
+              }
+            }
+          }
+          else {
+            myedegrees[k].ed -= adjwgt[j];
+            myedegrees[k].ned--;
+
+            /* Update the gv due to single 'ii' connection to 'from' */
+            if (myedegrees[k].ned == 1) {
+              /* find the vertex 'u' that 'ii' was connected into 'from' */
+              for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) {
+                u = adjncy[jj];
+                other = where[u];
+                orinfo = graph->vrinfo+u;
+                oedegrees = orinfo->edegrees;
+
+                if (other == from) {
+                  for (kk=0; kk<orinfo->ndegrees; kk++) 
+                    oedegrees[kk].gv += vsize[ii];
+                  break;  
+                }
+              }
+            }
+          }
+
+          break; 
+        }
+      }
+    }
+
+    /* Add the edgeweight to the 'pid == to' entry of the vertex */
+    if (me != to) {
+      for (k=0; k<myrinfo->ndegrees; k++) {
+        if (myedegrees[k].pid == to) {
+          myedegrees[k].ed += adjwgt[j];
+          myedegrees[k].ned++;
+
+          /* Update the gv due to non-single 'ii' connection to 'to' */
+          if (myedegrees[k].ned == 2) {
+            /* find the vertex 'u' that 'ii' was connected into 'to' */
+            for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) {
+              u = adjncy[jj];
+              other = where[u];
+              orinfo = graph->vrinfo+u;
+              oedegrees = orinfo->edegrees;
+
+              if (u != v && other == to) {
+                for (kk=0; kk<orinfo->ndegrees; kk++) 
+                  oedegrees[kk].gv -= vsize[ii];
+                break;  
+              }
+            }
+          }
+          break;
+        }
+      }
+
+      if (k == myrinfo->ndegrees) {
+        myedegrees[myrinfo->ndegrees].pid = to;
+        myedegrees[myrinfo->ndegrees].ed = adjwgt[j];
+        myedegrees[myrinfo->ndegrees++].ned = 1;
+        marker[ii] = 1;  /* You do a complete .gv calculation */
+
+        /* All vertices adjacent to 'ii' need to be updated */
+        for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) {
+          u = adjncy[jj];
+          other = where[u];
+          orinfo = graph->vrinfo+u;
+          oedegrees = orinfo->edegrees;
+
+          for (kk=0; kk<orinfo->ndegrees; kk++) {
+            if (oedegrees[kk].pid == to) {
+              oedegrees[kk].gv += vsize[ii];
+              if (!marker[u]) { /* Need to update boundary etc */
+                marker[u] = 2;
+                updind[nupd++] = u;
+              }
+              break;
+            }
+          }
+        }
+      }
+    }
+
+    ASSERT(myrinfo->ndegrees <= xadj[ii+1]-xadj[ii]);
+  }
+
+  /*======================================================================
+   * Add the contributions on the volume gain due to 'v'
+   *=====================================================================*/
+  myrinfo = graph->vrinfo+v;
+  myedegrees = myrinfo->edegrees;
+  for (k=0; k<myrinfo->ndegrees; k++)
+    phtable[myedegrees[k].pid] = k;
+  phtable[to] = k;
+
+  for (j=xadj[v]; j<xadj[v+1]; j++) {
+    ii = adjncy[j];
+    other = where[ii];
+    orinfo = graph->vrinfo+ii;
+    oedegrees = orinfo->edegrees;
+
+    if (other == to) {
+      for (k=0; k<orinfo->ndegrees; k++) {
+        if (phtable[oedegrees[k].pid] == -1) 
+          oedegrees[k].gv -= vsize[v];
+      }
+    }
+    else {
+      ASSERT(phtable[other] != -1);
+
+      if (myedegrees[phtable[other]].ned > 1) {
+        for (k=0; k<orinfo->ndegrees; k++) {
+          if (phtable[oedegrees[k].pid] == -1) 
+            oedegrees[k].gv -= vsize[v];
+        }
+      }
+      else { /* There is only one connection */
+        for (k=0; k<orinfo->ndegrees; k++) {
+          if (phtable[oedegrees[k].pid] != -1) 
+            oedegrees[k].gv += vsize[v];
+        }
+      }
+    }
+  }
+  for (k=0; k<myrinfo->ndegrees; k++)
+    phtable[myedegrees[k].pid] = -1;
+  phtable[to] = -1;
+
+
+  /*======================================================================
+   * Recompute the volume information of the 'hard' nodes, and update the
+   * max volume gain for all the update vertices
+   *=====================================================================*/
+  ComputeKWayVolume(graph, nupd, updind, marker, phtable);
+
+
+  /*======================================================================
+   * Maintain a consistent boundary
+   *=====================================================================*/
+  for (j=0; j<nupd; j++) {
+    k = updind[j];
+    marker[k] = 0;
+    myrinfo = graph->vrinfo+k;
+
+    if ((myrinfo->gv >= 0 || myrinfo->ed-myrinfo->id >= 0) && graph->bndptr[k] == -1)
+      BNDInsert(graph->nbnd, graph->bndind, graph->bndptr, k);
+
+    if (myrinfo->gv < 0 && myrinfo->ed-myrinfo->id < 0 && graph->bndptr[k] != -1)
+      BNDDelete(graph->nbnd, graph->bndind, graph->bndptr, k);
+  }
+
+}
+
+
+
+
+/*************************************************************************
+* This function computes the initial id/ed 
+**************************************************************************/
+void ComputeKWayVolume(GraphType *graph, int nupd, idxtype *updind, idxtype *marker, idxtype *phtable)
+{
+  int ii, iii, i, j, k, kk, l, nvtxs, me, other, pid;
+  idxtype *xadj, *vsize, *adjncy, *adjwgt, *where;
+  VRInfoType *rinfo, *myrinfo, *orinfo;
+  VEDegreeType *myedegrees, *oedegrees;
+
+  nvtxs = graph->nvtxs;
+  xadj = graph->xadj;
+  vsize = graph->vsize;
+  adjncy = graph->adjncy;
+  adjwgt = graph->adjwgt;
+  where = graph->where;
+  rinfo = graph->vrinfo;
+
+
+  /*------------------------------------------------------------
+  / Compute now the iv/ev degrees
+  /------------------------------------------------------------*/
+  for (iii=0; iii<nupd; iii++) {
+    i = updind[iii];
+    me = where[i];
+
+    myrinfo = rinfo+i;
+    myedegrees = myrinfo->edegrees;
+
+    if (marker[i] == 1) {  /* Only complete gain updates go through */
+      for (k=0; k<myrinfo->ndegrees; k++) 
+        myedegrees[k].gv = 0;
+
+      for (j=xadj[i]; j<xadj[i+1]; j++) {
+        ii = adjncy[j];
+        other = where[ii];
+        orinfo = rinfo+ii;
+        oedegrees = orinfo->edegrees;
+
+        for (kk=0; kk<orinfo->ndegrees; kk++) 
+          phtable[oedegrees[kk].pid] = kk;
+        phtable[other] = 1;
+
+        if (me == other) {
+          /* Find which domains 'i' is connected and 'ii' is not and update their gain */
+          for (k=0; k<myrinfo->ndegrees; k++) {
+            if (phtable[myedegrees[k].pid] == -1)
+              myedegrees[k].gv -= vsize[ii];
+          }
+        }
+        else {
+          ASSERT(phtable[me] != -1);
+
+          /* I'm the only connection of 'ii' in 'me' */
+          if (oedegrees[phtable[me]].ned == 1) { 
+            /* Increase the gains for all the common domains between 'i' and 'ii' */
+            for (k=0; k<myrinfo->ndegrees; k++) {
+              if (phtable[myedegrees[k].pid] != -1) 
+                myedegrees[k].gv += vsize[ii];
+            }
+          }
+          else {
+            /* Find which domains 'i' is connected and 'ii' is not and update their gain */
+            for (k=0; k<myrinfo->ndegrees; k++) {
+              if (phtable[myedegrees[k].pid] == -1) 
+                myedegrees[k].gv -= vsize[ii];
+            }
+          }
+        }
+
+        for (kk=0; kk<orinfo->ndegrees; kk++) 
+          phtable[oedegrees[kk].pid] = -1;
+        phtable[other] = -1;
+  
+      }
+    }
+
+    myrinfo->gv = -MAXIDX;
+    for (k=0; k<myrinfo->ndegrees; k++) {
+      if (myedegrees[k].gv > myrinfo->gv)
+        myrinfo->gv = myedegrees[k].gv;
+    }
+    if (myrinfo->ed > 0 && myrinfo->id == 0)
+      myrinfo->gv += vsize[i];
+
+  }
+
+}
+
+
+
+/*************************************************************************
+* This function computes the total volume
+**************************************************************************/
+int ComputeVolume(GraphType *graph, idxtype *where)
+{
+  int i, j, k, me, nvtxs, nparts, totalv;
+  idxtype *xadj, *adjncy, *vsize, *marker;
+
+
+  nvtxs = graph->nvtxs;
+  xadj = graph->xadj;
+  adjncy = graph->adjncy;
+  vsize = (graph->vsize == NULL ? graph->vwgt : graph->vsize);
+
+  nparts = where[idxamax(nvtxs, where)]+1;
+  marker = idxsmalloc(nparts, -1, "ComputeVolume: marker");
+
+  totalv = 0;
+
+  for (i=0; i<nvtxs; i++) {
+    marker[where[i]] = i;
+    for (j=xadj[i]; j<xadj[i+1]; j++) {
+      k = where[adjncy[j]];
+      if (marker[k] != i) {
+        marker[k] = i;
+        totalv += vsize[i];
+      }
+    }
+  }
+
+  free(marker);
+
+  return totalv;
+}
+
+
+
+
+
+/*************************************************************************
+* This function computes the initial id/ed 
+**************************************************************************/
+void CheckVolKWayPartitionParams(CtrlType *ctrl, GraphType *graph, int nparts)
+{
+  int i, ii, j, k, kk, l, nvtxs, nbnd, mincut, minvol, me, other, pid;
+  idxtype *xadj, *vsize, *adjncy, *adjwgt, *pwgts, *where, *bndind, *bndptr;
+  VRInfoType *rinfo, *myrinfo, *orinfo, tmprinfo;
+  VEDegreeType *myedegrees, *oedegrees, *tmpdegrees;
+
+  nvtxs = graph->nvtxs;
+  xadj = graph->xadj;
+  vsize = graph->vsize;
+  adjncy = graph->adjncy;
+  adjwgt = graph->adjwgt;
+  where = graph->where;
+  rinfo = graph->vrinfo;
+
+  tmpdegrees = (VEDegreeType *)GKmalloc(nparts*sizeof(VEDegreeType), "CheckVolKWayPartitionParams: tmpdegrees");
+
+  /*------------------------------------------------------------
+  / Compute now the iv/ev degrees
+  /------------------------------------------------------------*/
+  for (i=0; i<nvtxs; i++) {
+    me = where[i];
+
+    myrinfo = rinfo+i;
+    myedegrees = myrinfo->edegrees;
+
+    for (k=0; k<myrinfo->ndegrees; k++)
+      tmpdegrees[k] = myedegrees[k];
+
+    tmprinfo.ndegrees = myrinfo->ndegrees;
+    tmprinfo.id = myrinfo->id;
+    tmprinfo.ed = myrinfo->ed;
+
+    myrinfo = &tmprinfo;
+    myedegrees = tmpdegrees;
+
+
+    for (k=0; k<myrinfo->ndegrees; k++)
+      myedegrees[k].gv = 0;
+
+    for (j=xadj[i]; j<xadj[i+1]; j++) {
+      ii = adjncy[j];
+      other = where[ii];
+      orinfo = rinfo+ii;
+      oedegrees = orinfo->edegrees;
+
+      if (me == other) {
+        /* Find which domains 'i' is connected and 'ii' is not and update their gain */
+        for (k=0; k<myrinfo->ndegrees; k++) {
+          pid = myedegrees[k].pid;
+          for (kk=0; kk<orinfo->ndegrees; kk++) {
+            if (oedegrees[kk].pid == pid)
+              break;
+          }
+          if (kk == orinfo->ndegrees) 
+            myedegrees[k].gv -= vsize[ii];
+        }
+      }
+      else {
+        /* Find the orinfo[me].ed and see if I'm the only connection */
+        for (k=0; k<orinfo->ndegrees; k++) {
+          if (oedegrees[k].pid == me)
+            break;
+        }
+
+        if (oedegrees[k].ned == 1) { /* I'm the only connection of 'ii' in 'me' */
+          for (k=0; k<myrinfo->ndegrees; k++) {
+            if (myedegrees[k].pid == other) {
+              myedegrees[k].gv += vsize[ii];
+              break;
+            }
+          }
+
+          /* Increase the gains for all the common domains between 'i' and 'ii' */
+          for (k=0; k<myrinfo->ndegrees; k++) {
+            if ((pid = myedegrees[k].pid) == other)
+              continue;
+            for (kk=0; kk<orinfo->ndegrees; kk++) {
+              if (oedegrees[kk].pid == pid) {
+                myedegrees[k].gv += vsize[ii];
+                break;
+              }
+            }
+          }
+
+        }
+        else {
+          /* Find which domains 'i' is connected and 'ii' is not and update their gain */
+          for (k=0; k<myrinfo->ndegrees; k++) {
+            if ((pid = myedegrees[k].pid) == other)
+              continue;
+            for (kk=0; kk<orinfo->ndegrees; kk++) {
+              if (oedegrees[kk].pid == pid)
+                break;
+            }
+            if (kk == orinfo->ndegrees) 
+              myedegrees[k].gv -= vsize[ii];
+          }
+        }
+      }
+    }
+
+    myrinfo = rinfo+i;
+    myedegrees = myrinfo->edegrees;
+
+    for (k=0; k<myrinfo->ndegrees; k++) {
+      pid = myedegrees[k].pid;
+      for (kk=0; kk<tmprinfo.ndegrees; kk++) {
+        if (tmpdegrees[kk].pid == pid) {
+          if (tmpdegrees[kk].gv != myedegrees[k].gv)
+            printf("[%d %d %d %d]\n", i, pid, myedegrees[k].gv, tmpdegrees[kk].gv);
+          break;
+        }
+      }
+    }
+
+  }
+
+  free(tmpdegrees);
+
+}
+
+
+/*************************************************************************
+* This function computes the subdomain graph
+**************************************************************************/
+void ComputeVolSubDomainGraph(GraphType *graph, int nparts, idxtype *pmat, idxtype *ndoms)
+{
+  int i, j, k, me, nvtxs, ndegrees;
+  idxtype *xadj, *adjncy, *adjwgt, *where;
+  VRInfoType *rinfo;
+  VEDegreeType *edegrees;
+
+  nvtxs = graph->nvtxs;
+  xadj = graph->xadj;
+  adjncy = graph->adjncy;
+  adjwgt = graph->adjwgt;
+  where = graph->where;
+  rinfo = graph->vrinfo;
+
+  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 computes the subdomain graph
+**************************************************************************/
+void EliminateVolSubDomainEdges(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 = idxset(nparts, 0, graph->pwgts);
+
+  maxpwgt = idxwspacemalloc(ctrl, nparts);
+  ndoms = idxwspacemalloc(ctrl, nparts);
+  otherpmat = idxwspacemalloc(ctrl, nparts);
+  ind = idxwspacemalloc(ctrl, nvtxs);
+  pmat = idxset(nparts*nparts, 0, ctrl->wspace.pmat);
+
+  cand = (KeyValueType *)GKmalloc(nparts*sizeof(KeyValueType), "EliminateSubDomainEdges: cand");
+  cand2 = (KeyValueType *)GKmalloc(nparts*sizeof(KeyValueType), "EliminateSubDomainEdges: cand");
+
+  /* Compute the pmat matrix */
+  for (i=0; i<nvtxs; i++) {
+    me = where[i];
+    pwgts[me] += vwgt[i];
+    for (j=xadj[i]; j<xadj[i+1]; j++) {
+      k = adjncy[j];
+      if (where[k] != me) 
+        pmat[me*nparts+where[k]] += adjwgt[j];
+    }
+  }
+
+  /* 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;
+
+  /* Determine the domain connectivity */
+  for (i=0; i<nparts; i++) {
+    for (k=0, j=0; j<nparts; j++) {
+      if (pmat[i*nparts+j] > 0)
+        k++;
+    }
+    ndoms[i] = k;
+  }
+
+  /* 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\n", total, max, avg); */
+
+    if (max < 1.5*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++) {
+          k = adjncy[j];
+          if (where[k] != other) 
+            otherpmat[where[k]] += adjwgt[j];
+        }
+      }
+
+      for (ncand=0, i=0; i<nparts; i++) {
+        if (otherpmat[i] > 0) {
+          cand[ncand].key = -otherpmat[i];
+          cand[ncand++].val = i;
+        }
+      }
+      ikeysort(ncand, cand);
+
+      /* 
+       * Go through and the select the first domain that is common with 'me', and
+       * does not increase the ndoms[target] higher than my ndoms, subject to the
+       * maxpwgt constraint. Traversal is done from the mostly connected to the least.
+       */
+      target = target2 = -1;
+      for (i=0; i<ncand; i++) {
+        k = cand[i].val;
+
+        if (mypmat[k] > 0) {
+          if (pwgts[k] + cpwgt > maxpwgt[k])  /* Check if balance will go off */
+            continue;
+
+          for (j=0; j<nparts; j++) {
+            if (otherpmat[j] > 0 && ndoms[j] >= ndoms[me]-1 && pmat[nparts*j+k] == 0)
+              break;
+          }
+          if (j == nparts) { /* No bad second level effects */
+            for (nadd=0, j=0; j<nparts; j++) {
+              if (otherpmat[j] > 0 && pmat[nparts*k+j] == 0)
+                nadd++;
+            }
+
+            /*printf("\t\tto=%d, nadd=%d, %d\n", k, nadd, ndoms[k]);*/
+            if (target2 == -1 && ndoms[k]+nadd < ndoms[me]) {
+              target2 = k;
+            }
+            if (nadd == 0) {
+              target = k;
+              break;
+            }
+          }
+        }
+      }
+      if (target == -1 && target2 != -1)
+        target = target2;
+
+      if (target == -1) {
+        /* printf("\t\tCould not make the move\n");*/
+        continue;
+      }
+
+      /*printf("\t\tMoving to %d\n", target);*/
+
+      /* Update the partition weights */
+      INC_DEC(pwgts[target], pwgts[other], cpwgt);
+
+      /* Set all nind vertices to belong to 'target' */
+      for (ii=0; ii<nind; ii++) {
+        i = ind[ii];
+        where[i] = target;
+
+        /* First remove any contribution that this vertex may have made */
+        for (j=xadj[i]; j<xadj[i+1]; j++) {
+          k = adjncy[j];
+          if (where[k] != other) {
+            if (pmat[nparts*other + where[k]] == 0)
+              printf("Something wrong\n");
+            pmat[nparts*other + where[k]] -= adjwgt[j];
+            if (pmat[nparts*other + where[k]] == 0)
+              ndoms[other]--;
+
+            if (pmat[nparts*where[k] + other] == 0)
+              printf("Something wrong\n");
+            pmat[nparts*where[k] + other] -= adjwgt[j];
+            if (pmat[nparts*where[k] + other] == 0)
+              ndoms[where[k]]--;
+          }
+        }
+
+        /* Next add the new contributions as a result of the move */
+        for (j=xadj[i]; j<xadj[i+1]; j++) {
+          k = adjncy[j];
+          if (where[k] != target) {
+            if (pmat[nparts*target + where[k]] == 0)
+              ndoms[target]++;
+            pmat[nparts*target + where[k]] += adjwgt[j];
+  
+            if (pmat[nparts*where[k] + target] == 0)
+              ndoms[where[k]]++;
+            pmat[nparts*where[k] + target] += adjwgt[j];
+          }
+        }
+      }
+
+      move = 1;
+      break;
+    }
+
+    if (move == 0)
+      break;
+  }
+
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nvtxs);
+
+  GKfree(&cand, &cand2, 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 EliminateVolComponents(CtrlType *ctrl, GraphType *graph, int nparts, float *tpwgts, float ubfactor)
+{
+  int i, ii, j, jj, k, me, nvtxs, tvwgt, first, last, nleft, ncmps, cwgt, ncand, other, target, deltawgt;
+  idxtype *xadj, *adjncy, *vwgt, *adjwgt, *where, *pwgts, *maxpwgt;
+  idxtype *cpvec, *touched, *perm, *todo, *cind, *cptr, *npcmps;
+  KeyValueType *cand;
+  int recompute=0;
+
+  nvtxs = graph->nvtxs;
+  xadj = graph->xadj;
+  adjncy = graph->adjncy;
+  vwgt = graph->vwgt;
+  adjwgt = graph->adjwgt;
+
+  where = graph->where;
+  pwgts = idxset(nparts, 0, graph->pwgts);
+
+  touched = idxset(nvtxs, 0, idxwspacemalloc(ctrl, nvtxs));
+  cptr    = idxwspacemalloc(ctrl, nvtxs+1);
+  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 */
+    cand = (KeyValueType *)GKmalloc(nparts*sizeof(KeyValueType), "EliminateSubDomainEdges: cand");
+
+    /* First determine the partition sizes and max allowed load imbalance */
+    for (i=0; i<nvtxs; i++) 
+      pwgts[where[i]] += vwgt[i];
+    tvwgt = idxsum(nparts, pwgts);
+    for (i=0; i<nparts; i++)
+      maxpwgt[i] = ubfactor*tpwgts[i]*tvwgt;
+
+    deltawgt = tvwgt/(100*nparts);
+    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 connectivity */
+      idxset(nparts, 0, cpvec);
+      for (cwgt=0, j=cptr[i]; j<cptr[i+1]; j++) {
+        ii = cind[j];
+        cwgt += vwgt[ii];
+        for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) {
+          other = where[adjncy[jj]];
+          if (me != other)
+            cpvec[other] += adjwgt[jj];
+        }
+      }
+
+      /*printf("\tCmp weight: %d\n", cwgt);*/
+
+      if (cwgt > .30*pwgts[me])
+        continue;  /* Skip the component if it is over 30% of the weight */
+
+      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;
+
+      ikeysort(ncand, cand);
+
+      target = -1;
+      for (j=0; j<ncand; j++) {
+        k = cand[j].val;
+        if (cwgt < deltawgt || pwgts[k] + cwgt < maxpwgt[k]) {
+          target = k;
+          break;
+        }
+      }
+
+      /*printf("\tMoving it to %d [%d]\n", target, cpvec[target]);*/
+
+      if (target != -1) {
+        /* Assign all the vertices of 'me' to 'target' and update data structures */
+        pwgts[me] -= cwgt;
+        pwgts[target] += cwgt;
+        npcmps[me]--;
+
+        for (j=cptr[i]; j<cptr[i+1]; j++) 
+          where[cind[j]] = target;
+
+        graph->mincut -= cpvec[target];
+        recompute = 1;
+      }
+    }
+
+    free(cand);
+  }
+
+  if (recompute) {
+    int ttlv;
+    idxtype *marker;
+
+    marker = idxset(nparts, -1, cpvec);
+    for (ttlv=0, i=0; i<nvtxs; i++) {
+      marker[where[i]] = i;
+      for (j=xadj[i]; j<xadj[i+1]; j++) {
+        if (marker[where[adjncy[j]]] != i) {
+          ttlv += graph->vsize[i];
+          marker[where[adjncy[j]]] = i;
+        }
+      }
+    }
+    graph->minvol = ttlv;
+  }
+
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nparts);
+  idxwspacefree(ctrl, nvtxs);
+  idxwspacefree(ctrl, nvtxs);
+  idxwspacefree(ctrl, nvtxs);
+  idxwspacefree(ctrl, nvtxs);
+  idxwspacefree(ctrl, nvtxs+1);
+
+}
+