1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpqmd.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,584 @@
     1.4 +/* glpqmd.c (quotient minimum degree algorithm) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  THIS CODE IS THE RESULT OF TRANSLATION OF THE FORTRAN SUBROUTINES
    1.10 +*  GENQMD, QMDRCH, QMDQT, QMDUPD, AND QMDMRG FROM THE BOOK:
    1.11 +*
    1.12 +*  ALAN GEORGE, JOSEPH W-H LIU. COMPUTER SOLUTION OF LARGE SPARSE
    1.13 +*  POSITIVE DEFINITE SYSTEMS. PRENTICE-HALL, 1981.
    1.14 +*
    1.15 +*  THE TRANSLATION HAS BEEN DONE WITH THE PERMISSION OF THE AUTHORS
    1.16 +*  OF THE ORIGINAL FORTRAN SUBROUTINES: ALAN GEORGE AND JOSEPH LIU,
    1.17 +*  UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO, CANADA.
    1.18 +*
    1.19 +*  The translation was made by Andrew Makhorin <mao@gnu.org>.
    1.20 +*
    1.21 +*  GLPK is free software: you can redistribute it and/or modify it
    1.22 +*  under the terms of the GNU General Public License as published by
    1.23 +*  the Free Software Foundation, either version 3 of the License, or
    1.24 +*  (at your option) any later version.
    1.25 +*
    1.26 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.27 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.28 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.29 +*  License for more details.
    1.30 +*
    1.31 +*  You should have received a copy of the GNU General Public License
    1.32 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.33 +***********************************************************************/
    1.34 +
    1.35 +#include "glpqmd.h"
    1.36 +
    1.37 +/***********************************************************************
    1.38 +*  NAME
    1.39 +*
    1.40 +*  genqmd - GENeral Quotient Minimum Degree algorithm
    1.41 +*
    1.42 +*  SYNOPSIS
    1.43 +*
    1.44 +*  #include "glpqmd.h"
    1.45 +*  void genqmd(int *neqns, int xadj[], int adjncy[], int perm[],
    1.46 +*     int invp[], int deg[], int marker[], int rchset[], int nbrhd[],
    1.47 +*     int qsize[], int qlink[], int *nofsub);
    1.48 +*
    1.49 +*  PURPOSE
    1.50 +*
    1.51 +*  This routine implements the minimum degree algorithm. It makes use
    1.52 +*  of the implicit representation of the elimination graph by quotient
    1.53 +*  graphs, and the notion of indistinguishable nodes.
    1.54 +*
    1.55 +*  CAUTION
    1.56 +*
    1.57 +*  The adjancy vector adjncy will be destroyed.
    1.58 +*
    1.59 +*  INPUT PARAMETERS
    1.60 +*
    1.61 +*  neqns  - number of equations;
    1.62 +*  (xadj, adjncy) -
    1.63 +*           the adjancy structure.
    1.64 +*
    1.65 +*  OUTPUT PARAMETERS
    1.66 +*
    1.67 +*  perm   - the minimum degree ordering;
    1.68 +*  invp   - the inverse of perm.
    1.69 +*
    1.70 +*  WORKING PARAMETERS
    1.71 +*
    1.72 +*  deg    - the degree vector. deg[i] is negative means node i has been
    1.73 +*           numbered;
    1.74 +*  marker - a marker vector, where marker[i] is negative means node i
    1.75 +*           has been merged with another nodeand thus can be ignored;
    1.76 +*  rchset - vector used for the reachable set;
    1.77 +*  nbrhd  - vector used for neighborhood set;
    1.78 +*  qsize  - vector used to store the size of indistinguishable
    1.79 +*           supernodes;
    1.80 +*  qlink  - vector used to store indistinguishable nodes, i, qlink[i],
    1.81 +*           qlink[qlink[i]], ... are the members of the supernode
    1.82 +*           represented by i.
    1.83 +*
    1.84 +*  PROGRAM SUBROUTINES
    1.85 +*
    1.86 +*  qmdrch, qmdqt, qmdupd.
    1.87 +***********************************************************************/
    1.88 +
    1.89 +void genqmd(int *_neqns, int xadj[], int adjncy[], int perm[],
    1.90 +      int invp[], int deg[], int marker[], int rchset[], int nbrhd[],
    1.91 +      int qsize[], int qlink[], int *_nofsub)
    1.92 +{     int inode, ip, irch, j, mindeg, ndeg, nhdsze, node, np, num,
    1.93 +         nump1, nxnode, rchsze, search, thresh;
    1.94 +#     define neqns  (*_neqns)
    1.95 +#     define nofsub (*_nofsub)
    1.96 +      /* Initialize degree vector and other working variables. */
    1.97 +      mindeg = neqns;
    1.98 +      nofsub = 0;
    1.99 +      for (node = 1; node <= neqns; node++)
   1.100 +      {  perm[node] = node;
   1.101 +         invp[node] = node;
   1.102 +         marker[node] = 0;
   1.103 +         qsize[node] = 1;
   1.104 +         qlink[node] = 0;
   1.105 +         ndeg = xadj[node+1] - xadj[node];
   1.106 +         deg[node] = ndeg;
   1.107 +         if (ndeg < mindeg) mindeg = ndeg;
   1.108 +      }
   1.109 +      num = 0;
   1.110 +      /* Perform threshold search to get a node of min degree.
   1.111 +         Variable search point to where search should start. */
   1.112 +s200: search = 1;
   1.113 +      thresh = mindeg;
   1.114 +      mindeg = neqns;
   1.115 +s300: nump1 = num + 1;
   1.116 +      if (nump1 > search) search = nump1;
   1.117 +      for (j = search; j <= neqns; j++)
   1.118 +      {  node = perm[j];
   1.119 +         if (marker[node] >= 0)
   1.120 +         {  ndeg = deg[node];
   1.121 +            if (ndeg <= thresh) goto s500;
   1.122 +            if (ndeg < mindeg) mindeg = ndeg;
   1.123 +         }
   1.124 +      }
   1.125 +      goto s200;
   1.126 +      /* Node has minimum degree. Find its reachable sets by calling
   1.127 +         qmdrch. */
   1.128 +s500: search = j;
   1.129 +      nofsub += deg[node];
   1.130 +      marker[node] = 1;
   1.131 +      qmdrch(&node, xadj, adjncy, deg, marker, &rchsze, rchset, &nhdsze,
   1.132 +         nbrhd);
   1.133 +      /* Eliminate all nodes indistinguishable from node. They are given
   1.134 +         by node, qlink[node], ... . */
   1.135 +      nxnode = node;
   1.136 +s600: num++;
   1.137 +      np = invp[nxnode];
   1.138 +      ip = perm[num];
   1.139 +      perm[np] = ip;
   1.140 +      invp[ip] = np;
   1.141 +      perm[num] = nxnode;
   1.142 +      invp[nxnode] = num;
   1.143 +      deg[nxnode] = -1;
   1.144 +      nxnode = qlink[nxnode];
   1.145 +      if (nxnode > 0) goto s600;
   1.146 +      if (rchsze > 0)
   1.147 +      {  /* Update the degrees of the nodes in the reachable set and
   1.148 +            identify indistinguishable nodes. */
   1.149 +         qmdupd(xadj, adjncy, &rchsze, rchset, deg, qsize, qlink,
   1.150 +            marker, &rchset[rchsze+1], &nbrhd[nhdsze+1]);
   1.151 +         /* Reset marker value of nodes in reach set. Update threshold
   1.152 +            value for cyclic search. Also call qmdqt to form new
   1.153 +            quotient graph. */
   1.154 +         marker[node] = 0;
   1.155 +         for (irch = 1; irch <= rchsze; irch++)
   1.156 +         {  inode = rchset[irch];
   1.157 +            if (marker[inode] >= 0)
   1.158 +            {  marker[inode] = 0;
   1.159 +               ndeg = deg[inode];
   1.160 +               if (ndeg < mindeg) mindeg = ndeg;
   1.161 +               if (ndeg <= thresh)
   1.162 +               {  mindeg = thresh;
   1.163 +                  thresh = ndeg;
   1.164 +                  search = invp[inode];
   1.165 +               }
   1.166 +            }
   1.167 +         }
   1.168 +         if (nhdsze > 0)
   1.169 +            qmdqt(&node, xadj, adjncy, marker, &rchsze, rchset, nbrhd);
   1.170 +      }
   1.171 +      if (num < neqns) goto s300;
   1.172 +      return;
   1.173 +#     undef neqns
   1.174 +#     undef nofsub
   1.175 +}
   1.176 +
   1.177 +/***********************************************************************
   1.178 +*  NAME
   1.179 +*
   1.180 +*  qmdrch - Quotient MD ReaCHable set
   1.181 +*
   1.182 +*  SYNOPSIS
   1.183 +*
   1.184 +*  #include "glpqmd.h"
   1.185 +*  void qmdrch(int *root, int xadj[], int adjncy[], int deg[],
   1.186 +*     int marker[], int *rchsze, int rchset[], int *nhdsze,
   1.187 +*     int nbrhd[]);
   1.188 +*
   1.189 +*  PURPOSE
   1.190 +*
   1.191 +*  This subroutine determines the reachable set of a node through a
   1.192 +*  given subset. The adjancy structure is assumed to be stored in a
   1.193 +*  quotient graph format.
   1.194 +* 
   1.195 +*  INPUT PARAMETERS
   1.196 +*
   1.197 +*  root   - the given node not in the subset;
   1.198 +*  (xadj, adjncy) -
   1.199 +*           the adjancy structure pair;
   1.200 +*  deg    - the degree vector. deg[i] < 0 means the node belongs to the
   1.201 +*           given subset.
   1.202 +*
   1.203 +*  OUTPUT PARAMETERS
   1.204 +*
   1.205 +*  (rchsze, rchset) -
   1.206 +*           the reachable set;
   1.207 +*  (nhdsze, nbrhd) -
   1.208 +*           the neighborhood set.
   1.209 +*
   1.210 +*  UPDATED PARAMETERS
   1.211 +*
   1.212 +*  marker - the marker vector for reach and nbrhd sets. > 0 means the
   1.213 +*           node is in reach set. < 0 means the node has been merged
   1.214 +*           with others in the quotient or it is in nbrhd set.
   1.215 +***********************************************************************/
   1.216 +
   1.217 +void qmdrch(int *_root, int xadj[], int adjncy[], int deg[],
   1.218 +      int marker[], int *_rchsze, int rchset[], int *_nhdsze,
   1.219 +      int nbrhd[])
   1.220 +{     int i, istop, istrt, j, jstop, jstrt, nabor, node;
   1.221 +#     define root   (*_root)
   1.222 +#     define rchsze (*_rchsze)
   1.223 +#     define nhdsze (*_nhdsze)
   1.224 +      /* Loop through the neighbors of root in the quotient graph. */
   1.225 +      nhdsze = 0;
   1.226 +      rchsze = 0;
   1.227 +      istrt = xadj[root];
   1.228 +      istop = xadj[root+1] - 1;
   1.229 +      if (istop < istrt) return;
   1.230 +      for (i = istrt; i <= istop; i++)
   1.231 +      {  nabor = adjncy[i];
   1.232 +         if (nabor == 0) return;
   1.233 +         if (marker[nabor] == 0)
   1.234 +         {  if (deg[nabor] >= 0)
   1.235 +            {  /* Include nabor into the reachable set. */
   1.236 +               rchsze++;
   1.237 +               rchset[rchsze] = nabor;
   1.238 +               marker[nabor] = 1;
   1.239 +               goto s600;
   1.240 +            }
   1.241 +            /* nabor has been eliminated. Find nodes reachable from
   1.242 +               it. */
   1.243 +            marker[nabor] = -1;
   1.244 +            nhdsze++;
   1.245 +            nbrhd[nhdsze] = nabor;
   1.246 +s300:       jstrt = xadj[nabor];
   1.247 +            jstop = xadj[nabor+1] - 1;
   1.248 +            for (j = jstrt; j <= jstop; j++)
   1.249 +            {  node = adjncy[j];
   1.250 +               nabor = - node;
   1.251 +               if (node < 0) goto s300;
   1.252 +               if (node == 0) goto s600;
   1.253 +               if (marker[node] == 0)
   1.254 +               {  rchsze++;
   1.255 +                  rchset[rchsze] = node;
   1.256 +                  marker[node] = 1;
   1.257 +               }
   1.258 +            }
   1.259 +         }
   1.260 +s600:    ;
   1.261 +      }
   1.262 +      return;
   1.263 +#     undef root
   1.264 +#     undef rchsze
   1.265 +#     undef nhdsze
   1.266 +}
   1.267 +
   1.268 +/***********************************************************************
   1.269 +*  NAME
   1.270 +*
   1.271 +*  qmdqt - Quotient MD Quotient graph Transformation
   1.272 +*
   1.273 +*  SYNOPSIS
   1.274 +*
   1.275 +*  #include "glpqmd.h"
   1.276 +*  void qmdqt(int *root, int xadj[], int adjncy[], int marker[],
   1.277 +*     int *rchsze, int rchset[], int nbrhd[]);
   1.278 +*
   1.279 +*  PURPOSE
   1.280 +*
   1.281 +*  This subroutine performs the quotient graph transformation after a
   1.282 +*  node has been eliminated.
   1.283 +*
   1.284 +*  INPUT PARAMETERS
   1.285 +*
   1.286 +*  root   - the node just eliminated. It becomes the representative of
   1.287 +*           the new supernode;
   1.288 +*  (xadj, adjncy) -
   1.289 +*           the adjancy structure;
   1.290 +*  (rchsze, rchset) -
   1.291 +*           the reachable set of root in the old quotient graph;
   1.292 +*  nbrhd  - the neighborhood set which will be merged with root to form
   1.293 +*           the new supernode;
   1.294 +*  marker - the marker vector.
   1.295 +*
   1.296 +*  UPDATED PARAMETERS
   1.297 +*
   1.298 +*  adjncy - becomes the adjncy of the quotient graph.
   1.299 +***********************************************************************/
   1.300 +
   1.301 +void qmdqt(int *_root, int xadj[], int adjncy[], int marker[],
   1.302 +      int *_rchsze, int rchset[], int nbrhd[])
   1.303 +{     int inhd, irch, j, jstop, jstrt, link, nabor, node;
   1.304 +#     define root   (*_root)
   1.305 +#     define rchsze (*_rchsze)
   1.306 +      irch = 0;
   1.307 +      inhd = 0;
   1.308 +      node = root;
   1.309 +s100: jstrt = xadj[node];
   1.310 +      jstop = xadj[node+1] - 2;
   1.311 +      if (jstop >= jstrt)
   1.312 +      {  /* Place reach nodes into the adjacent list of node. */
   1.313 +         for (j = jstrt; j <= jstop; j++)
   1.314 +         {  irch++;
   1.315 +            adjncy[j] = rchset[irch];
   1.316 +            if (irch >= rchsze) goto s400;
   1.317 +         }
   1.318 +      }
   1.319 +      /* Link to other space provided by the nbrhd set. */
   1.320 +      link = adjncy[jstop+1];
   1.321 +      node = - link;
   1.322 +      if (link >= 0)
   1.323 +      {  inhd++;
   1.324 +         node = nbrhd[inhd];
   1.325 +         adjncy[jstop+1] = - node;
   1.326 +      }
   1.327 +      goto s100;
   1.328 +      /* All reachable nodes have been saved. End the adjacent list.
   1.329 +         Add root to the neighborhood list of each node in the reach
   1.330 +         set. */
   1.331 +s400: adjncy[j+1] = 0;
   1.332 +      for (irch = 1; irch <= rchsze; irch++)
   1.333 +      {  node = rchset[irch];
   1.334 +         if (marker[node] >= 0)
   1.335 +         {  jstrt = xadj[node];
   1.336 +            jstop = xadj[node+1] - 1;
   1.337 +            for (j = jstrt; j <= jstop; j++)
   1.338 +            {  nabor = adjncy[j];
   1.339 +               if (marker[nabor] < 0)
   1.340 +               {  adjncy[j] = root;
   1.341 +                  goto s600;
   1.342 +               }
   1.343 +            }
   1.344 +         }
   1.345 +s600:    ;
   1.346 +      }
   1.347 +      return;
   1.348 +#     undef root
   1.349 +#     undef rchsze
   1.350 +}
   1.351 +
   1.352 +/***********************************************************************
   1.353 +*  NAME
   1.354 +*
   1.355 +*  qmdupd - Quotient MD UPDate
   1.356 +*
   1.357 +*  SYNOPSIS
   1.358 +*
   1.359 +*  #include "glpqmd.h"
   1.360 +*  void qmdupd(int xadj[], int adjncy[], int *nlist, int list[],
   1.361 +*     int deg[], int qsize[], int qlink[], int marker[], int rchset[],
   1.362 +*     int nbrhd[]);
   1.363 +*
   1.364 +*  PURPOSE
   1.365 +*
   1.366 +*  This routine performs degree update for a set of nodes in the minimum
   1.367 +*  degree algorithm.
   1.368 +*
   1.369 +*  INPUT PARAMETERS
   1.370 +*
   1.371 +*  (xadj, adjncy) -
   1.372 +*           the adjancy structure;
   1.373 +*  (nlist, list) -
   1.374 +*           the list of nodes whose degree has to be updated.
   1.375 +*
   1.376 +*  UPDATED PARAMETERS
   1.377 +*
   1.378 +*  deg    - the degree vector;
   1.379 +*  qsize  - size of indistinguishable supernodes;
   1.380 +*  qlink  - linked list for indistinguishable nodes;
   1.381 +*  marker - used to mark those nodes in reach/nbrhd sets.
   1.382 +*
   1.383 +*  WORKING PARAMETERS
   1.384 +*
   1.385 +*  rchset - the reachable set;
   1.386 +*  nbrhd  - the neighborhood set.
   1.387 +*
   1.388 +*  PROGRAM SUBROUTINES
   1.389 +*
   1.390 +*  qmdmrg.
   1.391 +***********************************************************************/
   1.392 +
   1.393 +void qmdupd(int xadj[], int adjncy[], int *_nlist, int list[],
   1.394 +      int deg[], int qsize[], int qlink[], int marker[], int rchset[],
   1.395 +      int nbrhd[])
   1.396 +{     int deg0, deg1, il, inhd, inode, irch, j, jstop, jstrt, mark,
   1.397 +         nabor, nhdsze, node, rchsze;
   1.398 +#     define nlist  (*_nlist)
   1.399 +      /* Find all eliminated supernodes that are adjacent to some nodes
   1.400 +         in the given list. Put them into (nhdsze, nbrhd). deg0 contains
   1.401 +         the number of nodes in the list. */
   1.402 +      if (nlist <= 0) return;
   1.403 +      deg0 = 0;
   1.404 +      nhdsze = 0;
   1.405 +      for (il = 1; il <= nlist; il++)
   1.406 +      {  node = list[il];
   1.407 +         deg0 += qsize[node];
   1.408 +         jstrt = xadj[node];
   1.409 +         jstop = xadj[node+1] - 1;
   1.410 +         for (j = jstrt; j <= jstop; j++)
   1.411 +         {  nabor = adjncy[j];
   1.412 +            if (marker[nabor] == 0 && deg[nabor] < 0)
   1.413 +            {  marker[nabor] = -1;
   1.414 +               nhdsze++;
   1.415 +               nbrhd[nhdsze] = nabor;
   1.416 +            }
   1.417 +         }
   1.418 +      }
   1.419 +      /* Merge indistinguishable nodes in the list by calling the
   1.420 +         subroutine qmdmrg. */
   1.421 +      if (nhdsze > 0)
   1.422 +         qmdmrg(xadj, adjncy, deg, qsize, qlink, marker, &deg0, &nhdsze,
   1.423 +            nbrhd, rchset, &nbrhd[nhdsze+1]);
   1.424 +      /* Find the new degrees of the nodes that have not been merged. */
   1.425 +      for (il = 1; il <= nlist; il++)
   1.426 +      {  node = list[il];
   1.427 +         mark = marker[node];
   1.428 +         if (mark == 0 || mark == 1)
   1.429 +         {  marker[node] = 2;
   1.430 +            qmdrch(&node, xadj, adjncy, deg, marker, &rchsze, rchset,
   1.431 +               &nhdsze, nbrhd);
   1.432 +            deg1 = deg0;
   1.433 +            if (rchsze > 0)
   1.434 +            {  for (irch = 1; irch <= rchsze; irch++)
   1.435 +               {  inode = rchset[irch];
   1.436 +                  deg1 += qsize[inode];
   1.437 +                  marker[inode] = 0;
   1.438 +               }
   1.439 +            }
   1.440 +            deg[node] = deg1 - 1;
   1.441 +            if (nhdsze > 0)
   1.442 +            {  for (inhd = 1; inhd <= nhdsze; inhd++)
   1.443 +               {  inode = nbrhd[inhd];
   1.444 +                  marker[inode] = 0;
   1.445 +               }
   1.446 +            }
   1.447 +         }
   1.448 +      }
   1.449 +      return;
   1.450 +#     undef nlist
   1.451 +}
   1.452 +
   1.453 +/***********************************************************************
   1.454 +*  NAME
   1.455 +*
   1.456 +*  qmdmrg - Quotient MD MeRGe
   1.457 +*
   1.458 +*  SYNOPSIS
   1.459 +*
   1.460 +*  #include "qmdmrg.h"
   1.461 +*  void qmdmrg(int xadj[], int adjncy[], int deg[], int qsize[],
   1.462 +*     int qlink[], int marker[], int *deg0, int *nhdsze, int nbrhd[],
   1.463 +*     int rchset[], int ovrlp[]);
   1.464 +*
   1.465 +*  PURPOSE
   1.466 +*
   1.467 +*  This routine merges indistinguishable nodes in the minimum degree
   1.468 +*  ordering algorithm. It also computes the new degrees of these new
   1.469 +*  supernodes.
   1.470 +*
   1.471 +*  INPUT PARAMETERS
   1.472 +*
   1.473 +*  (xadj, adjncy) -
   1.474 +*           the adjancy structure;
   1.475 +*  deg0   - the number of nodes in the given set;
   1.476 +*  (nhdsze, nbrhd) -
   1.477 +*           the set of eliminated supernodes adjacent to some nodes in
   1.478 +*           the set.
   1.479 +*
   1.480 +*  UPDATED PARAMETERS
   1.481 +*
   1.482 +*  deg    - the degree vector;
   1.483 +*  qsize  - size of indistinguishable nodes;
   1.484 +*  qlink  - linked list for indistinguishable nodes;
   1.485 +*  marker - the given set is given by those nodes with marker value set
   1.486 +*           to 1. Those nodes with degree updated will have marker value
   1.487 +*           set to 2.
   1.488 +*
   1.489 +*  WORKING PARAMETERS
   1.490 +*
   1.491 +*  rchset - the reachable set;
   1.492 +*  ovrlp  - temp vector to store the intersection of two reachable sets.
   1.493 +***********************************************************************/
   1.494 +
   1.495 +void qmdmrg(int xadj[], int adjncy[], int deg[], int qsize[],
   1.496 +      int qlink[], int marker[], int *_deg0, int *_nhdsze, int nbrhd[],
   1.497 +      int rchset[], int ovrlp[])
   1.498 +{     int deg1, head, inhd, iov, irch, j, jstop, jstrt, link, lnode,
   1.499 +         mark, mrgsze, nabor, node, novrlp, rchsze, root;
   1.500 +#     define deg0   (*_deg0)
   1.501 +#     define nhdsze (*_nhdsze)
   1.502 +      /* Initialization. */
   1.503 +      if (nhdsze <= 0) return;
   1.504 +      for (inhd = 1; inhd <= nhdsze; inhd++)
   1.505 +      {  root = nbrhd[inhd];
   1.506 +         marker[root] = 0;
   1.507 +      }
   1.508 +      /* Loop through each eliminated supernode in the set
   1.509 +         (nhdsze, nbrhd). */
   1.510 +      for (inhd = 1; inhd <= nhdsze; inhd++)
   1.511 +      {  root = nbrhd[inhd];
   1.512 +         marker[root] = -1;
   1.513 +         rchsze = 0;
   1.514 +         novrlp = 0;
   1.515 +         deg1 = 0;
   1.516 +s200:    jstrt = xadj[root];
   1.517 +         jstop = xadj[root+1] - 1;
   1.518 +         /* Determine the reachable set and its intersection with the
   1.519 +            input reachable set. */
   1.520 +         for (j = jstrt; j <= jstop; j++)
   1.521 +         {  nabor = adjncy[j];
   1.522 +            root = - nabor;
   1.523 +            if (nabor < 0) goto s200;
   1.524 +            if (nabor == 0) break;
   1.525 +            mark = marker[nabor];
   1.526 +            if (mark == 0)
   1.527 +            {  rchsze++;
   1.528 +               rchset[rchsze] = nabor;
   1.529 +               deg1 += qsize[nabor];
   1.530 +               marker[nabor] = 1;
   1.531 +            }
   1.532 +            else if (mark == 1)
   1.533 +            {  novrlp++;
   1.534 +               ovrlp[novrlp] = nabor;
   1.535 +               marker[nabor] = 2;
   1.536 +            }
   1.537 +         }
   1.538 +         /* From the overlapped set, determine the nodes that can be
   1.539 +            merged together. */
   1.540 +         head = 0;
   1.541 +         mrgsze = 0;
   1.542 +         for (iov = 1; iov <= novrlp; iov++)
   1.543 +         {  node = ovrlp[iov];
   1.544 +            jstrt = xadj[node];
   1.545 +            jstop = xadj[node+1] - 1;
   1.546 +            for (j = jstrt; j <= jstop; j++)
   1.547 +            {  nabor = adjncy[j];
   1.548 +               if (marker[nabor] == 0)
   1.549 +               {  marker[node] = 1;
   1.550 +                  goto s1100;
   1.551 +               }
   1.552 +            }
   1.553 +            /* Node belongs to the new merged supernode. Update the
   1.554 +               vectors qlink and qsize. */
   1.555 +            mrgsze += qsize[node];
   1.556 +            marker[node] = -1;
   1.557 +            lnode = node;
   1.558 +s900:       link = qlink[lnode];
   1.559 +            if (link > 0)
   1.560 +            {  lnode = link;
   1.561 +               goto s900;
   1.562 +            }
   1.563 +            qlink[lnode] = head;
   1.564 +            head = node;
   1.565 +s1100:      ;
   1.566 +         }
   1.567 +         if (head > 0)
   1.568 +         {  qsize[head] = mrgsze;
   1.569 +            deg[head] = deg0 + deg1 - 1;
   1.570 +            marker[head] = 2;
   1.571 +         }
   1.572 +         /* Reset marker values. */
   1.573 +         root = nbrhd[inhd];
   1.574 +         marker[root] = 0;
   1.575 +         if (rchsze > 0)
   1.576 +         {  for (irch = 1; irch <= rchsze; irch++)
   1.577 +            {  node = rchset[irch];
   1.578 +               marker[node] = 0;
   1.579 +            }
   1.580 +         }
   1.581 +      }
   1.582 +      return;
   1.583 +#     undef deg0
   1.584 +#     undef nhdsze
   1.585 +}
   1.586 +
   1.587 +/* eof */