kwayvolfm.c

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

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

}