kwayvolfm.c

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

    }

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

  }