/* * mmd.c * * ************************************************************** * The following C function was developed from a FORTRAN subroutine * in SPARSPAK written by Eleanor Chu, Alan George, Joseph Liu * and Esmond Ng. * * The FORTRAN-to-C transformation and modifications such as dynamic * memory allocation and deallocation were performed by Chunguang * Sun. * ************************************************************** * * Taken from SMMS, George 12/13/94 * * The meaning of invperm, and perm vectors is different from that * in genqmd_ of SparsPak * * $Id: mmd.c 5993 2009-01-07 02:09:57Z karypis $ */ #include "metislib.h" /************************************************************************* * genmmd -- multiple minimum external degree * purpose -- this routine implements the minimum degree * algorithm. it makes use of the implicit representation * of elimination graphs by quotient graphs, and the notion * of indistinguishable nodes. It also implements the modifications * by multiple elimination and minimum external degree. * Caution -- the adjacency vector adjncy will be destroyed. * Input parameters -- * neqns -- number of equations. * (xadj, adjncy) -- the adjacency structure. * delta -- tolerance value for multiple elimination. * maxint -- maximum machine representable (short) integer * (any smaller estimate will do) for marking nodes. * Output parameters -- * perm -- the minimum degree ordering. * invp -- the inverse of perm. * *ncsub -- an upper bound on the number of nonzero subscripts * for the compressed storage scheme. * Working parameters -- * head -- vector for head of degree lists. * invp -- used temporarily for degree forward link. * perm -- used temporarily for degree backward link. * qsize -- vector for size of supernodes. * list -- vector for temporary linked lists. * marker -- a temporary marker vector. * Subroutines used -- mmdelm, mmdint, mmdnum, mmdupd. **************************************************************************/ 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) { idx_t ehead, i, mdeg, mdlmt, mdeg_node, nextmd, num, tag; if (neqns <= 0) return; /* Adjust from C to Fortran */ xadj--; adjncy--; invp--; perm--; head--; qsize--; list--; marker--; /* initialization for the minimum degree algorithm. */ *ncsub = 0; mmdint(neqns, xadj, adjncy, head, invp, perm, qsize, list, marker); /* 'num' counts the number of ordered nodes plus 1. */ num = 1; /* eliminate all isolated nodes. */ nextmd = head[1]; while (nextmd > 0) { mdeg_node = nextmd; nextmd = invp[mdeg_node]; marker[mdeg_node] = maxint; invp[mdeg_node] = -num; num = num + 1; } /* search for node of the minimum degree. 'mdeg' is the current */ /* minimum degree; 'tag' is used to facilitate marking nodes. */ if (num > neqns) goto n1000; tag = 1; head[1] = 0; mdeg = 2; /* infinite loop here ! */ while (1) { while (head[mdeg] <= 0) mdeg++; /* use value of 'delta' to set up 'mdlmt', which governs */ /* when a degree update is to be performed. */ mdlmt = mdeg + delta; ehead = 0; n500: mdeg_node = head[mdeg]; while (mdeg_node <= 0) { mdeg++; if (mdeg > mdlmt) goto n900; mdeg_node = head[mdeg]; }; /* remove 'mdeg_node' from the degree structure. */ nextmd = invp[mdeg_node]; head[mdeg] = nextmd; if (nextmd > 0) perm[nextmd] = -mdeg; invp[mdeg_node] = -num; *ncsub += mdeg + qsize[mdeg_node] - 2; if ((num+qsize[mdeg_node]) > neqns) goto n1000; /* eliminate 'mdeg_node' and perform quotient graph */ /* transformation. reset 'tag' value if necessary. */ tag++; if (tag >= maxint) { tag = 1; for (i = 1; i <= neqns; i++) if (marker[i] < maxint) marker[i] = 0; }; mmdelm(mdeg_node, xadj, adjncy, head, invp, perm, qsize, list, marker, maxint, tag); num += qsize[mdeg_node]; list[mdeg_node] = ehead; ehead = mdeg_node; if (delta >= 0) goto n500; n900: /* update degrees of the nodes involved in the */ /* minimum degree nodes elimination. */ if (num > neqns) goto n1000; mmdupd( ehead, neqns, xadj, adjncy, delta, &mdeg, head, invp, perm, qsize, list, marker, maxint, &tag); }; /* end of -- while ( 1 ) -- */ n1000: mmdnum( neqns, perm, invp, qsize ); /* Adjust from Fortran back to C*/ xadj++; adjncy++; invp++; perm++; head++; qsize++; list++; marker++; } /************************************************************************** * mmdelm ...... multiple minimum degree elimination * Purpose -- This routine eliminates the node mdeg_node of minimum degree * from the adjacency structure, which is stored in the quotient * graph format. It also transforms the quotient graph representation * of the elimination graph. * Input parameters -- * mdeg_node -- node of minimum degree. * maxint -- estimate of maximum representable (short) integer. * tag -- tag value. * Updated parameters -- * (xadj, adjncy) -- updated adjacency structure. * (head, forward, backward) -- degree doubly linked structure. * qsize -- size of supernode. * marker -- marker vector. * list -- temporary linked list of eliminated nabors. ***************************************************************************/ 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) { idx_t element, i, istop, istart, j, jstop, jstart, link, nabor, node, npv, nqnbrs, nxnode, pvnode, rlmt, rloc, rnode, xqnbr; /* find the reachable set of 'mdeg_node' and */ /* place it in the data structure. */ marker[mdeg_node] = tag; istart = xadj[mdeg_node]; istop = xadj[mdeg_node+1] - 1; /* 'element' points to the beginning of the list of */ /* eliminated nabors of 'mdeg_node', and 'rloc' gives the */ /* storage location for the next reachable node. */ element = 0; rloc = istart; rlmt = istop; for ( i = istart; i <= istop; i++ ) { nabor = adjncy[i]; if ( nabor == 0 ) break; if ( marker[nabor] < tag ) { marker[nabor] = tag; if ( forward[nabor] < 0 ) { list[nabor] = element; element = nabor; } else { adjncy[rloc] = nabor; rloc++; }; }; /* end of -- if -- */ }; /* end of -- for -- */ /* merge with reachable nodes from generalized elements. */ while ( element > 0 ) { adjncy[rlmt] = -element; link = element; n400: jstart = xadj[link]; jstop = xadj[link+1] - 1; for ( j = jstart; j <= jstop; j++ ) { node = adjncy[j]; link = -node; if ( node < 0 ) goto n400; if ( node == 0 ) break; if ((marker[node]=0)) { marker[node] = tag; /*use storage from eliminated nodes if necessary.*/ while ( rloc >= rlmt ) { link = -adjncy[rlmt]; rloc = xadj[link]; rlmt = xadj[link+1] - 1; }; adjncy[rloc] = node; rloc++; }; }; /* end of -- for ( j = jstart; -- */ element = list[element]; }; /* end of -- while ( element > 0 ) -- */ if ( rloc <= rlmt ) adjncy[rloc] = 0; /* for each node in the reachable set, do the following. */ link = mdeg_node; n1100: istart = xadj[link]; istop = xadj[link+1] - 1; for ( i = istart; i <= istop; i++ ) { rnode = adjncy[i]; link = -rnode; if ( rnode < 0 ) goto n1100; if ( rnode == 0 ) return; /* 'rnode' is in the degree list structure. */ pvnode = backward[rnode]; if (( pvnode != 0 ) && ( pvnode != (-maxint) )) { /* then remove 'rnode' from the structure. */ nxnode = forward[rnode]; if ( nxnode > 0 ) backward[nxnode] = pvnode; if ( pvnode > 0 ) forward[pvnode] = nxnode; npv = -pvnode; if ( pvnode < 0 ) head[npv] = nxnode; }; /* purge inactive quotient nabors of 'rnode'. */ jstart = xadj[rnode]; jstop = xadj[rnode+1] - 1; xqnbr = jstart; for ( j = jstart; j <= jstop; j++ ) { nabor = adjncy[j]; if ( nabor == 0 ) break; if ( marker[nabor] < tag ) { adjncy[xqnbr] = nabor; xqnbr++; }; }; /* no active nabor after the purging. */ nqnbrs = xqnbr - jstart; if ( nqnbrs <= 0 ) { /* merge 'rnode' with 'mdeg_node'. */ qsize[mdeg_node] += qsize[rnode]; qsize[rnode] = 0; marker[rnode] = maxint; forward[rnode] = -mdeg_node; backward[rnode] = -maxint; } else { /* flag 'rnode' for degree update, and */ /* add 'mdeg_node' as a nabor of 'rnode'. */ forward[rnode] = nqnbrs + 1; backward[rnode] = 0; adjncy[xqnbr] = mdeg_node; xqnbr++; if ( xqnbr <= jstop ) adjncy[xqnbr] = 0; }; }; /* end of -- for ( i = istart; -- */ return; } /*************************************************************************** * mmdint ---- mult minimum degree initialization * purpose -- this routine performs initialization for the * multiple elimination version of the minimum degree algorithm. * input parameters -- * neqns -- number of equations. * (xadj, adjncy) -- adjacency structure. * output parameters -- * (head, dfrow, backward) -- degree doubly linked structure. * qsize -- size of supernode ( initialized to one). * list -- linked list. * marker -- marker vector. ****************************************************************************/ 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) { idx_t fnode, ndeg, node; for ( node = 1; node <= neqns; node++ ) { head[node] = 0; qsize[node] = 1; marker[node] = 0; list[node] = 0; }; /* initialize the degree doubly linked lists. */ for ( node = 1; node <= neqns; node++ ) { ndeg = xadj[node+1] - xadj[node]/* + 1*/; /* george */ if (ndeg == 0) ndeg = 1; fnode = head[ndeg]; forward[node] = fnode; head[ndeg] = node; if ( fnode > 0 ) backward[fnode] = node; backward[node] = -ndeg; }; return 0; } /**************************************************************************** * mmdnum --- multi minimum degree numbering * purpose -- this routine performs the final step in producing * the permutation and inverse permutation vectors in the * multiple elimination version of the minimum degree * ordering algorithm. * input parameters -- * neqns -- number of equations. * qsize -- size of supernodes at elimination. * updated parameters -- * invp -- inverse permutation vector. on input, * if qsize[node] = 0, then node has been merged * into the node -invp[node]; otherwise, * -invp[node] is its inverse labelling. * output parameters -- * perm -- the permutation vector. ****************************************************************************/ void mmdnum(idx_t neqns, idx_t *perm, idx_t *invp, idx_t *qsize) { idx_t father, nextf, node, nqsize, num, root; for ( node = 1; node <= neqns; node++ ) { nqsize = qsize[node]; if ( nqsize <= 0 ) perm[node] = invp[node]; if ( nqsize > 0 ) perm[node] = -invp[node]; }; /* for each node which has been merged, do the following. */ for ( node = 1; node <= neqns; node++ ) { if ( perm[node] <= 0 ) { /* trace the merged tree until one which has not */ /* been merged, call it root. */ father = node; while ( perm[father] <= 0 ) father = - perm[father]; /* number node after root. */ root = father; num = perm[root] + 1; invp[node] = -num; perm[root] = num; /* shorten the merged tree. */ father = node; nextf = - perm[father]; while ( nextf > 0 ) { perm[father] = -root; father = nextf; nextf = -perm[father]; }; }; /* end of -- if ( perm[node] <= 0 ) -- */ }; /* end of -- for ( node = 1; -- */ /* ready to compute perm. */ for ( node = 1; node <= neqns; node++ ) { num = -invp[node]; invp[node] = num; perm[num] = node; }; return; } /**************************************************************************** * mmdupd ---- multiple minimum degree update * purpose -- this routine updates the degrees of nodes after a * multiple elimination step. * input parameters -- * ehead -- the beginning of the list of eliminated nodes * (i.e., newly formed elements). * neqns -- number of equations. * (xadj, adjncy) -- adjacency structure. * delta -- tolerance value for multiple elimination. * maxint -- maximum machine representable (short) integer. * updated parameters -- * mdeg -- new minimum degree after degree update. * (head, forward, backward) -- degree doubly linked structure. * qsize -- size of supernode. * list -- marker vector for degree update. * *tag -- tag value. ****************************************************************************/ 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) { idx_t deg, deg0, element, enode, fnode, i, iq2, istop, istart, j, jstop, jstart, link, mdeg0, mtag, nabor, node, q2head, qxhead; mdeg0 = *mdeg + delta; element = ehead; n100: if ( element <= 0 ) return; /* for each of the newly formed element, do the following. */ /* reset tag value if necessary. */ mtag = *tag + mdeg0; if ( mtag >= maxint ) { *tag = 1; for ( i = 1; i <= neqns; i++ ) if ( marker[i] < maxint ) marker[i] = 0; mtag = *tag + mdeg0; }; /* create two linked lists from nodes associated with 'element': */ /* one with two nabors (q2head) in the adjacency structure, and the*/ /* other with more than two nabors (qxhead). also compute 'deg0',*/ /* number of nodes in this element. */ q2head = 0; qxhead = 0; deg0 = 0; link =element; n400: istart = xadj[link]; istop = xadj[link+1] - 1; for ( i = istart; i <= istop; i++ ) { enode = adjncy[i]; link = -enode; if ( enode < 0 ) goto n400; if ( enode == 0 ) break; if ( qsize[enode] != 0 ) { deg0 += qsize[enode]; marker[enode] = mtag; /*'enode' requires a degree update*/ if ( backward[enode] == 0 ) { /* place either in qxhead or q2head list. */ if ( forward[enode] != 2 ) { list[enode] = qxhead; qxhead = enode; } else { list[enode] = q2head; q2head = enode; }; }; }; /* enf of -- if ( qsize[enode] != 0 ) -- */ }; /* end of -- for ( i = istart; -- */ /* for each node in q2 list, do the following. */ enode = q2head; iq2 = 1; n900: if ( enode <= 0 ) goto n1500; if ( backward[enode] != 0 ) goto n2200; (*tag)++; deg = deg0; /* identify the other adjacent element nabor. */ istart = xadj[enode]; nabor = adjncy[istart]; if ( nabor == element ) nabor = adjncy[istart+1]; link = nabor; if ( forward[nabor] >= 0 ) { /* nabor is uneliminated, increase degree count. */ deg += qsize[nabor]; goto n2100; }; /* the nabor is eliminated. for each node in the 2nd element */ /* do the following. */ n1000: istart = xadj[link]; istop = xadj[link+1] - 1; for ( i = istart; i <= istop; i++ ) { node = adjncy[i]; link = -node; if ( node != enode ) { if ( node < 0 ) goto n1000; if ( node == 0 ) goto n2100; if ( qsize[node] != 0 ) { if ( marker[node] < *tag ) { /* 'node' is not yet considered. */ marker[node] = *tag; deg += qsize[node]; } else { if ( backward[node] == 0 ) { if ( forward[node] == 2 ) { /* 'node' is indistinguishable from 'enode'.*/ /* merge them into a new supernode. */ qsize[enode] += qsize[node]; qsize[node] = 0; marker[node] = maxint; forward[node] = -enode; backward[node] = -maxint; } else { /* 'node' is outmacthed by 'enode' */ if (backward[node]==0) backward[node] = -maxint; }; }; /* end of -- if ( backward[node] == 0 ) -- */ }; /* end of -- if ( marker[node] < *tag ) -- */ }; /* end of -- if ( qsize[node] != 0 ) -- */ }; /* end of -- if ( node != enode ) -- */ }; /* end of -- for ( i = istart; -- */ goto n2100; n1500: /* for each 'enode' in the 'qx' list, do the following. */ enode = qxhead; iq2 = 0; n1600: if ( enode <= 0 ) goto n2300; if ( backward[enode] != 0 ) goto n2200; (*tag)++; deg = deg0; /*for each unmarked nabor of 'enode', do the following.*/ istart = xadj[enode]; istop = xadj[enode+1] - 1; for ( i = istart; i <= istop; i++ ) { nabor = adjncy[i]; if ( nabor == 0 ) break; if ( marker[nabor] < *tag ) { marker[nabor] = *tag; link = nabor; if ( forward[nabor] >= 0 ) /*if uneliminated, include it in deg count.*/ deg += qsize[nabor]; else { n1700: /* if eliminated, include unmarked nodes in this*/ /* element into the degree count. */ jstart = xadj[link]; jstop = xadj[link+1] - 1; for ( j = jstart; j <= jstop; j++ ) { node = adjncy[j]; link = -node; if ( node < 0 ) goto n1700; if ( node == 0 ) break; if ( marker[node] < *tag ) { marker[node] = *tag; deg += qsize[node]; }; }; /* end of -- for ( j = jstart; -- */ }; /* end of -- if ( forward[nabor] >= 0 ) -- */ }; /* end of -- if ( marker[nabor] < *tag ) -- */ }; /* end of -- for ( i = istart; -- */ n2100: /* update external degree of 'enode' in degree structure, */ /* and '*mdeg' if necessary. */ deg = deg - qsize[enode] + 1; fnode = head[deg]; forward[enode] = fnode; backward[enode] = -deg; if ( fnode > 0 ) backward[fnode] = enode; head[deg] = enode; if ( deg < *mdeg ) *mdeg = deg; n2200: /* get next enode in current element. */ enode = list[enode]; if ( iq2 == 1 ) goto n900; goto n1600; n2300: /* get next element in the list. */ *tag = mtag; element = list[element]; goto n100; }