serial.c

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

/*
 * serial.c
 *
 * This file contains code that implements k-way refinement
 *
 * Started 7/28/97
 * George
 *
 * $Id: serial.c 2501 2007-11-20 02:33:29Z benkirk $
 *
 */

#include <parmetislib.h>


/*************************************************************************
* This function performs k-way refinement
**************************************************************************/
void Moc_SerialKWayAdaptRefine(GraphType *graph, int nparts, idxtype *home,
     float *orgubvec, int npasses)
{
  int i, ii, iii, j, k;
  int nvtxs, ncon, pass, nmoves, myndegrees;
  int from, me, myhome, to, oldcut, gain, tmp;
  idxtype *xadj, *adjncy, *adjwgt;
  idxtype *where;
  EdgeType *mydegrees;
  RInfoType *rinfo, *myrinfo;
  float *npwgts, *nvwgt, *minwgt, *maxwgt, ubvec[MAXNCON];
  int gain_is_greater, gain_is_same, fit_in_to, fit_in_from, going_home;
  int zero_gain, better_balance_ft, better_balance_tt;
  KeyValueType *cand;
int mype;
MPI_Comm_rank(MPI_COMM_WORLD, &mype);

  nvtxs = graph->nvtxs;
  ncon = graph->ncon;
  xadj = graph->xadj;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  where = graph->where;
  rinfo = graph->rinfo;
  npwgts = graph->gnpwgts;
  
  /* Setup the weight intervals of the various subdomains */
  cand = (KeyValueType *)GKmalloc(nvtxs*sizeof(KeyValueType), "cand");
  minwgt =  fmalloc(nparts*ncon, "minwgt");
  maxwgt = fmalloc(nparts*ncon, "maxwgt");

  ComputeHKWayLoadImbalance(ncon, nparts, npwgts, ubvec);
  for (i=0; i<ncon; i++)
    ubvec[i] = amax(ubvec[i], orgubvec[i]);

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

  for (pass=0; pass<npasses; pass++) {
    oldcut = graph->mincut;

    for (i=0; i<nvtxs; i++) {
      cand[i].key = rinfo[i].id-rinfo[i].ed;
      cand[i].val = i;
    }
    ikeysort(nvtxs, cand);

    nmoves = 0;
    for (iii=0; iii<nvtxs; iii++) {
      i = cand[iii].val;

      myrinfo = rinfo+i;

      if (myrinfo->ed >= myrinfo->id) {
        from = where[i];
        myhome = home[i];
        nvwgt = graph->nvwgt+i*ncon;

        if (myrinfo->id > 0 &&
        AreAllHVwgtsBelow(ncon, 1.0, npwgts+from*ncon, -1.0, nvwgt, minwgt+from*ncon)) 
          continue;

        mydegrees = myrinfo->degrees;
        myndegrees = myrinfo->ndegrees;

        for (k=0; k<myndegrees; k++) {
          to = mydegrees[k].edge;
          gain = mydegrees[k].ewgt - myrinfo->id; 
          if (gain >= 0 && 
             (AreAllHVwgtsBelow(ncon, 1.0, npwgts+to*ncon, 1.0, nvwgt, maxwgt+to*ncon) ||
             IsHBalanceBetterFT(ncon,npwgts+from*ncon,npwgts+to*ncon,nvwgt,ubvec))) {
            break;
          }
        }

        /* break out if you did not find a candidate */
        if (k == myndegrees)
          continue;

        for (j=k+1; j<myndegrees; j++) {
          to = mydegrees[j].edge;
          going_home = (myhome == to);
          gain_is_same = (mydegrees[j].ewgt == mydegrees[k].ewgt);
          gain_is_greater = (mydegrees[j].ewgt > mydegrees[k].ewgt);
          fit_in_to = AreAllHVwgtsBelow(ncon,1.0,npwgts+to*ncon,1.0,nvwgt,maxwgt+to*ncon);
          better_balance_ft = IsHBalanceBetterFT(ncon,npwgts+from*ncon,
                              npwgts+to*ncon,nvwgt,ubvec);
          better_balance_tt = IsHBalanceBetterTT(ncon,npwgts+mydegrees[k].edge*ncon,
                              npwgts+to*ncon,nvwgt,ubvec);

          if (
               (gain_is_greater &&
                 (fit_in_to ||
                  better_balance_ft)
               )
            ||
               (gain_is_same &&
                 (
                   (fit_in_to &&
                    going_home)
                ||
                    better_balance_tt
                 )
               )
             ) {
            k = j;
          }
        }

        to = mydegrees[k].edge;
        going_home = (myhome == to);
        zero_gain = (mydegrees[k].ewgt == myrinfo->id);

        fit_in_from = AreAllHVwgtsBelow(ncon,1.0,npwgts+from*ncon,0.0,npwgts+from*ncon,
                      maxwgt+from*ncon);
        better_balance_ft = IsHBalanceBetterFT(ncon,npwgts+from*ncon,
                            npwgts+to*ncon,nvwgt,ubvec);

        if (zero_gain &&
            !going_home &&
            !better_balance_ft &&
            fit_in_from)
          continue;

        /*=====================================================================
        * If we got here, we can now move the vertex from 'from' to 'to' 
        *======================================================================*/
        graph->mincut -= mydegrees[k].ewgt-myrinfo->id;

        /* Update where, weight, and ID/ED information of the vertex you moved */
        saxpy2(ncon, 1.0, nvwgt, 1, npwgts+to*ncon, 1);
        saxpy2(ncon, -1.0, nvwgt, 1, npwgts+from*ncon, 1);
        where[i] = to;
        myrinfo->ed += myrinfo->id-mydegrees[k].ewgt;
        SWAP(myrinfo->id, mydegrees[k].ewgt, tmp);

        if (mydegrees[k].ewgt == 0) {
          myrinfo->ndegrees--;
          mydegrees[k].edge = mydegrees[myrinfo->ndegrees].edge;
          mydegrees[k].ewgt = mydegrees[myrinfo->ndegrees].ewgt;
        }
        else
          mydegrees[k].edge = from;

        /* Update the degrees of adjacent vertices */
        for (j=xadj[i]; j<xadj[i+1]; j++) {
          ii = adjncy[j];
          me = where[ii];

          myrinfo = rinfo+ii;
          mydegrees = myrinfo->degrees;

          if (me == from) {
            INC_DEC(myrinfo->ed, myrinfo->id, adjwgt[j]);
          }
          else {
            if (me == to) {
              INC_DEC(myrinfo->id, myrinfo->ed, adjwgt[j]);
            }
          }

          /* Remove contribution of the ed from 'from' */
          if (me != from) {
            for (k=0; k<myrinfo->ndegrees; k++) {
              if (mydegrees[k].edge == from) {
                if (mydegrees[k].ewgt == adjwgt[j]) {
                  myrinfo->ndegrees--;
                  mydegrees[k].edge = mydegrees[myrinfo->ndegrees].edge;
                  mydegrees[k].ewgt = mydegrees[myrinfo->ndegrees].ewgt;
                }
                else
                  mydegrees[k].ewgt -= adjwgt[j];
                break;
              }
            }
          }

          /* Add contribution of the ed to 'to' */
          if (me != to) {
            for (k=0; k<myrinfo->ndegrees; k++) {
              if (mydegrees[k].edge == to) {
                mydegrees[k].ewgt += adjwgt[j];
                break;
              }
            }
            if (k == myrinfo->ndegrees) {
              mydegrees[myrinfo->ndegrees].edge = to;
              mydegrees[myrinfo->ndegrees++].ewgt = adjwgt[j];
            }
          }

        }
        nmoves++;
      }
    }

    if (graph->mincut == oldcut)
      break;
  }

  GKfree((void **)&minwgt, (void **)&maxwgt, (void **)&cand, LTERM);

  return;
}


/*************************************************************************
* This function computes the initial id/ed
**************************************************************************/
void Moc_ComputeSerialPartitionParams(GraphType *graph, int nparts,
     EdgeType *degrees)
{
  int i, j, k;
  int nvtxs, nedges, ncon, mincut, me, other;
  idxtype *xadj, *adjncy, *adjwgt, *where;
  RInfoType *rinfo, *myrinfo;
  EdgeType *mydegrees;
  float *nvwgt, *npwgts;
int mype;
MPI_Comm_rank(MPI_COMM_WORLD, &mype);


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

  npwgts = sset(ncon*nparts, 0.0, graph->gnpwgts);

  /*------------------------------------------------------------
  / Compute now the id/ed degrees
  /------------------------------------------------------------*/
  nedges = mincut = 0;
  for (i=0; i<nvtxs; i++) {
    me = where[i];
    saxpy2(ncon, 1.0, nvwgt+i*ncon, 1, npwgts+me*ncon, 1);

    myrinfo = rinfo+i;
    myrinfo->id = myrinfo->ed = myrinfo->ndegrees = 0;
    myrinfo->degrees = degrees + nedges;
    nedges += xadj[i+1]-xadj[i];

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

    mincut += myrinfo->ed;

    /* Time to compute the particular external degrees */
    if (myrinfo->ed > 0) {
      mydegrees = myrinfo->degrees;

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

  graph->mincut = mincut/2;

  return;
}


/*************************************************************************
* This function checks if the vertex weights of two vertices are below 
* a given set of values
**************************************************************************/
int AreAllHVwgtsBelow(int ncon, float alpha, float *vwgt1, float beta, float *vwgt2, float *limit)
{
  int i;

  for (i=0; i<ncon; i++)
    if (alpha*vwgt1[i] + beta*vwgt2[i] > limit[i])
      return 0;

  return 1;
}


/*************************************************************************
* This function computes the load imbalance over all the constrains
* For now assume that we just want balanced partitionings
**************************************************************************/ 
void ComputeHKWayLoadImbalance(int ncon, int nparts, float *npwgts, float *lbvec)
{
  int i, j;
  float max;

  for (i=0; i<ncon; i++) {
    max = 0.0;
    for (j=0; j<nparts; j++) {
      if (npwgts[j*ncon+i] > max)
        max = npwgts[j*ncon+i];
    }

    lbvec[i] = max*nparts;
  }
}


/**************************************************************
*  This subroutine remaps a partitioning on a single processor
**************************************************************/
void SerialRemap(GraphType *graph, int nparts, idxtype *base, idxtype *scratch,
     idxtype *remap, float *tpwgts)
{
  int i, ii, j, k;
  int nvtxs, nmapped, max_mult;
  int from, to, current_from, smallcount, bigcount;
  KeyValueType *flowto, *bestflow;
  KeyKeyValueType *sortvtx;
  idxtype *vsize, *htable, *map, *rowmap;

  nvtxs = graph->nvtxs;
  vsize = graph->vsize;
  max_mult = amin(MAX_NPARTS_MULTIPLIER, nparts);

  sortvtx = (KeyKeyValueType *)GKmalloc(nvtxs*sizeof(KeyKeyValueType), "sortvtx");
  flowto = (KeyValueType *)GKmalloc((nparts*max_mult+nparts)*sizeof(KeyValueType), "flowto");
  bestflow = flowto+nparts;
  map = htable = idxsmalloc(nparts*2, -1, "htable");
  rowmap = map+nparts;

  for (i=0; i<nvtxs; i++) {
    sortvtx[i].key1 = base[i];
    sortvtx[i].key2 = vsize[i];
    sortvtx[i].val = i;
  }

  qsort((void *)sortvtx, (size_t)nvtxs, (size_t)sizeof(KeyKeyValueType), SSMIncKeyCmp);

  for (j=0; j<nparts; j++) {
    flowto[j].key = 0;
    flowto[j].val = j;
  }

  /* this step has nparts*nparts*log(nparts) computational complexity */
  bigcount = smallcount = current_from = 0;
  for (ii=0; ii<nvtxs; ii++) {
    i = sortvtx[ii].val;
    from = base[i];
    to = scratch[i];

    if (from > current_from) {
      /* reset the hash table */
      for (j=0; j<smallcount; j++)
        htable[flowto[j].val] = -1;
      ASSERTS(idxsum(nparts, htable) == -nparts);

      ikeysort(smallcount, flowto);

      for (j=0; j<amin(smallcount, max_mult); j++, bigcount++) {
        bestflow[bigcount].key = flowto[j].key;
        bestflow[bigcount].val = current_from*nparts+flowto[j].val;
      }

      smallcount = 0;
      current_from = from;
    }

    if (htable[to] == -1) {
      htable[to] = smallcount;
      flowto[smallcount].key = -vsize[i];
      flowto[smallcount].val = to;
      smallcount++;
    }
    else {
      flowto[htable[to]].key += -vsize[i];
    }
  }

  /* reset the hash table */
  for (j=0; j<smallcount; j++)
    htable[flowto[j].val] = -1;