/* glpapi16.c (graph and network analysis routines) */

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

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

 *

 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

 *  E-mail: <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 "glpapi.h"

 #include "glpnet.h"

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

 *  NAME

 *

 *  glp_weak_comp - find all weakly connected components of graph

 *

 *  SYNOPSIS

 *

 *  int glp_weak_comp(glp_graph *G, int v_num);

 *

 *  DESCRIPTION

 *

 *  The routine glp_weak_comp finds all weakly connected components of

 *  the specified graph.

 *

 *  The parameter v_num specifies an offset of the field of type int

 *  in the vertex data block, to which the routine stores the number of

 *  a (weakly) connected component containing that vertex. If v_num < 0,

 *  no component numbers are stored.

 *

 *  The components are numbered in arbitrary order from 1 to nc, where

 *  nc is the total number of components found, 0 <= nc <= |V|.

 *

 *  RETURNS

 *

 *  The routine returns nc, the total number of components found. */

 int glp_weak_comp(glp_graph *G, int v_num)

 {     glp_vertex *v;

       glp_arc *a;

       int f, i, j, nc, nv, pos1, pos2, *prev, *next, *list;

       if (v_num >= 0 && v_num > G->v_size - (int)sizeof(int))

          xerror("glp_weak_comp: v_num = %d; invalid offset\n", v_num);

       nv = G->nv;

       if (nv == 0)

       {  nc = 0;

          goto done;

       }

       /* allocate working arrays */

       prev = xcalloc(1+nv, sizeof(int));

       next = xcalloc(1+nv, sizeof(int));

       list = xcalloc(1+nv, sizeof(int));

       /* if vertex i is unlabelled, prev[i] is the index of previous

          unlabelled vertex, and next[i] is the index of next unlabelled

          vertex; if vertex i is labelled, then prev[i] < 0, and next[i]

          is the connected component number */

       /* initially all vertices are unlabelled */

       f = 1;

       for (i = 1; i <= nv; i++)

          prev[i] = i - 1, next[i] = i + 1;

       next[nv] = 0;

       /* main loop (until all vertices have been labelled) */

       nc = 0;

       while (f != 0)

       {  /* take an unlabelled vertex */

          i = f;

          /* and remove it from the list of unlabelled vertices */

          f = next[i];

          if (f != 0) prev[f] = 0;

          /* label the vertex; it begins a new component */

          prev[i] = -1, next[i] = ++nc;

          /* breadth first search */

          list[1] = i, pos1 = pos2 = 1;

          while (pos1 <= pos2)

          {  /* dequeue vertex i */

             i = list[pos1++];

             /* consider all arcs incoming to vertex i */

             for (a = G->v[i]->in; a != NULL; a = a->h_next)

             {  /* vertex j is adjacent to vertex i */

                j = a->tail->i;

                if (prev[j] >= 0)

                {  /* vertex j is unlabelled */

                   /* remove it from the list of unlabelled vertices */

                   if (prev[j] == 0)

                      f = next[j];

                   else

                      next[prev[j]] = next[j];

                   if (next[j] == 0)

                      ;

                   else

                      prev[next[j]] = prev[j];

                   /* label the vertex */

                   prev[j] = -1, next[j] = nc;

                   /* and enqueue it for further consideration */

                   list[++pos2] = j;

                }

             }

             /* consider all arcs outgoing from vertex i */

             for (a = G->v[i]->out; a != NULL; a = a->t_next)

             {  /* vertex j is adjacent to vertex i */

                j = a->head->i;

                if (prev[j] >= 0)

                {  /* vertex j is unlabelled */

                   /* remove it from the list of unlabelled vertices */

                   if (prev[j] == 0)

                      f = next[j];

                   else

                      next[prev[j]] = next[j];

                   if (next[j] == 0)

                      ;

                   else

                      prev[next[j]] = prev[j];

                   /* label the vertex */

                   prev[j] = -1, next[j] = nc;

                   /* and enqueue it for further consideration */

                   list[++pos2] = j;

                }

             }

          }

       }

       /* store component numbers */

       if (v_num >= 0)

       {  for (i = 1; i <= nv; i++)

          {  v = G->v[i];

             memcpy((char *)v->data + v_num, &next[i], sizeof(int));

          }

       }

       /* free working arrays */

       xfree(prev);

       xfree(next);

       xfree(list);

 done: return nc;

 }

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

 *  NAME

 *

 *  glp_strong_comp - find all strongly connected components of graph

 *

 *  SYNOPSIS

 *

 *  int glp_strong_comp(glp_graph *G, int v_num);

 *

 *  DESCRIPTION

 *

 *  The routine glp_strong_comp finds all strongly connected components

 *  of the specified graph.

 *

 *  The parameter v_num specifies an offset of the field of type int

 *  in the vertex data block, to which the routine stores the number of

 *  a strongly connected component containing that vertex. If v_num < 0,

 *  no component numbers are stored.

 *

 *  The components are numbered in arbitrary order from 1 to nc, where

 *  nc is the total number of components found, 0 <= nc <= |V|. However,

 *  the component numbering has the property that for every arc (i->j)

 *  in the graph the condition num(i) >= num(j) holds.

 *

 *  RETURNS

 *

 *  The routine returns nc, the total number of components found. */

 int glp_strong_comp(glp_graph *G, int v_num)

 {     glp_vertex *v;

       glp_arc *a;

       int i, k, last, n, na, nc, *icn, *ip, *lenr, *ior, *ib, *lowl,

          *numb, *prev;

       if (v_num >= 0 && v_num > G->v_size - (int)sizeof(int))

          xerror("glp_strong_comp: v_num = %d; invalid offset\n",

             v_num);

       n = G->nv;

       if (n == 0)

       {  nc = 0;

          goto done;

       }

       na = G->na;

       icn = xcalloc(1+na, sizeof(int));

       ip = xcalloc(1+n, sizeof(int));

       lenr = xcalloc(1+n, sizeof(int));

       ior = xcalloc(1+n, sizeof(int));

       ib = xcalloc(1+n, sizeof(int));

       lowl = xcalloc(1+n, sizeof(int));

       numb = xcalloc(1+n, sizeof(int));

       prev = xcalloc(1+n, sizeof(int));

       k = 1;

       for (i = 1; i <= n; i++)

       {  v = G->v[i];

          ip[i] = k;

          for (a = v->out; a != NULL; a = a->t_next)

             icn[k++] = a->head->i;

          lenr[i] = k - ip[i];

       }

       xassert(na == k-1);

       nc = mc13d(n, icn, ip, lenr, ior, ib, lowl, numb, prev);

       if (v_num >= 0)

       {  xassert(ib[1] == 1);

          for (k = 1; k <= nc; k++)

          {  last = (k < nc ? ib[k+1] : n+1);

             xassert(ib[k] < last);

             for (i = ib[k]; i < last; i++)

             {  v = G->v[ior[i]];

                memcpy((char *)v->data + v_num, &k, sizeof(int));

             }

          }

       }

       xfree(icn);

       xfree(ip);

       xfree(lenr);

       xfree(ior);

       xfree(ib);

       xfree(lowl);

       xfree(numb);

       xfree(prev);

 done: return nc;

 }

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

 *  NAME

 *

 *  glp_top_sort - topological sorting of acyclic digraph

 *

 *  SYNOPSIS

 *

 *  int glp_top_sort(glp_graph *G, int v_num);

 *

 *  DESCRIPTION

 *

 *  The routine glp_top_sort performs topological sorting of vertices of

 *  the specified acyclic digraph.

 *

 *  The parameter v_num specifies an offset of the field of type int in

 *  the vertex data block, to which the routine stores the vertex number

 *  assigned. If v_num < 0, vertex numbers are not stored.

 *

 *  The vertices are numbered from 1 to n, where n is the total number

 *  of vertices in the graph. The vertex numbering has the property that

 *  for every arc (i->j) in the graph the condition num(i) < num(j)

 *  holds. Special case num(i) = 0 means that vertex i is not assigned a

 *  number, because the graph is *not* acyclic.

 *

 *  RETURNS

 *

 *  If the graph is acyclic and therefore all the vertices have been

 *  assigned numbers, the routine glp_top_sort returns zero. Otherwise,

 *  if the graph is not acyclic, the routine returns the number of

 *  vertices which have not been numbered, i.e. for which num(i) = 0. */

 static int top_sort(glp_graph *G, int num[])

 {     glp_arc *a;

       int i, j, cnt, top, *stack, *indeg;

       /* allocate working arrays */

       indeg = xcalloc(1+G->nv, sizeof(int));

       stack = xcalloc(1+G->nv, sizeof(int));

       /* determine initial indegree of each vertex; push into the stack

          the vertices having zero indegree */

       top = 0;

       for (i = 1; i <= G->nv; i++)

       {  num[i] = indeg[i] = 0;

          for (a = G->v[i]->in; a != NULL; a = a->h_next)

             indeg[i]++;

          if (indeg[i] == 0)

             stack[++top] = i;

       }

       /* assign numbers to vertices in the sorted order */

       cnt = 0;

       while (top > 0)

       {  /* pull vertex i from the stack */

          i = stack[top--];

          /* it has zero indegree in the current graph */

          xassert(indeg[i] == 0);

          /* so assign it a next number */

          xassert(num[i] == 0);

          num[i] = ++cnt;

          /* remove vertex i from the current graph, update indegree of

             its adjacent vertices, and push into the stack new vertices

             whose indegree becomes zero */

          for (a = G->v[i]->out; a != NULL; a = a->t_next)

          {  j = a->head->i;

             /* there exists arc (i->j) in the graph */

             xassert(indeg[j] > 0);

             indeg[j]--;

             if (indeg[j] == 0)

                stack[++top] = j;

          }

       }

       /* free working arrays */

       xfree(indeg);

       xfree(stack);

       return G->nv - cnt;

 }

 int glp_top_sort(glp_graph *G, int v_num)

 {     glp_vertex *v;

       int i, cnt, *num;

       if (v_num >= 0 && v_num > G->v_size - (int)sizeof(int))

          xerror("glp_top_sort: v_num = %d; invalid offset\n", v_num);

       if (G->nv == 0)

       {  cnt = 0;

          goto done;

       }

       num = xcalloc(1+G->nv, sizeof(int));

       cnt = top_sort(G, num);

       if (v_num >= 0)

       {  for (i = 1; i <= G->nv; i++)

          {  v = G->v[i];

             memcpy((char *)v->data + v_num, &num[i], sizeof(int));

          }

       }

       xfree(num);

 done: return cnt;

 }

 /* eof */