1 /*
2 * mmd.c
3 *
4 * **************************************************************
5 * The following C function was developed from a FORTRAN subroutine
6 * in SPARSPAK written by Eleanor Chu, Alan George, Joseph Liu
7 * and Esmond Ng.
8 *
9 * The FORTRAN-to-C transformation and modifications such as dynamic
10 * memory allocation and deallocation were performed by Chunguang
11 * Sun.
12 * **************************************************************
13 *
14 * Taken from SMMS, George 12/13/94
15 *
16 * The meaning of invperm, and perm vectors is different from that
17 * in genqmd_ of SparsPak
18 *
19 * $Id: mmd.c,v 1.1 1998/11/27 17:59:25 karypis Exp $
20 */
21
22 #include "metis.h"
23
24
25 /*************************************************************************
26 * genmmd -- multiple minimum external degree
27 * purpose -- this routine implements the minimum degree
28 * algorithm. it makes use of the implicit representation
29 * of elimination graphs by quotient graphs, and the notion
30 * of indistinguishable nodes. It also implements the modifications
31 * by multiple elimination and minimum external degree.
32 * Caution -- the adjacency vector adjncy will be destroyed.
33 * Input parameters --
34 * neqns -- number of equations.
35 * (xadj, adjncy) -- the adjacency structure.
36 * delta -- tolerance value for multiple elimination.
37 * maxint -- maximum machine representable (short) integer
38 * (any smaller estimate will do) for marking nodes.
39 * Output parameters --
40 * perm -- the minimum degree ordering.
41 * invp -- the inverse of perm.
42 * *ncsub -- an upper bound on the number of nonzero subscripts
43 * for the compressed storage scheme.
44 * Working parameters --
45 * head -- vector for head of degree lists.
46 * invp -- used temporarily for degree forward link.
47 * perm -- used temporarily for degree backward link.
48 * qsize -- vector for size of supernodes.
49 * list -- vector for temporary linked lists.
50 * marker -- a temporary marker vector.
51 * Subroutines used -- mmdelm, mmdint, mmdnum, mmdupd.
52 **************************************************************************/
genmmd(int neqns,idxtype * xadj,idxtype * adjncy,idxtype * invp,idxtype * perm,int delta,idxtype * head,idxtype * qsize,idxtype * list,idxtype * marker,int maxint,int * ncsub)53 void genmmd(int neqns, idxtype *xadj, idxtype *adjncy, idxtype *invp, idxtype *perm,
54 int delta, idxtype *head, idxtype *qsize, idxtype *list, idxtype *marker,
55 int maxint, int *ncsub)
56 {
57 int ehead, i, mdeg, mdlmt, mdeg_node, nextmd, num, tag;
58
59 if (neqns <= 0)
60 return;
61
62 /* Adjust from C to Fortran */
63 xadj--; adjncy--; invp--; perm--; head--; qsize--; list--; marker--;
64
65 /* initialization for the minimum degree algorithm. */
66 *ncsub = 0;
67 mmdint(neqns, xadj, adjncy, head, invp, perm, qsize, list, marker);
68
69 /* 'num' counts the number of ordered nodes plus 1. */
70 num = 1;
71
72 /* eliminate all isolated nodes. */
73 nextmd = head[1];
74 while (nextmd > 0) {
75 mdeg_node = nextmd;
76 nextmd = invp[mdeg_node];
77 marker[mdeg_node] = maxint;
78 invp[mdeg_node] = -num;
79 num = num + 1;
80 }
81
82 /* search for node of the minimum degree. 'mdeg' is the current */
83 /* minimum degree; 'tag' is used to facilitate marking nodes. */
84 if (num > neqns)
85 goto n1000;
86 tag = 1;
87 head[1] = 0;
88 mdeg = 2;
89
90 /* infinite loop here ! */
91 while (1) {
92 while (head[mdeg] <= 0)
93 mdeg++;
94
95 /* use value of 'delta' to set up 'mdlmt', which governs */
96 /* when a degree update is to be performed. */
97 mdlmt = mdeg + delta;
98 ehead = 0;
99
100 n500:
101 mdeg_node = head[mdeg];
102 while (mdeg_node <= 0) {
103 mdeg++;
104
105 if (mdeg > mdlmt)
106 goto n900;
107 mdeg_node = head[mdeg];
108 };
109
110 /* remove 'mdeg_node' from the degree structure. */
111 nextmd = invp[mdeg_node];
112 head[mdeg] = nextmd;
113 if (nextmd > 0)
114 perm[nextmd] = -mdeg;
115 invp[mdeg_node] = -num;
116 *ncsub += mdeg + qsize[mdeg_node] - 2;
117 if ((num+qsize[mdeg_node]) > neqns)
118 goto n1000;
119
120 /* eliminate 'mdeg_node' and perform quotient graph */
121 /* transformation. reset 'tag' value if necessary. */
122 tag++;
123 if (tag >= maxint) {
124 tag = 1;
125 for (i = 1; i <= neqns; i++)
126 if (marker[i] < maxint)
127 marker[i] = 0;
128 };
129
130 mmdelm(mdeg_node, xadj, adjncy, head, invp, perm, qsize, list, marker, maxint, tag);
131
132 num += qsize[mdeg_node];
133 list[mdeg_node] = ehead;
134 ehead = mdeg_node;
135 if (delta >= 0)
136 goto n500;
137
138 n900:
139 /* update degrees of the nodes involved in the */
140 /* minimum degree nodes elimination. */
141 if (num > neqns)
142 goto n1000;
143 mmdupd( ehead, neqns, xadj, adjncy, delta, &mdeg, head, invp, perm, qsize, list, marker, maxint, &tag);
144 }; /* end of -- while ( 1 ) -- */
145
146 n1000:
147 mmdnum( neqns, perm, invp, qsize );
148
149 /* Adjust from Fortran back to C*/
150 xadj++; adjncy++; invp++; perm++; head++; qsize++; list++; marker++;
151 }
152
153
154 /**************************************************************************
155 * mmdelm ...... multiple minimum degree elimination
156 * Purpose -- This routine eliminates the node mdeg_node of minimum degree
157 * from the adjacency structure, which is stored in the quotient
158 * graph format. It also transforms the quotient graph representation
159 * of the elimination graph.
160 * Input parameters --
161 * mdeg_node -- node of minimum degree.
162 * maxint -- estimate of maximum representable (short) integer.
163 * tag -- tag value.
164 * Updated parameters --
165 * (xadj, adjncy) -- updated adjacency structure.
166 * (head, forward, backward) -- degree doubly linked structure.
167 * qsize -- size of supernode.
168 * marker -- marker vector.
169 * list -- temporary linked list of eliminated nabors.
170 ***************************************************************************/
mmdelm(int mdeg_node,idxtype * xadj,idxtype * adjncy,idxtype * head,idxtype * forward,idxtype * backward,idxtype * qsize,idxtype * list,idxtype * marker,int maxint,int tag)171 void mmdelm(int mdeg_node, idxtype *xadj, idxtype *adjncy, idxtype *head, idxtype *forward,
172 idxtype *backward, idxtype *qsize, idxtype *list, idxtype *marker, int maxint,int tag)
173 {
174 int element, i, istop, istart, j,
175 jstop, jstart, link,
176 nabor, node, npv, nqnbrs, nxnode,
177 pvnode, rlmt, rloc, rnode, xqnbr;
178
179 /* find the reachable set of 'mdeg_node' and */
180 /* place it in the data structure. */
181 marker[mdeg_node] = tag;
182 istart = xadj[mdeg_node];
183 istop = xadj[mdeg_node+1] - 1;
184
185 /* 'element' points to the beginning of the list of */
186 /* eliminated nabors of 'mdeg_node', and 'rloc' gives the */
187 /* storage location for the next reachable node. */
188 element = 0;
189 rloc = istart;
190 rlmt = istop;
191 for ( i = istart; i <= istop; i++ ) {
192 nabor = adjncy[i];
193 if ( nabor == 0 ) break;
194 if ( marker[nabor] < tag ) {
195 marker[nabor] = tag;
196 if ( forward[nabor] < 0 ) {
197 list[nabor] = element;
198 element = nabor;
199 } else {
200 adjncy[rloc] = nabor;
201 rloc++;
202 };
203 }; /* end of -- if -- */
204 }; /* end of -- for -- */
205
206 /* merge with reachable nodes from generalized elements. */
207 while ( element > 0 ) {
208 adjncy[rlmt] = -element;
209 link = element;
210
211 n400:
212 jstart = xadj[link];
213 jstop = xadj[link+1] - 1;
214 for ( j = jstart; j <= jstop; j++ ) {
215 node = adjncy[j];
216 link = -node;
217 if ( node < 0 ) goto n400;
218 if ( node == 0 ) break;
219 if ((marker[node]<tag)&&(forward[node]>=0)) {
220 marker[node] = tag;
221 /*use storage from eliminated nodes if necessary.*/
222 while ( rloc >= rlmt ) {
223 link = -adjncy[rlmt];
224 rloc = xadj[link];
225 rlmt = xadj[link+1] - 1;
226 };
227 adjncy[rloc] = node;
228 rloc++;
229 };
230 }; /* end of -- for ( j = jstart; -- */
231 element = list[element];
232 }; /* end of -- while ( element > 0 ) -- */
233 if ( rloc <= rlmt ) adjncy[rloc] = 0;
234 /* for each node in the reachable set, do the following. */
235 link = mdeg_node;
236
237 n1100:
238 istart = xadj[link];
239 istop = xadj[link+1] - 1;
240 for ( i = istart; i <= istop; i++ ) {
241 rnode = adjncy[i];
242 link = -rnode;
243 if ( rnode < 0 ) goto n1100;
244 if ( rnode == 0 ) return;
245
246 /* 'rnode' is in the degree list structure. */
247 pvnode = backward[rnode];
248 if (( pvnode != 0 ) && ( pvnode != (-maxint) )) {
249 /* then remove 'rnode' from the structure. */
250 nxnode = forward[rnode];
251 if ( nxnode > 0 ) backward[nxnode] = pvnode;
252 if ( pvnode > 0 ) forward[pvnode] = nxnode;
253 npv = -pvnode;
254 if ( pvnode < 0 ) head[npv] = nxnode;
255 };
256
257 /* purge inactive quotient nabors of 'rnode'. */
258 jstart = xadj[rnode];
259 jstop = xadj[rnode+1] - 1;
260 xqnbr = jstart;
261 for ( j = jstart; j <= jstop; j++ ) {
262 nabor = adjncy[j];
263 if ( nabor == 0 ) break;
264 if ( marker[nabor] < tag ) {
265 adjncy[xqnbr] = nabor;
266 xqnbr++;
267 };
268 };
269
270 /* no active nabor after the purging. */
271 nqnbrs = xqnbr - jstart;
272 if ( nqnbrs <= 0 ) {
273 /* merge 'rnode' with 'mdeg_node'. */
274 qsize[mdeg_node] += qsize[rnode];
275 qsize[rnode] = 0;
276 marker[rnode] = maxint;
277 forward[rnode] = -mdeg_node;
278 backward[rnode] = -maxint;
279 } else {
280 /* flag 'rnode' for degree update, and */
281 /* add 'mdeg_node' as a nabor of 'rnode'. */
282 forward[rnode] = nqnbrs + 1;
283 backward[rnode] = 0;
284 adjncy[xqnbr] = mdeg_node;
285 xqnbr++;
286 if ( xqnbr <= jstop ) adjncy[xqnbr] = 0;
287 };
288 }; /* end of -- for ( i = istart; -- */
289 return;
290 }
291
292 /***************************************************************************
293 * mmdint ---- mult minimum degree initialization
294 * purpose -- this routine performs initialization for the
295 * multiple elimination version of the minimum degree algorithm.
296 * input parameters --
297 * neqns -- number of equations.
298 * (xadj, adjncy) -- adjacency structure.
299 * output parameters --
300 * (head, dfrow, backward) -- degree doubly linked structure.
301 * qsize -- size of supernode ( initialized to one).
302 * list -- linked list.
303 * marker -- marker vector.
304 ****************************************************************************/
mmdint(int neqns,idxtype * xadj,idxtype * adjncy,idxtype * head,idxtype * forward,idxtype * backward,idxtype * qsize,idxtype * list,idxtype * marker)305 int mmdint(int neqns, idxtype *xadj, idxtype *adjncy, idxtype *head, idxtype *forward,
306 idxtype *backward, idxtype *qsize, idxtype *list, idxtype *marker)
307 {
308 int fnode, ndeg, node;
309
310 for ( node = 1; node <= neqns; node++ ) {
311 head[node] = 0;
312 qsize[node] = 1;
313 marker[node] = 0;
314 list[node] = 0;
315 };
316
317 /* initialize the degree doubly linked lists. */
318 for ( node = 1; node <= neqns; node++ ) {
319 ndeg = xadj[node+1] - xadj[node]/* + 1*/; /* george */
320 if (ndeg == 0)
321 ndeg = 1;
322 fnode = head[ndeg];
323 forward[node] = fnode;
324 head[ndeg] = node;
325 if ( fnode > 0 ) backward[fnode] = node;
326 backward[node] = -ndeg;
327 };
328 return 0;
329 }
330
331 /****************************************************************************
332 * mmdnum --- multi minimum degree numbering
333 * purpose -- this routine performs the final step in producing
334 * the permutation and inverse permutation vectors in the
335 * multiple elimination version of the minimum degree
336 * ordering algorithm.
337 * input parameters --
338 * neqns -- number of equations.
339 * qsize -- size of supernodes at elimination.
340 * updated parameters --
341 * invp -- inverse permutation vector. on input,
342 * if qsize[node] = 0, then node has been merged
343 * into the node -invp[node]; otherwise,
344 * -invp[node] is its inverse labelling.
345 * output parameters --
346 * perm -- the permutation vector.
347 ****************************************************************************/
mmdnum(int neqns,idxtype * perm,idxtype * invp,idxtype * qsize)348 void mmdnum(int neqns, idxtype *perm, idxtype *invp, idxtype *qsize)
349 {
350 int father, nextf, node, nqsize, num, root;
351
352 for ( node = 1; node <= neqns; node++ ) {
353 nqsize = qsize[node];
354 if ( nqsize <= 0 ) perm[node] = invp[node];
355 if ( nqsize > 0 ) perm[node] = -invp[node];
356 };
357
358 /* for each node which has been merged, do the following. */
359 for ( node = 1; node <= neqns; node++ ) {
360 if ( perm[node] <= 0 ) {
361
362 /* trace the merged tree until one which has not */
363 /* been merged, call it root. */
364 father = node;
365 while ( perm[father] <= 0 )
366 father = - perm[father];
367
368 /* number node after root. */
369 root = father;
370 num = perm[root] + 1;
371 invp[node] = -num;
372 perm[root] = num;
373
374 /* shorten the merged tree. */
375 father = node;
376 nextf = - perm[father];
377 while ( nextf > 0 ) {
378 perm[father] = -root;
379 father = nextf;
380 nextf = -perm[father];
381 };
382 }; /* end of -- if ( perm[node] <= 0 ) -- */
383 }; /* end of -- for ( node = 1; -- */
384
385 /* ready to compute perm. */
386 for ( node = 1; node <= neqns; node++ ) {
387 num = -invp[node];
388 invp[node] = num;
389 perm[num] = node;
390 };
391 return;
392 }
393
394 /****************************************************************************
395 * mmdupd ---- multiple minimum degree update
396 * purpose -- this routine updates the degrees of nodes after a
397 * multiple elimination step.
398 * input parameters --
399 * ehead -- the beginning of the list of eliminated nodes
400 * (i.e., newly formed elements).
401 * neqns -- number of equations.
402 * (xadj, adjncy) -- adjacency structure.
403 * delta -- tolerance value for multiple elimination.
404 * maxint -- maximum machine representable (short) integer.
405 * updated parameters --
406 * mdeg -- new minimum degree after degree update.
407 * (head, forward, backward) -- degree doubly linked structure.
408 * qsize -- size of supernode.
409 * list -- marker vector for degree update.
410 * *tag -- tag value.
411 ****************************************************************************/
mmdupd(int ehead,int neqns,idxtype * xadj,idxtype * adjncy,int delta,int * mdeg,idxtype * head,idxtype * forward,idxtype * backward,idxtype * qsize,idxtype * list,idxtype * marker,int maxint,int * tag)412 void mmdupd(int ehead, int neqns, idxtype *xadj, idxtype *adjncy, int delta, int *mdeg,
413 idxtype *head, idxtype *forward, idxtype *backward, idxtype *qsize, idxtype *list,
414 idxtype *marker, int maxint,int *tag)
415 {
416 int deg, deg0, element, enode, fnode, i, iq2, istop,
417 istart, j, jstop, jstart, link, mdeg0, mtag, nabor,
418 node, q2head, qxhead;
419
420 mdeg0 = *mdeg + delta;
421 element = ehead;
422
423 n100:
424 if ( element <= 0 ) return;
425
426 /* for each of the newly formed element, do the following. */
427 /* reset tag value if necessary. */
428 mtag = *tag + mdeg0;
429 if ( mtag >= maxint ) {
430 *tag = 1;
431 for ( i = 1; i <= neqns; i++ )
432 if ( marker[i] < maxint ) marker[i] = 0;
433 mtag = *tag + mdeg0;
434 };
435
436 /* create two linked lists from nodes associated with 'element': */
437 /* one with two nabors (q2head) in the adjacency structure, and the*/
438 /* other with more than two nabors (qxhead). also compute 'deg0',*/
439 /* number of nodes in this element. */
440 q2head = 0;
441 qxhead = 0;
442 deg0 = 0;
443 link =element;
444
445 n400:
446 istart = xadj[link];
447 istop = xadj[link+1] - 1;
448 for ( i = istart; i <= istop; i++ ) {
449 enode = adjncy[i];
450 link = -enode;
451 if ( enode < 0 ) goto n400;
452 if ( enode == 0 ) break;
453 if ( qsize[enode] != 0 ) {
454 deg0 += qsize[enode];
455 marker[enode] = mtag;
456
457 /*'enode' requires a degree update*/
458 if ( backward[enode] == 0 ) {
459 /* place either in qxhead or q2head list. */
460 if ( forward[enode] != 2 ) {
461 list[enode] = qxhead;
462 qxhead = enode;
463 } else {
464 list[enode] = q2head;
465 q2head = enode;
466 };
467 };
468 }; /* enf of -- if ( qsize[enode] != 0 ) -- */
469 }; /* end of -- for ( i = istart; -- */
470
471 /* for each node in q2 list, do the following. */
472 enode = q2head;
473 iq2 = 1;
474
475 n900:
476 if ( enode <= 0 ) goto n1500;
477 if ( backward[enode] != 0 ) goto n2200;
478 (*tag)++;
479 deg = deg0;
480
481 /* identify the other adjacent element nabor. */
482 istart = xadj[enode];
483 nabor = adjncy[istart];
484 if ( nabor == element ) nabor = adjncy[istart+1];
485 link = nabor;
486 if ( forward[nabor] >= 0 ) {
487 /* nabor is uneliminated, increase degree count. */
488 deg += qsize[nabor];
489 goto n2100;
490 };
491
492 /* the nabor is eliminated. for each node in the 2nd element */
493 /* do the following. */
494 n1000:
495 istart = xadj[link];
496 istop = xadj[link+1] - 1;
497 for ( i = istart; i <= istop; i++ ) {
498 node = adjncy[i];
499 link = -node;
500 if ( node != enode ) {
501 if ( node < 0 ) goto n1000;
502 if ( node == 0 ) goto n2100;
503 if ( qsize[node] != 0 ) {
504 if ( marker[node] < *tag ) {
505 /* 'node' is not yet considered. */
506 marker[node] = *tag;
507 deg += qsize[node];
508 } else {
509 if ( backward[node] == 0 ) {
510 if ( forward[node] == 2 ) {
511 /* 'node' is indistinguishable from 'enode'.*/
512 /* merge them into a new supernode. */
513 qsize[enode] += qsize[node];
514 qsize[node] = 0;
515 marker[node] = maxint;
516 forward[node] = -enode;
517 backward[node] = -maxint;
518 } else {
519 /* 'node' is outmacthed by 'enode' */
520 if (backward[node]==0) backward[node] = -maxint;
521 };
522 }; /* end of -- if ( backward[node] == 0 ) -- */
523 }; /* end of -- if ( marker[node] < *tag ) -- */
524 }; /* end of -- if ( qsize[node] != 0 ) -- */
525 }; /* end of -- if ( node != enode ) -- */
526 }; /* end of -- for ( i = istart; -- */
527 goto n2100;
528
529 n1500:
530 /* for each 'enode' in the 'qx' list, do the following. */
531 enode = qxhead;
532 iq2 = 0;
533
534 n1600: if ( enode <= 0 ) goto n2300;
535 if ( backward[enode] != 0 ) goto n2200;
536 (*tag)++;
537 deg = deg0;
538
539 /*for each unmarked nabor of 'enode', do the following.*/
540 istart = xadj[enode];
541 istop = xadj[enode+1] - 1;
542 for ( i = istart; i <= istop; i++ ) {
543 nabor = adjncy[i];
544 if ( nabor == 0 ) break;
545 if ( marker[nabor] < *tag ) {
546 marker[nabor] = *tag;
547 link = nabor;
548 if ( forward[nabor] >= 0 )
549 /*if uneliminated, include it in deg count.*/
550 deg += qsize[nabor];
551 else {
552 n1700:
553 /* if eliminated, include unmarked nodes in this*/
554 /* element into the degree count. */
555 jstart = xadj[link];
556 jstop = xadj[link+1] - 1;
557 for ( j = jstart; j <= jstop; j++ ) {
558 node = adjncy[j];
559 link = -node;
560 if ( node < 0 ) goto n1700;
561 if ( node == 0 ) break;
562 if ( marker[node] < *tag ) {
563 marker[node] = *tag;
564 deg += qsize[node];
565 };
566 }; /* end of -- for ( j = jstart; -- */
567 }; /* end of -- if ( forward[nabor] >= 0 ) -- */
568 }; /* end of -- if ( marker[nabor] < *tag ) -- */
569 }; /* end of -- for ( i = istart; -- */
570
571 n2100:
572 /* update external degree of 'enode' in degree structure, */
573 /* and '*mdeg' if necessary. */
574 deg = deg - qsize[enode] + 1;
575 fnode = head[deg];
576 forward[enode] = fnode;
577 backward[enode] = -deg;
578 if ( fnode > 0 ) backward[fnode] = enode;
579 head[deg] = enode;
580 if ( deg < *mdeg ) *mdeg = deg;
581
582 n2200:
583 /* get next enode in current element. */
584 enode = list[enode];
585 if ( iq2 == 1 ) goto n900;
586 goto n1600;
587
588 n2300:
589 /* get next element in the list. */
590 *tag = mtag;
591 element = list[element];
592 goto n100;
593 }
594