void METIS_NodeRefine(int nvtxs, idxtype *xadj, idxtype *vwgt, idxtype *adjncy, 
           idxtype *adjwgt, idxtype *where, idxtype *hmarker, float ubfactor)
{
  GraphType *graph;
  CtrlType ctrl;

  ctrl.dbglvl    = ONMETIS_DBGLVL;
  ctrl.optype    = OP_ONMETIS;

  graph = CreateGraph();
  SetUpGraph(graph, OP_ONMETIS, nvtxs, 1, xadj, adjncy, vwgt, adjwgt, 3);

  AllocateWorkSpace(&ctrl, graph, 2);

  Allocate2WayNodePartitionMemory(&ctrl, graph);
  idxcopy(nvtxs, where, graph->where);

  Compute2WayNodePartitionParams(&ctrl, graph);

  FM_2WayNodeRefine_OneSidedP(&ctrl, graph, hmarker, ubfactor, 10); 
  /* FM_2WayNodeRefine_TwoSidedP(&ctrl, graph, hmarker, ubfactor, 10); */

  FreeWorkSpace(&ctrl, graph);

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

  FreeGraph(graph);

}

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

Matrix* setupCanonicalMatrix(int nvtxs, int nedges, idxtype* xadj,
	idxtype* adjncy, idxtype* adjwgt, int ncutify)
{
	int i,j;
	Matrix* ret;

	if ( ncutify )
		ret=allocMatrix(nvtxs,nedges,1,0,0);
	else
		ret=allocMatrix(nvtxs,nedges,0,0,0);
		
	idxcopy(nvtxs+1, xadj, ret->xadj); 
	idxcopy(nedges, adjncy, ret->adjncy);
	if ( adjwgt != NULL )
	{
		if ( ncutify )
		{
			for(i=0;i<ret->nvtxs;i++)
			{
				ret->adjwgtsum[i]=0;
				for(j=ret->xadj[i];j<ret->xadj[i+1];j++)
				{

					ret->adjwgt[j]=(wgttype)adjwgt[j];
				
					ret->adjwgtsum[i]+=ret->adjwgt[j];
				}
			}
			//ncutifyWeights(ret,1,ncutify);  //YK removed
		}
		else
		{
			for(i=0;i<nedges;i++)
				ret->adjwgt[i]=(wgttype)adjwgt[i];
		}
		

		normalizeColumns(ret,1,0);
	}
	else
	{
		if ( ncutify )
			ncutifyWeights(ret,0,ncutify);
		normalizeColumns(ret,0,0);
	}

	// sort each column in ascending order. This is necessary for
	// getDprAdjMatrix. 
	for(i=0;i<nvtxs;i++)
	{
		ParallelQSort(ret->adjncy,ret->adjwgt,ret->xadj[i],ret->xadj[i+1]-1);
	}
	return ret;
}

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 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 re-adjusts the amount of memory that was allocated if
* it will lead to significant savings
**************************************************************************/
void ReAdjustMemory(GraphType *graph, GraphType *cgraph, int dovsize) 
{

  if (cgraph->nedges > 100000 && graph->nedges < 0.7*graph->nedges) {
    idxcopy(cgraph->nedges, cgraph->adjwgt, cgraph->adjncy+cgraph->nedges);

    if (graph->ncon == 1) {
      if (dovsize) {
        cgraph->gdata = realloc(cgraph->gdata, (5*cgraph->nvtxs+1 + 2*cgraph->nedges)*sizeof(idxtype));

        /* Do this, in case everything was copied into new space */
        cgraph->xadj 		= cgraph->gdata;
        cgraph->vwgt 		= cgraph->gdata + cgraph->nvtxs+1;
        cgraph->vsize 		= cgraph->gdata + 2*cgraph->nvtxs+1;
        cgraph->adjwgtsum	= cgraph->gdata + 3*cgraph->nvtxs+1;
        cgraph->cmap 		= cgraph->gdata + 4*cgraph->nvtxs+1;
        cgraph->adjncy 		= cgraph->gdata + 5*cgraph->nvtxs+1;
        cgraph->adjwgt 		= cgraph->gdata + 5*cgraph->nvtxs+1 + cgraph->nedges;
      }
      else {
        cgraph->gdata = realloc(cgraph->gdata, (4*cgraph->nvtxs+1 + 2*cgraph->nedges)*sizeof(idxtype));

        /* Do this, in case everything was copied into new space */
        cgraph->xadj 	= cgraph->gdata;
        cgraph->vwgt 	= cgraph->gdata + cgraph->nvtxs+1;
        cgraph->adjwgtsum	= cgraph->gdata + 2*cgraph->nvtxs+1;
        cgraph->cmap 	= cgraph->gdata + 3*cgraph->nvtxs+1;
        cgraph->adjncy 	= cgraph->gdata + 4*cgraph->nvtxs+1;
        cgraph->adjwgt 	= cgraph->gdata + 4*cgraph->nvtxs+1 + cgraph->nedges;
      }
    }
    else {
      if (dovsize) {
        cgraph->gdata = realloc(cgraph->gdata, (4*cgraph->nvtxs+1 + 2*cgraph->nedges)*sizeof(idxtype));

        /* Do this, in case everything was copied into new space */
        cgraph->xadj 		= cgraph->gdata;
        cgraph->vsize		= cgraph->gdata + cgraph->nvtxs+1;
        cgraph->adjwgtsum	= cgraph->gdata + 2*cgraph->nvtxs+1;
        cgraph->cmap 		= cgraph->gdata + 3*cgraph->nvtxs+1;
        cgraph->adjncy 		= cgraph->gdata + 4*cgraph->nvtxs+1;
        cgraph->adjwgt 		= cgraph->gdata + 4*cgraph->nvtxs+1 + cgraph->nedges;
      }
      else {
        cgraph->gdata = realloc(cgraph->gdata, (3*cgraph->nvtxs+1 + 2*cgraph->nedges)*sizeof(idxtype));

        /* Do this, in case everything was copied into new space */
        cgraph->xadj 		= cgraph->gdata;
        cgraph->adjwgtsum	= cgraph->gdata + cgraph->nvtxs+1;
        cgraph->cmap 		= cgraph->gdata + 2*cgraph->nvtxs+1;
        cgraph->adjncy 		= cgraph->gdata + 3*cgraph->nvtxs+1;
        cgraph->adjwgt 		= cgraph->gdata + 3*cgraph->nvtxs+1 + cgraph->nedges;
      }
    }
  }

}

/*************************************************************************
* This function takes a graph and produces a bisection by using a region
* growing algorithm. The resulting partition is returned in
* graph->where
**************************************************************************/
void MocGrowBisection(CtrlType *ctrl, GraphType *graph, float *tpwgts, float ubfactor)
{
  int i, j, k, nvtxs, ncon, from, bestcut, mincut, nbfs;
  idxtype *bestwhere, *where;

  nvtxs = graph->nvtxs;

  MocAllocate2WayPartitionMemory(ctrl, graph);
  where = graph->where;

  bestwhere = idxmalloc(nvtxs, "BisectGraph: bestwhere");
  nbfs = 2*(nvtxs <= ctrl->CoarsenTo ? SMALLNIPARTS : LARGENIPARTS);
  bestcut = idxsum(graph->nedges, graph->adjwgt);  

  for (; nbfs>0; nbfs--) {
    idxset(nvtxs, 1, where);
    where[RandomInRange(nvtxs)] = 0;

    MocCompute2WayPartitionParams(ctrl, graph);

    MocInit2WayBalance(ctrl, graph, tpwgts);

    MocFM_2WayEdgeRefine(ctrl, graph, tpwgts, 4); 

    MocBalance2Way(ctrl, graph, tpwgts, 1.02);
    MocFM_2WayEdgeRefine(ctrl, graph, tpwgts, 4); 

    if (bestcut >= graph->mincut) {
      bestcut = graph->mincut;
      idxcopy(nvtxs, where, bestwhere);
      if (bestcut == 0)
        break;
    }
  }

  graph->mincut = bestcut;
  idxcopy(nvtxs, bestwhere, where);

  /*GKfree(&bestwhere, LTERM);*/
  GKfree1((void**)&bestwhere);
}

/*************************************************************************
* This function takes a graph and produces a bisection by using a region
* growing algorithm. The resulting partition is returned in
* graph->where
**************************************************************************/
void MocGrowBisection2(CtrlType *ctrl, GraphType *graph, float *tpwgts, float *ubvec)
{
  int /*i, j, k,*/ nvtxs, /*ncon, from,*/ bestcut, /*mincut,*/ nbfs;
  idxtype *bestwhere, *where;

  nvtxs = graph->nvtxs;

  MocAllocate2WayPartitionMemory(ctrl, graph);
  where = graph->where;

  bestwhere = idxmalloc(nvtxs, "BisectGraph: bestwhere");
  nbfs = 2*(nvtxs <= ctrl->CoarsenTo ? SMALLNIPARTS : LARGENIPARTS);
  bestcut = idxsum(graph->nedges, graph->adjwgt);  

  for (; nbfs>0; nbfs--) {
    idxset(nvtxs, 1, where);
    where[RandomInRange(nvtxs)] = 0;

    MocCompute2WayPartitionParams(ctrl, graph);

    MocBalance2Way2(ctrl, graph, tpwgts, ubvec);

    MocFM_2WayEdgeRefine2(ctrl, graph, tpwgts, ubvec, 4); 

    MocBalance2Way2(ctrl, graph, tpwgts, ubvec);
    MocFM_2WayEdgeRefine2(ctrl, graph, tpwgts, ubvec, 4); 

    if (bestcut > graph->mincut) {
      bestcut = graph->mincut;
      idxcopy(nvtxs, where, bestwhere);
      if (bestcut == 0)
        break;
    }
  }

  graph->mincut = bestcut;
  idxcopy(nvtxs, bestwhere, where);

  GKfree((void**)&bestwhere, 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 ONWMETIS. It requires weights on the
* vertices. It is for the case that the matrix has been pre-compressed.
**************************************************************************/
void METIS_NodeComputeSeparator(idxtype *nvtxs, idxtype *xadj, idxtype *adjncy, idxtype *vwgt, 
           idxtype *adjwgt, idxtype *options, idxtype *sepsize, idxtype *part) 
{
  idxtype 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, 1);

  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  = 0;
  ctrl.pfactor = 0;
  ctrl.nseps = 3;
  ctrl.optype = OP_ONMETIS;
  ctrl.CoarsenTo = amin(100, *nvtxs-1);
  ctrl.maxvwgt = 1.5*tvwgt/ctrl.CoarsenTo;

  InitRandom(options[7]);

  AllocateWorkSpace(&ctrl, &graph, 2);

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

  MlevelNodeBisectionMultiple(&ctrl, &graph, tpwgts, 1.02);

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

  FreeGraph(&graph, 0);

  FreeWorkSpace(&ctrl, &graph);

}

/*************************************************************************
* This function takes a graph and produces a bisection by using a region
* growing algorithm. The resulting partition is returned in
* graph->where
**************************************************************************/
void MocGrowBisectionNew2(CtrlType *ctrl, GraphType *graph, float *tpwgts, float *ubvec)
{
  idxtype i, j, k, nvtxs, ncon, from, bestcut, mincut, nbfs, inbfs;
  idxtype *bestwhere, *where;

  nvtxs = graph->nvtxs;

  MocAllocate2WayPartitionMemory(ctrl, graph);
  where = graph->where;

  bestwhere = idxmalloc(nvtxs, "BisectGraph: bestwhere");
  nbfs = 2*(nvtxs <= ctrl->CoarsenTo ? SMALLNIPARTS : LARGENIPARTS);

  for (inbfs=0; inbfs<nbfs; inbfs++) {
    idxset(nvtxs, 1, where);
    where[RandomInRange(nvtxs)] = 0;

    MocCompute2WayPartitionParams(ctrl, graph);

    MocInit2WayBalance2(ctrl, graph, tpwgts, ubvec);

    MocFM_2WayEdgeRefine2(ctrl, graph, tpwgts, ubvec, 4); 

    if (inbfs == 0 || bestcut > graph->mincut) {
      bestcut = graph->mincut;
      idxcopy(nvtxs, where, bestwhere);
      if (bestcut == 0)
        break;
    }
  }

  graph->mincut = bestcut;
  idxcopy(nvtxs, bestwhere, where);

  gk_free((void **)&bestwhere, 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 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 creates the fused-element-graph and returns the partition
************************************************************************************/
void ParMETIS_FusedElementGraph(idxtype *vtxdist, idxtype *xadj, realtype *vvol,
              realtype *vsurf, idxtype *adjncy, idxtype *vwgt, realtype *adjwgt,
              int *wgtflag, int *numflag, int *nparts, int *options,
              idxtype *part, MPI_Comm *comm)
{
  int npes, mype, nvtxs;
  CtrlType ctrl;
  WorkSpaceType wspace;
  GraphType *graph;

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

  nvtxs = vtxdist[mype+1]-vtxdist[mype];

  /* IFSET(options[OPTION_DBGLVL], DBG_TRACK, printf("%d ParMETIS_FEG npes=%d\n",mype, npes)); */

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

  graph = SetUpGraph(&ctrl, vtxdist, xadj, vwgt, adjncy, adjwgt, *wgtflag);

  graph->where = part;

  PreAllocateMemory(&ctrl, graph, &wspace);

  IFSET(ctrl.dbglvl, DBG_TIME, InitTimers(&ctrl));
  IFSET(ctrl.dbglvl, DBG_TIME, MPI_Barrier(ctrl.gcomm));
  IFSET(ctrl.dbglvl, DBG_TIME, starttimer(ctrl.TotalTmr));

  CreateFusedElementGraph(&ctrl, graph, &wspace, numflag);

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

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

  if (((*wgtflag)&2) == 0)
    IMfree((void**)&graph->vwgt, LTERM);
  IMfree((void**)&graph->lperm, &graph->peind, &graph->pexadj, &graph->peadjncy,
         &graph->peadjloc, &graph->recvptr, &graph->recvind, &graph->sendptr,
         &graph->imap, &graph->sendind, &graph, LTERM);
  FreeWSpace(&wspace);
  FreeCtrl(&ctrl);
}

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

}

/******************************************************************************
* This function takes a partition vector that is distributed and reads in
* the original graph and computes the edgecut
*******************************************************************************/
int ComputeRealCut(idxtype *vtxdist, idxtype *part, char *filename, MPI_Comm comm)
{
  int i, j, nvtxs, mype, npes, cut;
  idxtype *xadj, *adjncy, *gpart;
  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);
  }
  else {  /* Processor 0 does all the rest */
    gpart = idxmalloc(vtxdist[npes], "ComputeRealCut: gpart");
    idxcopy(vtxdist[1], part, gpart);

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

    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 (gpart[i] != gpart[adjncy[j]])
          cut++;
      }
    }
    cut = cut/2;

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

    return cut;
  }
  return 0;
}

/*************************************************************************
* This function performs k-way refinement
**************************************************************************/
void Moc_KWayFM(CtrlType *ctrl, GraphType *graph, WorkSpaceType *wspace, int npasses)
{
  int h, i, ii, iii, j, k, c;
  int pass, nvtxs, nedges, ncon;
  int nmoves, nmoved, nswaps, nzgswaps;
/*  int gnswaps, gnzgswaps; */
  int me, firstvtx, lastvtx, yourlastvtx;
  int from, to = -1, oldto, oldcut, mydomain, yourdomain, imbalanced, overweight;
  int npes = ctrl->npes, mype = ctrl->mype, nparts = ctrl->nparts;
  int nlupd, nsupd, nnbrs, nchanged;
  idxtype *xadj, *ladjncy, *adjwgt, *vtxdist;
  idxtype *where, *tmp_where, *moved;
  floattype *lnpwgts, *gnpwgts, *ognpwgts, *pgnpwgts, *movewgts, *overfill;
  idxtype *update, *supdate, *rupdate, *pe_updates;
  idxtype *changed, *perm, *pperm, *htable;
  idxtype *peind, *recvptr, *sendptr;
  KeyValueType *swchanges, *rwchanges;
  RInfoType *rinfo, *myrinfo, *tmp_myrinfo, *tmp_rinfo;
  EdgeType *tmp_edegrees, *my_edegrees, *your_edegrees;
  floattype lbvec[MAXNCON], *nvwgt, *badmaxpwgt, *ubvec, *tpwgts, lbavg, ubavg;
  int *nupds_pe;

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

  /*************************/
  /* set up common aliases */
  /*************************/
  nvtxs = graph->nvtxs;
  nedges = graph->nedges;
  ncon = graph->ncon;

  vtxdist = graph->vtxdist;
  xadj = graph->xadj;
  ladjncy = graph->adjncy;
  adjwgt = graph->adjwgt;

  firstvtx = vtxdist[mype];
  lastvtx = vtxdist[mype+1];

  where   = graph->where;
  rinfo   = graph->rinfo;
  lnpwgts = graph->lnpwgts;
  gnpwgts = graph->gnpwgts;
  ubvec   = ctrl->ubvec;
  tpwgts  = ctrl->tpwgts;

  nnbrs = graph->nnbrs;
  peind = graph->peind;
  recvptr = graph->recvptr;
  sendptr = graph->sendptr;

  changed = idxmalloc(nvtxs, "KWR: changed");
  rwchanges = wspace->pairs;
  swchanges = rwchanges + recvptr[nnbrs];

  /************************************/
  /* set up important data structures */
  /************************************/
  perm = idxmalloc(nvtxs, "KWR: perm");
  pperm = idxmalloc(nparts, "KWR: pperm");

  update = idxmalloc(nvtxs, "KWR: update");
  supdate = wspace->indices;
  rupdate = supdate + recvptr[nnbrs];
  nupds_pe = imalloc(npes, "KWR: nupds_pe");
  htable = idxsmalloc(nvtxs+graph->nrecv, 0, "KWR: lhtable");
  badmaxpwgt = fmalloc(nparts*ncon, "badmaxpwgt");

  for (i=0; i<nparts; i++) {
    for (h=0; h<ncon; h++) {
      badmaxpwgt[i*ncon+h] = ubvec[h]*tpwgts[i*ncon+h];
    }
  }

  movewgts = fmalloc(nparts*ncon, "KWR: movewgts");
  ognpwgts = fmalloc(nparts*ncon, "KWR: ognpwgts");
  pgnpwgts = fmalloc(nparts*ncon, "KWR: pgnpwgts");
  overfill = fmalloc(nparts*ncon, "KWR: overfill");
  moved = idxmalloc(nvtxs, "KWR: moved");
  tmp_where = idxmalloc(nvtxs+graph->nrecv, "KWR: tmp_where");
  tmp_rinfo = (RInfoType *)GKmalloc(sizeof(RInfoType)*nvtxs, "KWR: tmp_rinfo");
  tmp_edegrees = (EdgeType *)GKmalloc(sizeof(EdgeType)*nedges, "KWR: tmp_edegrees");

  idxcopy(nvtxs+graph->nrecv, where, tmp_where);
  for (i=0; i<nvtxs; i++) {
    tmp_rinfo[i].id = rinfo[i].id;
    tmp_rinfo[i].ed = rinfo[i].ed;
    tmp_rinfo[i].ndegrees = rinfo[i].ndegrees;
    tmp_rinfo[i].degrees = tmp_edegrees+xadj[i];

    for (j=0; j<rinfo[i].ndegrees; j++) {
      tmp_rinfo[i].degrees[j].edge = rinfo[i].degrees[j].edge;
      tmp_rinfo[i].degrees[j].ewgt = rinfo[i].degrees[j].ewgt;
    }
  }

  nswaps = nzgswaps = 0;
//.........这里部分代码省略.........

/*************************************************************************
* This function projects a partition, and at the same time computes the
* parameters for refinement.
**************************************************************************/
void Project2WayPartition(CtrlType *ctrl, GraphType *graph)
{
  int i, j, k, nvtxs, nbnd, me;
  idxtype *xadj, *adjncy, *adjwgt, *adjwgtsum;
  idxtype *cmap, *where, *id, *ed, *bndptr, *bndind;
  idxtype *cwhere, *cid, *ced, *cbndptr;
  GraphType *cgraph;

  cgraph = graph->coarser;
  cwhere = cgraph->where;
  cid = cgraph->id;
  ced = cgraph->ed;
  cbndptr = cgraph->bndptr;

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

  Allocate2WayPartitionMemory(ctrl, graph);

  where = graph->where;
  id = idxset(nvtxs, 0, graph->id);
  ed = idxset(nvtxs, 0, graph->ed);
  bndptr = idxset(nvtxs, -1, graph->bndptr);
  bndind = graph->bndind;


  /* Go through and project partition and compute id/ed for the nodes */
  for (i=0; i<nvtxs; i++) {
    k = cmap[i];
    where[i] = cwhere[k];
    cmap[i] = cbndptr[k];
  }

  for (nbnd=0, i=0; i<nvtxs; i++) {
    me = where[i];

    id[i] = adjwgtsum[i];

    if (xadj[i] == xadj[i+1]) {
      bndptr[i] = nbnd;
      bndind[nbnd++] = i;
    }
    else {
      if (cmap[i] != -1) { /* If it is an interface node. Note that cmap[i] = cbndptr[cmap[i]] */
        for (j=xadj[i]; j<xadj[i+1]; j++) {
          if (me != where[adjncy[j]])
            ed[i] += adjwgt[j];
        }
        id[i] -= ed[i];

        if (ed[i] > 0 || xadj[i] == xadj[i+1]) {
          bndptr[i] = nbnd;
          bndind[nbnd++] = i;
        }
      }
    }
  }

  graph->mincut = cgraph->mincut;
  graph->nbnd = nbnd;
  idxcopy(2, cgraph->pwgts, graph->pwgts);

  FreeGraph(graph->coarser);
  graph->coarser = NULL;

}