source: branches/numpy/pymetis/metis-4.0/Lib/mmatch.c @ 6971

Last change on this file since 6971 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: 13.0 KB
Line 
1/*
2 * Copyright 1997, Regents of the University of Minnesota
3 *
4 * mmatch.c
5 *
6 * This file contains the code that computes matchings and creates the next
7 * level coarse graph.
8 *
9 * Started 7/23/97
10 * George
11 *
12 * $Id: mmatch.c,v 1.3 1998/11/30 14:50:44 karypis Exp $
13 *
14 */
15
16#include <metis.h>
17
18
19/*************************************************************************
20* This function finds a matching using the HEM heuristic
21**************************************************************************/
22void MCMatch_RM(CtrlType *ctrl, GraphType *graph)
23{
24  int i, ii, j, k, nvtxs, ncon, cnvtxs, maxidx;
25  idxtype *xadj, *adjncy, *adjwgt;
26  idxtype *match, *cmap, *perm;
27  float *nvwgt;
28
29  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->MatchTmr));
30
31  nvtxs = graph->nvtxs;
32  ncon = graph->ncon;
33  xadj = graph->xadj;
34  nvwgt = graph->nvwgt;
35  adjncy = graph->adjncy;
36  adjwgt = graph->adjwgt;
37
38  cmap = graph->cmap;
39  match = idxset(nvtxs, UNMATCHED, idxwspacemalloc(ctrl, nvtxs));
40
41  perm = idxwspacemalloc(ctrl, nvtxs);
42  RandomPermute(nvtxs, perm, 1);
43
44  cnvtxs = 0;
45  for (ii=0; ii<nvtxs; ii++) {
46    i = perm[ii];
47
48    if (match[i] == UNMATCHED) {  /* Unmatched */
49      maxidx = i;
50
51      /* Find a random matching, subject to maxvwgt constraints */
52      for (j=xadj[i]; j<xadj[i+1]; j++) {
53        k = adjncy[j];
54        if (match[k] == UNMATCHED && AreAllVwgtsBelowFast(ncon, nvwgt+i*ncon, nvwgt+k*ncon, ctrl->nmaxvwgt)) {
55          maxidx = k;
56          break;
57        }
58      }
59
60      cmap[i] = cmap[maxidx] = cnvtxs++;
61      match[i] = maxidx;
62      match[maxidx] = i;
63    }
64  }
65
66  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->MatchTmr));
67
68  CreateCoarseGraph(ctrl, graph, cnvtxs, match, perm);
69
70  idxwspacefree(ctrl, nvtxs);
71  idxwspacefree(ctrl, nvtxs);
72}
73
74
75
76/*************************************************************************
77* This function finds a matching using the HEM heuristic
78**************************************************************************/
79void MCMatch_HEM(CtrlType *ctrl, GraphType *graph)
80{
81  int i, ii, j, k, l, nvtxs, cnvtxs, ncon, maxidx, maxwgt;
82  idxtype *xadj, *adjncy, *adjwgt;
83  idxtype *match, *cmap, *perm;
84  float *nvwgt;
85
86  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->MatchTmr));
87
88  nvtxs = graph->nvtxs;
89  ncon = graph->ncon;
90  xadj = graph->xadj;
91  nvwgt = graph->nvwgt;
92  adjncy = graph->adjncy;
93  adjwgt = graph->adjwgt;
94
95  cmap = graph->cmap;
96  match = idxset(nvtxs, UNMATCHED, idxwspacemalloc(ctrl, nvtxs));
97
98  perm = idxwspacemalloc(ctrl, nvtxs);
99  RandomPermute(nvtxs, perm, 1);
100
101  cnvtxs = 0;
102  for (ii=0; ii<nvtxs; ii++) {
103    i = perm[ii];
104
105    if (match[i] == UNMATCHED) {  /* Unmatched */
106      maxidx = i;
107      maxwgt = 0;
108
109      /* Find a heavy-edge matching, subject to maxvwgt constraints */
110      for (j=xadj[i]; j<xadj[i+1]; j++) {
111        k = adjncy[j];
112        if (match[k] == UNMATCHED && maxwgt <= adjwgt[j] &&
113               AreAllVwgtsBelowFast(ncon, nvwgt+i*ncon, nvwgt+k*ncon, ctrl->nmaxvwgt)) {
114          maxwgt = adjwgt[j];
115          maxidx = adjncy[j];
116        }
117      }
118
119      cmap[i] = cmap[maxidx] = cnvtxs++;
120      match[i] = maxidx;
121      match[maxidx] = i;
122    }
123  }
124
125  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->MatchTmr));
126
127  CreateCoarseGraph(ctrl, graph, cnvtxs, match, perm);
128
129  idxwspacefree(ctrl, nvtxs);
130  idxwspacefree(ctrl, nvtxs);
131}
132
133
134
135/*************************************************************************
136* This function finds a matching using the HEM heuristic
137**************************************************************************/
138void MCMatch_SHEM(CtrlType *ctrl, GraphType *graph)
139{
140  int i, ii, j, k, nvtxs, cnvtxs, ncon, maxidx, maxwgt, avgdegree;
141  idxtype *xadj, *adjncy, *adjwgt;
142  idxtype *match, *cmap, *degrees, *perm, *tperm;
143  float *nvwgt;
144
145  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->MatchTmr));
146
147  nvtxs = graph->nvtxs;
148  ncon = graph->ncon;
149  xadj = graph->xadj;
150  nvwgt = graph->nvwgt;
151  adjncy = graph->adjncy;
152  adjwgt = graph->adjwgt;
153
154  cmap = graph->cmap;
155  match = idxset(nvtxs, UNMATCHED, idxwspacemalloc(ctrl, nvtxs));
156
157  perm = idxwspacemalloc(ctrl, nvtxs);
158  tperm = idxwspacemalloc(ctrl, nvtxs);
159  degrees = idxwspacemalloc(ctrl, nvtxs);
160
161  RandomPermute(nvtxs, tperm, 1);
162  avgdegree = 0.7*(xadj[nvtxs]/nvtxs);
163  for (i=0; i<nvtxs; i++) 
164    degrees[i] = (xadj[i+1]-xadj[i] > avgdegree ? avgdegree : xadj[i+1]-xadj[i]);
165  BucketSortKeysInc(nvtxs, avgdegree, degrees, tperm, perm);
166
167  cnvtxs = 0;
168
169  /* Take care any islands. Islands are matched with non-islands due to coarsening */
170  for (ii=0; ii<nvtxs; ii++) {
171    i = perm[ii];
172
173    if (match[i] == UNMATCHED) {  /* Unmatched */
174      if (xadj[i] < xadj[i+1])
175        break;
176
177      maxidx = i;
178      for (j=nvtxs-1; j>ii; j--) {
179        k = perm[j];
180        if (match[k] == UNMATCHED && xadj[k] < xadj[k+1]) {
181          maxidx = k;
182          break;
183        }
184      }
185
186      cmap[i] = cmap[maxidx] = cnvtxs++;
187      match[i] = maxidx;
188      match[maxidx] = i;
189    }
190  }
191
192  /* Continue with normal matching */
193  for (; ii<nvtxs; ii++) {
194    i = perm[ii];
195
196    if (match[i] == UNMATCHED) {  /* Unmatched */
197      maxidx = i;
198      maxwgt = 0;
199
200      /* Find a heavy-edge matching, subject to maxvwgt constraints */
201      for (j=xadj[i]; j<xadj[i+1]; j++) {
202        k = adjncy[j];
203        if (match[k] == UNMATCHED && maxwgt <= adjwgt[j] &&
204               AreAllVwgtsBelowFast(ncon, nvwgt+i*ncon, nvwgt+k*ncon, ctrl->nmaxvwgt)) {
205          maxwgt = adjwgt[j];
206          maxidx = adjncy[j];
207        }
208      }
209
210      cmap[i] = cmap[maxidx] = cnvtxs++;
211      match[i] = maxidx;
212      match[maxidx] = i;
213    }
214  }
215
216  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->MatchTmr));
217
218  idxwspacefree(ctrl, nvtxs);  /* degrees */
219  idxwspacefree(ctrl, nvtxs);  /* tperm */
220
221  CreateCoarseGraph(ctrl, graph, cnvtxs, match, perm);
222
223  idxwspacefree(ctrl, nvtxs);
224  idxwspacefree(ctrl, nvtxs);
225}
226
227
228
229/*************************************************************************
230* This function finds a matching using the HEM heuristic
231**************************************************************************/
232void MCMatch_SHEBM(CtrlType *ctrl, GraphType *graph, int norm)
233{
234  int i, ii, j, k, nvtxs, cnvtxs, ncon, maxidx, maxwgt, avgdegree;
235  idxtype *xadj, *adjncy, *adjwgt;
236  idxtype *match, *cmap, *degrees, *perm, *tperm;
237  float *nvwgt;
238
239  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->MatchTmr));
240
241  nvtxs = graph->nvtxs;
242  ncon = graph->ncon;
243  xadj = graph->xadj;
244  nvwgt = graph->nvwgt;
245  adjncy = graph->adjncy;
246  adjwgt = graph->adjwgt;
247
248  cmap = graph->cmap;
249  match = idxset(nvtxs, UNMATCHED, idxwspacemalloc(ctrl, nvtxs));
250
251  perm = idxwspacemalloc(ctrl, nvtxs);
252  tperm = idxwspacemalloc(ctrl, nvtxs);
253  degrees = idxwspacemalloc(ctrl, nvtxs);
254
255  RandomPermute(nvtxs, tperm, 1);
256  avgdegree = 0.7*(xadj[nvtxs]/nvtxs);
257  for (i=0; i<nvtxs; i++) 
258    degrees[i] = (xadj[i+1]-xadj[i] > avgdegree ? avgdegree : xadj[i+1]-xadj[i]);
259  BucketSortKeysInc(nvtxs, avgdegree, degrees, tperm, perm);
260
261  cnvtxs = 0;
262
263  /* Take care any islands. Islands are matched with non-islands due to coarsening */
264  for (ii=0; ii<nvtxs; ii++) {
265    i = perm[ii];
266
267    if (match[i] == UNMATCHED) {  /* Unmatched */
268      if (xadj[i] < xadj[i+1])
269        break;
270
271      maxidx = i;
272      for (j=nvtxs-1; j>ii; j--) {
273        k = perm[j];
274        if (match[k] == UNMATCHED && xadj[k] < xadj[k+1]) {
275          maxidx = k;
276          break;
277        }
278      }
279
280      cmap[i] = cmap[maxidx] = cnvtxs++;
281      match[i] = maxidx;
282      match[maxidx] = i;
283    }
284  }
285
286  /* Continue with normal matching */
287  for (; ii<nvtxs; ii++) {
288    i = perm[ii];
289
290    if (match[i] == UNMATCHED) {  /* Unmatched */
291      maxidx = i;
292      maxwgt = -1;
293
294      /* Find a heavy-edge matching, subject to maxvwgt constraints */
295      for (j=xadj[i]; j<xadj[i+1]; j++) {
296        k = adjncy[j];
297
298        if (match[k] == UNMATCHED && 
299            AreAllVwgtsBelowFast(ncon, nvwgt+i*ncon, nvwgt+k*ncon, ctrl->nmaxvwgt) &&
300            (maxwgt < adjwgt[j] || 
301              (maxwgt == adjwgt[j] && 
302               BetterVBalance(ncon, norm, nvwgt+i*ncon, nvwgt+maxidx*ncon, nvwgt+k*ncon) >= 0
303              )
304            )
305           ) {
306          maxwgt = adjwgt[j];
307          maxidx = k;
308        }
309      }
310
311      cmap[i] = cmap[maxidx] = cnvtxs++;
312      match[i] = maxidx;
313      match[maxidx] = i;
314    }
315  }
316
317  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->MatchTmr));
318
319  idxwspacefree(ctrl, nvtxs);  /* degrees */
320  idxwspacefree(ctrl, nvtxs);  /* tperm */
321
322  CreateCoarseGraph(ctrl, graph, cnvtxs, match, perm);
323
324  idxwspacefree(ctrl, nvtxs);
325  idxwspacefree(ctrl, nvtxs);
326}
327
328
329
330/*************************************************************************
331* This function finds a matching using the HEM heuristic
332**************************************************************************/
333void MCMatch_SBHEM(CtrlType *ctrl, GraphType *graph, int norm)
334{
335  int i, ii, j, k, nvtxs, cnvtxs, ncon, maxidx, maxwgt, avgdegree;
336  idxtype *xadj, *adjncy, *adjwgt;
337  idxtype *match, *cmap, *degrees, *perm, *tperm;
338  float *nvwgt, vbal;
339
340  IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->MatchTmr));
341
342  nvtxs = graph->nvtxs;
343  ncon = graph->ncon;
344  xadj = graph->xadj;
345  nvwgt = graph->nvwgt;
346  adjncy = graph->adjncy;
347  adjwgt = graph->adjwgt;
348
349  cmap = graph->cmap;
350  match = idxset(nvtxs, UNMATCHED, idxwspacemalloc(ctrl, nvtxs));
351
352  perm = idxwspacemalloc(ctrl, nvtxs);
353  tperm = idxwspacemalloc(ctrl, nvtxs);
354  degrees = idxwspacemalloc(ctrl, nvtxs);
355
356  RandomPermute(nvtxs, tperm, 1);
357  avgdegree = 0.7*(xadj[nvtxs]/nvtxs);
358  for (i=0; i<nvtxs; i++) 
359    degrees[i] = (xadj[i+1]-xadj[i] > avgdegree ? avgdegree : xadj[i+1]-xadj[i]);
360  BucketSortKeysInc(nvtxs, avgdegree, degrees, tperm, perm);
361
362  cnvtxs = 0;
363
364  /* Take care any islands. Islands are matched with non-islands due to coarsening */
365  for (ii=0; ii<nvtxs; ii++) {
366    i = perm[ii];
367
368    if (match[i] == UNMATCHED) {  /* Unmatched */
369      if (xadj[i] < xadj[i+1])
370        break;
371
372      maxidx = i;
373      for (j=nvtxs-1; j>ii; j--) {
374        k = perm[j];
375        if (match[k] == UNMATCHED && xadj[k] < xadj[k+1]) {
376          maxidx = k;
377          break;
378        }
379      }
380
381      cmap[i] = cmap[maxidx] = cnvtxs++;
382      match[i] = maxidx;
383      match[maxidx] = i;
384    }
385  }
386
387  /* Continue with normal matching */
388  for (; ii<nvtxs; ii++) {
389    i = perm[ii];
390
391    if (match[i] == UNMATCHED) {  /* Unmatched */
392      maxidx = i;
393      maxwgt = -1;
394      vbal = 0.0;
395
396      /* Find a heavy-edge matching, subject to maxvwgt constraints */
397      for (j=xadj[i]; j<xadj[i+1]; j++) {
398        k = adjncy[j];
399        if (match[k] == UNMATCHED && AreAllVwgtsBelowFast(ncon, nvwgt+i*ncon, nvwgt+k*ncon, ctrl->nmaxvwgt)) {
400          if (maxidx != i)
401            vbal = BetterVBalance(ncon, norm, nvwgt+i*ncon, nvwgt+maxidx*ncon, nvwgt+k*ncon);
402
403          if (vbal > 0 || (vbal > -.01 && maxwgt < adjwgt[j])) {
404            maxwgt = adjwgt[j];
405            maxidx = k;
406          }
407        }
408      }
409
410      cmap[i] = cmap[maxidx] = cnvtxs++;
411      match[i] = maxidx;
412      match[maxidx] = i;
413    }
414  }
415
416  IFSET(ctrl->dbglvl, DBG_TIME, stoptimer(ctrl->MatchTmr));
417
418  idxwspacefree(ctrl, nvtxs);  /* degrees */
419  idxwspacefree(ctrl, nvtxs);  /* tperm */
420
421  CreateCoarseGraph(ctrl, graph, cnvtxs, match, perm);
422
423  idxwspacefree(ctrl, nvtxs);
424  idxwspacefree(ctrl, nvtxs);
425}
426
427
428
429
430
431/*************************************************************************
432* This function checks if v+u2 provides a better balance in the weight
433* vector that v+u1
434**************************************************************************/
435float BetterVBalance(int ncon, int norm, float *vwgt, float *u1wgt, float *u2wgt)
436{
437  int i;
438  float sum1, sum2, max1, max2, min1, min2, diff1, diff2;
439
440  if (norm == -1) {
441    max1 = min1 = vwgt[0]+u1wgt[0];
442    max2 = min2 = vwgt[0]+u2wgt[0];
443    sum1 = vwgt[0]+u1wgt[0];
444    sum2 = vwgt[0]+u2wgt[0];
445
446    for (i=1; i<ncon; i++) {
447      if (max1 < vwgt[i]+u1wgt[i])
448        max1 = vwgt[i]+u1wgt[i];
449      if (min1 > vwgt[i]+u1wgt[i])
450        min1 = vwgt[i]+u1wgt[i];
451
452      if (max2 < vwgt[i]+u2wgt[i])
453        max2 = vwgt[i]+u2wgt[i];
454      if (min2 > vwgt[i]+u2wgt[i])
455        min2 = vwgt[i]+u2wgt[i];
456
457      sum1 += vwgt[i]+u1wgt[i];
458      sum2 += vwgt[i]+u2wgt[i];
459    }
460
461    if (sum1 == 0.0)
462      return 1;
463    else if (sum2 == 0.0)
464      return -1;
465    else
466      return ((max1-min1)/sum1) - ((max2-min2)/sum2);
467  }
468  else if (norm == 1) {
469    sum1 = sum2 = 0.0;
470    for (i=0; i<ncon; i++) {
471      sum1 += vwgt[i]+u1wgt[i];
472      sum2 += vwgt[i]+u2wgt[i];
473    }
474    sum1 = sum1/(1.0*ncon);
475    sum2 = sum2/(1.0*ncon);
476
477    diff1 = diff2 = 0.0;
478    for (i=0; i<ncon; i++) {
479      diff1 += fabs(sum1 - (vwgt[i]+u1wgt[i]));
480      diff2 += fabs(sum2 - (vwgt[i]+u2wgt[i]));
481    }
482
483    return diff1 - diff2;
484  }
485  else {
486    errexit("Unknown norm: %d\n", norm);
487  }
488  return 0.0;
489}
490
491
492/*************************************************************************
493* This function checks if the vertex weights of two vertices are below
494* a given set of values
495**************************************************************************/
496int AreAllVwgtsBelowFast(int ncon, float *vwgt1, float *vwgt2, float limit)
497{
498  int i;
499
500  for (i=0; i<ncon; i++)
501    if (vwgt1[i] + vwgt2[i] > limit)
502      return 0;
503
504  return 1;
505}
506
Note: See TracBrowser for help on using the repository browser.