kwayvolfm.c

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

/*
 * kwayvolfm.c
 *
 * This file contains code that implements the multilevel k-way refinement
 *
 * Started 7/8/98
 * George
 *
 * $Id: kwayvolfm.c 2501 2007-11-20 02:33:29Z benkirk $
 *
 */

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