/*************************************************************************
* This function finds a matching using the HEM heuristic
**************************************************************************/
void Match_HEM(CtrlType *ctrl, GraphType *graph)
{
  int i, ii, j, k, nvtxs, cnvtxs, maxidx, dim;
  idxtype *xadj, *vwgt, *adjncy;
  idxtype *match, *cmap, *perm, *tperm;
  realtype curwgt, maxwgt;
  realtype *vvol, *vsurf, *adjwgt, *adjwgtsum;

  dim       = ctrl->dim;
  nvtxs     = graph->nvtxs;
  xadj      = graph->xadj;
  vwgt      = graph->vwgt;
  vvol      = graph->vvol;
  vsurf     = graph->vsurf;
  adjncy    = graph->adjncy;
  adjwgt    = graph->adjwgt;
  adjwgtsum = graph->adjwgtsum;

  cmap = graph->cmap = idxsmalloc(nvtxs, -1, "cmap");
  match = idxsmalloc(nvtxs, -1, "match");

  perm = idxmalloc(nvtxs, "perm");
  tperm = idxmalloc(nvtxs, "tperm");

  RandomPermute(nvtxs, tperm, 1);
  BucketSortKeysInc(nvtxs, vwgt[iamax(nvtxs, vwgt)], vwgt, tperm, perm);
  /* RandomPermute(nvtxs, perm, 1);  */

  cnvtxs = 0;

  /* Compute a heavy-edge style matching giving preferance to small vertices */
  for (ii=0; ii<nvtxs; ii++) {
     i = perm[ii];

     if (match[i] == UNMATCHED) {
       maxidx = i;
       maxwgt = 0.0;

       /* Find a heavy-edge matching, subject to maxvwgt constraints */
       for (j=xadj[i]; j<xadj[i+1]; j++) {
          k = adjncy[j];
          curwgt = 1.0/ARATIO2(dim, vsurf[i]+vsurf[k]+adjwgtsum[i]+adjwgtsum[k]-
                   2.0*adjwgt[j], vvol[i]+vvol[k]);
          if (match[k] == UNMATCHED && vwgt[i]+vwgt[k] <= ctrl->maxsize &&
              curwgt > maxwgt) {
            maxwgt = curwgt;
            maxidx = k;
          }
       }

       cmap[i] = cmap[maxidx] = cnvtxs++;
       match[i] = maxidx;
       match[maxidx] = i;
     }
  }

  CreateCoarseGraph(graph, cnvtxs, match, perm);

  IMfree((void**)&tperm, &perm, &match, LTERM);
}

/*************************************************************************
* This function computes movement statistics for adaptive refinement
* schemes
**************************************************************************/
void ComputeMoveStatistics(CtrlType *ctrl, GraphType *graph, int *nmoved, int *maxin, int *maxout)
{
    int i, j, nvtxs;
    idxtype *vwgt, *where;
    idxtype *lpvtxs, *gpvtxs;

    nvtxs = graph->nvtxs;
    vwgt = graph->vwgt;
    where = graph->where;

    lpvtxs = idxsmalloc(ctrl->nparts, 0, "ComputeMoveStatistics: lpvtxs");
    gpvtxs = idxsmalloc(ctrl->nparts, 0, "ComputeMoveStatistics: gpvtxs");

    for (j=i=0; i<nvtxs; i++) {
        lpvtxs[where[i]]++;
        if (where[i] != ctrl->mype)
            j++;
    }

    /* PrintVector(ctrl, ctrl->npes, 0, lpvtxs, "Lpvtxs: "); */

    MPI_Allreduce((void *)lpvtxs, (void *)gpvtxs, ctrl->nparts, IDX_DATATYPE, MPI_SUM, ctrl->comm);

    *nmoved = GlobalSESum(ctrl, j);
    *maxout = GlobalSEMax(ctrl, j);
    *maxin = GlobalSEMax(ctrl, gpvtxs[ctrl->mype]-(nvtxs-j));

    GKfree((void **)&lpvtxs, (void **)&gpvtxs, LTERM);
}

/*****************************************************************************
* This function creates the nodal graph of a finite element mesh
******************************************************************************/
void QUADNODALMETIS(int nelmnts, int nvtxs, idxtype *elmnts, idxtype *dxadj, idxtype *dadjncy)
{
    int i, j, jj, k, kk, /*kkk, l, m, n,*/ nedges;
    idxtype *nptr, *nind;
    idxtype *mark;
    int table[4][2] = {{1, 3},
        {0, 2},
        {1, 3},
        {0, 2}
    };

    /* Construct the node-element list first */
    nptr = idxsmalloc(nvtxs+1, 0, "QUADNODALMETIS: nptr");
    for (j=4*nelmnts, i=0; i<j; i++)
        nptr[elmnts[i]]++;
    MAKECSR(i, nvtxs, nptr);

    nind = idxmalloc(nptr[nvtxs], "QUADNODALMETIS: nind");
    for (k=i=0; i<nelmnts; i++) {
        for (j=0; j<4; j++, k++)
            nind[nptr[elmnts[k]]++] = i;
    }
    for (i=nvtxs; i>0; i--)
        nptr[i] = nptr[i-1];
    nptr[0] = 0;


    mark = idxsmalloc(nvtxs, -1, "QUADNODALMETIS: mark");

    nedges = dxadj[0] = 0;
    for (i=0; i<nvtxs; i++) {
        mark[i] = i;
        for (j=nptr[i]; j<nptr[i+1]; j++) {
            jj=4*nind[j];
            for (k=0; k<4; k++) {
                if (elmnts[jj+k] == i)
                    break;
            }
            ASSERT(k != 4);

            /* You found the index, now go and put the 2 neighbors */
            kk = elmnts[jj+table[k][0]];
            if (mark[kk] != i) {
                mark[kk] = i;
                dadjncy[nedges++] = kk;
            }
            kk = elmnts[jj+table[k][1]];
            if (mark[kk] != i) {
                mark[kk] = i;
                dadjncy[nedges++] = kk;
            }
        }
        dxadj[i+1] = nedges;
    }

    free(mark);
    free(nptr);
    free(nind);

}

/*************************************************************************
* This function setsup the CtrlType structure
**************************************************************************/
GraphType *Moc_SetUpGraph(CtrlType *ctrl, int ncon, idxtype *vtxdist, idxtype *xadj,
                          idxtype *vwgt, idxtype *adjncy, idxtype *adjwgt, int *wgtflag)
{
    int i, j;
    GraphType *graph;
    int ltvwgts[MAXNCON];

    graph = CreateGraph();
    graph->level   = 0;
    graph->gnvtxs  = vtxdist[ctrl->npes];
    graph->nvtxs   = vtxdist[ctrl->mype+1]-vtxdist[ctrl->mype];
    graph->ncon    = ncon;
    graph->nedges  = xadj[graph->nvtxs];
    graph->xadj    = xadj;
    graph->vwgt    = vwgt;
    graph->adjncy  = adjncy;
    graph->adjwgt  = adjwgt;
    graph->vtxdist = vtxdist;


    if (((*wgtflag)&2) == 0)
        graph->vwgt = idxsmalloc(graph->nvtxs*ncon, 1, "Par_KMetis: vwgt");

    if (((*wgtflag)&1) == 0)
        graph->adjwgt = idxsmalloc(graph->nedges, 1, "Par_KMetis: adjwgt");

    /* compute tvwgts */
    for (j=0; j<ncon; j++)
        ltvwgts[j] = 0;

    for (i=0; i<graph->nvtxs; i++)
        for (j=0; j<ncon; j++)
            ltvwgts[j] += graph->vwgt[i*ncon+j];

    for (j=0; j<ncon; j++)
        ctrl->tvwgts[j] = GlobalSESum(ctrl, ltvwgts[j]);

    /* check for zero wgt constraints */
    for (i=0; i<ncon; i++) {
        /* ADD: take care of the case in which tvwgts is zero */
        if (ctrl->tvwgts[i] == 0) {
            rprintf(ctrl, "ERROR: sum weight for constraint %d is zero\n", i);
            MPI_Finalize();
            exit(-1);
        }
    }

    /* compute nvwgts */
    graph->nvwgt = fmalloc(graph->nvtxs*ncon, "graph->nvwgt");
    for (i=0; i<graph->nvtxs; i++) {
        for (j=0; j<ncon; j++)
            graph->nvwgt[i*ncon+j] = (floattype)(graph->vwgt[i*ncon+j]) / (floattype)(ctrl->tvwgts[j]);
    }

    srand(ctrl->seed);

    return graph;
}

/*************************************************************************
* This function computes the normalized cut given the graph and a where vector
**************************************************************************/
float ComputeNCutVector(GraphType *graph, idxtype *where, int
npart,float* ncutVector)
{
  int i, j, cm, nvtxs;
  idxtype *ncut, *degree, *xadj, *adjncy;
  float result;
  idxtype * adjwgt;

  ncut = idxsmalloc(npart, 0, "ComputeNCut: ncut");
  degree = idxsmalloc(npart, 0, "ComputeNCut: degree");
  if ( ncutVector == NULL )
  {
  	ncutVector=(float*)malloc(sizeof(float)*npart);
  }
  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;

  if (graph->adjwgt == NULL) {
    for (i=0; i<nvtxs; i++) {
      cm = where[i];
      for (j=xadj[i]; j<xadj[i+1]; j++){
	  	if ( adjncy[j] != i )
       		degree[cm] ++; 
        if (cm != where[adjncy[j]])
          ncut[cm] ++;
      }
    }
  }
  else {
    for (i=0; i<nvtxs; i++) {
      cm = where[i];
      for (j=xadj[i]; j<xadj[i+1]; j++){
	  	if ( adjncy[j] != i )
			degree[cm] += adjwgt[j];
        if (cm != where[adjncy[j]])
          ncut[cm] += adjwgt[j];
      }
    }
  }
  int empty = 0;
  result =0;
  for (i=0; i<npart; i++){
    if (degree[i] == 0)
      empty++;
    if (degree[i] >0)
	{
	  ncutVector[i] =ncut[i] *1.0/ degree[i]; 
      result += ncutVector[i];
	}
  }
  //printf("Empty clusters: %d\n", empty);
  free(ncut);
  free(degree);
  return result+empty;
}

/*************************************************************************
* This function checks whether a graph is contigous or not
**************************************************************************/
int IsConnected(CtrlType *ctrl, GraphType *graph, int report)
{
  int i, j, k, nvtxs, first, last;
  idxtype *xadj, *adjncy, *touched, *queue;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  adjncy = graph->adjncy;

  touched = idxsmalloc(nvtxs, 0, "IsConnected: touched");
  queue = idxmalloc(nvtxs, "IsConnected: queue");

  touched[0] = 1;
  queue[0] = 0;
  first = 0; last = 1;

  while (first < last) {
    i = queue[first++];
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      if (!touched[k]) {
        queue[last++] = k;
        touched[k] = 1;
      }
    }
  }

  if (first != nvtxs && report)
    printf("The graph is not connected. It has %d disconnected vertices!\n", nvtxs-first);

  return (first == nvtxs ? 1 : 0);
}

/*************************************************************************
* This function computes the balance of the partitioning
**************************************************************************/
void ComputePartitionBalance(GraphType *graph, int nparts, idxtype *where, float *ubvec)
{
  int i, j, nvtxs, ncon;
  idxtype *kpwgts, *vwgt;
  /*float balance;*/

  nvtxs = graph->nvtxs;
  ncon = graph->ncon;
  vwgt = graph->vwgt;

  kpwgts = idxsmalloc(nparts, 0, "ComputePartitionInfo: kpwgts");

  if (vwgt == NULL) {
    for (i=0; i<nvtxs; i++)
      kpwgts[where[i]]++;
    ubvec[0] = 1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*nvtxs);
  }
  else {
    for (j=0; j<ncon; j++) {
      idxset(nparts, 0, kpwgts);
      for (i=0; i<graph->nvtxs; i++)
        kpwgts[where[i]] += vwgt[i*ncon+j];

      ubvec[j] = 1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts));
    }
  }

  free(kpwgts);

}

/*************************************************************************
* This function returns the number of connected components in cptr,cind
* The separator of the graph is used to split it and then find its components.
**************************************************************************/
int FindComponents(CtrlType *ctrl, GraphType *graph, idxtype *cptr, idxtype *cind)
{
  int i, j, k, nvtxs, first, last, nleft, ncmps, wgt;
  idxtype *xadj, *adjncy, *where, *touched, *queue;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  adjncy = graph->adjncy;
  where = graph->where;

  touched = idxsmalloc(nvtxs, 0, "IsConnected: queue");

  for (i=0; i<graph->nbnd; i++)
    touched[graph->bndind[i]] = 1;

  queue = cind;

  nleft = 0;
  for (i=0; i<nvtxs; i++) {
    if (where[i] != 2) 
      nleft++;
  }

  for (i=0; i<nvtxs; i++) {
    if (where[i] != 2)
      break;
  }

  touched[i] = 1;
  queue[0] = i;
  first = 0; last = 1;

  cptr[0] = 0;  /* This actually points to queue */
  ncmps = 0;
  while (first != nleft) {
    if (first == last) { /* Find another starting vertex */
      cptr[++ncmps] = first;
      for (i=0; i<nvtxs; i++) {
        if (!touched[i])
          break;
      }
      queue[last++] = i;
      touched[i] = 1;
    }

    i = queue[first++];
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      if (!touched[k]) {
        queue[last++] = k;
        touched[k] = 1;
      }
    }
  }
  cptr[++ncmps] = first;

  free(touched);

  return ncmps;
}

/*************************************************************************
* This function computes the cut given the graph and a where vector
**************************************************************************/
idxtype ComputeMaxCut(GraphType *graph, idxtype nparts, idxtype *where)
{
  idxtype i, j, maxcut;
  idxtype *cuts;

  cuts = idxsmalloc(nparts, 0, "ComputeMaxCut: cuts");

  if (graph->adjwgt == NULL) {
    for (i=0; i<graph->nvtxs; i++) {
      for (j=graph->xadj[i]; j<graph->xadj[i+1]; j++)
        if (where[i] != where[graph->adjncy[j]]) 
          cuts[where[i]]++;
    }
  }
  else {
    for (i=0; i<graph->nvtxs; i++) {
      for (j=graph->xadj[i]; j<graph->xadj[i+1]; j++)
        if (where[i] != where[graph->adjncy[j]])
          cuts[where[i]] += graph->adjwgt[j];
    }
  }

  maxcut = cuts[idxargmax(nparts, cuts)];

  mprintf("%D => %D\n", idxargmax(nparts, cuts), maxcut);

  gk_free((void **)&cuts, LTERM);

  return maxcut;
}

/*************************************************************************
* This function checks whether or not partition pid is contigous
**************************************************************************/
int IsConnected2(GraphType *graph, int report)
{
  int i, j, k, nvtxs, first, last, nleft, ncmps, wgt;
  idxtype *xadj, *adjncy, *where, *touched, *queue;
  idxtype *cptr;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  adjncy = graph->adjncy;
  where = graph->where;

  touched = idxsmalloc(nvtxs, 0, "IsConnected: touched");
  queue = idxmalloc(nvtxs, "IsConnected: queue");
  cptr = idxmalloc(nvtxs, "IsConnected: cptr");

  nleft = nvtxs;
  touched[0] = 1;
  queue[0] = 0;
  first = 0; last = 1;

  cptr[0] = 0;  /* This actually points to queue */
  ncmps = 0;
  while (first != nleft) {
    if (first == last) { /* Find another starting vertex */
      cptr[++ncmps] = first;
      for (i=0; i<nvtxs; i++) {
        if (!touched[i])
          break;
      }
      queue[last++] = i;
      touched[i] = 1;
    }

    i = queue[first++];
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      if (!touched[k]) {
        queue[last++] = k;
        touched[k] = 1;
      }
    }
  }
  cptr[++ncmps] = first;

  if (ncmps > 1 && report) {
    printf("%d connected components:\t", ncmps);
    for (i=0; i<ncmps; i++) {
      if (cptr[i+1]-cptr[i] > 200)
        printf("[%5d] ", cptr[i+1]-cptr[i]);
    }
    printf("\n");
  }

  GKfree(&touched, &queue, &cptr, LTERM);

  return (ncmps == 1 ? 1 : 0);
}

/*************************************************************************
* This function finds a matching using the HEM heuristic
**************************************************************************/
void Match_RM(CtrlType *ctrl, GraphType *graph)
{
  int i, ii, j, k, nvtxs, cnvtxs, maxidx;
  idxtype *xadj, *vwgt, *adjncy;
  idxtype *match, *cmap, *perm;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  vwgt = graph->vwgt;
  adjncy = graph->adjncy;

  cmap = graph->cmap = idxsmalloc(nvtxs, -1, "graph->cmap");
  match = idxsmalloc(nvtxs, -1, "match");
  perm = idxmalloc(nvtxs, "perm");

  RandomPermute(nvtxs, perm, 1);

  cnvtxs = 0;
  for (ii=0; ii<nvtxs; ii++) {
     i = perm[ii];

     if (match[i] == UNMATCHED) {
       maxidx = i;

       /* Find a random matching, subject to maxvwgt constraints */
       for (j=xadj[i]; j<xadj[i+1]; j++) {
          k = adjncy[j];
          if (match[k] == UNMATCHED && vwgt[i]+vwgt[k] <= ctrl->maxsize) {
            maxidx = k;
            break;
          }
       }

       cmap[i] = cmap[maxidx] = cnvtxs++;
       match[i] = maxidx;
       match[maxidx] = i;
     }
  }

  CreateCoarseGraph(graph, cnvtxs, match, perm);

  IMfree((void**)&match, &perm, LTERM);
}

/*************************************************************************
* This function computes the ratio assoc. given the graph and a where vector
**************************************************************************/
float ComputeRAsso(GraphType *graph, idxtype *where, int npart)
{
  int i, j, cm, nvtxs;
  idxtype *rasso, *clusterSize, *xadj, *adjncy;
  float result;
  idxtype * adjwgt;

  rasso = idxsmalloc(npart, 0, "ComputeNCut: ncut");
  clusterSize = idxsmalloc(npart, 0, "ComputeNCut: degree");
  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;

  for (i=0; i<nvtxs; i++)
    clusterSize[where[i]] ++;
  
  if (graph->adjwgt == NULL) {
    for (i=0; i<nvtxs; i++) {
      cm = where[i];
      for (j=xadj[i]; j<xadj[i+1]; j++)
	if (cm == where[adjncy[j]])
	  rasso[where[adjncy[j]]] ++;
    }
  }
  else {
    for (i=0; i<nvtxs; i++){
      cm = where[i];
      for (j=xadj[i]; j<xadj[i+1]; j++)
	if (cm == where[adjncy[j]])
	  rasso[where[adjncy[j]]] += adjwgt[j];
    }
  }
    
  result =0;
  for (i=0; i<npart; i++){
    if (clusterSize[i] >0)
      result +=  rasso[i] *1.0/ clusterSize[i];
  }
  free(rasso);
  free(clusterSize);
  return result;
}

/******************************************************************************
* This function takes a partition vector that is distributed and reads in
* the original graph and computes the edgecut
*******************************************************************************/
int ComputeRealCut2(idxtype *vtxdist, idxtype *mvtxdist, idxtype *part, idxtype *mpart, char *filename, MPI_Comm comm)
{
  int i, j, nvtxs, mype, npes, cut;
  idxtype *xadj, *adjncy, *gpart, *gmpart, *perm, *sizes;
  MPI_Status status;


  MPI_Comm_size(comm, &npes);
  MPI_Comm_rank(comm, &mype);

  if (mype != 0) {
    MPI_Send((void *)part, vtxdist[mype+1]-vtxdist[mype], IDX_DATATYPE, 0, 1, comm);
    MPI_Send((void *)mpart, mvtxdist[mype+1]-mvtxdist[mype], IDX_DATATYPE, 0, 1, comm);
  }
  else {  /* Processor 0 does all the rest */
    gpart = idxmalloc(vtxdist[npes], "ComputeRealCut: gpart");
    idxcopy(vtxdist[1], part, gpart);
    gmpart = idxmalloc(mvtxdist[npes], "ComputeRealCut: gmpart");
    idxcopy(mvtxdist[1], mpart, gmpart);

    for (i=1; i<npes; i++) {
      MPI_Recv((void *)(gpart+vtxdist[i]), vtxdist[i+1]-vtxdist[i], IDX_DATATYPE, i, 1, comm, &status);
      MPI_Recv((void *)(gmpart+mvtxdist[i]), mvtxdist[i+1]-mvtxdist[i], IDX_DATATYPE, i, 1, comm, &status);
    }

    /* OK, now go and reconstruct the permutation to go from the graph to mgraph */
    perm = idxmalloc(vtxdist[npes], "ComputeRealCut: perm");
    sizes = idxsmalloc(npes+1, 0, "ComputeRealCut: sizes");

    for (i=0; i<vtxdist[npes]; i++)
      sizes[gpart[i]]++;
    MAKECSR(i, npes, sizes);
    for (i=0; i<vtxdist[npes]; i++)
      perm[i] = sizes[gpart[i]]++;

    /* Ok, now read the graph from the file */
    ReadMetisGraph(filename, &nvtxs, &xadj, &adjncy);

    /* OK, now compute the cut */
    for (cut=0, i=0; i<nvtxs; i++) {
      for (j=xadj[i]; j<xadj[i+1]; j++) {
        if (gmpart[perm[i]] != gmpart[perm[adjncy[j]]])
          cut++;
      }
    }
    cut = cut/2;

    GKfree(&gpart, &gmpart, &perm, &sizes, &xadj, &adjncy, LTERM);

    return cut;
  }

  return 0;
}

/*************************************************************************
* This function is the entry point for detecting contacts between 
* bounding boxes and surface nodes
**************************************************************************/
void METIS_FindContacts(void *raw_cinfo, idxtype *nboxes, double *boxcoords, idxtype *nparts, 
               idxtype **r_cntptr, idxtype **r_cntind)
{
  idxtype i, ncnts, tncnts, maxtncnts;
  idxtype *cntptr, *cntind, *auxcntind, *stack, *marker;
  ContactInfoType *cinfo;

  cinfo = (ContactInfoType *)raw_cinfo;

  maxtncnts = 6*(*nboxes);
  cntptr    = idxsmalloc(*nboxes+1, 0, "METIS_FindContacts: cntptr");
  cntind    = idxmalloc(maxtncnts, "METIS_FindContacts: cntind");
  auxcntind = idxmalloc(*nparts, "METIS_FindContacts: auxcntind");
  stack     = idxmalloc(cinfo->nnodes, "METIS_FindContacts: stack");
  marker    = idxsmalloc(*nparts, 0, "METIS_FindContacts: marker");
  

  /* Go through each box and determine its contacting partitions */
  for (tncnts=0, i=0; i<*nboxes; i++) {
    ncnts = FindBoxContacts(cinfo, boxcoords+i*6, stack, auxcntind, marker);

    if (ncnts == 0)
      mprintf("CSearchError: Box has no contacts!\n");
  
    if (ncnts + tncnts >= maxtncnts) {
      maxtncnts += (tncnts+ncnts)*(*nboxes-i)/i;
      if ((cntind = (idxtype *)realloc(cntind, maxtncnts*sizeof(idxtype))) == NULL)
        errexit("Realloc failed! of %d words!\n", maxtncnts);
    }
    cntptr[i] = ncnts;
    idxcopy(ncnts, auxcntind, cntind+tncnts);
    tncnts += ncnts;
  }
  MAKECSR(i, *nboxes, cntptr); 

  *r_cntptr = cntptr;
  *r_cntind = cntind;

  gk_free((void **)&auxcntind, &stack, &marker, LTERM);

}

idxtype* getDegreeHistogram(GraphType* graph, int* maxDegree, int
	logScale)
{
	int i;
	*maxDegree=0;
	int maxLogDegree;
	for ( i=0; i<graph->nvtxs; i++ )
	{
		int k;
		if ( (k=(graph->xadj[i+1] - graph->xadj[i])) > *maxDegree )
		{
			*maxDegree = k;
			maxLogDegree = getLogBin(k);
		}
	}

	idxtype* hist;
	if ( logScale > 0 )
	{
		hist = idxsmalloc(maxLogDegree+1, 0,
					"getDegreeHistogram:hist");
	}
	else
	{
		hist = idxsmalloc(*maxDegree+1, 0,
							"getDegreeHistogram:hist");
	}

	for ( i=0; i<graph->nvtxs; i++ )
	{
		int l = graph->xadj[i+1]-graph->xadj[i];
		if ( logScale > 0 )
		{
			l = getLogBin(l);
		}
		hist[l]++;
	}
	
	return hist;
}

idxtype* getWeightsHistogram(GraphType* graph, int* maxWeight, int
	logScale)
{
	int i;
	*maxWeight=0;
	int maxLogWeight;
	for ( i=0; i<graph->xadj[graph->nvtxs]; i++ )
	{
		if ( graph->adjwgt[i] > *maxWeight )
		{
			*maxWeight = graph->adjwgt[i];
			maxLogWeight = getLogBin(graph->adjwgt[i]);
		}
	}

	idxtype* hist;
	if ( logScale > 0 )
	{
		hist = idxsmalloc(maxLogWeight+1, 0,
					"getDegreeHistogram:hist");
	}
	else
	{
		hist = idxsmalloc(*maxWeight+1, 0,
							"getDegreeHistogram:hist");
	}

	for ( i=0; i<graph->xadj[graph->nvtxs]; i++ )
	{
		int l = graph->adjwgt[i];
		if ( logScale > 0 )
		{
			l = getLogBin(l);
		}
		hist[l]++;
	}
	
	return hist;
}

/***********************************************************************************
* This function is the entry point of the parallel multilevel local diffusion
* algorithm. It uses parallel undirected diffusion followed by adaptive k-way 
* refinement. This function utilizes local coarsening.
************************************************************************************/
void ParMETIS_RepartLDiffusion(idxtype *vtxdist, idxtype *xadj, idxtype *adjncy, 
       idxtype *vwgt, realtype *adjwgt, int *wgtflag, int *numflag, int *options,
       int *edgecut, idxtype *part, MPI_Comm *comm)
{
  int npes, mype;
  CtrlType ctrl;
  WorkSpaceType wspace;
  GraphType *graph;

  MPI_Comm_size(*comm, &npes);
  MPI_Comm_rank(*comm, &mype);

  if (npes == 1) { /* Take care the npes = 1 case */
    idxset(vtxdist[1], 0, part);
    *edgecut = 0;
    return;
  }

  if (*numflag == 1) 
    ChangeNumbering(vtxdist, xadj, adjncy, part, npes, mype, 1);

  SetUpCtrl(&ctrl, npes, options, *comm);
  ctrl.CoarsenTo = amin(vtxdist[npes]+1, 70*npes);

  graph = SetUpGraph(&ctrl, vtxdist, xadj, vwgt, adjncy, adjwgt, *wgtflag);
  graph->vsize = idxsmalloc(graph->nvtxs, 1, "Par_KMetis: vsize");

  PreAllocateMemory(&ctrl, graph, &wspace);

  IFSET(ctrl.dbglvl, DBG_TRACK, printf("%d ParMETIS_RepartLDiffusion about to call AdaptiveUndirected_Partition\n",mype));
  AdaptiveUndirected_Partition(&ctrl, graph, &wspace);

  IFSET(ctrl.dbglvl, DBG_TRACK, printf("%d ParMETIS_RepartLDiffusion about to call ReMapGraph\n",mype));
  ReMapGraph(&ctrl, graph, 0, &wspace);

  idxcopy(graph->nvtxs, graph->where, part);
  *edgecut = graph->mincut;

  IMfree((void**)&graph->vsize, LTERM);
  FreeInitialGraphAndRemap(graph, *wgtflag);
  FreeWSpace(&wspace);
  FreeCtrl(&ctrl);

  if (*numflag == 1)
    ChangeNumbering(vtxdist, xadj, adjncy, part, npes, mype, 0);
}

/*************************************************************************
* This function computes the size of the coarse graph
**************************************************************************/
int ComputeCoarseGraphSize(int nvtxs, idxtype *xadj, idxtype *adjncy, int cnvtxs, idxtype *cmap, idxtype *match, idxtype *perm)
{
  int i, j, k, istart, iend, nedges, cnedges, v, u;
  idxtype *htable;

  htable = idxsmalloc(cnvtxs, -1, "htable");

  cnvtxs = cnedges = 0;
  for (i=0; i<nvtxs; i++) {
    v = perm[i];
    if (cmap[v] != cnvtxs) 
      continue;

    htable[cnvtxs] = cnvtxs;

    u = match[v];

    istart = xadj[v];
    iend = xadj[v+1];
    for (j=istart; j<iend; j++) {
      k = cmap[adjncy[j]];
      if (htable[k] != cnvtxs) {
        htable[k] = cnvtxs;
        cnedges++;
      }
    }

    if (v != u) { 
      istart = xadj[u];
      iend = xadj[u+1];
      for (j=istart; j<iend; j++) {
        k = cmap[adjncy[j]];
        if (htable[k] != cnvtxs) {
          htable[k] = cnvtxs;
          cnedges++;
        }
      }
    }
    cnvtxs++;
  }

  GKfree(&htable, LTERM);

  return cnedges;
}

/*************************************************************************
* This function computes the subdomain graph
**************************************************************************/
void PrintSubDomainGraph(GraphType *graph, int nparts, idxtype *where)
{
  int i, j, k, me, nvtxs, total, max;
  idxtype *xadj, *adjncy, *adjwgt, *pmat;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;

  pmat = idxsmalloc(nparts*nparts, 0, "ComputeSubDomainGraph: pmat");

  for (i=0; i<nvtxs; i++) {
    me = where[i];
    for (j=xadj[i]; j<xadj[i+1]; j++) {
      k = adjncy[j];
      if (where[k] != me) 
        pmat[me*nparts+where[k]] += adjwgt[j];
    }
  }

  /* printf("Subdomain Info\n"); */
  total = max = 0;
  for (i=0; i<nparts; i++) {
    for (k=0, j=0; j<nparts; j++) {
      if (pmat[i*nparts+j] > 0)
        k++;
    }
    total += k;

    if (k > max)
      max = k;
/*
    printf("%2d -> %2d  ", i, k);
    for (j=0; j<nparts; j++) {
      if (pmat[i*nparts+j] > 0)
        printf("[%2d %4d] ", j, pmat[i*nparts+j]);
    }
    printf("\n");
*/
  }
  printf("Total adjacent subdomains: %d, Max: %d\n", total, max);

  free(pmat);
}

/*************************************************************************
* This function computes the balance of the element partitioning
**************************************************************************/
float ComputeElementBalance(int ne, int nparts, idxtype *where)
{
  int i;
  idxtype *kpwgts;
  float balance;

  kpwgts = idxsmalloc(nparts, 0, "ComputeElementBalance: kpwgts");

  for (i=0; i<ne; i++)
    kpwgts[where[i]]++;

  balance = 1.0*nparts*kpwgts[idxamax(nparts, kpwgts)]/(1.0*idxsum(nparts, kpwgts));

  free(kpwgts);

  return balance;

}