source: branches/numpy/pymetis/metis-4.0/Lib/mmd.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: 20.7 KB
Line 
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**************************************************************************/
53void 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
100n500:
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
146n1000:
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***************************************************************************/
171void 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
211n400:
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
237n1100:
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****************************************************************************/
305int  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****************************************************************************/
348void 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****************************************************************************/
412void 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
423n100:
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
445n400:
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
475n900:
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.                                         */
494n1000:
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
529n1500:
530          /* for each 'enode' in the 'qx' list, do the following. */
531          enode = qxhead;
532          iq2 = 0;
533
534n1600:    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 {
552n1700:
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
571n2100:
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
582n2200:
583          /* get next enode in current element. */
584          enode = list[enode];
585          if ( iq2 == 1 ) goto n900;
586          goto n1600;
587
588n2300:
589          /* get next element in the list. */
590          *tag = mtag;
591          element = list[element];
592          goto n100;
593    }
Note: See TracBrowser for help on using the repository browser.