/*************************************************************************
* This function takes a graph and produces a bisection of it
**************************************************************************/
int MlevelRecursiveBisection(CtrlType *ctrl, GraphType *graph, int nparts, idxtype *part, floattype *tpwgts, floattype ubfactor, int fpart)
{
  int i, j, nvtxs, cut, tvwgt, tpwgts2[2];
  idxtype *label, *where;
  GraphType lgraph, rgraph;
  floattype wsum;

  nvtxs = graph->nvtxs;
  if (nvtxs == 0) {
    printf("\t***Cannot bisect a graph with 0 vertices!\n\t***You are trying to partition a graph into too many parts!\n");
    return 0;
  }

  /* Determine the weights of the partitions */
  tvwgt = idxsum(nvtxs, graph->vwgt);
  tpwgts2[0] = tvwgt*ssum(nparts/2, tpwgts);
  tpwgts2[1] = tvwgt-tpwgts2[0];

  MlevelEdgeBisection(ctrl, graph, tpwgts2, ubfactor);
  cut = graph->mincut;

  /* printf("%5d %5d %5d [%5d %f]\n", tpwgts2[0], tpwgts2[1], cut, tvwgt, ssum(nparts/2, tpwgts));*/

  label = graph->label;
  where = graph->where;
  for (i=0; i<nvtxs; i++)
    part[label[i]] = where[i] + fpart;

  if (nparts > 2) {
    SplitGraphPart(ctrl, graph, &lgraph, &rgraph);
    /* printf("%d %d\n", lgraph.nvtxs, rgraph.nvtxs); */
  }


  /* Free the memory of the top level graph */
  GKfree(&graph->gdata, &graph->rdata, &graph->label, LTERM);

  /* Scale the fractions in the tpwgts according to the true weight */
  wsum = ssum(nparts/2, tpwgts);
  sscale(nparts/2, 1.0/wsum, tpwgts);
  sscale(nparts-nparts/2, 1.0/(1.0-wsum), tpwgts+nparts/2);
  /*
  for (i=0; i<nparts; i++)
    printf("%5.3f ", tpwgts[i]);
  printf("[%5.3f]\n", wsum);
  */

  /* Do the recursive call */
  if (nparts > 3) {
    cut += MlevelRecursiveBisection(ctrl, &lgraph, nparts/2, part, tpwgts, ubfactor, fpart);
    cut += MlevelRecursiveBisection(ctrl, &rgraph, nparts-nparts/2, part, tpwgts+nparts/2, ubfactor, fpart+nparts/2);
  }
  else if (nparts == 3) {
    cut += MlevelRecursiveBisection(ctrl, &rgraph, nparts-nparts/2, part, tpwgts+nparts/2, ubfactor, fpart+nparts/2);
    GKfree(&lgraph.gdata, &lgraph.label, LTERM);
  }

  return cut;

}

/*************************************************************************
* This function takes a graph and produces a bisection of it
**************************************************************************/
void MlevelNestedDissectionCC(CtrlType *ctrl, GraphType *graph, idxtype *order, float ubfactor, int lastvtx)
{
  int i, j, nvtxs, nbnd, tvwgt, tpwgts2[2], nsgraphs, ncmps, rnvtxs;
  idxtype *label, *bndind;
  idxtype *cptr, *cind;
  GraphType *sgraphs;

  nvtxs = graph->nvtxs;

  /* Determine the weights of the partitions */
  tvwgt = idxsum(nvtxs, graph->vwgt);
  tpwgts2[0] = tvwgt/2;
  tpwgts2[1] = tvwgt-tpwgts2[0];

  MlevelNodeBisectionMultiple(ctrl, graph, tpwgts2, ubfactor);
  IFSET(ctrl->dbglvl, DBG_SEPINFO, printf("Nvtxs: %6d, [%6d %6d %6d]\n", graph->nvtxs, graph->pwgts[0], graph->pwgts[1], graph->pwgts[2]));

  /* Order the nodes in the separator */
  nbnd = graph->nbnd;
  bndind = graph->bndind;
  label = graph->label;
  for (i=0; i<nbnd; i++) 
    order[label[bndind[i]]] = --lastvtx;

  cptr = idxmalloc(nvtxs, "MlevelNestedDissectionCC: cptr");
  cind = idxmalloc(nvtxs, "MlevelNestedDissectionCC: cind");
  ncmps = FindComponents(ctrl, graph, cptr, cind);

/*
  if (ncmps > 2)
    printf("[%5d] has %3d components\n", nvtxs, ncmps);
*/

  sgraphs = (GraphType *)GKmalloc(ncmps*sizeof(GraphType), "MlevelNestedDissectionCC: sgraphs");

  nsgraphs = SplitGraphOrderCC(ctrl, graph, sgraphs, ncmps, cptr, cind);

  GKfree(&cptr, &cind, LTERM);

  /* Free the memory of the top level graph */
  GKfree(&graph->gdata, &graph->rdata, &graph->label, LTERM);

  /* Go and process the subgraphs */
  for (rnvtxs=i=0; i<nsgraphs; i++) {
    if (sgraphs[i].adjwgt == NULL) {
      MMDOrder(ctrl, sgraphs+i, order, lastvtx-rnvtxs);
      GKfree(&sgraphs[i].gdata, &sgraphs[i].label, LTERM);
    }
    else {
      MlevelNestedDissectionCC(ctrl, sgraphs+i, order, ubfactor, lastvtx-rnvtxs);
    }
    rnvtxs += sgraphs[i].nvtxs;
  }

  free(sgraphs);
}

/*************************************************************************
* This function takes a graph and produces a bisection of it
**************************************************************************/
void MlevelNestedDissection(CtrlType *ctrl, GraphType *graph, idxtype *order, float ubfactor, int lastvtx)
{
  int i, j, nvtxs, nbnd, tvwgt, tpwgts2[2];
  idxtype *label, *bndind;
  GraphType lgraph, rgraph;

  nvtxs = graph->nvtxs;

  /* Determine the weights of the partitions */
  tvwgt = idxsum(nvtxs, graph->vwgt);
  tpwgts2[0] = tvwgt/2;
  tpwgts2[1] = tvwgt-tpwgts2[0];

  switch (ctrl->optype) {
    case OP_OEMETIS:
      MlevelEdgeBisection(ctrl, graph, tpwgts2, ubfactor);

      IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->SepTmr));
      ConstructMinCoverSeparator(ctrl, graph, ubfactor);
      IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->SepTmr));

      break;
    case OP_ONMETIS:
      MlevelNodeBisectionMultiple(ctrl, graph, tpwgts2, ubfactor);

      IFSET(ctrl->dbglvl, DBG_SEPINFO, printf("Nvtxs: %6d, [%6d %6d %6d]\n", graph->nvtxs, graph->pwgts[0], graph->pwgts[1], graph->pwgts[2]));

      break;
  }

  /* Order the nodes in the separator */
  nbnd = graph->nbnd;
  bndind = graph->bndind;
  label = graph->label;
  for (i=0; i<nbnd; i++) 
    order[label[bndind[i]]] = --lastvtx;

  SplitGraphOrder(ctrl, graph, &lgraph, &rgraph);

  /* Free the memory of the top level graph */
  GKfree(&graph->gdata, &graph->rdata, &graph->label, LTERM);

  if (rgraph.nvtxs > MMDSWITCH) 
    MlevelNestedDissection(ctrl, &rgraph, order, ubfactor, lastvtx);
  else {
    MMDOrder(ctrl, &rgraph, order, lastvtx); 
    GKfree(&rgraph.gdata, &rgraph.rdata, &rgraph.label, LTERM);
  }
  if (lgraph.nvtxs > MMDSWITCH) 
    MlevelNestedDissection(ctrl, &lgraph, order, ubfactor, lastvtx-rgraph.nvtxs);
  else {
    MMDOrder(ctrl, &lgraph, order, lastvtx-rgraph.nvtxs); 
    GKfree(&lgraph.gdata, &lgraph.rdata, &lgraph.label, LTERM);
  }
}

/*************************************************************************
* This function deallocates any memory stored in a graph
**************************************************************************/
void FreeInitialGraphAndRemap(GraphType *graph, int wgtflag, int freevsize)
{
    int i, nedges;
    idxtype *adjncy, *imap;

    nedges = graph->nedges;
    adjncy = graph->adjncy;
    imap   = graph->imap;

    if (imap != NULL) {
        for (i=0; i<nedges; i++)
            adjncy[i] = imap[adjncy[i]];  /* Apply local to global transformation */
    }

    /* Free Metis's things */
    GKfree((void **)&graph->match,
           (void **)&graph->cmap,
           (void **)&graph->lperm,
           (void **)&graph->where,
           (void **)&graph->label,
           (void **)&graph->rinfo,
           (void **)&graph->nrinfo,
           (void **)&graph->nvwgt,
           (void **)&graph->lpwgts,
           (void **)&graph->gpwgts,
           (void **)&graph->lnpwgts,
           (void **)&graph->gnpwgts,
           (void **)&graph->sepind,
           (void **)&graph->peind,
           (void **)&graph->sendptr,
           (void **)&graph->sendind,
           (void **)&graph->recvptr,
           (void **)&graph->recvind,
           (void **)&graph->imap,
           (void **)&graph->rlens,
           (void **)&graph->slens,
           (void **)&graph->rcand,
           (void **)&graph->pexadj,
           (void **)&graph->peadjncy,
           (void **)&graph->peadjloc,
           LTERM);

    if (freevsize)
        GKfree((void **)&graph->vsize, LTERM);
    if ((wgtflag&2) == 0)
        GKfree((void **)&graph->vwgt, LTERM);
    if ((wgtflag&1) == 0)
        GKfree((void **)&graph->adjwgt, LTERM);

    free(graph);
}

/*************************************************************************
* This function takes a graph and produces a bisection of it
**************************************************************************/
int MlevelKWayPartitioning(CtrlType *ctrl, GraphType *graph, int nparts, idxtype *part, float *tpwgts, float ubfactor)
{
  int i, j, nvtxs, tvwgt, tpwgts2[2];
  GraphType *cgraph;
  int wgtflag=3, numflag=0, options[10], edgecut;

  cgraph = Coarsen2Way(ctrl, graph);

  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->InitPartTmr));
  AllocateKWayPartitionMemory(ctrl, cgraph, nparts);

  options[0] = 1; 
  options[OPTION_CTYPE] = MATCH_SHEMKWAY;
  options[OPTION_ITYPE] = IPART_GGPKL;
  options[OPTION_RTYPE] = RTYPE_FM;
  options[OPTION_DBGLVL] = 0;

  METIS_WPartGraphRecursive(&cgraph->nvtxs, cgraph->xadj, cgraph->adjncy, cgraph->vwgt, 
                            cgraph->adjwgt, &wgtflag, &numflag, &nparts, tpwgts, options, 
                            &edgecut, cgraph->where);

  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->InitPartTmr));
  IFSET(ctrl->dbglvl, DBG_IPART, printf("Initial %d-way partitioning cut: %d\n", nparts, edgecut));

  IFSET(ctrl->dbglvl, DBG_KWAYPINFO, ComputePartitionInfo(cgraph, nparts, cgraph->where));

  RefineKWay(ctrl, graph, cgraph, nparts, tpwgts, ubfactor);

  idxcopy(graph->nvtxs, graph->where, part);

  GKfree(&graph->gdata, &graph->rdata, LTERM);

  return graph->mincut;

}

/***********************************************************************************
* This function is the entry point of the parallel ordering algorithm.
* This function assumes that the graph is already nice partitioned among the
* processors and then proceeds to perform recursive bisection.
************************************************************************************/
void ParMETIS_V3_PartGeom(idxtype *vtxdist, int *ndims, float *xyz, idxtype *part, MPI_Comm *comm)
{
  int i, npes, mype, nvtxs, firstvtx, dbglvl;
  idxtype *xadj, *adjncy;
  CtrlType ctrl;
  WorkSpaceType wspace;
  GraphType *graph;
  int zeroflg = 0;

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

  if (npes == 1) {
    idxset(vtxdist[mype+1]-vtxdist[mype], 0, part);
    return;
  }

  /* Setup a fake graph to allow the rest of the code to work unchanged */
  dbglvl = 0;

  nvtxs = vtxdist[mype+1]-vtxdist[mype];
  firstvtx = vtxdist[mype];
  xadj = idxmalloc(nvtxs+1, "ParMETIS_PartGeom: xadj");
  adjncy = idxmalloc(nvtxs, "ParMETIS_PartGeom: adjncy");
  for (i=0; i<nvtxs; i++) {
    xadj[i] = i;
    adjncy[i] = firstvtx + (i+1)%nvtxs;
  }
  xadj[nvtxs] = nvtxs;

  /* Proceed with the rest of the code */
  SetUpCtrl(&ctrl, npes, dbglvl, *comm);
  ctrl.seed      = mype;
  ctrl.CoarsenTo = amin(vtxdist[npes]+1, 25*npes);

  graph = Moc_SetUpGraph(&ctrl, 1, vtxdist, xadj, NULL, adjncy, NULL, &zeroflg);

  PreAllocateMemory(&ctrl, graph, &wspace);

  /*=======================================================
   * Compute the initial geometric partitioning
   =======================================================*/
  IFSET(ctrl.dbglvl, DBG_TIME, InitTimers(&ctrl));
  IFSET(ctrl.dbglvl, DBG_TIME, MPI_Barrier(ctrl.gcomm));
  IFSET(ctrl.dbglvl, DBG_TIME, starttimer(ctrl.TotalTmr));

  Coordinate_Partition(&ctrl, graph, *ndims, xyz, 0, &wspace);

  idxcopy(graph->nvtxs, graph->where, part);

  IFSET(ctrl.dbglvl, DBG_TIME, MPI_Barrier(ctrl.gcomm));
  IFSET(ctrl.dbglvl, DBG_TIME, stoptimer(ctrl.TotalTmr));
  IFSET(ctrl.dbglvl, DBG_TIME, PrintTimingInfo(&ctrl));

  FreeInitialGraphAndRemap(graph, 0);
  FreeWSpace(&wspace);
  FreeCtrl(&ctrl);

  GKfree((void **)&xadj, (void **)&adjncy, 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);
}

void AllocateNodePartitionParams(CtrlType *ctrl, GraphType *graph, WorkSpaceType *wspace)
{
  int nparts, nvtxs;
  idxtype *vwgt;
  NRInfoType *rinfo, *myrinfo;

  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->KWayInitTmr));

  nvtxs  = graph->nvtxs;
  nparts = ctrl->nparts;

  graph->nrinfo  = (NRInfoType *)GKmalloc(sizeof(NRInfoType)*nvtxs, "AllocateNodePartitionParams: rinfo");
  graph->lpwgts  = idxmalloc(2*nparts, "AllocateNodePartitionParams: lpwgts");
  graph->gpwgts  = idxmalloc(2*nparts, "AllocateNodePartitionParams: gpwgts");
  graph->sepind  = idxmalloc(nvtxs, "AllocateNodePartitionParams: sepind");
  graph->hmarker = idxmalloc(nvtxs, "AllocateNodePartitionParams: hmarker");

  /* Allocate additional memory for graph->vwgt in order to store the weights
     of the remote vertices */
  vwgt        = graph->vwgt;
  graph->vwgt = idxmalloc(nvtxs+graph->nrecv, "AllocateNodePartitionParams: graph->vwgt");
  idxcopy(nvtxs, vwgt, graph->vwgt);
  GKfree((void **)&vwgt, LTERM);

  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->KWayInitTmr));
}

/*****************************************************************************
*  This function computes a partitioning using coordinate data.
*****************************************************************************/
void ParMETIS_PartGeomKway(idxtype *vtxdist, idxtype *xadj, idxtype *adjncy, idxtype *vwgt,
       idxtype *adjwgt, int *wgtflag, int *numflag, int *ndims, float *xyz, int *nparts,
       int *options, int *edgecut, idxtype *part, MPI_Comm *comm)
{
  int i;
  int ncon = 1;
  float *tpwgts, ubvec[MAXNCON];
  int myoptions[10];

  tpwgts = fmalloc(*nparts*ncon, "tpwgts");
  for (i=0; i<*nparts*ncon; i++)
    tpwgts[i] = 1.0/(float)(*nparts);
  for (i=0; i<ncon; i++)
    ubvec[i] = UNBALANCE_FRACTION;

  if (options[0] == 0) {
    myoptions[0] = 0;
  }
  else {
    myoptions[0] = 1;
    myoptions[PMV3_OPTION_DBGLVL] = options[OPTION_DBGLVL];
    myoptions[PMV3_OPTION_SEED] = GLOBAL_SEED;
  }

  ParMETIS_V3_PartGeomKway(vtxdist, xadj, adjncy, vwgt, adjwgt, wgtflag, numflag, ndims, xyz,
  &ncon, nparts, tpwgts, ubvec, myoptions, edgecut, part, comm);

  GKfree((void **)&tpwgts, LTERM);
  return;
}

/*****************************************************************************
*  This function computes a repartitioning by LMSR scratch-remap.
*****************************************************************************/
void ParMETIS_RepartMLRemap(idxtype *vtxdist, idxtype *xadj, idxtype *adjncy,
       idxtype *vwgt, idxtype *adjwgt, int *wgtflag, int *numflag, int *options,
       int *edgecut, idxtype *part, MPI_Comm *comm)
{
  int i;
  int nparts;
  int ncon = 1;
  float *tpwgts, ubvec[MAXNCON];
  float ipc_factor = 1000.0;
  int myoptions[10];

  MPI_Comm_size(*comm, &nparts);
  tpwgts = fmalloc(nparts*ncon, "tpwgts");
  for (i=0; i<nparts*ncon; i++)
    tpwgts[i] = 1.0/(float)(nparts);
  for (i=0; i<ncon; i++)
    ubvec[i] = UNBALANCE_FRACTION;

  if (options[0] == 0) {
    myoptions[0] = 0;
  }
  else {
    myoptions[0] = 1;
    myoptions[PMV3_OPTION_DBGLVL] = options[OPTION_DBGLVL];
    myoptions[PMV3_OPTION_SEED]   = GLOBAL_SEED;
    myoptions[PMV3_OPTION_PSR]    = PARMETIS_PSR_COUPLED;
  }

  ParMETIS_V3_AdaptiveRepart(vtxdist, xadj, adjncy, vwgt, NULL, adjwgt, wgtflag, numflag,
  &ncon, &nparts, tpwgts, ubvec, &ipc_factor, myoptions, edgecut, part, comm);

  GKfree((void **)&tpwgts, LTERM);
}

/*************************************************************************
* This function takes a bisection and constructs a minimum weight vertex 
* separator out of it. It uses the node-based separator refinement for it.
**************************************************************************/
void ConstructSeparator(CtrlType *ctrl, GraphType *graph, float ubfactor)
{
  int i, j, k, nvtxs, nbnd;
  idxtype *xadj, *where, *bndind;

  nvtxs = graph->nvtxs;
  xadj = graph->xadj;
  nbnd = graph->nbnd;
  bndind = graph->bndind;

  where = idxcopy(nvtxs, graph->where, idxwspacemalloc(ctrl, nvtxs));

  /* Put the nodes in the boundary into the separator */
  for (i=0; i<nbnd; i++) {
    j = bndind[i];
    if (xadj[j+1]-xadj[j] > 0)  /* Ignore islands */
      where[j] = 2;
  }

  GKfree(&graph->rdata, LTERM);
  Allocate2WayNodePartitionMemory(ctrl, graph);
  idxcopy(nvtxs, where, graph->where);
  idxwspacefree(ctrl, nvtxs);

  ASSERT(IsSeparable(graph));

  Compute2WayNodePartitionParams(ctrl, graph);

  ASSERT(CheckNodePartitionParams(graph));

  FM_2WayNodeRefine(ctrl, graph, ubfactor, 8); 

  ASSERT(IsSeparable(graph));
}

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

void FreeWSpace(WorkSpaceType *wspace)
{

    GKfree((void **)&wspace->core,
           (void **)&wspace->pv1,
           (void **)&wspace->pv2,
           (void **)&wspace->pv3,
           (void **)&wspace->pv4,
           (void **)&wspace->pepairs1,
           (void **)&wspace->pepairs2,
           LTERM);
}

/*************************************************************************
* This function is the entry point for ONWMETIS. It requires weights on the
* vertices. It is for the case that the matrix has been pre-compressed.
**************************************************************************/
void METIS_NodeComputeSeparator(int *nvtxs, idxtype *xadj, idxtype *adjncy, idxtype *vwgt, 
           idxtype *adjwgt, float *ubfactor, int *options, int *sepsize, idxtype *part) 
{
  int i, j, tvwgt, tpwgts[2];
  GraphType graph;
  CtrlType ctrl;

  SetUpGraph(&graph, OP_ONMETIS, *nvtxs, 1, xadj, adjncy, vwgt, adjwgt, 3);
  tvwgt = idxsum(*nvtxs, graph.vwgt);

  if (options[0] == 0) {  /* Use the default parameters */
    ctrl.CType  = ONMETIS_CTYPE;
    ctrl.IType  = ONMETIS_ITYPE;
    ctrl.RType  = ONMETIS_RTYPE;
    ctrl.dbglvl = ONMETIS_DBGLVL;
  }
  else {
    ctrl.CType  = options[OPTION_CTYPE];
    ctrl.IType  = options[OPTION_ITYPE];
    ctrl.RType  = options[OPTION_RTYPE];
    ctrl.dbglvl = options[OPTION_DBGLVL];
  }

  ctrl.oflags    = OFLAG_COMPRESS; /* For by-passing the pre-coarsening for multiple runs */
  ctrl.RType     = 2;  /* Standard 1-sided node refinement code */
  ctrl.pfactor   = 0;
  ctrl.nseps     = 5;  /* This should match NUM_INIT_MSECTIONS in ParMETISLib/defs.h */
  ctrl.optype    = OP_ONMETIS;

  InitRandom(options[7]);

  AllocateWorkSpace(&ctrl, &graph, 2);

  /*============================================================
   * Perform the bisection
   *============================================================*/ 
  tpwgts[0] = tvwgt/2;
  tpwgts[1] = tvwgt-tpwgts[0];

  MlevelNodeBisectionMultiple(&ctrl, &graph, tpwgts, *ubfactor*.95);

  *sepsize = graph.pwgts[2];
  idxcopy(*nvtxs, graph.where, part);

  GKfree((void **)&graph.gdata, &graph.rdata, &graph.label, LTERM);


  FreeWorkSpace(&ctrl, &graph);

}

/*************************************************************************
* This function is the entry point for KMETIS with seed specification
* in options[7] 
**************************************************************************/
void METIS_PartGraphKway2(int *nvtxs, idxtype *xadj, idxtype *adjncy, idxtype *vwgt, 
                         idxtype *adjwgt, int *wgtflag, int *numflag, int *nparts, 
                         int *options, int *edgecut, idxtype *part)
{
  int i;
  float *tpwgts;

  tpwgts = fmalloc(*nparts, "KMETIS: tpwgts");
  for (i=0; i<*nparts; i++) 
    tpwgts[i] = 1.0/(1.0*(*nparts));

  METIS_WPartGraphKway2(nvtxs, xadj, adjncy, vwgt, adjwgt, wgtflag, numflag, nparts, 
                       tpwgts, options, edgecut, part);
  GKfree((void **)&tpwgts, LTERM);
}

/*************************************************************************
* Let the game begin
**************************************************************************/
main(int argc, char *argv[])
{
  int i, j, ne, nn, etype, numflag=0;
  idxtype *elmnts, *xadj, *adjncy;
  timer IOTmr, DUALTmr;
  char fileout[256], etypestr[4][5] = {"TRI", "TET", "HEX", "QUAD"};

  if (argc != 2) {
    printf("Usage: %s <meshfile>\n",argv[0]);
    exit(0);
  }

  cleartimer(IOTmr);
  cleartimer(DUALTmr);

  starttimer(IOTmr);
  elmnts = ReadMesh(argv[1], &ne, &nn, &etype);
  stoptimer(IOTmr);

  printf("**********************************************************************\n");
  printf("%s", METISTITLE);
  printf("Mesh Information ----------------------------------------------------\n");
  printf("  Name: %s, #Elements: %d, #Nodes: %d, Etype: %s\n\n", argv[1], ne, nn, etypestr[etype-1]);
  printf("Forming Dual Graph... -----------------------------------------------\n");

  xadj = idxmalloc(ne+1, "main: xadj");
  adjncy = idxmalloc(10*ne, "main: adjncy");

  starttimer(DUALTmr);
  METIS_MeshToDual(&ne, &nn, elmnts, &etype, &numflag, xadj, adjncy);
  stoptimer(DUALTmr);

  printf("  Dual Information: #Vertices: %d, #Edges: %d\n", ne, xadj[ne]/2);

  sprintf(fileout, "%s.dgraph", argv[1]);
  starttimer(IOTmr);
  WriteGraph(fileout, ne, xadj, adjncy);
  stoptimer(IOTmr);


  printf("\nTiming Information --------------------------------------------------\n");
  printf("  I/O:          \t\t %7.3f\n", gettimer(IOTmr));
  printf("  Dual Creation:\t\t %7.3f\n", gettimer(DUALTmr));
  printf("**********************************************************************\n");

  GKfree(&elmnts, &xadj, &adjncy, LTERM);

}

/*************************************************************************
* This function takes a graph and produces a bisection of it
**************************************************************************/
int MCMlevelKWayPartitioning(CtrlType *ctrl, GraphType *graph, int nparts, idxtype *part, 
      float *rubvec)
{
  int i, j, nvtxs;
  GraphType *cgraph;
  int options[10], edgecut;

  cgraph = MCCoarsen2Way(ctrl, graph);

  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->InitPartTmr));
  MocAllocateKWayPartitionMemory(ctrl, cgraph, nparts);

  options[0] = 1; 
  options[OPTION_CTYPE] = MATCH_SBHEM_INFNORM;
  options[OPTION_ITYPE] = IPART_RANDOM;
  options[OPTION_RTYPE] = RTYPE_FM;
  options[OPTION_DBGLVL] = 0;

  /* Determine what you will use as the initial partitioner, based on tolerances */
  for (i=0; i<graph->ncon; i++) {
    if (rubvec[i] > 1.2)
      break;
  }
  if (i == graph->ncon)
    METIS_mCPartGraphRecursiveInternal(&cgraph->nvtxs, &cgraph->ncon, 
          cgraph->xadj, cgraph->adjncy, cgraph->nvwgt, cgraph->adjwgt, &nparts, 
          options, &edgecut, cgraph->where);
  else
    METIS_mCHPartGraphRecursiveInternal(&cgraph->nvtxs, &cgraph->ncon, 
          cgraph->xadj, cgraph->adjncy, cgraph->nvwgt, cgraph->adjwgt, &nparts, 
          rubvec, options, &edgecut, cgraph->where);


  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->InitPartTmr));
  IFSET(ctrl->dbglvl, DBG_IPART, printf("Initial %d-way partitioning cut: %d\n", nparts, edgecut));

  IFSET(ctrl->dbglvl, DBG_KWAYPINFO, ComputePartitionInfo(cgraph, nparts, cgraph->where));

  MocRefineKWayHorizontal(ctrl, graph, cgraph, nparts, rubvec);

  idxcopy(graph->nvtxs, graph->where, part);

  GKfree(&graph->nvwgt, &graph->npwgts, &graph->gdata, &graph->rdata, LTERM);

  return graph->mincut;

}

/*************************************************************************
* 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 actually undoes the permutation
Matrix* permuteRowsAndColumns(const Matrix *M, const idxtype*
rowPerm, const idxtype* colPerm)
{
	Matrix *ret = allocMatrix(M->nvtxs, M->nnz, 0, 0, 0);
	int i;
	idxtype* revRowPerm = idxmalloc(M->nvtxs,
						"permuteRowsAndColumns:revRowPerm");
//	idxtype* revColPerm = idxmalloc(M->nvtxs, 
//						"permuteRowsAndColumns:revColPerm");

	for ( i=0; i<M->nvtxs; i++ )
	{
		revRowPerm[rowPerm[i]] = i;
//		revColPerm[colPerm[i]] = i;
	}

	ret->xadj[0]=0;
	for ( i=0; i<M->nvtxs; i++ )
	{
		int orgI = revRowPerm[i];
		int j;
		ret->xadj[i+1] = ret->xadj[i] + ( M->xadj[orgI+1] -
								M->xadj[orgI] );
		for ( j=ret->xadj[i]; j<ret->xadj[i+1]; j++ )
		{
			int orgJ = M->xadj[orgI] + j - ret->xadj[i];
			ret->adjncy[j] = colPerm[M->adjncy[orgJ]];
			ret->adjwgt[j] = M->adjwgt[orgJ];
		}
								
		ParallelQSort( ret->adjncy, ret->adjwgt, ret->xadj[i],
		ret->xadj[i+1]-1 );						
	}

	ret->nnz = M->nnz;

//	GKfree ( (void**)&revRowPerm, (void**)&revColPerm, LTERM );
	GKfree ( (void**)&revRowPerm, LTERM );

	return ret;
	
}