/*
* 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]<tag)&&(forward[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;
}