diff options
Diffstat (limited to 'benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/ParMETISLib/coarsen.c')
| -rw-r--r-- | benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/ParMETISLib/coarsen.c | 485 |
1 files changed, 485 insertions, 0 deletions
diff --git a/benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/ParMETISLib/coarsen.c b/benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/ParMETISLib/coarsen.c new file mode 100644 index 0000000..70f48c2 --- /dev/null +++ b/benchmarks/CUDA/DG/3rdParty/ParMetis-3.1/ParMETISLib/coarsen.c @@ -0,0 +1,485 @@ +/* + * Copyright 1997, Regents of the University of Minnesota + * + * mcoarsen.c + * + * This file contains code that performs graph coarsening + * + * Started 2/22/96 + * George + * + * $Id: coarsen.c,v 1.2 2003/07/21 17:18:48 karypis Exp $ + * + */ + +#include <parmetislib.h> + + +/************************************************************************* +* This function creates the coarser graph +**************************************************************************/ +void Moc_Global_CreateCoarseGraph(CtrlType *ctrl, GraphType *graph, + WorkSpaceType *wspace, int cnvtxs) +{ + int h, i, j, k, l, ii, jj, ll, nnbrs, nvtxs, nedges, ncon; + int firstvtx, lastvtx, cfirstvtx, clastvtx, otherlastvtx; + int npes=ctrl->npes, mype=ctrl->mype; + int cnedges, nsend, nrecv, nkeepsize, nrecvsize, nsendsize, v, u; + idxtype *xadj, *ladjncy, *adjwgt, *vwgt, *vsize, *vtxdist, *home; + idxtype *match, *cmap, *rcmap, *scmap; + idxtype *cxadj, *cadjncy, *cadjwgt, *cvwgt, *cvsize = NULL, *chome = NULL, *cvtxdist; + idxtype *rsizes, *ssizes, *rlens, *slens, *rgraph, *sgraph, *perm; + idxtype *peind, *recvptr, *recvind; + float *nvwgt, *cnvwgt; + GraphType *cgraph; + KeyValueType *scand, *rcand; + int mask=(1<<13)-1, htable[8192], htableidx[8192]; + + nvtxs = graph->nvtxs; + ncon = graph->ncon; + + vtxdist = graph->vtxdist; + xadj = graph->xadj; + vwgt = graph->vwgt; + vsize = graph->vsize; + nvwgt = graph->nvwgt; + home = graph->home; + ladjncy = graph->adjncy; + adjwgt = graph->adjwgt; + + match = graph->match; + + firstvtx = vtxdist[mype]; + lastvtx = vtxdist[mype+1]; + + cmap = graph->cmap = idxmalloc(nvtxs+graph->nrecv, "CreateCoarseGraph: cmap"); + + nnbrs = graph->nnbrs; + peind = graph->peind; + recvind = graph->recvind; + recvptr = graph->recvptr; + + /* Use wspace->indices as the tmp space for map of the boundary + * vertices that are sent and received */ + scmap = wspace->indices; + rcmap = cmap + nvtxs; + + + /* Initialize the coarser graph */ + cgraph = CreateGraph(); + cgraph->nvtxs = cnvtxs; + cgraph->ncon = ncon; + cgraph->level = graph->level+1; + cgraph->finer = graph; + graph->coarser = cgraph; + + + + /************************************************************* + * Obtain the vtxdist of the coarser graph + **************************************************************/ + cvtxdist = cgraph->vtxdist = idxmalloc(npes+1, "CreateCoarseGraph: cvtxdist"); + cvtxdist[npes] = cnvtxs; /* Use last position in the cvtxdist as a temp buffer */ + + MPI_Allgather((void *)(cvtxdist+npes), 1, IDX_DATATYPE, (void *)cvtxdist, 1, IDX_DATATYPE, ctrl->comm); + + MAKECSR(i, npes, cvtxdist); + + cgraph->gnvtxs = cvtxdist[npes]; + +#ifdef DEBUG_CONTRACT + PrintVector(ctrl, npes+1, 0, cvtxdist, "cvtxdist"); +#endif + + + /************************************************************* + * Construct the cmap vector + **************************************************************/ + cfirstvtx = cvtxdist[mype]; + clastvtx = cvtxdist[mype+1]; + + /* Create the cmap of what you know so far locally */ + cnvtxs = 0; + for (i=0; i<nvtxs; i++) { + if (match[i] >= KEEP_BIT) { + k = match[i] - KEEP_BIT; + if (k>=firstvtx && k<firstvtx+i) + continue; /* Both (i,k) are local and i has been matched via the (k,i) side */ + + cmap[i] = cfirstvtx + cnvtxs++; + if (k != firstvtx+i && (k>=firstvtx && k<lastvtx)) { /* I'm matched locally */ + cmap[k-firstvtx] = cmap[i]; + match[k-firstvtx] += KEEP_BIT; /* Add the KEEP_BIT to simplify coding */ + } + } + } + ASSERT(ctrl, cnvtxs == clastvtx-cfirstvtx); + + CommInterfaceData(ctrl, graph, cmap, scmap, rcmap); + + /* Update the cmap of the locally stored vertices that will go away. + * The remote processor assigned cmap for them */ + for (i=0; i<nvtxs; i++) { + if (match[i] < KEEP_BIT) { /* Only vertices that go away satisfy this*/ + cmap[i] = rcmap[BSearch(graph->nrecv, recvind, match[i])]; + } + } + + CommInterfaceData(ctrl, graph, cmap, scmap, rcmap); + + +#ifdef DEBUG_CONTRACT + PrintVector(ctrl, nvtxs, firstvtx, cmap, "Cmap"); +#endif + + + /************************************************************* + * Determine how many adjcency lists you need to send/receive. + **************************************************************/ + /* Use wspace->pairs as the tmp space for the boundary vertices that are sent and received */ + scand = wspace->pairs; + rcand = graph->rcand = (KeyValueType *)GKmalloc(recvptr[nnbrs]*sizeof(KeyValueType), "CreateCoarseGraph: rcand"); + + nkeepsize = nsend = nrecv = 0; + for (i=0; i<nvtxs; i++) { + if (match[i] < KEEP_BIT) { /* This is going away */ + scand[nsend].key = match[i]; + scand[nsend].val = i; + nsend++; + } + else { + nkeepsize += (xadj[i+1]-xadj[i]); + + k = match[i]-KEEP_BIT; + if (k<firstvtx || k>=lastvtx) { /* This is comming from afar */ + rcand[nrecv].key = k; + rcand[nrecv].val = cmap[i] - cfirstvtx; /* Set it for use during the partition projection */ + ASSERT(ctrl, rcand[nrecv].val>=0 && rcand[nrecv].val<cnvtxs); + nrecv++; + } + } + } + + +#ifdef DEBUG_CONTRACT + PrintPairs(ctrl, nsend, scand, "scand"); + PrintPairs(ctrl, nrecv, rcand, "rcand"); +#endif + + /*************************************************************** + * Determine how many lists and their sizes you will send and + * received for each of the neighboring PEs + ****************************************************************/ + rsizes = wspace->pv1; + ssizes = wspace->pv2; + idxset(nnbrs, 0, ssizes); + idxset(nnbrs, 0, rsizes); + rlens = graph->rlens = idxmalloc(nnbrs+1, "CreateCoarseGraph: graph->rlens"); + slens = graph->slens = idxmalloc(nnbrs+1, "CreateCoarseGraph: graph->slens"); + + /* Take care the sending data first */ + ikeyvalsort(nsend, scand); + slens[0] = 0; + for (k=i=0; i<nnbrs; i++) { + otherlastvtx = vtxdist[peind[i]+1]; + for (; k<nsend && scand[k].key < otherlastvtx; k++) + ssizes[i] += (xadj[scand[k].val+1]-xadj[scand[k].val]); + slens[i+1] = k; + } + + /* Take care the receiving data next. You cannot yet determine the rsizes[] */ + ikeyvalsort(nrecv, rcand); + rlens[0] = 0; + for (k=i=0; i<nnbrs; i++) { + otherlastvtx = vtxdist[peind[i]+1]; + for (; k<nrecv && rcand[k].key < otherlastvtx; k++); + rlens[i+1] = k; + } + +#ifdef DEBUG_CONTRACT + PrintVector(ctrl, nnbrs+1, 0, slens, "slens"); + PrintVector(ctrl, nnbrs+1, 0, rlens, "rlens"); +#endif + + /*************************************************************** + * Exchange size information + ****************************************************************/ + /* Issue the receives first. */ + for (i=0; i<nnbrs; i++) { + if (rlens[i+1]-rlens[i] > 0) /* Issue a receive only if you are getting something */ + MPI_Irecv((void *)(rsizes+i), 1, IDX_DATATYPE, peind[i], 1, ctrl->comm, ctrl->rreq+i); + } + + /* Take care the sending data next */ + for (i=0; i<nnbrs; i++) { + if (slens[i+1]-slens[i] > 0) /* Issue a send only if you are sending something */ + MPI_Isend((void *)(ssizes+i), 1, IDX_DATATYPE, peind[i], 1, ctrl->comm, ctrl->sreq+i); + } + + /* OK, now get into the loop waiting for the operations to finish */ + for (i=0; i<nnbrs; i++) { + if (rlens[i+1]-rlens[i] > 0) + MPI_Wait(ctrl->rreq+i, &ctrl->status); + } + for (i=0; i<nnbrs; i++) { + if (slens[i+1]-slens[i] > 0) + MPI_Wait(ctrl->sreq+i, &ctrl->status); + } + + +#ifdef DEBUG_CONTRACT + PrintVector(ctrl, nnbrs, 0, rsizes, "rsizes"); + PrintVector(ctrl, nnbrs, 0, ssizes, "ssizes"); +#endif + + /************************************************************* + * Allocate memory for received/sent graphs and start sending + * and receiving data. + * rgraph and sgraph is a different data structure than CSR + * to facilitate single message exchange. + **************************************************************/ + nrecvsize = idxsum(nnbrs, rsizes); + nsendsize = idxsum(nnbrs, ssizes); + if ((4+ncon)*(nrecv+nsend) + 2*(nrecvsize+nsendsize) <= wspace->nlarge) { + rgraph = (idxtype *)wspace->degrees; + sgraph = rgraph + (4+ncon)*nrecv+2*nrecvsize; + } + else { + rgraph = idxmalloc((4+ncon)*nrecv+2*nrecvsize, "CreateCoarseGraph: rgraph"); + sgraph = idxmalloc((4+ncon)*nsend+2*nsendsize, "CreateCoarseGraph: sgraph"); + } + + /* Deal with the received portion first */ + for (l=i=0; i<nnbrs; i++) { + /* Issue a receive only if you are getting something */ + if (rlens[i+1]-rlens[i] > 0) { + MPI_Irecv((void *)(rgraph+l), (4+ncon)*(rlens[i+1]-rlens[i])+2*rsizes[i], IDX_DATATYPE, peind[i], 1, ctrl->comm, ctrl->rreq+i); + l += (4+ncon)*(rlens[i+1]-rlens[i])+2*rsizes[i]; + } + } + + + /* Deal with the sent portion now */ + for (ll=l=i=0; i<nnbrs; i++) { + if (slens[i+1]-slens[i] > 0) { /* Issue a send only if you are sending something */ + for (k=slens[i]; k<slens[i+1]; k++) { + ii = scand[k].val; + sgraph[ll++] = firstvtx+ii; + sgraph[ll++] = xadj[ii+1]-xadj[ii]; + for (h=0; h<ncon; h++) + sgraph[ll++] = vwgt[ii*ncon+h]; + sgraph[ll++] = (ctrl->partType == STATIC_PARTITION) ? -1 : vsize[ii]; + sgraph[ll++] = (ctrl->partType == STATIC_PARTITION) ? -1 : home[ii]; + for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) { + sgraph[ll++] = cmap[ladjncy[jj]]; + sgraph[ll++] = adjwgt[jj]; + } + } + + ASSERT(ctrl, ll-l == (4+ncon)*(slens[i+1]-slens[i])+2*ssizes[i]); + + /* myprintf(ctrl, "Sending to pe:%d, %d lists of size %d\n", peind[i], slens[i+1]-slens[i], ssizes[i]); */ + MPI_Isend((void *)(sgraph+l), ll-l, IDX_DATATYPE, peind[i], 1, ctrl->comm, ctrl->sreq+i); + l = ll; + } + } + + /* OK, now get into the loop waiting for the operations to finish */ + for (i=0; i<nnbrs; i++) { + if (rlens[i+1]-rlens[i] > 0) + MPI_Wait(ctrl->rreq+i, &ctrl->status); + } + for (i=0; i<nnbrs; i++) { + if (slens[i+1]-slens[i] > 0) + MPI_Wait(ctrl->sreq+i, &ctrl->status); + } + + +#ifdef DEBUG_CONTRACT + rprintf(ctrl, "Graphs were sent!\n"); + PrintTransferedGraphs(ctrl, nnbrs, peind, slens, rlens, sgraph, rgraph); +#endif + + /************************************************************* + * Setup the mapping from indices returned by BSearch to + * those that are actually stored + **************************************************************/ + perm = idxsmalloc(recvptr[nnbrs], -1, "CreateCoarseGraph: perm"); + for (j=i=0; i<nrecv; i++) { + /* myprintf(ctrl, "For received vertex %d, set perm[%d]=%d\n", rgraph[j], BSearch(graph->nrecv, recvind, rgraph[j]), j+ncon); */ + perm[BSearch(graph->nrecv, recvind, rgraph[j])] = j+1; + j += (4+ncon)+2*rgraph[j+1]; + } + + /************************************************************* + * Finally, create the coarser graph + **************************************************************/ + /* Allocate memory for the coarser graph, and fire up coarsening */ + cxadj = cgraph->xadj = idxmalloc(cnvtxs+1, "CreateCoarserGraph: cxadj"); + cvwgt = cgraph->vwgt = idxmalloc(cnvtxs*ncon, "CreateCoarserGraph: cvwgt"); + if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) { + cvsize = cgraph->vsize = idxmalloc(cnvtxs, "CreateCoarserGraph: cvsize"); + chome = cgraph->home = idxmalloc(cnvtxs, "CreateCoarserGraph: chome"); + } + cnvwgt = cgraph->nvwgt = fmalloc(cnvtxs*ncon, "CreateCoarserGraph: cnvwgt"); + cadjncy = idxmalloc(2*(nkeepsize+nrecvsize), "CreateCoarserGraph: cadjncy"); + cadjwgt = cadjncy + nkeepsize+nrecvsize; + + iset(8192, -1, htable); + + cxadj[0] = cnvtxs = cnedges = 0; + for (i=0; i<nvtxs; i++) { + if (match[i] >= KEEP_BIT) { + v = firstvtx+i; + u = match[i]-KEEP_BIT; + + if (u>=firstvtx && u<lastvtx && v > u) + continue; /* I have already collapsed it as (u,v) */ + + /* Collapse the v vertex first, which you know is local */ + for (h=0; h<ncon; h++) + cvwgt[cnvtxs*ncon+h] = vwgt[i*ncon+h]; + if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) { + cvsize[cnvtxs] = vsize[i]; + chome[cnvtxs] = home[i]; + } + nedges = 0; + + for (j=xadj[i]; j<xadj[i+1]; j++) { + k = cmap[ladjncy[j]]; + if (k != cfirstvtx+cnvtxs) { /* If this is not an internal edge */ + l = k&mask; + if (htable[l] == -1) { /* Seeing this for first time */ + htable[l] = k; + htableidx[l] = cnedges+nedges; + cadjncy[cnedges+nedges] = k; + cadjwgt[cnedges+nedges++] = adjwgt[j]; + } + else if (htable[l] == k) { + cadjwgt[htableidx[l]] += adjwgt[j]; + } + else { /* Now you have to go and do a search. Expensive case */ + for (l=0; l<nedges; l++) { + if (cadjncy[cnedges+l] == k) + break; + } + if (l < nedges) { + cadjwgt[cnedges+l] += adjwgt[j]; + } + else { + cadjncy[cnedges+nedges] = k; + cadjwgt[cnedges+nedges++] = adjwgt[j]; + } + } + } + } + + /* Collapse the u vertex next */ + if (v != u) { + if (u>=firstvtx && u<lastvtx) { /* Local vertex */ + u -= firstvtx; + for (h=0; h<ncon; h++) + cvwgt[cnvtxs*ncon+h] += vwgt[u*ncon+h]; + if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) { + cvsize[cnvtxs] += vsize[u]; + /* chome[cnvtxs] = home[u]; */ + } + + for (j=xadj[u]; j<xadj[u+1]; j++) { + k = cmap[ladjncy[j]]; + if (k != cfirstvtx+cnvtxs) { /* If this is not an internal edge */ + l = k&mask; + if (htable[l] == -1) { /* Seeing this for first time */ + htable[l] = k; + htableidx[l] = cnedges+nedges; + cadjncy[cnedges+nedges] = k; + cadjwgt[cnedges+nedges++] = adjwgt[j]; + } + else if (htable[l] == k) { + cadjwgt[htableidx[l]] += adjwgt[j]; + } + else { /* Now you have to go and do a search. Expensive case */ + for (l=0; l<nedges; l++) { + if (cadjncy[cnedges+l] == k) + break; + } + if (l < nedges) { + cadjwgt[cnedges+l] += adjwgt[j]; + } + else { + cadjncy[cnedges+nedges] = k; + cadjwgt[cnedges+nedges++] = adjwgt[j]; + } + } + } + } + } + else { /* Remote vertex */ + u = perm[BSearch(graph->nrecv, recvind, u)]; + for (h=0; h<ncon; h++) + /* Remember that the +1 stores the vertex weight */ + cvwgt[cnvtxs*ncon+h] += rgraph[(u+1)+h]; + if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) { + cvsize[cnvtxs] += rgraph[u+1+ncon]; + chome[cnvtxs] = rgraph[u+2+ncon]; + } + for (j=0; j<rgraph[u]; j++) { + k = rgraph[u+3+ncon+2*j]; + if (k != cfirstvtx+cnvtxs) { /* If this is not an internal edge */ + l = k&mask; + if (htable[l] == -1) { /* Seeing this for first time */ + htable[l] = k; + htableidx[l] = cnedges+nedges; + cadjncy[cnedges+nedges] = k; + cadjwgt[cnedges+nedges++] = rgraph[u+3+ncon+2*j+1]; + } + else if (htable[l] == k) { + cadjwgt[htableidx[l]] += rgraph[u+3+ncon+2*j+1]; + } + else { /* Now you have to go and do a search. Expensive case */ + for (l=0; l<nedges; l++) { + if (cadjncy[cnedges+l] == k) + break; + } + if (l < nedges) { + cadjwgt[cnedges+l] += rgraph[u+3+ncon+2*j+1]; + } + else { + cadjncy[cnedges+nedges] = k; + cadjwgt[cnedges+nedges++] = rgraph[u+3+ncon+2*j+1]; + } + } + } + } + } + } + + cnedges += nedges; + for (j=cxadj[cnvtxs]; j<cnedges; j++) + htable[cadjncy[j]&mask] = -1; /* reset the htable */ + cxadj[++cnvtxs] = cnedges; + } + } + + cgraph->nedges = cnedges; + + /* ADD: In order to keep from having to change this too much */ + /* ADD: I kept vwgt array and recomputed nvwgt for each coarser graph */ + for (j=0; j<cnvtxs; j++) + for (h=0; h<ncon; h++) + cgraph->nvwgt[j*ncon+h] = (float)(cvwgt[j*ncon+h])/(float)(ctrl->tvwgts[h]); + + cgraph->adjncy = idxmalloc(cnedges, "CreateCoarserGraph: cadjncy"); + cgraph->adjwgt = idxmalloc(cnedges, "CreateCoarserGraph: cadjwgt"); + idxcopy(cnedges, cadjncy, cgraph->adjncy); + idxcopy(cnedges, cadjwgt, cgraph->adjwgt); + free(cadjncy); + + free(perm); + + if (rgraph != (idxtype *)wspace->degrees) + GKfree((void **)&rgraph, (void **)&sgraph, LTERM); + +} + + |
