alpar@9: /* glpapi16.c (graph and network analysis routines) */
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@9: *
alpar@9: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9: *  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9: *  E-mail: <mao@gnu.org>.
alpar@9: *
alpar@9: *  GLPK is free software: you can redistribute it and/or modify it
alpar@9: *  under the terms of the GNU General Public License as published by
alpar@9: *  the Free Software Foundation, either version 3 of the License, or
alpar@9: *  (at your option) any later version.
alpar@9: *
alpar@9: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9: *  License for more details.
alpar@9: *
alpar@9: *  You should have received a copy of the GNU General Public License
alpar@9: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9: ***********************************************************************/
alpar@9: 
alpar@9: #include "glpapi.h"
alpar@9: #include "glpnet.h"
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  NAME
alpar@9: *
alpar@9: *  glp_weak_comp - find all weakly connected components of graph
alpar@9: *
alpar@9: *  SYNOPSIS
alpar@9: *
alpar@9: *  int glp_weak_comp(glp_graph *G, int v_num);
alpar@9: *
alpar@9: *  DESCRIPTION
alpar@9: *
alpar@9: *  The routine glp_weak_comp finds all weakly connected components of
alpar@9: *  the specified graph.
alpar@9: *
alpar@9: *  The parameter v_num specifies an offset of the field of type int
alpar@9: *  in the vertex data block, to which the routine stores the number of
alpar@9: *  a (weakly) connected component containing that vertex. If v_num < 0,
alpar@9: *  no component numbers are stored.
alpar@9: *
alpar@9: *  The components are numbered in arbitrary order from 1 to nc, where
alpar@9: *  nc is the total number of components found, 0 <= nc <= |V|.
alpar@9: *
alpar@9: *  RETURNS
alpar@9: *
alpar@9: *  The routine returns nc, the total number of components found. */
alpar@9: 
alpar@9: int glp_weak_comp(glp_graph *G, int v_num)
alpar@9: {     glp_vertex *v;
alpar@9:       glp_arc *a;
alpar@9:       int f, i, j, nc, nv, pos1, pos2, *prev, *next, *list;
alpar@9:       if (v_num >= 0 && v_num > G->v_size - (int)sizeof(int))
alpar@9:          xerror("glp_weak_comp: v_num = %d; invalid offset\n", v_num);
alpar@9:       nv = G->nv;
alpar@9:       if (nv == 0)
alpar@9:       {  nc = 0;
alpar@9:          goto done;
alpar@9:       }
alpar@9:       /* allocate working arrays */
alpar@9:       prev = xcalloc(1+nv, sizeof(int));
alpar@9:       next = xcalloc(1+nv, sizeof(int));
alpar@9:       list = xcalloc(1+nv, sizeof(int));
alpar@9:       /* if vertex i is unlabelled, prev[i] is the index of previous
alpar@9:          unlabelled vertex, and next[i] is the index of next unlabelled
alpar@9:          vertex; if vertex i is labelled, then prev[i] < 0, and next[i]
alpar@9:          is the connected component number */
alpar@9:       /* initially all vertices are unlabelled */
alpar@9:       f = 1;
alpar@9:       for (i = 1; i <= nv; i++)
alpar@9:          prev[i] = i - 1, next[i] = i + 1;
alpar@9:       next[nv] = 0;
alpar@9:       /* main loop (until all vertices have been labelled) */
alpar@9:       nc = 0;
alpar@9:       while (f != 0)
alpar@9:       {  /* take an unlabelled vertex */
alpar@9:          i = f;
alpar@9:          /* and remove it from the list of unlabelled vertices */
alpar@9:          f = next[i];
alpar@9:          if (f != 0) prev[f] = 0;
alpar@9:          /* label the vertex; it begins a new component */
alpar@9:          prev[i] = -1, next[i] = ++nc;
alpar@9:          /* breadth first search */
alpar@9:          list[1] = i, pos1 = pos2 = 1;
alpar@9:          while (pos1 <= pos2)
alpar@9:          {  /* dequeue vertex i */
alpar@9:             i = list[pos1++];
alpar@9:             /* consider all arcs incoming to vertex i */
alpar@9:             for (a = G->v[i]->in; a != NULL; a = a->h_next)
alpar@9:             {  /* vertex j is adjacent to vertex i */
alpar@9:                j = a->tail->i;
alpar@9:                if (prev[j] >= 0)
alpar@9:                {  /* vertex j is unlabelled */
alpar@9:                   /* remove it from the list of unlabelled vertices */
alpar@9:                   if (prev[j] == 0)
alpar@9:                      f = next[j];
alpar@9:                   else
alpar@9:                      next[prev[j]] = next[j];
alpar@9:                   if (next[j] == 0)
alpar@9:                      ;
alpar@9:                   else
alpar@9:                      prev[next[j]] = prev[j];
alpar@9:                   /* label the vertex */
alpar@9:                   prev[j] = -1, next[j] = nc;
alpar@9:                   /* and enqueue it for further consideration */
alpar@9:                   list[++pos2] = j;
alpar@9:                }
alpar@9:             }
alpar@9:             /* consider all arcs outgoing from vertex i */
alpar@9:             for (a = G->v[i]->out; a != NULL; a = a->t_next)
alpar@9:             {  /* vertex j is adjacent to vertex i */
alpar@9:                j = a->head->i;
alpar@9:                if (prev[j] >= 0)
alpar@9:                {  /* vertex j is unlabelled */
alpar@9:                   /* remove it from the list of unlabelled vertices */
alpar@9:                   if (prev[j] == 0)
alpar@9:                      f = next[j];
alpar@9:                   else
alpar@9:                      next[prev[j]] = next[j];
alpar@9:                   if (next[j] == 0)
alpar@9:                      ;
alpar@9:                   else
alpar@9:                      prev[next[j]] = prev[j];
alpar@9:                   /* label the vertex */
alpar@9:                   prev[j] = -1, next[j] = nc;
alpar@9:                   /* and enqueue it for further consideration */
alpar@9:                   list[++pos2] = j;
alpar@9:                }
alpar@9:             }
alpar@9:          }
alpar@9:       }
alpar@9:       /* store component numbers */
alpar@9:       if (v_num >= 0)
alpar@9:       {  for (i = 1; i <= nv; i++)
alpar@9:          {  v = G->v[i];
alpar@9:             memcpy((char *)v->data + v_num, &next[i], sizeof(int));
alpar@9:          }
alpar@9:       }
alpar@9:       /* free working arrays */
alpar@9:       xfree(prev);
alpar@9:       xfree(next);
alpar@9:       xfree(list);
alpar@9: done: return nc;
alpar@9: }
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  NAME
alpar@9: *
alpar@9: *  glp_strong_comp - find all strongly connected components of graph
alpar@9: *
alpar@9: *  SYNOPSIS
alpar@9: *
alpar@9: *  int glp_strong_comp(glp_graph *G, int v_num);
alpar@9: *
alpar@9: *  DESCRIPTION
alpar@9: *
alpar@9: *  The routine glp_strong_comp finds all strongly connected components
alpar@9: *  of the specified graph.
alpar@9: *
alpar@9: *  The parameter v_num specifies an offset of the field of type int
alpar@9: *  in the vertex data block, to which the routine stores the number of
alpar@9: *  a strongly connected component containing that vertex. If v_num < 0,
alpar@9: *  no component numbers are stored.
alpar@9: *
alpar@9: *  The components are numbered in arbitrary order from 1 to nc, where
alpar@9: *  nc is the total number of components found, 0 <= nc <= |V|. However,
alpar@9: *  the component numbering has the property that for every arc (i->j)
alpar@9: *  in the graph the condition num(i) >= num(j) holds.
alpar@9: *
alpar@9: *  RETURNS
alpar@9: *
alpar@9: *  The routine returns nc, the total number of components found. */
alpar@9: 
alpar@9: int glp_strong_comp(glp_graph *G, int v_num)
alpar@9: {     glp_vertex *v;
alpar@9:       glp_arc *a;
alpar@9:       int i, k, last, n, na, nc, *icn, *ip, *lenr, *ior, *ib, *lowl,
alpar@9:          *numb, *prev;
alpar@9:       if (v_num >= 0 && v_num > G->v_size - (int)sizeof(int))
alpar@9:          xerror("glp_strong_comp: v_num = %d; invalid offset\n",
alpar@9:             v_num);
alpar@9:       n = G->nv;
alpar@9:       if (n == 0)
alpar@9:       {  nc = 0;
alpar@9:          goto done;
alpar@9:       }
alpar@9:       na = G->na;
alpar@9:       icn = xcalloc(1+na, sizeof(int));
alpar@9:       ip = xcalloc(1+n, sizeof(int));
alpar@9:       lenr = xcalloc(1+n, sizeof(int));
alpar@9:       ior = xcalloc(1+n, sizeof(int));
alpar@9:       ib = xcalloc(1+n, sizeof(int));
alpar@9:       lowl = xcalloc(1+n, sizeof(int));
alpar@9:       numb = xcalloc(1+n, sizeof(int));
alpar@9:       prev = xcalloc(1+n, sizeof(int));
alpar@9:       k = 1;
alpar@9:       for (i = 1; i <= n; i++)
alpar@9:       {  v = G->v[i];
alpar@9:          ip[i] = k;
alpar@9:          for (a = v->out; a != NULL; a = a->t_next)
alpar@9:             icn[k++] = a->head->i;
alpar@9:          lenr[i] = k - ip[i];
alpar@9:       }
alpar@9:       xassert(na == k-1);
alpar@9:       nc = mc13d(n, icn, ip, lenr, ior, ib, lowl, numb, prev);
alpar@9:       if (v_num >= 0)
alpar@9:       {  xassert(ib[1] == 1);
alpar@9:          for (k = 1; k <= nc; k++)
alpar@9:          {  last = (k < nc ? ib[k+1] : n+1);
alpar@9:             xassert(ib[k] < last);
alpar@9:             for (i = ib[k]; i < last; i++)
alpar@9:             {  v = G->v[ior[i]];
alpar@9:                memcpy((char *)v->data + v_num, &k, sizeof(int));
alpar@9:             }
alpar@9:          }
alpar@9:       }
alpar@9:       xfree(icn);
alpar@9:       xfree(ip);
alpar@9:       xfree(lenr);
alpar@9:       xfree(ior);
alpar@9:       xfree(ib);
alpar@9:       xfree(lowl);
alpar@9:       xfree(numb);
alpar@9:       xfree(prev);
alpar@9: done: return nc;
alpar@9: }
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  NAME
alpar@9: *
alpar@9: *  glp_top_sort - topological sorting of acyclic digraph
alpar@9: *
alpar@9: *  SYNOPSIS
alpar@9: *
alpar@9: *  int glp_top_sort(glp_graph *G, int v_num);
alpar@9: *
alpar@9: *  DESCRIPTION
alpar@9: *
alpar@9: *  The routine glp_top_sort performs topological sorting of vertices of
alpar@9: *  the specified acyclic digraph.
alpar@9: *
alpar@9: *  The parameter v_num specifies an offset of the field of type int in
alpar@9: *  the vertex data block, to which the routine stores the vertex number
alpar@9: *  assigned. If v_num < 0, vertex numbers are not stored.
alpar@9: *
alpar@9: *  The vertices are numbered from 1 to n, where n is the total number
alpar@9: *  of vertices in the graph. The vertex numbering has the property that
alpar@9: *  for every arc (i->j) in the graph the condition num(i) < num(j)
alpar@9: *  holds. Special case num(i) = 0 means that vertex i is not assigned a
alpar@9: *  number, because the graph is *not* acyclic.
alpar@9: *
alpar@9: *  RETURNS
alpar@9: *
alpar@9: *  If the graph is acyclic and therefore all the vertices have been
alpar@9: *  assigned numbers, the routine glp_top_sort returns zero. Otherwise,
alpar@9: *  if the graph is not acyclic, the routine returns the number of
alpar@9: *  vertices which have not been numbered, i.e. for which num(i) = 0. */
alpar@9: 
alpar@9: static int top_sort(glp_graph *G, int num[])
alpar@9: {     glp_arc *a;
alpar@9:       int i, j, cnt, top, *stack, *indeg;
alpar@9:       /* allocate working arrays */
alpar@9:       indeg = xcalloc(1+G->nv, sizeof(int));
alpar@9:       stack = xcalloc(1+G->nv, sizeof(int));
alpar@9:       /* determine initial indegree of each vertex; push into the stack
alpar@9:          the vertices having zero indegree */
alpar@9:       top = 0;
alpar@9:       for (i = 1; i <= G->nv; i++)
alpar@9:       {  num[i] = indeg[i] = 0;
alpar@9:          for (a = G->v[i]->in; a != NULL; a = a->h_next)
alpar@9:             indeg[i]++;
alpar@9:          if (indeg[i] == 0)
alpar@9:             stack[++top] = i;
alpar@9:       }
alpar@9:       /* assign numbers to vertices in the sorted order */
alpar@9:       cnt = 0;
alpar@9:       while (top > 0)
alpar@9:       {  /* pull vertex i from the stack */
alpar@9:          i = stack[top--];
alpar@9:          /* it has zero indegree in the current graph */
alpar@9:          xassert(indeg[i] == 0);
alpar@9:          /* so assign it a next number */
alpar@9:          xassert(num[i] == 0);
alpar@9:          num[i] = ++cnt;
alpar@9:          /* remove vertex i from the current graph, update indegree of
alpar@9:             its adjacent vertices, and push into the stack new vertices
alpar@9:             whose indegree becomes zero */
alpar@9:          for (a = G->v[i]->out; a != NULL; a = a->t_next)
alpar@9:          {  j = a->head->i;
alpar@9:             /* there exists arc (i->j) in the graph */
alpar@9:             xassert(indeg[j] > 0);
alpar@9:             indeg[j]--;
alpar@9:             if (indeg[j] == 0)
alpar@9:                stack[++top] = j;
alpar@9:          }
alpar@9:       }
alpar@9:       /* free working arrays */
alpar@9:       xfree(indeg);
alpar@9:       xfree(stack);
alpar@9:       return G->nv - cnt;
alpar@9: }
alpar@9: 
alpar@9: int glp_top_sort(glp_graph *G, int v_num)
alpar@9: {     glp_vertex *v;
alpar@9:       int i, cnt, *num;
alpar@9:       if (v_num >= 0 && v_num > G->v_size - (int)sizeof(int))
alpar@9:          xerror("glp_top_sort: v_num = %d; invalid offset\n", v_num);
alpar@9:       if (G->nv == 0)
alpar@9:       {  cnt = 0;
alpar@9:          goto done;
alpar@9:       }
alpar@9:       num = xcalloc(1+G->nv, sizeof(int));
alpar@9:       cnt = top_sort(G, num);
alpar@9:       if (v_num >= 0)
alpar@9:       {  for (i = 1; i <= G->nv; i++)
alpar@9:          {  v = G->v[i];
alpar@9:             memcpy((char *)v->data + v_num, &num[i], sizeof(int));
alpar@9:          }
alpar@9:       }
alpar@9:       xfree(num);
alpar@9: done: return cnt;
alpar@9: }
alpar@9: 
alpar@9: /* eof */