kwayvolfm.c

一个用来实现偏微分方程中网格的计算库

  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;