Charm++: Charm Source Code Documentation

00001 /*
00002  * Copyright 1997, Regents of the University of Minnesota
00003  *
00004  * mcoarsen.c
00005  *
00006  * This file contains code that performs graph coarsening
00007  *
00008  * Started 2/22/96
00009  * George
00010  *
00011  * $Id$
00012  *
00013  */
00014 
00015 #include <parmetislib.h>
00016 
00017 
00018 /*************************************************************************
00019 * This function creates the coarser graph
00020 **************************************************************************/
00021 void Moc_Global_CreateCoarseGraph(CtrlType *ctrl, GraphType *graph,
00022         WorkSpaceType *wspace, int cnvtxs)
00023 {
00024   int h, i, j, k, l, ii, jj, ll, nnbrs, nvtxs, nedges, ncon;
00025   int firstvtx, lastvtx, cfirstvtx, clastvtx, otherlastvtx;
00026   int npes=ctrl->npes, mype=ctrl->mype;
00027   int cnedges, nsend, nrecv, nkeepsize, nrecvsize, nsendsize, v, u;
00028   idxtype *xadj, *ladjncy, *adjwgt, *vwgt, *vsize, *vtxdist, *home;
00029   idxtype *match, *cmap, *rcmap, *scmap;
00030   idxtype *cxadj, *cadjncy, *cadjwgt, *cvwgt, *cvsize = NULL, *chome = NULL, *cvtxdist;
00031   idxtype *rsizes, *ssizes, *rlens, *slens, *rgraph, *sgraph, *perm;
00032   idxtype *peind, *recvptr, *recvind;
00033   floattype *nvwgt, *cnvwgt;
00034   GraphType *cgraph;
00035   KeyValueType *scand, *rcand;
00036   int mask=(1<<13)-1, htable[8192], htableidx[8192];
00037 
00038   nvtxs = graph->nvtxs;
00039   ncon = graph->ncon;
00040 
00041   vtxdist = graph->vtxdist;
00042   xadj = graph->xadj;
00043   vwgt = graph->vwgt;
00044   vsize = graph->vsize;
00045   nvwgt = graph->nvwgt;
00046   home = graph->home;
00047   ladjncy = graph->adjncy;
00048   adjwgt = graph->adjwgt;
00049 
00050   match = graph->match;
00051 
00052   firstvtx = vtxdist[mype];
00053   lastvtx = vtxdist[mype+1];
00054 
00055   cmap = graph->cmap = idxmalloc(nvtxs+graph->nrecv, "CreateCoarseGraph: cmap");
00056 
00057   nnbrs = graph->nnbrs;
00058   peind = graph->peind;
00059   recvind = graph->recvind;
00060   recvptr = graph->recvptr;
00061 
00062   /* Use wspace->indices as the tmp space for map of the boundary
00063    * vertices that are sent and received */
00064   scmap = wspace->indices;
00065   rcmap = cmap + nvtxs;
00066 
00067 
00068   /* Initialize the coarser graph */
00069   cgraph = CreateGraph();
00070   cgraph->nvtxs = cnvtxs;
00071   cgraph->ncon = ncon;
00072   cgraph->level = graph->level+1;
00073   cgraph->finer = graph;
00074   graph->coarser = cgraph;
00075 
00076 
00077 
00078   /*************************************************************
00079   * Obtain the vtxdist of the coarser graph 
00080   **************************************************************/
00081   cvtxdist = cgraph->vtxdist = idxmalloc(npes+1, "CreateCoarseGraph: cvtxdist");
00082   cvtxdist[npes] = cnvtxs;  /* Use last position in the cvtxdist as a temp buffer */
00083 
00084   MPI_Allgather((void *)(cvtxdist+npes), 1, IDX_DATATYPE, (void *)cvtxdist, 1, IDX_DATATYPE, ctrl->comm);
00085 
00086   MAKECSR(i, npes, cvtxdist);
00087 
00088   cgraph->gnvtxs = cvtxdist[npes];
00089 
00090 #ifdef DEBUG_CONTRACT
00091   PrintVector(ctrl, npes+1, 0, cvtxdist, "cvtxdist");
00092 #endif
00093 
00094 
00095   /*************************************************************
00096   * Construct the cmap vector 
00097   **************************************************************/
00098   cfirstvtx = cvtxdist[mype];
00099   clastvtx = cvtxdist[mype+1];
00100 
00101   /* Create the cmap of what you know so far locally */
00102   cnvtxs = 0;
00103   for (i=0; i<nvtxs; i++) {
00104     if (match[i] >= KEEP_BIT) {
00105       k = match[i] - KEEP_BIT;
00106       if (k>=firstvtx && k<firstvtx+i)
00107         continue;  /* Both (i,k) are local and i has been matched via the (k,i) side */
00108 
00109       cmap[i] = cfirstvtx + cnvtxs++;
00110       if (k != firstvtx+i && (k>=firstvtx && k<lastvtx)) { /* I'm matched locally */
00111         cmap[k-firstvtx] = cmap[i];
00112         match[k-firstvtx] += KEEP_BIT;  /* Add the KEEP_BIT to simplify coding */
00113       }
00114     }
00115   }
00116   ASSERT(ctrl, cnvtxs == clastvtx-cfirstvtx);
00117 
00118   CommInterfaceData(ctrl, graph, cmap, scmap, rcmap);
00119 
00120   /* Update the cmap of the locally stored vertices that will go away. 
00121    * The remote processor assigned cmap for them */
00122   for (i=0; i<nvtxs; i++) {
00123     if (match[i] < KEEP_BIT) { /* Only vertices that go away satisfy this*/
00124       cmap[i] = rcmap[BSearch(graph->nrecv, recvind, match[i])];
00125     }
00126   }
00127 
00128   CommInterfaceData(ctrl, graph, cmap, scmap, rcmap);
00129 
00130 
00131 #ifdef DEBUG_CONTRACT
00132   PrintVector(ctrl, nvtxs, firstvtx, cmap, "Cmap");
00133 #endif
00134 
00135 
00136   /*************************************************************
00137   * Determine how many adjcency lists you need to send/receive.
00138   **************************************************************/
00139   /* Use wspace->pairs as the tmp space for the boundary vertices that are sent and received */
00140   scand = wspace->pairs;
00141   rcand = graph->rcand = (KeyValueType *)GKmalloc(recvptr[nnbrs]*sizeof(KeyValueType), "CreateCoarseGraph: rcand");
00142 
00143   nkeepsize = nsend = nrecv = 0;
00144   for (i=0; i<nvtxs; i++) {
00145     if (match[i] < KEEP_BIT) { /* This is going away */
00146       scand[nsend].key = match[i];
00147       scand[nsend].val = i;
00148       nsend++;
00149     }
00150     else {
00151       nkeepsize += (xadj[i+1]-xadj[i]);
00152 
00153       k = match[i]-KEEP_BIT;
00154       if (k<firstvtx || k>=lastvtx) { /* This is comming from afar */
00155         rcand[nrecv].key = k;
00156         rcand[nrecv].val = cmap[i] - cfirstvtx;  /* Set it for use during the partition projection */
00157         ASSERT(ctrl, rcand[nrecv].val>=0 && rcand[nrecv].val<cnvtxs);
00158         nrecv++;
00159       }
00160     }
00161   }
00162 
00163 
00164 #ifdef DEBUG_CONTRACT
00165   PrintPairs(ctrl, nsend, scand, "scand");
00166   PrintPairs(ctrl, nrecv, rcand, "rcand");
00167 #endif
00168 
00169   /***************************************************************
00170   * Determine how many lists and their sizes  you will send and 
00171   * received for each of the neighboring PEs
00172   ****************************************************************/
00173   rsizes = wspace->pv1;
00174   ssizes = wspace->pv2;
00175   idxset(nnbrs, 0, ssizes);
00176   idxset(nnbrs, 0, rsizes);
00177   rlens = graph->rlens = idxmalloc(nnbrs+1, "CreateCoarseGraph: graph->rlens");
00178   slens = graph->slens = idxmalloc(nnbrs+1, "CreateCoarseGraph: graph->slens");
00179 
00180   /* Take care the sending data first */
00181   ikeyvalsort(nsend, scand);
00182   slens[0] = 0;
00183   for (k=i=0; i<nnbrs; i++) {
00184     otherlastvtx = vtxdist[peind[i]+1];
00185     for (; k<nsend && scand[k].key < otherlastvtx; k++)
00186       ssizes[i] += (xadj[scand[k].val+1]-xadj[scand[k].val]);
00187     slens[i+1] = k;
00188   }
00189 
00190   /* Take care the receiving data next. You cannot yet determine the rsizes[] */
00191   ikeyvalsort(nrecv, rcand);
00192   rlens[0] = 0;
00193   for (k=i=0; i<nnbrs; i++) {
00194     otherlastvtx = vtxdist[peind[i]+1];
00195     for (; k<nrecv && rcand[k].key < otherlastvtx; k++);
00196     rlens[i+1] = k;
00197   }
00198 
00199 #ifdef DEBUG_CONTRACT
00200   PrintVector(ctrl, nnbrs+1, 0, slens, "slens");
00201   PrintVector(ctrl, nnbrs+1, 0, rlens, "rlens");
00202 #endif
00203 
00204   /***************************************************************
00205   * Exchange size information
00206   ****************************************************************/
00207   /* Issue the receives first. */
00208   for (i=0; i<nnbrs; i++) {
00209     if (rlens[i+1]-rlens[i] > 0)  /* Issue a receive only if you are getting something */
00210       MPI_Irecv((void *)(rsizes+i), 1, IDX_DATATYPE, peind[i], 1, ctrl->comm, ctrl->rreq+i);
00211   }
00212 
00213   /* Take care the sending data next */
00214   for (i=0; i<nnbrs; i++) {
00215     if (slens[i+1]-slens[i] > 0)  /* Issue a send only if you are sending something */
00216       MPI_Isend((void *)(ssizes+i), 1, IDX_DATATYPE, peind[i], 1, ctrl->comm, ctrl->sreq+i);
00217   }
00218 
00219   /* OK, now get into the loop waiting for the operations to finish */
00220   for (i=0; i<nnbrs; i++) {
00221     if (rlens[i+1]-rlens[i] > 0)  
00222       MPI_Wait(ctrl->rreq+i, &ctrl->status);
00223   }
00224   for (i=0; i<nnbrs; i++) {
00225     if (slens[i+1]-slens[i] > 0)  
00226       MPI_Wait(ctrl->sreq+i, &ctrl->status);
00227   }
00228 
00229 
00230 #ifdef DEBUG_CONTRACT
00231   PrintVector(ctrl, nnbrs, 0, rsizes, "rsizes");
00232   PrintVector(ctrl, nnbrs, 0, ssizes, "ssizes");
00233 #endif
00234 
00235   /*************************************************************
00236   * Allocate memory for received/sent graphs and start sending 
00237   * and receiving data.
00238   * rgraph and sgraph is a different data structure than CSR
00239   * to facilitate single message exchange.
00240   **************************************************************/
00241   nrecvsize = idxsum(nnbrs, rsizes);
00242   nsendsize = idxsum(nnbrs, ssizes);
00243   if ((4+ncon)*(nrecv+nsend) + 2*(nrecvsize+nsendsize) <= wspace->nlarge) {  
00244     rgraph = (idxtype *)wspace->degrees;
00245     sgraph = rgraph + (4+ncon)*nrecv+2*nrecvsize;
00246   }
00247   else {
00248     rgraph = idxmalloc((4+ncon)*nrecv+2*nrecvsize, "CreateCoarseGraph: rgraph");
00249     sgraph = idxmalloc((4+ncon)*nsend+2*nsendsize, "CreateCoarseGraph: sgraph");
00250   }
00251 
00252   /* Deal with the received portion first */
00253   for (l=i=0; i<nnbrs; i++) {
00254     /* Issue a receive only if you are getting something */
00255     if (rlens[i+1]-rlens[i] > 0) {
00256       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);
00257       l += (4+ncon)*(rlens[i+1]-rlens[i])+2*rsizes[i];
00258     }
00259   }
00260 
00261 
00262   /* Deal with the sent portion now */
00263   for (ll=l=i=0; i<nnbrs; i++) {
00264     if (slens[i+1]-slens[i] > 0) {  /* Issue a send only if you are sending something */
00265       for (k=slens[i]; k<slens[i+1]; k++) {
00266         ii = scand[k].val;
00267         sgraph[ll++] = firstvtx+ii;
00268         sgraph[ll++] = xadj[ii+1]-xadj[ii];
00269         for (h=0; h<ncon; h++)
00270           sgraph[ll++] = vwgt[ii*ncon+h];
00271         sgraph[ll++] = (ctrl->partType == STATIC_PARTITION) ? -1 : vsize[ii];
00272         sgraph[ll++] = (ctrl->partType == STATIC_PARTITION) ? -1 : home[ii];
00273         for (jj=xadj[ii]; jj<xadj[ii+1]; jj++) {
00274           sgraph[ll++] = cmap[ladjncy[jj]];
00275           sgraph[ll++] = adjwgt[jj];
00276         }
00277       }
00278 
00279       ASSERT(ctrl, ll-l == (4+ncon)*(slens[i+1]-slens[i])+2*ssizes[i]);
00280 
00281       /* myprintf(ctrl, "Sending to pe:%d, %d lists of size %d\n", peind[i], slens[i+1]-slens[i], ssizes[i]); */
00282       MPI_Isend((void *)(sgraph+l), ll-l, IDX_DATATYPE, peind[i], 1, ctrl->comm, ctrl->sreq+i);
00283       l = ll;
00284     }
00285   }
00286 
00287   /* OK, now get into the loop waiting for the operations to finish */
00288   for (i=0; i<nnbrs; i++) {
00289     if (rlens[i+1]-rlens[i] > 0)  
00290       MPI_Wait(ctrl->rreq+i, &ctrl->status);
00291   }
00292   for (i=0; i<nnbrs; i++) {
00293     if (slens[i+1]-slens[i] > 0)  
00294       MPI_Wait(ctrl->sreq+i, &ctrl->status);
00295   }
00296 
00297 
00298 #ifdef DEBUG_CONTRACT
00299   rprintf(ctrl, "Graphs were sent!\n");
00300   PrintTransferedGraphs(ctrl, nnbrs, peind, slens, rlens, sgraph, rgraph);
00301 #endif
00302 
00303   /*************************************************************
00304   * Setup the mapping from indices returned by BSearch to 
00305   * those that are actually stored 
00306   **************************************************************/
00307   perm = idxsmalloc(recvptr[nnbrs], -1, "CreateCoarseGraph: perm");
00308   for (j=i=0; i<nrecv; i++) {
00309   /*   myprintf(ctrl, "For received vertex %d, set perm[%d]=%d\n", rgraph[j], BSearch(graph->nrecv, recvind, rgraph[j]), j+ncon); */
00310     perm[BSearch(graph->nrecv, recvind, rgraph[j])] = j+1;
00311     j += (4+ncon)+2*rgraph[j+1];
00312   }
00313 
00314   /*************************************************************
00315   * Finally, create the coarser graph
00316   **************************************************************/
00317   /* Allocate memory for the coarser graph, and fire up coarsening */
00318   cxadj = cgraph->xadj = idxmalloc(cnvtxs+1, "CreateCoarserGraph: cxadj");
00319   cvwgt = cgraph->vwgt = idxmalloc(cnvtxs*ncon, "CreateCoarserGraph: cvwgt");
00320   if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) {
00321     cvsize = cgraph->vsize = idxmalloc(cnvtxs, "CreateCoarserGraph: cvsize");
00322     chome = cgraph->home = idxmalloc(cnvtxs, "CreateCoarserGraph: chome");
00323   }
00324   cnvwgt = cgraph->nvwgt = fmalloc(cnvtxs*ncon, "CreateCoarserGraph: cnvwgt");
00325   cadjncy = idxmalloc(2*(nkeepsize+nrecvsize), "CreateCoarserGraph: cadjncy");
00326   cadjwgt = cadjncy + nkeepsize+nrecvsize;
00327 
00328   iset(8192, -1, htable);
00329 
00330   cxadj[0] = cnvtxs = cnedges = 0;
00331   for (i=0; i<nvtxs; i++) {
00332     if (match[i] >= KEEP_BIT) {
00333       v = firstvtx+i; 
00334       u = match[i]-KEEP_BIT;
00335 
00336       if (u>=firstvtx && u<lastvtx && v > u) 
00337         continue;  /* I have already collapsed it as (u,v) */
00338 
00339       /* Collapse the v vertex first, which you know is local */
00340       for (h=0; h<ncon; h++)
00341         cvwgt[cnvtxs*ncon+h] = vwgt[i*ncon+h];
00342       if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) {
00343         cvsize[cnvtxs] = vsize[i];
00344         chome[cnvtxs] = home[i];
00345       }
00346       nedges = 0;
00347 
00348       for (j=xadj[i]; j<xadj[i+1]; j++) {
00349         k = cmap[ladjncy[j]];
00350         if (k != cfirstvtx+cnvtxs) {  /* If this is not an internal edge */
00351           l = k&mask;
00352           if (htable[l] == -1) { /* Seeing this for first time */
00353             htable[l] = k;
00354             htableidx[l] = cnedges+nedges;
00355             cadjncy[cnedges+nedges] = k; 
00356             cadjwgt[cnedges+nedges++] = adjwgt[j];
00357           }
00358           else if (htable[l] == k) {
00359             cadjwgt[htableidx[l]] += adjwgt[j];
00360           }
00361           else { /* Now you have to go and do a search. Expensive case */
00362             for (l=0; l<nedges; l++) {
00363               if (cadjncy[cnedges+l] == k)
00364                 break;
00365             }
00366             if (l < nedges) {
00367               cadjwgt[cnedges+l] += adjwgt[j];
00368             }
00369             else {
00370               cadjncy[cnedges+nedges] = k; 
00371               cadjwgt[cnedges+nedges++] = adjwgt[j];
00372             }
00373           }
00374         }
00375       }
00376 
00377       /* Collapse the u vertex next */
00378       if (v != u) { 
00379         if (u>=firstvtx && u<lastvtx) { /* Local vertex */
00380           u -= firstvtx;
00381           for (h=0; h<ncon; h++)
00382             cvwgt[cnvtxs*ncon+h] += vwgt[u*ncon+h];
00383           if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) {
00384             cvsize[cnvtxs] += vsize[u];
00385             /* chome[cnvtxs] = home[u]; */
00386           }
00387 
00388           for (j=xadj[u]; j<xadj[u+1]; j++) {
00389             k = cmap[ladjncy[j]];
00390             if (k != cfirstvtx+cnvtxs) {  /* If this is not an internal edge */
00391               l = k&mask;
00392               if (htable[l] == -1) { /* Seeing this for first time */
00393                 htable[l] = k;
00394                 htableidx[l] = cnedges+nedges;
00395                 cadjncy[cnedges+nedges] = k; 
00396                 cadjwgt[cnedges+nedges++] = adjwgt[j];
00397               }
00398               else if (htable[l] == k) {
00399                 cadjwgt[htableidx[l]] += adjwgt[j];
00400               }
00401               else { /* Now you have to go and do a search. Expensive case */
00402                 for (l=0; l<nedges; l++) {
00403                   if (cadjncy[cnedges+l] == k)
00404                     break;
00405                 }
00406                 if (l < nedges) {
00407                   cadjwgt[cnedges+l] += adjwgt[j];
00408                 }
00409                 else {
00410                   cadjncy[cnedges+nedges] = k; 
00411                   cadjwgt[cnedges+nedges++] = adjwgt[j];
00412                 }
00413               }
00414             }
00415           }
00416         }
00417         else { /* Remote vertex */
00418           u = perm[BSearch(graph->nrecv, recvind, u)];
00419           for (h=0; h<ncon; h++)
00420             /* Remember that the +1 stores the vertex weight */
00421             cvwgt[cnvtxs*ncon+h] += rgraph[(u+1)+h];
00422             if (ctrl->partType == ADAPTIVE_PARTITION || ctrl->partType == REFINE_PARTITION) {
00423               cvsize[cnvtxs] += rgraph[u+1+ncon];
00424               chome[cnvtxs] = rgraph[u+2+ncon];
00425             }
00426           for (j=0; j<rgraph[u]; j++) {
00427             k = rgraph[u+3+ncon+2*j];
00428             if (k != cfirstvtx+cnvtxs) {  /* If this is not an internal edge */
00429               l = k&mask;
00430               if (htable[l] == -1) { /* Seeing this for first time */
00431                 htable[l] = k;
00432                 htableidx[l] = cnedges+nedges;
00433                 cadjncy[cnedges+nedges] = k; 
00434                 cadjwgt[cnedges+nedges++] = rgraph[u+3+ncon+2*j+1];
00435               }
00436               else if (htable[l] == k) {
00437                 cadjwgt[htableidx[l]] += rgraph[u+3+ncon+2*j+1];
00438               }
00439               else { /* Now you have to go and do a search. Expensive case */
00440                 for (l=0; l<nedges; l++) {
00441                   if (cadjncy[cnedges+l] == k)
00442                     break;
00443                 }
00444                 if (l < nedges) {
00445                   cadjwgt[cnedges+l] += rgraph[u+3+ncon+2*j+1];
00446                 }
00447                 else {
00448                   cadjncy[cnedges+nedges] = k; 
00449                   cadjwgt[cnedges+nedges++] = rgraph[u+3+ncon+2*j+1];
00450                 }
00451               }
00452             }
00453           }
00454         }
00455       }
00456 
00457       cnedges += nedges;
00458       for (j=cxadj[cnvtxs]; j<cnedges; j++)
00459         htable[cadjncy[j]&mask] = -1;  /* reset the htable */
00460       cxadj[++cnvtxs] = cnedges;
00461     }
00462   }
00463 
00464   cgraph->nedges = cnedges;
00465 
00466   /* ADD:  In order to keep from having to change this too much */
00467   /* ADD:  I kept vwgt array and recomputed nvwgt for each coarser graph */
00468   for (j=0; j<cnvtxs; j++)
00469     for (h=0; h<ncon; h++)
00470       cgraph->nvwgt[j*ncon+h] = (floattype)(cvwgt[j*ncon+h])/(floattype)(ctrl->tvwgts[h]);
00471 
00472   cgraph->adjncy = idxmalloc(cnedges, "CreateCoarserGraph: cadjncy");
00473   cgraph->adjwgt = idxmalloc(cnedges, "CreateCoarserGraph: cadjwgt");
00474   idxcopy(cnedges, cadjncy, cgraph->adjncy);
00475   idxcopy(cnedges, cadjwgt, cgraph->adjwgt);
00476   free(cadjncy);
00477 
00478   free(perm);
00479 
00480   if (rgraph != (idxtype *)wspace->degrees) 
00481     GKfree((void **)&rgraph, (void **)&sgraph, LTERM);
00482 
00483 }
00484 
00485
libs/ck-libs/parmetis/ParMETISLib/coarsen.c
