source: inundation/pymetis/metis-4.0/Lib/mfm.c @ 3381

Last change on this file since 3381 was 2051, checked in by jack, 19 years ago

Python interface to metis. Currently provides only the
METIS_PartMeshNodal function, since that is what is currently needed for partitioning.
Module name is metis.

File size: 10.8 KB
Line 
1/*
2 * Copyright 1997, Regents of the University of Minnesota
3 *
4 * mfm.c
5 *
6 * This file contains code that implements the edge-based FM refinement
7 *
8 * Started 7/23/97
9 * George
10 *
11 * $Id: mfm.c,v 1.3 1998/11/30 14:50:44 karypis Exp $
12 */
13
14#include <metis.h>
15
16
17/*************************************************************************
18* This function performs an edge-based FM refinement
19**************************************************************************/
20void MocFM_2WayEdgeRefine(CtrlType *ctrl, GraphType *graph, float *tpwgts, int npasses)
21{
22  int i, ii, j, k, l, kwgt, nvtxs, ncon, nbnd, nswaps, from, to, pass, me, limit, tmp, cnum;
23  idxtype *xadj, *adjncy, *adjwgt, *where, *id, *ed, *bndptr, *bndind;
24  idxtype *moved, *swaps, *perm, *qnum;
25  float *nvwgt, *npwgts, mindiff[MAXNCON], origbal, minbal, newbal;
26  PQueueType parts[MAXNCON][2];
27  int higain, oldgain, mincut, initcut, newcut, mincutorder;
28  float rtpwgts[2];
29
30  nvtxs = graph->nvtxs;
31  ncon = graph->ncon;
32  xadj = graph->xadj;
33  nvwgt = graph->nvwgt;
34  adjncy = graph->adjncy;
35  adjwgt = graph->adjwgt;
36  where = graph->where;
37  id = graph->id;
38  ed = graph->ed;
39  npwgts = graph->npwgts;
40  bndptr = graph->bndptr;
41  bndind = graph->bndind;
42
43  moved = idxwspacemalloc(ctrl, nvtxs);
44  swaps = idxwspacemalloc(ctrl, nvtxs);
45  perm = idxwspacemalloc(ctrl, nvtxs);
46  qnum = idxwspacemalloc(ctrl, nvtxs);
47
48  limit = amin(amax(0.01*nvtxs, 25), 150);
49
50  /* Initialize the queues */
51  for (i=0; i<ncon; i++) {
52    PQueueInit(ctrl, &parts[i][0], nvtxs, PLUS_GAINSPAN+1);
53    PQueueInit(ctrl, &parts[i][1], nvtxs, PLUS_GAINSPAN+1);
54  }
55  for (i=0; i<nvtxs; i++)
56    qnum[i] = samax(ncon, nvwgt+i*ncon);
57
58  origbal = Compute2WayHLoadImbalance(ncon, npwgts, tpwgts);
59
60  rtpwgts[0] = origbal*tpwgts[0];
61  rtpwgts[1] = origbal*tpwgts[1];
62
63
64  if (ctrl->dbglvl&DBG_REFINE) {
65    printf("Parts: [");
66    for (l=0; l<ncon; l++)
67      printf("(%.3f, %.3f) ", npwgts[l], npwgts[ncon+l]);
68    printf("] T[%.3f %.3f], Nv-Nb[%5d, %5d]. ICut: %6d, LB: %.3f\n", tpwgts[0], tpwgts[1], graph->nvtxs, graph->nbnd, graph->mincut, origbal);
69  }
70
71  idxset(nvtxs, -1, moved);
72  for (pass=0; pass<npasses; pass++) { /* Do a number of passes */
73    for (i=0; i<ncon; i++) { 
74      PQueueReset(&parts[i][0]);
75      PQueueReset(&parts[i][1]);
76    }
77
78    mincutorder = -1;
79    newcut = mincut = initcut = graph->mincut;
80    for (i=0; i<ncon; i++)
81      mindiff[i] = fabs(tpwgts[0]-npwgts[i]);
82    minbal = Compute2WayHLoadImbalance(ncon, npwgts, tpwgts);
83
84    ASSERT(ComputeCut(graph, where) == graph->mincut);
85    ASSERT(CheckBnd(graph));
86
87    /* Insert boundary nodes in the priority queues */
88    nbnd = graph->nbnd;
89    RandomPermute(nbnd, perm, 1);
90    for (ii=0; ii<nbnd; ii++) {
91      i = bndind[perm[ii]];
92      ASSERT(ed[i] > 0 || id[i] == 0);
93      ASSERT(bndptr[i] != -1);
94      PQueueInsert(&parts[qnum[i]][where[i]], i, ed[i]-id[i]);
95    }
96
97    for (nswaps=0; nswaps<nvtxs; nswaps++) {
98      SelectQueue(ncon, npwgts, rtpwgts, &from, &cnum, parts);
99      to = (from+1)%2;
100
101      if (from == -1 || (higain = PQueueGetMax(&parts[cnum][from])) == -1)
102        break;
103      ASSERT(bndptr[higain] != -1);
104
105      saxpy(ncon, 1.0, nvwgt+higain*ncon, 1, npwgts+to*ncon, 1);
106      saxpy(ncon, -1.0, nvwgt+higain*ncon, 1, npwgts+from*ncon, 1);
107
108      newcut -= (ed[higain]-id[higain]);
109      newbal = Compute2WayHLoadImbalance(ncon, npwgts, tpwgts);
110
111      if ((newcut < mincut && newbal-origbal <= .00001) || 
112          (newcut == mincut && (newbal < minbal || 
113                                (newbal == minbal && BetterBalance(ncon, npwgts, tpwgts, mindiff))))) {
114        mincut = newcut;
115        minbal = newbal;
116        mincutorder = nswaps;
117        for (i=0; i<ncon; i++)
118          mindiff[i] = fabs(tpwgts[0]-npwgts[i]);
119      }
120      else if (nswaps-mincutorder > limit) { /* We hit the limit, undo last move */
121        newcut += (ed[higain]-id[higain]);
122        saxpy(ncon, 1.0, nvwgt+higain*ncon, 1, npwgts+from*ncon, 1);
123        saxpy(ncon, -1.0, nvwgt+higain*ncon, 1, npwgts+to*ncon, 1);
124        break;
125      }
126
127      where[higain] = to;
128      moved[higain] = nswaps;
129      swaps[nswaps] = higain;
130
131      if (ctrl->dbglvl&DBG_MOVEINFO) {
132        printf("Moved %6d from %d(%d). Gain: %5d, Cut: %5d, NPwgts: ", higain, from, cnum, ed[higain]-id[higain], newcut);
133        for (l=0; l<ncon; l++) 
134          printf("(%.3f, %.3f) ", npwgts[l], npwgts[ncon+l]);
135        printf(", %.3f LB: %.3f\n", minbal, newbal);
136      }
137
138
139      /**************************************************************
140      * Update the id[i]/ed[i] values of the affected nodes
141      ***************************************************************/
142      SWAP(id[higain], ed[higain], tmp);
143      if (ed[higain] == 0 && xadj[higain] < xadj[higain+1]) 
144        BNDDelete(nbnd, bndind,  bndptr, higain);
145
146      for (j=xadj[higain]; j<xadj[higain+1]; j++) {
147        k = adjncy[j];
148        oldgain = ed[k]-id[k];
149
150        kwgt = (to == where[k] ? adjwgt[j] : -adjwgt[j]);
151        INC_DEC(id[k], ed[k], kwgt);
152
153        /* Update its boundary information and queue position */
154        if (bndptr[k] != -1) { /* If k was a boundary vertex */
155          if (ed[k] == 0) { /* Not a boundary vertex any more */
156            BNDDelete(nbnd, bndind, bndptr, k);
157            if (moved[k] == -1)  /* Remove it if in the queues */
158              PQueueDelete(&parts[qnum[k]][where[k]], k, oldgain);
159          }
160          else { /* If it has not been moved, update its position in the queue */
161            if (moved[k] == -1)
162              PQueueUpdate(&parts[qnum[k]][where[k]], k, oldgain, ed[k]-id[k]);
163          }
164        }
165        else {
166          if (ed[k] > 0) {  /* It will now become a boundary vertex */
167            BNDInsert(nbnd, bndind, bndptr, k);
168            if (moved[k] == -1) 
169              PQueueInsert(&parts[qnum[k]][where[k]], k, ed[k]-id[k]);
170          }
171        }
172      }
173
174    }
175
176
177    /****************************************************************
178    * Roll back computations
179    *****************************************************************/
180    for (i=0; i<nswaps; i++)
181      moved[swaps[i]] = -1;  /* reset moved array */
182    for (nswaps--; nswaps>mincutorder; nswaps--) {
183      higain = swaps[nswaps];
184
185      to = where[higain] = (where[higain]+1)%2;
186      SWAP(id[higain], ed[higain], tmp);
187      if (ed[higain] == 0 && bndptr[higain] != -1 && xadj[higain] < xadj[higain+1])
188        BNDDelete(nbnd, bndind,  bndptr, higain);
189      else if (ed[higain] > 0 && bndptr[higain] == -1)
190        BNDInsert(nbnd, bndind,  bndptr, higain);
191
192      saxpy(ncon, 1.0, nvwgt+higain*ncon, 1, npwgts+to*ncon, 1);
193      saxpy(ncon, -1.0, nvwgt+higain*ncon, 1, npwgts+((to+1)%2)*ncon, 1);
194      for (j=xadj[higain]; j<xadj[higain+1]; j++) {
195        k = adjncy[j];
196
197        kwgt = (to == where[k] ? adjwgt[j] : -adjwgt[j]);
198        INC_DEC(id[k], ed[k], kwgt);
199
200        if (bndptr[k] != -1 && ed[k] == 0)
201          BNDDelete(nbnd, bndind, bndptr, k);
202        if (bndptr[k] == -1 && ed[k] > 0)
203          BNDInsert(nbnd, bndind, bndptr, k);
204      }
205    }
206
207    if (ctrl->dbglvl&DBG_REFINE) {
208      printf("\tMincut: %6d at %5d, NBND: %6d, NPwgts: [", mincut, mincutorder, nbnd);
209      for (l=0; l<ncon; l++)
210        printf("(%.3f, %.3f) ", npwgts[l], npwgts[ncon+l]);
211      printf("], LB: %.3f\n", Compute2WayHLoadImbalance(ncon, npwgts, tpwgts));
212    }
213
214    graph->mincut = mincut;
215    graph->nbnd = nbnd;
216
217    if (mincutorder == -1 || mincut == initcut)
218      break;
219  }
220
221  for (i=0; i<ncon; i++) {
222    PQueueFree(ctrl, &parts[i][0]);
223    PQueueFree(ctrl, &parts[i][1]);
224  }
225
226  idxwspacefree(ctrl, nvtxs);
227  idxwspacefree(ctrl, nvtxs);
228  idxwspacefree(ctrl, nvtxs);
229  idxwspacefree(ctrl, nvtxs);
230
231}
232
233
234/*************************************************************************
235* This function selects the partition number and the queue from which
236* we will move vertices out
237**************************************************************************/ 
238void SelectQueue(int ncon, float *npwgts, float *tpwgts, int *from, int *cnum, PQueueType queues[MAXNCON][2])
239{
240  int i, part, maxgain=0;
241  float max, maxdiff=0.0;
242
243  *from = -1;
244  *cnum = -1;
245
246  /* First determine the side and the queue, irrespective of the presence of nodes */
247  for (part=0; part<2; part++) {
248    for (i=0; i<ncon; i++) {
249      if (npwgts[part*ncon+i]-tpwgts[part] >= maxdiff) {
250        maxdiff = npwgts[part*ncon+i]-tpwgts[part];
251        *from = part;
252        *cnum = i;
253      }
254    }
255  }
256
257  /* printf("Selected1 %d(%d) -> %d [%5f]\n", *from, *cnum, PQueueGetSize(&queues[*cnum][*from]), maxdiff); */
258
259  if (*from != -1 && PQueueGetSize(&queues[*cnum][*from]) == 0) {
260    /* The desired queue is empty, select a node from that side anyway */
261    for (i=0; i<ncon; i++) {
262      if (PQueueGetSize(&queues[i][*from]) > 0) {
263        max = npwgts[(*from)*ncon + i];
264        *cnum = i;
265        break;
266      }
267    }
268
269    for (i++; i<ncon; i++) {
270      if (npwgts[(*from)*ncon + i] > max && PQueueGetSize(&queues[i][*from]) > 0) {
271        max = npwgts[(*from)*ncon + i];
272        *cnum = i;
273      }
274    }
275  }
276
277  /* Check to see if you can focus on the cut */
278  if (maxdiff <= 0.0 || *from == -1) {
279    maxgain = -100000;
280
281    for (part=0; part<2; part++) {
282      for (i=0; i<ncon; i++) {
283        if (PQueueGetSize(&queues[i][part]) > 0 && PQueueGetKey(&queues[i][part]) > maxgain) {
284          maxgain = PQueueGetKey(&queues[i][part]); 
285          *from = part;
286          *cnum = i;
287        }
288      }
289    }
290  }
291
292  /* printf("Selected2 %d(%d) -> %d\n", *from, *cnum, PQueueGetSize(&queues[*cnum][*from])); */
293}
294
295
296
297
298
299/*************************************************************************
300* This function checks if the balance achieved is better than the diff
301* For now, it uses a 2-norm measure
302**************************************************************************/ 
303int BetterBalance(int ncon, float *npwgts, float *tpwgts, float *diff)
304{
305  int i;
306  float ndiff[MAXNCON];
307
308  for (i=0; i<ncon; i++)
309    ndiff[i] = fabs(tpwgts[0]-npwgts[i]);
310   
311  return snorm2(ncon, ndiff) < snorm2(ncon, diff);
312}
313
314
315
316/*************************************************************************
317* This function computes the load imbalance over all the constrains
318**************************************************************************/ 
319float Compute2WayHLoadImbalance(int ncon, float *npwgts, float *tpwgts)
320{
321  int i;
322  float max=0.0, temp;
323
324  for (i=0; i<ncon; i++) {
325    /* temp = amax(npwgts[i]/tpwgts[0], npwgts[ncon+i]/tpwgts[1]); */
326    temp = fabs(tpwgts[0]-npwgts[i])/tpwgts[0];
327    max = (max < temp ? temp : max);
328  }
329  return 1.0+max;
330}
331
332
333/*************************************************************************
334* This function computes the load imbalance over all the constrains
335* For now assume that we just want balanced partitionings
336**************************************************************************/ 
337void Compute2WayHLoadImbalanceVec(int ncon, float *npwgts, float *tpwgts, float *lbvec)
338{
339  int i;
340
341  for (i=0; i<ncon; i++) 
342    lbvec[i] = 1.0 + fabs(tpwgts[0]-npwgts[i])/tpwgts[0];
343}
344
Note: See TracBrowser for help on using the repository browser.