glpqmd.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

/* 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 <[e-mail]>.
*
*  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