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

/*************************************************************************
* Setup the various arrays for the splitted graph
**************************************************************************/
void SetUpSplitGraph(GraphType *graph, GraphType *sgraph, int snvtxs, int snedges)
{
  InitGraph(sgraph);
  sgraph->nvtxs = snvtxs;
  sgraph->nedges = snedges;
  sgraph->ncon = graph->ncon;

  /* Allocate memory for the splitted graph */
  if (graph->ncon == 1) {
    sgraph->gdata = idxmalloc(4*snvtxs+1 + 2*snedges, "SetUpSplitGraph: gdata");

    sgraph->xadj        = sgraph->gdata;
    sgraph->vwgt        = sgraph->gdata + snvtxs+1;
    sgraph->adjwgtsum   = sgraph->gdata + 2*snvtxs+1;
    sgraph->cmap        = sgraph->gdata + 3*snvtxs+1;
    sgraph->adjncy      = sgraph->gdata + 4*snvtxs+1;
    sgraph->adjwgt      = sgraph->gdata + 4*snvtxs+1 + snedges;
  }
  else {
    sgraph->gdata = idxmalloc(3*snvtxs+1 + 2*snedges, "SetUpSplitGraph: gdata");

    sgraph->xadj        = sgraph->gdata;
    sgraph->adjwgtsum   = sgraph->gdata + snvtxs+1;
    sgraph->cmap        = sgraph->gdata + 2*snvtxs+1;
    sgraph->adjncy      = sgraph->gdata + 3*snvtxs+1;
    sgraph->adjwgt      = sgraph->gdata + 3*snvtxs+1 + snedges;

    sgraph->nvwgt       = fmalloc(graph->ncon*snvtxs, "SetUpSplitGraph: nvwgt");
  }

  sgraph->label	= idxmalloc(snvtxs, "SetUpSplitGraph: sgraph->label");
}

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 sets up the graph from the user input
**************************************************************************/
void SetUpGraph2(GraphType *graph, int nvtxs, int ncon, idxtype *xadj, 
       idxtype *adjncy, float *nvwgt, idxtype *adjwgt)
{
  int i, j, sum;

  InitGraph(graph);

  graph->nvtxs = nvtxs;
  graph->nedges = xadj[nvtxs];
  graph->ncon = ncon;
  graph->xadj = xadj;
  graph->adjncy = adjncy;
  graph->adjwgt = adjwgt;

  graph->nvwgt = fmalloc(nvtxs*ncon, "SetUpGraph2: graph->nvwgt");
  scopy(nvtxs*ncon, nvwgt, graph->nvwgt);

  graph->gdata = idxmalloc(2*nvtxs, "SetUpGraph: gdata");

  /* Compute the initial values of the adjwgtsum */
  graph->adjwgtsum = graph->gdata;
  for (i=0; i<nvtxs; i++) {
    sum = 0;
    for (j=xadj[i]; j<xadj[i+1]; j++)
      sum += adjwgt[j];
    graph->adjwgtsum[i] = sum;
  }

  graph->cmap = graph->gdata+nvtxs;

  graph->label = idxmalloc(nvtxs, "SetUpGraph: label");
  for (i=0; i<nvtxs; i++)
    graph->label[i] = i;

}

/***********************************************************************************
* 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 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);*/
  GKfree2((void **)&cptr, (void **)&cind);

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

  /* 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);*/
      GKfree2((void**)&sgraphs[i].gdata, (void**)&sgraphs[i].label);
    }
    else {
      MlevelNestedDissectionCC(ctrl, sgraphs+i, order, ubfactor, lastvtx-rnvtxs);
    }
    rnvtxs += sgraphs[i].nvtxs;
  }

  free(sgraphs);
}

/*************************************************************************
* Setup the various arrays for the coarse graph
**************************************************************************/
GraphType *SetUpCoarseGraph(GraphType *graph, int cnvtxs, int dovsize)
{
  GraphType *cgraph;

  cgraph = CreateGraph();
  cgraph->nvtxs = cnvtxs;
  cgraph->ncon = graph->ncon;

  cgraph->finer = graph;
  graph->coarser = cgraph;


  /* Allocate memory for the coarser graph */
  if (graph->ncon == 1) {
    if (dovsize) {
      cgraph->gdata = idxmalloc(5*cnvtxs+1 + 2*graph->nedges, "SetUpCoarseGraph: gdata");
      cgraph->xadj 		= cgraph->gdata;
      cgraph->vwgt 		= cgraph->gdata + cnvtxs+1;
      cgraph->vsize 		= cgraph->gdata + 2*cnvtxs+1;
      cgraph->adjwgtsum 	= cgraph->gdata + 3*cnvtxs+1;
      cgraph->cmap 		= cgraph->gdata + 4*cnvtxs+1;
      cgraph->adjncy 		= cgraph->gdata + 5*cnvtxs+1;
      cgraph->adjwgt 		= cgraph->gdata + 5*cnvtxs+1 + graph->nedges;
    }
    else {
      cgraph->gdata = idxmalloc(4*cnvtxs+1 + 2*graph->nedges, "SetUpCoarseGraph: gdata");
      cgraph->xadj 		= cgraph->gdata;
      cgraph->vwgt 		= cgraph->gdata + cnvtxs+1;
      cgraph->adjwgtsum 	= cgraph->gdata + 2*cnvtxs+1;
      cgraph->cmap 		= cgraph->gdata + 3*cnvtxs+1;
      cgraph->adjncy 		= cgraph->gdata + 4*cnvtxs+1;
      cgraph->adjwgt 		= cgraph->gdata + 4*cnvtxs+1 + graph->nedges;
    }
  }
  else {
    if (dovsize) {
      cgraph->gdata = idxmalloc(4*cnvtxs+1 + 2*graph->nedges, "SetUpCoarseGraph: gdata");
      cgraph->xadj 		= cgraph->gdata;
      cgraph->vsize 		= cgraph->gdata + cnvtxs+1;
      cgraph->adjwgtsum 	= cgraph->gdata + 2*cnvtxs+1;
      cgraph->cmap 		= cgraph->gdata + 3*cnvtxs+1;
      cgraph->adjncy 		= cgraph->gdata + 4*cnvtxs+1;
      cgraph->adjwgt 		= cgraph->gdata + 4*cnvtxs+1 + graph->nedges;
    }
    else {
      cgraph->gdata = idxmalloc(3*cnvtxs+1 + 2*graph->nedges, "SetUpCoarseGraph: gdata");
      cgraph->xadj 		= cgraph->gdata;
      cgraph->adjwgtsum 	= cgraph->gdata + cnvtxs+1;
      cgraph->cmap 		= cgraph->gdata + 2*cnvtxs+1;
      cgraph->adjncy 		= cgraph->gdata + 3*cnvtxs+1;
      cgraph->adjwgt 		= cgraph->gdata + 3*cnvtxs+1 + graph->nedges;
    }

    cgraph->nvwgt 	= fmalloc(graph->ncon*cnvtxs, "SetUpCoarseGraph: nvwgt");
  }

  return cgraph;
}

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

/*************************************************************************
* This function initializes the data structures of the priority queue
**************************************************************************/
void PQueueInit(CtrlType *ctrl, PQueueType *queue, int maxnodes, int maxgain)
{
  int i, j, ncore;

  queue->nnodes = 0;
  queue->maxnodes = maxnodes;

  queue->buckets = NULL;
  queue->nodes = NULL;
  queue->heap = NULL;
  queue->locator = NULL;

  if (maxgain > PLUS_GAINSPAN || maxnodes < 500)
    queue->type = 2;
  else
    queue->type = 1;

  if (queue->type == 1) {
    queue->pgainspan = amin(PLUS_GAINSPAN, maxgain);
    queue->ngainspan = amin(NEG_GAINSPAN, maxgain);

    j = queue->ngainspan+queue->pgainspan+1;

    ncore = 2 + (sizeof(ListNodeType)/sizeof(idxtype))*maxnodes + (sizeof(ListNodeType *)/sizeof(idxtype))*j;

    if (WspaceAvail(ctrl) > ncore) {
      queue->nodes = (ListNodeType *)idxwspacemalloc(ctrl, (sizeof(ListNodeType)/sizeof(idxtype))*maxnodes);
      queue->buckets = (ListNodeType **)idxwspacemalloc(ctrl, (sizeof(ListNodeType *)/sizeof(idxtype))*j);
      queue->mustfree = 0;
    }
    else { /* Not enough memory in the wspace, allocate it */
      queue->nodes = (ListNodeType *)idxmalloc((sizeof(ListNodeType)/sizeof(idxtype))*maxnodes, "PQueueInit: queue->nodes");
      queue->buckets = (ListNodeType **)idxmalloc((sizeof(ListNodeType *)/sizeof(idxtype))*j, "PQueueInit: queue->buckets");
      queue->mustfree = 1;
    }

    for (i=0; i<maxnodes; i++) 
      queue->nodes[i].id = i;

    for (i=0; i<j; i++)
      queue->buckets[i] = NULL;

    queue->buckets += queue->ngainspan;  /* Advance buckets by the ngainspan proper indexing */
    queue->maxgain = -queue->ngainspan;
  }
  else {
    queue->heap = (KeyValueType *)idxwspacemalloc(ctrl, (sizeof(KeyValueType)/sizeof(idxtype))*maxnodes);
    queue->locator = idxwspacemalloc(ctrl, maxnodes);
    idxset(maxnodes, -1, queue->locator);
  }

}

/*************************************************************************
* This function allocates memory for 2-way edge refinement
**************************************************************************/
void MocAllocate2WayPartitionMemory(CtrlType *ctrl, GraphType *graph)
{
  idxtype nvtxs, ncon;

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

  graph->npwgts	= gk_fmalloc(2*ncon, "MocAllocate2WayPartitionMemory: npwgts");
  graph->where	= idxmalloc(nvtxs, "MocAllocate2WayPartitionMemory: where");
  graph->id   	= idxmalloc(nvtxs, "MocAllocate2WayPartitionMemory: id");
  graph->ed   	= idxmalloc(nvtxs, "MocAllocate2WayPartitionMemory: ed");
  graph->bndptr	= idxmalloc(nvtxs, "MocAllocate2WayPartitionMemory: bndptr");
  graph->bndind	= idxmalloc(nvtxs, "MocAllocate2WayPartitionMemory: bndind");
}

/*************************************************************************
* 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 writes out a partition vector
**************************************************************************/
void WritePVector(char *gname, idxtype *vtxdist, idxtype *part, MPI_Comm comm)
{
  int i, j, k, l, rnvtxs, npes, mype, penum;
  FILE *fpin;
  idxtype *rpart;
  char partfile[256];
  MPI_Status status;

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

  if (mype == 0) {
    sprintf(partfile, "%s.part", gname);
    if ((fpin = fopen(partfile, "w")) == NULL) 
      errexit("Failed to open file %s", partfile);

    for (i=0; i<vtxdist[1]; i++)
      fprintf(fpin, "%d\n", part[i]);

    for (penum=1; penum<npes; penum++) {
      rnvtxs = vtxdist[penum+1]-vtxdist[penum];
      rpart = idxmalloc(rnvtxs, "rpart");
      MPI_Recv((void *)rpart, rnvtxs, IDX_DATATYPE, penum, 1, comm, &status);

      for (i=0; i<rnvtxs; i++)
        fprintf(fpin, "%d\n", rpart[i]);

      free(rpart);
    }
    fclose(fpin);
  }
  else
    MPI_Send((void *)part, vtxdist[mype+1]-vtxdist[mype], IDX_DATATYPE, 0, 1, comm); 

}

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

}

int mapPartition(idxtype* part, idxtype nvtxs)
{
	int i,j,htcounter;
	idxtype* hashtable = idxmalloc(nvtxs, "hashtable for mapPartition");  //hash tabke will map to cluster id
	for(i=0;i<nvtxs;i++)
		hashtable[i]=-1;
	
	for(i=0,htcounter=0;i<nvtxs;i++)
	{		
		if ( part[i] != -1  && hashtable[part[i]] == -1 ){
			hashtable[part[i]]=htcounter++;
			//printf("map attractor %d to cluster %d\n",part[i]+1,hashtable[part[i]]);
		}
	}

	for(i=0;i<nvtxs;i++){
		if(part[i]==-1)  //singleton cluster
			part[i] = htcounter++;
		else
			part[i]=hashtable[part[i]];
	}
//	GKfree(&hashtable);
	free(hashtable);

	return htcounter;
}

/*************************************************************************
* This function reads the element node array of a  Mixed mesh with weight
**************************************************************************/
idxtype *ReadMixedMeshWgt(char *filename, idxtype *ne, idxtype *nn, 
           idxtype  *etype, idxtype *vwgt)
{
  idxtype i, j, k, esize;
  idxtype *elmnts;
  FILE *fpin;
  idxtype sizes[]={-1,3,4,8,4,2};

  fpin = gk_fopen(filename, "r", __func__);

  mfscanf(fpin, "%D", ne);

  mfscanf(fpin, "%D", nn);

  elmnts = idxmalloc(8*(*ne), "ReadMixedMeshWgt: elmnts");

  for (j=0, i=0; i<*ne; i++) {
    mfscanf(fpin, "%D",etype+i);
    mfscanf(fpin, "%D",vwgt+i);

    for (k=0;k<sizes[etype[i]];k++) {
      mfscanf(fpin, "%D", elmnts+j);
      elmnts[j++]--;
    }
  }

  gk_fclose(fpin);

  *nn = elmnts[idxargmax(j, elmnts)]+1;

  return elmnts;
}

idxtype* lookForSingletons(GraphType* graph, int* noOfSingletons)
{
	idxtype* newIds = idxmalloc(graph->nvtxs, 
							"lookForSingletons:newIds");
	int i, newIdCounter=0;

	for ( i=0; i<graph->nvtxs; i++ )
	{
		newIds[i]=newIdCounter++;
		if ( graph->xadj[i+1]==graph->xadj[i] 
		 || ( graph->xadj[i+1]==graph->xadj[i]+1 
		   && graph->adjncy[graph->xadj[i]] == i ) )
		{
			newIds[i] = -1;
			newIdCounter--;
		}
	}

	*noOfSingletons = graph->nvtxs - newIdCounter;

	if ( *noOfSingletons > 0 )
		return newIds;
	else
	{
		free(newIds);
		return NULL;
	}

}

/*************************************************************************
* 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 reads the spd matrix
**************************************************************************/
void ReadMetisGraph(char *filename, int *r_nvtxs, idxtype **r_xadj, idxtype **r_adjncy)
{
  int i, k, edge, nvtxs, nedges;
  idxtype *xadj, *adjncy;
  char *line, *oldstr, *newstr;
  FILE *fpin;

  line = (char *)malloc(sizeof(char)*(8192+1));

  if ((fpin = fopen(filename, "r")) == NULL) {
    printf("Failed to open file %s\n", filename);
    exit(0);
  }

  fgets(line, 8192, fpin);
  sscanf(line, "%d %d", &nvtxs, &nedges);
  nedges *=2;

  xadj = idxmalloc(nvtxs+1, "ReadGraph: xadj");
  adjncy = idxmalloc(nedges, "ReadGraph: adjncy");

  /* Start reading the graph file */
  for (xadj[0]=0, k=0, i=0; i<nvtxs; i++) {
    fgets(line, 8192, fpin);
    oldstr = line;
    newstr = NULL;

    for (;;) {
      edge = (int)strtol(oldstr, &newstr, 10) -1;
      oldstr = newstr;

      if (edge < 0)
        break;

      adjncy[k++] = edge;
    } 
    xadj[i+1] = k;
  }

  fclose(fpin);

  free(line);

  *r_nvtxs = nvtxs;
  *r_xadj = xadj;
  *r_adjncy = adjncy;
}

/*************************************************************************
* Setup the various arrays for the coarse graph
**************************************************************************/
GraphType *SetUpCoarseGraph(GraphType *graph, idxtype cnvtxs, idxtype dovsize)
{
  GraphType *cgraph;

  cgraph = CreateGraph();

  cgraph->nvtxs = cnvtxs;
  cgraph->ncon  = graph->ncon;

  cgraph->finer  = graph;
  graph->coarser = cgraph;


  /* Allocate memory for the coarser graph */
  cgraph->xadj       = idxmalloc(cnvtxs+1, "SetUpCoarseGraph: xadj");
  cgraph->adjwgtsum  = idxmalloc(cnvtxs,   "SetUpCoarseGraph: adjwgtsum");
  cgraph->cmap       = idxmalloc(cnvtxs,   "SetUpCoarseGraph: cmap");
  cgraph->adjncy     = idxmalloc(graph->nedges,   "SetUpCoarseGraph: adjncy");
  cgraph->adjwgt     = idxmalloc(graph->nedges,   "SetUpCoarseGraph: adjwgt");

  if (graph->ncon == 1)
    cgraph->vwgt     = idxmalloc(cnvtxs,   "SetUpCoarseGraph: vwgt");
  else
    cgraph->nvwgt    = gk_fmalloc(graph->ncon*cnvtxs, "SetUpCoarseGraph: nvwgt");

  if (dovsize)
    cgraph->vsize    = idxmalloc(cnvtxs,   "SetUpCoarseGraph: vsize");

  return cgraph;
}