/* glpqmd.c (quotient minimum degree algorithm) */

/***********************************************************************

*  This code is part of GLPK (GNU Linear Programming Kit).

*

*  THIS CODE IS THE RESULT OF TRANSLATION OF THE FORTRAN SUBROUTINES

*  GENQMD, QMDRCH, QMDQT, QMDUPD, AND QMDMRG FROM THE BOOK:

*

*  ALAN GEORGE, JOSEPH W-H LIU. COMPUTER SOLUTION OF LARGE SPARSE

*  POSITIVE DEFINITE SYSTEMS. PRENTICE-HALL, 1981.

*

*  THE TRANSLATION HAS BEEN DONE WITH THE PERMISSION OF THE AUTHORS

*  OF THE ORIGINAL FORTRAN SUBROUTINES: ALAN GEORGE AND JOSEPH LIU,

*  UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO, CANADA.

*

*  The translation was made by Andrew Makhorin <mao@gnu.org>.

*

*  GLPK is free software: you can redistribute it and/or modify it

*  under the terms of the GNU General Public License as published by

*  the Free Software Foundation, either version 3 of the License, or

*  (at your option) any later version.

*

*  GLPK is distributed in the hope that it will be useful, but WITHOUT

*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

*  License for more details.

*

*  You should have received a copy of the GNU General Public License

*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

***********************************************************************/

#include "glpqmd.h"

/***********************************************************************

*  NAME

*

*  genqmd - GENeral Quotient Minimum Degree algorithm

*

*  SYNOPSIS

*

*  #include "glpqmd.h"

*  void genqmd(int *neqns, int xadj[], int adjncy[], int perm[],

*     int invp[], int deg[], int marker[], int rchset[], int nbrhd[],

*     int qsize[], int qlink[], int *nofsub);

*

*  PURPOSE

*

*  This routine implements the minimum degree algorithm. It makes use

*  of the implicit representation of the elimination graph by quotient

*  graphs, and the notion of indistinguishable nodes.

*

*  CAUTION

*

*  The adjancy vector adjncy will be destroyed.

*

*  INPUT PARAMETERS

*

*  neqns  - number of equations;

*  (xadj, adjncy) -

*           the adjancy structure.

*

*  OUTPUT PARAMETERS

*

*  perm   - the minimum degree ordering;

*  invp   - the inverse of perm.

*

*  WORKING PARAMETERS

*

*  deg    - the degree vector. deg[i] is negative means node i has been

*           numbered;

*  marker - a marker vector, where marker[i] is negative means node i

*           has been merged with another nodeand thus can be ignored;

*  rchset - vector used for the reachable set;

*  nbrhd  - vector used for neighborhood set;

*  qsize  - vector used to store the size of indistinguishable

*           supernodes;

*  qlink  - vector used to store indistinguishable nodes, i, qlink[i],

*           qlink[qlink[i]], ... are the members of the supernode

*           represented by i.

*

*  PROGRAM SUBROUTINES

*

*  qmdrch, qmdqt, qmdupd.

***********************************************************************/

void genqmd(int *_neqns, int xadj[], int adjncy[], int perm[],

      int invp[], int deg[], int marker[], int rchset[], int nbrhd[],

      int qsize[], int qlink[], int *_nofsub)

{     int inode, ip, irch, j, mindeg, ndeg, nhdsze, node, np, num,

         nump1, nxnode, rchsze, search, thresh;

#     define neqns  (*_neqns)

#     define nofsub (*_nofsub)

      /* Initialize degree vector and other working variables. */

      mindeg = neqns;

      nofsub = 0;

      for (node = 1; node <= neqns; node++)

      {  perm[node] = node;

         invp[node] = node;

         marker[node] = 0;

         qsize[node] = 1;

         qlink[node] = 0;

         ndeg = xadj[node+1] - xadj[node];

         deg[node] = ndeg;

         if (ndeg < mindeg) mindeg = ndeg;

      }

      num = 0;

      /* Perform threshold search to get a node of min degree.

         Variable search point to where search should start. */

s200: search = 1;

      thresh = mindeg;

      mindeg = neqns;

s300: nump1 = num + 1;

      if (nump1 > search) search = nump1;

      for (j = search; j <= neqns; j++)

      {  node = perm[j];

         if (marker[node] >= 0)

         {  ndeg = deg[node];

            if (ndeg <= thresh) goto s500;

            if (ndeg < mindeg) mindeg = ndeg;

         }

      }

      goto s200;

      /* Node has minimum degree. Find its reachable sets by calling

         qmdrch. */

s500: search = j;

      nofsub += deg[node];

      marker[node] = 1;

      qmdrch(&node, xadj, adjncy, deg, marker, &rchsze, rchset, &nhdsze,

         nbrhd);

      /* Eliminate all nodes indistinguishable from node. They are given

         by node, qlink[node], ... . */

      nxnode = node;

s600: num++;

      np = invp[nxnode];

      ip = perm[num];

      perm[np] = ip;

      invp[ip] = np;

      perm[num] = nxnode;

      invp[nxnode] = num;

      deg[nxnode] = -1;

      nxnode = qlink[nxnode];

      if (nxnode > 0) goto s600;

      if (rchsze > 0)

      {  /* Update the degrees of the nodes in the reachable set and

            identify indistinguishable nodes. */

         qmdupd(xadj, adjncy, &rchsze, rchset, deg, qsize, qlink,

            marker, &rchset[rchsze+1], &nbrhd[nhdsze+1]);

         /* Reset marker value of nodes in reach set. Update threshold

            value for cyclic search. Also call qmdqt to form new

            quotient graph. */

         marker[node] = 0;

         for (irch = 1; irch <= rchsze; irch++)

         {  inode = rchset[irch];

            if (marker[inode] >= 0)

            {  marker[inode] = 0;

               ndeg = deg[inode];

               if (ndeg < mindeg) mindeg = ndeg;

               if (ndeg <= thresh)

               {  mindeg = thresh;

                  thresh = ndeg;

                  search = invp[inode];

               }

            }

         }

         if (nhdsze > 0)

            qmdqt(&node, xadj, adjncy, marker, &rchsze, rchset, nbrhd);

      }

      if (num < neqns) goto s300;

      return;

#     undef neqns

#     undef nofsub

}

/***********************************************************************

*  NAME

*

*  qmdrch - Quotient MD ReaCHable set

*

*  SYNOPSIS

*

*  #include "glpqmd.h"

*  void qmdrch(int *root, int xadj[], int adjncy[], int deg[],

*     int marker[], int *rchsze, int rchset[], int *nhdsze,

*     int nbrhd[]);

*

*  PURPOSE

*

*  This subroutine determines the reachable set of a node through a

*  given subset. The adjancy structure is assumed to be stored in a

*  quotient graph format.

* 

*  INPUT PARAMETERS

*

*  root   - the given node not in the subset;

*  (xadj, adjncy) -

*           the adjancy structure pair;

*  deg    - the degree vector. deg[i] < 0 means the node belongs to the

*           given subset.

*

*  OUTPUT PARAMETERS

*

*  (rchsze, rchset) -

*           the reachable set;

*  (nhdsze, nbrhd) -

*           the neighborhood set.

*

*  UPDATED PARAMETERS

*

*  marker - the marker vector for reach and nbrhd sets. > 0 means the

*           node is in reach set. < 0 means the node has been merged

*           with others in the quotient or it is in nbrhd set.

***********************************************************************/

void qmdrch(int *_root, int xadj[], int adjncy[], int deg[],

      int marker[], int *_rchsze, int rchset[], int *_nhdsze,

      int nbrhd[])

{     int i, istop, istrt, j, jstop, jstrt, nabor, node;

#     define root   (*_root)

#     define rchsze (*_rchsze)

#     define nhdsze (*_nhdsze)

      /* Loop through the neighbors of root in the quotient graph. */

      nhdsze = 0;

      rchsze = 0;

      istrt = xadj[root];

      istop = xadj[root+1] - 1;

      if (istop < istrt) return;

      for (i = istrt; i <= istop; i++)

      {  nabor = adjncy[i];

         if (nabor == 0) return;

         if (marker[nabor] == 0)

         {  if (deg[nabor] >= 0)

            {  /* Include nabor into the reachable set. */

               rchsze++;

               rchset[rchsze] = nabor;

               marker[nabor] = 1;

               goto s600;

            }

            /* nabor has been eliminated. Find nodes reachable from

               it. */

            marker[nabor] = -1;

            nhdsze++;

            nbrhd[nhdsze] = nabor;

s300:       jstrt = xadj[nabor];

            jstop = xadj[nabor+1] - 1;

            for (j = jstrt; j <= jstop; j++)

            {  node = adjncy[j];

               nabor = - node;

               if (node < 0) goto s300;

               if (node == 0) goto s600;

               if (marker[node] == 0)

               {  rchsze++;

                  rchset[rchsze] = node;

                  marker[node] = 1;

               }

            }

         }

s600:    ;

      }

      return;

#     undef root

#     undef rchsze

#     undef nhdsze

}

/***********************************************************************

*  NAME

*

*  qmdqt - Quotient MD Quotient graph Transformation

*

*  SYNOPSIS

*

*  #include "glpqmd.h"

*  void qmdqt(int *root, int xadj[], int adjncy[], int marker[],

*     int *rchsze, int rchset[], int nbrhd[]);

*

*  PURPOSE

*

*  This subroutine performs the quotient graph transformation after a

*  node has been eliminated.

*

*  INPUT PARAMETERS

*

*  root   - the node just eliminated. It becomes the representative of

*           the new supernode;

*  (xadj, adjncy) -

*           the adjancy structure;

*  (rchsze, rchset) -

*           the reachable set of root in the old quotient graph;

*  nbrhd  - the neighborhood set which will be merged with root to form

*           the new supernode;

*  marker - the marker vector.

*

*  UPDATED PARAMETERS

*

*  adjncy - becomes the adjncy of the quotient graph.

***********************************************************************/

void qmdqt(int *_root, int xadj[], int adjncy[], int marker[],

      int *_rchsze, int rchset[], int nbrhd[])

{     int inhd, irch, j, jstop, jstrt, link, nabor, node;

#     define root   (*_root)

#     define rchsze (*_rchsze)

      irch = 0;

      inhd = 0;

      node = root;

s100: jstrt = xadj[node];

      jstop = xadj[node+1] - 2;

      if (jstop >= jstrt)

      {  /* Place reach nodes into the adjacent list of node. */

         for (j = jstrt; j <= jstop; j++)

         {  irch++;

            adjncy[j] = rchset[irch];

            if (irch >= rchsze) goto s400;

         }

      }

      /* Link to other space provided by the nbrhd set. */

      link = adjncy[jstop+1];

      node = - link;

      if (link >= 0)

      {  inhd++;

         node = nbrhd[inhd];

         adjncy[jstop+1] = - node;

      }

      goto s100;

      /* All reachable nodes have been saved. End the adjacent list.

         Add root to the neighborhood list of each node in the reach

         set. */

s400: adjncy[j+1] = 0;

      for (irch = 1; irch <= rchsze; irch++)

      {  node = rchset[irch];

         if (marker[node] >= 0)

         {  jstrt = xadj[node];

            jstop = xadj[node+1] - 1;

            for (j = jstrt; j <= jstop; j++)

            {  nabor = adjncy[j];

               if (marker[nabor] < 0)

               {  adjncy[j] = root;

                  goto s600;

               }

            }

         }

s600:    ;

      }

      return;

#     undef root

#     undef rchsze

}

/***********************************************************************

*  NAME

*

*  qmdupd - Quotient MD UPDate

*

*  SYNOPSIS

*

*  #include "glpqmd.h"

*  void qmdupd(int xadj[], int adjncy[], int *nlist, int list[],

*     int deg[], int qsize[], int qlink[], int marker[], int rchset[],

*     int nbrhd[]);

*

*  PURPOSE

*

*  This routine performs degree update for a set of nodes in the minimum

*  degree algorithm.

*

*  INPUT PARAMETERS

*

*  (xadj, adjncy) -

*           the adjancy structure;

*  (nlist, list) -

*           the list of nodes whose degree has to be updated.

*

*  UPDATED PARAMETERS

*

*  deg    - the degree vector;

*  qsize  - size of indistinguishable supernodes;

*  qlink  - linked list for indistinguishable nodes;

*  marker - used to mark those nodes in reach/nbrhd sets.

*

*  WORKING PARAMETERS

*

*  rchset - the reachable set;

*  nbrhd  - the neighborhood set.

*

*  PROGRAM SUBROUTINES

*

*  qmdmrg.

***********************************************************************/

void qmdupd(int xadj[], int adjncy[], int *_nlist, int list[],

      int deg[], int qsize[], int qlink[], int marker[], int rchset[],

      int nbrhd[])

{     int deg0, deg1, il, inhd, inode, irch, j, jstop, jstrt, mark,

         nabor, nhdsze, node, rchsze;

#     define nlist  (*_nlist)

      /* Find all eliminated supernodes that are adjacent to some nodes

         in the given list. Put them into (nhdsze, nbrhd). deg0 contains

         the number of nodes in the list. */

      if (nlist <= 0) return;

      deg0 = 0;

      nhdsze = 0;

      for (il = 1; il <= nlist; il++)

      {  node = list[il];

         deg0 += qsize[node];

         jstrt = xadj[node];

         jstop = xadj[node+1] - 1;

         for (j = jstrt; j <= jstop; j++)

         {  nabor = adjncy[j];

            if (marker[nabor] == 0 && deg[nabor] < 0)

            {  marker[nabor] = -1;

               nhdsze++;

               nbrhd[nhdsze] = nabor;

            }

         }

      }

      /* Merge indistinguishable nodes in the list by calling the

         subroutine qmdmrg. */

      if (nhdsze > 0)

         qmdmrg(xadj, adjncy, deg, qsize, qlink, marker, &deg0, &nhdsze,

            nbrhd, rchset, &nbrhd[nhdsze+1]);

      /* Find the new degrees of the nodes that have not been merged. */

      for (il = 1; il <= nlist; il++)

      {  node = list[il];

         mark = marker[node];

         if (mark == 0 || mark == 1)

         {  marker[node] = 2;

            qmdrch(&node, xadj, adjncy, deg, marker, &rchsze, rchset,

               &nhdsze, nbrhd);

            deg1 = deg0;

            if (rchsze > 0)

            {  for (irch = 1; irch <= rchsze; irch++)

               {  inode = rchset[irch];

                  deg1 += qsize[inode];

                  marker[inode] = 0;

               }

            }

            deg[node] = deg1 - 1;

            if (nhdsze > 0)

            {  for (inhd = 1; inhd <= nhdsze; inhd++)

               {  inode = nbrhd[inhd];

                  marker[inode] = 0;

               }

            }

         }

      }

      return;

#     undef nlist

}

/***********************************************************************

*  NAME

*

*  qmdmrg - Quotient MD MeRGe

*

*  SYNOPSIS

*

*  #include "qmdmrg.h"

*  void qmdmrg(int xadj[], int adjncy[], int deg[], int qsize[],

*     int qlink[], int marker[], int *deg0, int *nhdsze, int nbrhd[],

*     int rchset[], int ovrlp[]);

*

*  PURPOSE

*

*  This routine merges indistinguishable nodes in the minimum degree

*  ordering algorithm. It also computes the new degrees of these new

*  supernodes.

*

*  INPUT PARAMETERS

*

*  (xadj, adjncy) -

*           the adjancy structure;

*  deg0   - the number of nodes in the given set;

*  (nhdsze, nbrhd) -

*           the set of eliminated supernodes adjacent to some nodes in

*           the set.

*

*  UPDATED PARAMETERS

*

*  deg    - the degree vector;

*  qsize  - size of indistinguishable nodes;

*  qlink  - linked list for indistinguishable nodes;

*  marker - the given set is given by those nodes with marker value set

*           to 1. Those nodes with degree updated will have marker value

*           set to 2.

*

*  WORKING PARAMETERS

*

*  rchset - the reachable set;

*  ovrlp  - temp vector to store the intersection of two reachable sets.

***********************************************************************/

void qmdmrg(int xadj[], int adjncy[], int deg[], int qsize[],

      int qlink[], int marker[], int *_deg0, int *_nhdsze, int nbrhd[],

      int rchset[], int ovrlp[])

{     int deg1, head, inhd, iov, irch, j, jstop, jstrt, link, lnode,

         mark, mrgsze, nabor, node, novrlp, rchsze, root;

#     define deg0   (*_deg0)

#     define nhdsze (*_nhdsze)

      /* Initialization. */

      if (nhdsze <= 0) return;

      for (inhd = 1; inhd <= nhdsze; inhd++)

      {  root = nbrhd[inhd];

         marker[root] = 0;

      }

      /* Loop through each eliminated supernode in the set

         (nhdsze, nbrhd). */

      for (inhd = 1; inhd <= nhdsze; inhd++)

      {  root = nbrhd[inhd];

         marker[root] = -1;

         rchsze = 0;

         novrlp = 0;

         deg1 = 0;

s200:    jstrt = xadj[root];

         jstop = xadj[root+1] - 1;

         /* Determine the reachable set and its intersection with the

            input reachable set. */

         for (j = jstrt; j <= jstop; j++)

         {  nabor = adjncy[j];

            root = - nabor;

            if (nabor < 0) goto s200;

            if (nabor == 0) break;

            mark = marker[nabor];

            if (mark == 0)

            {  rchsze++;

               rchset[rchsze] = nabor;

               deg1 += qsize[nabor];

               marker[nabor] = 1;

            }

            else if (mark == 1)

            {  novrlp++;

               ovrlp[novrlp] = nabor;

               marker[nabor] = 2;

            }

         }

         /* From the overlapped set, determine the nodes that can be

            merged together. */

         head = 0;

         mrgsze = 0;

         for (iov = 1; iov <= novrlp; iov++)

         {  node = ovrlp[iov];

            jstrt = xadj[node];

            jstop = xadj[node+1] - 1;

            for (j = jstrt; j <= jstop; j++)

            {  nabor = adjncy[j];

               if (marker[nabor] == 0)

               {  marker[node] = 1;

                  goto s1100;

               }

            }

            /* Node belongs to the new merged supernode. Update the

               vectors qlink and qsize. */

            mrgsze += qsize[node];

            marker[node] = -1;

            lnode = node;

s900:       link = qlink[lnode];

            if (link > 0)

            {  lnode = link;

               goto s900;

            }

            qlink[lnode] = head;

            head = node;

s1100:      ;

         }

         if (head > 0)

         {  qsize[head] = mrgsze;

            deg[head] = deg0 + deg1 - 1;

            marker[head] = 2;

         }

         /* Reset marker values. */

         root = nbrhd[inhd];

         marker[root] = 0;

         if (rchsze > 0)

         {  for (irch = 1; irch <= rchsze; irch++)

            {  node = rchset[irch];

               marker[node] = 0;

            }

         }

      }

      return;

#     undef deg0

#     undef nhdsze

}

/* eof */