/* 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 */