mmd.c
Go to the documentation of this file.
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 5993 2009-01-07 02:09:57Z karypis $
20  */
21 
22 #include "metislib.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 **************************************************************************/
53 void genmmd(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *invp, idx_t *perm,
54  idx_t delta, idx_t *head, idx_t *qsize, idx_t *list, idx_t *marker,
55  idx_t maxint, idx_t *ncsub)
56 {
57  idx_t 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 ***************************************************************************/
171 void mmdelm(idx_t mdeg_node, idx_t *xadj, idx_t *adjncy, idx_t *head, idx_t *forward,
172  idx_t *backward, idx_t *qsize, idx_t *list, idx_t *marker, idx_t maxint, idx_t tag)
173 {
174  idx_t 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 ****************************************************************************/
305 idx_t mmdint(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *head, idx_t *forward,
306  idx_t *backward, idx_t *qsize, idx_t *list, idx_t *marker)
307 {
308  idx_t 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 ****************************************************************************/
348 void mmdnum(idx_t neqns, idx_t *perm, idx_t *invp, idx_t *qsize)
349 {
350  idx_t 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 ****************************************************************************/
412 void mmdupd(idx_t ehead, idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t delta, idx_t *mdeg,
413  idx_t *head, idx_t *forward, idx_t *backward, idx_t *qsize, idx_t *list,
414  idx_t *marker, idx_t maxint, idx_t *tag)
415 {
416  idx_t 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  }
genmmd
void genmmd(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *invp, idx_t *perm, idx_t delta, idx_t *head, idx_t *qsize, idx_t *list, idx_t *marker, idx_t maxint, idx_t *ncsub)
Definition: mmd.c:53
xadj
idx_t idx_t * xadj
Definition: include/metis.h:197
mpfr::root
const mpreal root(const mpreal &x, unsigned long int k, mp_rnd_t r=mpreal::get_default_rnd())
Definition: mpreal.h:2194
idx_t
int32_t idx_t
Definition: include/metis.h:101
perm
idx_t idx_t idx_t idx_t idx_t * perm
Definition: include/metis.h:222
mmdelm
void mmdelm(idx_t mdeg_node, idx_t *xadj, idx_t *adjncy, idx_t *head, idx_t *forward, idx_t *backward, idx_t *qsize, idx_t *list, idx_t *marker, idx_t maxint, idx_t tag)
Definition: mmd.c:171
list
Definition: pytypes.h:1301
head
EIGEN_DEVICE_FUNC SegmentReturnType head(Index n)
This is the const version of head(Index).
Definition: BlockMethods.h:919
adjncy
idx_t idx_t idx_t * adjncy
Definition: include/metis.h:197
mmdnum
void mmdnum(idx_t neqns, idx_t *perm, idx_t *invp, idx_t *qsize)
Definition: mmd.c:348
mmdupd
void mmdupd(idx_t ehead, idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t delta, idx_t *mdeg, idx_t *head, idx_t *forward, idx_t *backward, idx_t *qsize, idx_t *list, idx_t *marker, idx_t maxint, idx_t *tag)
Definition: mmd.c:412
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
mmdint
idx_t mmdint(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *head, idx_t *forward, idx_t *backward, idx_t *qsize, idx_t *list, idx_t *marker)
Definition: mmd.c:305
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:01