1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnet07.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,222 @@
     1.4 +/* glpnet07.c (Ford-Fulkerson algorithm) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpenv.h"
    1.29 +#include "glpnet.h"
    1.30 +
    1.31 +/***********************************************************************
    1.32 +*  NAME
    1.33 +*
    1.34 +*  ffalg - Ford-Fulkerson algorithm
    1.35 +*
    1.36 +*  SYNOPSIS
    1.37 +*
    1.38 +*  #include "glpnet.h"
    1.39 +*  void ffalg(int nv, int na, const int tail[], const int head[],
    1.40 +*     int s, int t, const int cap[], int x[], char cut[]);
    1.41 +*
    1.42 +*  DESCRIPTION
    1.43 +*
    1.44 +*  The routine ffalg implements the Ford-Fulkerson algorithm to find a
    1.45 +*  maximal flow in the specified flow network.
    1.46 +*
    1.47 +*  INPUT PARAMETERS
    1.48 +*
    1.49 +*  nv is the number of nodes, nv >= 2.
    1.50 +*
    1.51 +*  na is the number of arcs, na >= 0.
    1.52 +*
    1.53 +*  tail[a], a = 1,...,na, is the index of tail node of arc a.
    1.54 +*
    1.55 +*  head[a], a = 1,...,na, is the index of head node of arc a.
    1.56 +*
    1.57 +*  s is the source node index, 1 <= s <= nv.
    1.58 +*
    1.59 +*  t is the sink node index, 1 <= t <= nv, t != s.
    1.60 +*
    1.61 +*  cap[a], a = 1,...,na, is the capacity of arc a, cap[a] >= 0.
    1.62 +*
    1.63 +*  NOTE: Multiple arcs are allowed, but self-loops are not allowed.
    1.64 +*
    1.65 +*  OUTPUT PARAMETERS
    1.66 +*
    1.67 +*  x[a], a = 1,...,na, is optimal value of the flow through arc a.
    1.68 +*
    1.69 +*  cut[i], i = 1,...,nv, is 1 if node i is labelled, and 0 otherwise.
    1.70 +*  The set of arcs, whose one endpoint is labelled and other is not,
    1.71 +*  defines the minimal cut corresponding to the maximal flow found.
    1.72 +*  If the parameter cut is NULL, the cut information are not stored.
    1.73 +*
    1.74 +*  REFERENCES
    1.75 +*
    1.76 +*  L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND
    1.77 +*  Corp., Report R-375-PR (August 1962), Chap. I "Static Maximal Flow,"
    1.78 +*  pp.30-33. */
    1.79 +
    1.80 +void ffalg(int nv, int na, const int tail[], const int head[],
    1.81 +      int s, int t, const int cap[], int x[], char cut[])
    1.82 +{     int a, delta, i, j, k, pos1, pos2, temp,
    1.83 +         *ptr, *arc, *link, *list;
    1.84 +      /* sanity checks */
    1.85 +      xassert(nv >= 2);
    1.86 +      xassert(na >= 0);
    1.87 +      xassert(1 <= s && s <= nv);
    1.88 +      xassert(1 <= t && t <= nv);
    1.89 +      xassert(s != t);
    1.90 +      for (a = 1; a <= na; a++)
    1.91 +      {  i = tail[a], j = head[a];
    1.92 +         xassert(1 <= i && i <= nv);
    1.93 +         xassert(1 <= j && j <= nv);
    1.94 +         xassert(i != j);
    1.95 +         xassert(cap[a] >= 0);
    1.96 +      }
    1.97 +      /* allocate working arrays */
    1.98 +      ptr = xcalloc(1+nv+1, sizeof(int));
    1.99 +      arc = xcalloc(1+na+na, sizeof(int));
   1.100 +      link = xcalloc(1+nv, sizeof(int));
   1.101 +      list = xcalloc(1+nv, sizeof(int));
   1.102 +      /* ptr[i] := (degree of node i) */
   1.103 +      for (i = 1; i <= nv; i++)
   1.104 +         ptr[i] = 0;
   1.105 +      for (a = 1; a <= na; a++)
   1.106 +      {  ptr[tail[a]]++;
   1.107 +         ptr[head[a]]++;
   1.108 +      }
   1.109 +      /* initialize arc pointers */
   1.110 +      ptr[1]++;
   1.111 +      for (i = 1; i < nv; i++)
   1.112 +         ptr[i+1] += ptr[i];
   1.113 +      ptr[nv+1] = ptr[nv];
   1.114 +      /* build arc lists */
   1.115 +      for (a = 1; a <= na; a++)
   1.116 +      {  arc[--ptr[tail[a]]] = a;
   1.117 +         arc[--ptr[head[a]]] = a;
   1.118 +      }
   1.119 +      xassert(ptr[1] == 1);
   1.120 +      xassert(ptr[nv+1] == na+na+1);
   1.121 +      /* now the indices of arcs incident to node i are stored in
   1.122 +         locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */
   1.123 +      /* initialize arc flows */
   1.124 +      for (a = 1; a <= na; a++)
   1.125 +         x[a] = 0;
   1.126 +loop: /* main loop starts here */
   1.127 +      /* build augmenting tree rooted at s */
   1.128 +      /* link[i] = 0 means that node i is not labelled yet;
   1.129 +         link[i] = a means that arc a immediately precedes node i */
   1.130 +      /* initially node s is labelled as the root */
   1.131 +      for (i = 1; i <= nv; i++)
   1.132 +         link[i] = 0;
   1.133 +      link[s] = -1, list[1] = s, pos1 = pos2 = 1;
   1.134 +      /* breadth first search */
   1.135 +      while (pos1 <= pos2)
   1.136 +      {  /* dequeue node i */
   1.137 +         i = list[pos1++];
   1.138 +         /* consider all arcs incident to node i */
   1.139 +         for (k = ptr[i]; k < ptr[i+1]; k++)
   1.140 +         {  a = arc[k];
   1.141 +            if (tail[a] == i)
   1.142 +            {  /* a = i->j is a forward arc from s to t */
   1.143 +               j = head[a];
   1.144 +               /* if node j has been labelled, skip the arc */
   1.145 +               if (link[j] != 0) continue;
   1.146 +               /* if the arc does not allow increasing the flow through
   1.147 +                  it, skip the arc */
   1.148 +               if (x[a] == cap[a]) continue;
   1.149 +            }
   1.150 +            else if (head[a] == i)
   1.151 +            {  /* a = i<-j is a backward arc from s to t */
   1.152 +               j = tail[a];
   1.153 +               /* if node j has been labelled, skip the arc */
   1.154 +               if (link[j] != 0) continue;
   1.155 +               /* if the arc does not allow decreasing the flow through
   1.156 +                  it, skip the arc */
   1.157 +               if (x[a] == 0) continue;
   1.158 +            }
   1.159 +            else
   1.160 +               xassert(a != a);
   1.161 +            /* label node j and enqueue it */
   1.162 +            link[j] = a, list[++pos2] = j;
   1.163 +            /* check for breakthrough */
   1.164 +            if (j == t) goto brkt;
   1.165 +         }
   1.166 +      }
   1.167 +      /* NONBREAKTHROUGH */
   1.168 +      /* no augmenting path exists; current flow is maximal */
   1.169 +      /* store minimal cut information, if necessary */
   1.170 +      if (cut != NULL)
   1.171 +      {  for (i = 1; i <= nv; i++)
   1.172 +            cut[i] = (char)(link[i] != 0);
   1.173 +      }
   1.174 +      goto done;
   1.175 +brkt: /* BREAKTHROUGH */
   1.176 +      /* walk through arcs of the augmenting path (s, ..., t) found in
   1.177 +         the reverse order and determine maximal change of the flow */
   1.178 +      delta = 0;
   1.179 +      for (j = t; j != s; j = i)
   1.180 +      {  /* arc a immediately precedes node j in the path */
   1.181 +         a = link[j];
   1.182 +         if (head[a] == j)
   1.183 +         {  /* a = i->j is a forward arc of the cycle */
   1.184 +            i = tail[a];
   1.185 +            /* x[a] may be increased until its upper bound */
   1.186 +            temp = cap[a] - x[a];
   1.187 +         }
   1.188 +         else if (tail[a] == j)
   1.189 +         {  /* a = i<-j is a backward arc of the cycle */
   1.190 +            i = head[a];
   1.191 +            /* x[a] may be decreased until its lower bound */
   1.192 +            temp = x[a];
   1.193 +         }
   1.194 +         else
   1.195 +            xassert(a != a);
   1.196 +         if (delta == 0 || delta > temp) delta = temp;
   1.197 +      }
   1.198 +      xassert(delta > 0);
   1.199 +      /* increase the flow along the path */
   1.200 +      for (j = t; j != s; j = i)
   1.201 +      {  /* arc a immediately precedes node j in the path */
   1.202 +         a = link[j];
   1.203 +         if (head[a] == j)
   1.204 +         {  /* a = i->j is a forward arc of the cycle */
   1.205 +            i = tail[a];
   1.206 +            x[a] += delta;
   1.207 +         }
   1.208 +         else if (tail[a] == j)
   1.209 +         {  /* a = i<-j is a backward arc of the cycle */
   1.210 +            i = head[a];
   1.211 +            x[a] -= delta;
   1.212 +         }
   1.213 +         else
   1.214 +            xassert(a != a);
   1.215 +      }
   1.216 +      goto loop;
   1.217 +done: /* free working arrays */
   1.218 +      xfree(ptr);
   1.219 +      xfree(arc);
   1.220 +      xfree(link);
   1.221 +      xfree(list);
   1.222 +      return;
   1.223 +}
   1.224 +
   1.225 +/* eof */