/* glpnet07.c (Ford-Fulkerson algorithm) */

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

 *  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 "glpenv.h"

 #include "glpnet.h"

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

 *  NAME

 *

 *  ffalg - Ford-Fulkerson algorithm

 *

 *  SYNOPSIS

 *

 *  #include "glpnet.h"

 *  void ffalg(int nv, int na, const int tail[], const int head[],

 *     int s, int t, const int cap[], int x[], char cut[]);

 *

 *  DESCRIPTION

 *

 *  The routine ffalg implements the Ford-Fulkerson algorithm to find a

 *  maximal flow in the specified flow network.

 *

 *  INPUT PARAMETERS

 *

 *  nv is the number of nodes, nv >= 2.

 *

 *  na is the number of arcs, na >= 0.

 *

 *  tail[a], a = 1,...,na, is the index of tail node of arc a.

 *

 *  head[a], a = 1,...,na, is the index of head node of arc a.

 *

 *  s is the source node index, 1 <= s <= nv.

 *

 *  t is the sink node index, 1 <= t <= nv, t != s.

 *

 *  cap[a], a = 1,...,na, is the capacity of arc a, cap[a] >= 0.

 *

 *  NOTE: Multiple arcs are allowed, but self-loops are not allowed.

 *

 *  OUTPUT PARAMETERS

 *

 *  x[a], a = 1,...,na, is optimal value of the flow through arc a.

 *

 *  cut[i], i = 1,...,nv, is 1 if node i is labelled, and 0 otherwise.

 *  The set of arcs, whose one endpoint is labelled and other is not,

 *  defines the minimal cut corresponding to the maximal flow found.

 *  If the parameter cut is NULL, the cut information are not stored.

 *

 *  REFERENCES

 *

 *  L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND

 *  Corp., Report R-375-PR (August 1962), Chap. I "Static Maximal Flow,"

 *  pp.30-33. */

 void ffalg(int nv, int na, const int tail[], const int head[],

       int s, int t, const int cap[], int x[], char cut[])

 {     int a, delta, i, j, k, pos1, pos2, temp,

          *ptr, *arc, *link, *list;

       /* sanity checks */

       xassert(nv >= 2);

       xassert(na >= 0);

       xassert(1 <= s && s <= nv);

       xassert(1 <= t && t <= nv);

       xassert(s != t);

       for (a = 1; a <= na; a++)

       {  i = tail[a], j = head[a];

          xassert(1 <= i && i <= nv);

          xassert(1 <= j && j <= nv);

          xassert(i != j);

          xassert(cap[a] >= 0);

       }

       /* allocate working arrays */

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

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

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

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

       /* ptr[i] := (degree of node i) */

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

          ptr[i] = 0;

       for (a = 1; a <= na; a++)

       {  ptr[tail[a]]++;

          ptr[head[a]]++;

       }

       /* initialize arc pointers */

       ptr[1]++;

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

          ptr[i+1] += ptr[i];

       ptr[nv+1] = ptr[nv];

       /* build arc lists */

       for (a = 1; a <= na; a++)

       {  arc[--ptr[tail[a]]] = a;

          arc[--ptr[head[a]]] = a;

       }

       xassert(ptr[1] == 1);

       xassert(ptr[nv+1] == na+na+1);

       /* now the indices of arcs incident to node i are stored in

          locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */

       /* initialize arc flows */

       for (a = 1; a <= na; a++)

          x[a] = 0;

 loop: /* main loop starts here */

       /* build augmenting tree rooted at s */

       /* link[i] = 0 means that node i is not labelled yet;

          link[i] = a means that arc a immediately precedes node i */

       /* initially node s is labelled as the root */

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

          link[i] = 0;

       link[s] = -1, list[1] = s, pos1 = pos2 = 1;

       /* breadth first search */

       while (pos1 <= pos2)

       {  /* dequeue node i */

          i = list[pos1++];

          /* consider all arcs incident to node i */

          for (k = ptr[i]; k < ptr[i+1]; k++)

          {  a = arc[k];

             if (tail[a] == i)

             {  /* a = i->j is a forward arc from s to t */

                j = head[a];

                /* if node j has been labelled, skip the arc */

                if (link[j] != 0) continue;

                /* if the arc does not allow increasing the flow through

                   it, skip the arc */

                if (x[a] == cap[a]) continue;

             }

             else if (head[a] == i)

             {  /* a = i<-j is a backward arc from s to t */

                j = tail[a];

                /* if node j has been labelled, skip the arc */

                if (link[j] != 0) continue;

                /* if the arc does not allow decreasing the flow through

                   it, skip the arc */

                if (x[a] == 0) continue;

             }

             else

                xassert(a != a);

             /* label node j and enqueue it */

             link[j] = a, list[++pos2] = j;

             /* check for breakthrough */

             if (j == t) goto brkt;

          }

       }

       /* NONBREAKTHROUGH */

       /* no augmenting path exists; current flow is maximal */

       /* store minimal cut information, if necessary */

       if (cut != NULL)

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

             cut[i] = (char)(link[i] != 0);

       }

       goto done;

 brkt: /* BREAKTHROUGH */

       /* walk through arcs of the augmenting path (s, ..., t) found in

          the reverse order and determine maximal change of the flow */

       delta = 0;

       for (j = t; j != s; j = i)

       {  /* arc a immediately precedes node j in the path */

          a = link[j];

          if (head[a] == j)

          {  /* a = i->j is a forward arc of the cycle */

             i = tail[a];

             /* x[a] may be increased until its upper bound */

             temp = cap[a] - x[a];

          }

          else if (tail[a] == j)

          {  /* a = i<-j is a backward arc of the cycle */

             i = head[a];

             /* x[a] may be decreased until its lower bound */

             temp = x[a];

          }

          else

             xassert(a != a);

          if (delta == 0 || delta > temp) delta = temp;

       }

       xassert(delta > 0);

       /* increase the flow along the path */

       for (j = t; j != s; j = i)

       {  /* arc a immediately precedes node j in the path */

          a = link[j];

          if (head[a] == j)

          {  /* a = i->j is a forward arc of the cycle */

             i = tail[a];

             x[a] += delta;

          }

          else if (tail[a] == j)

          {  /* a = i<-j is a backward arc of the cycle */

             i = head[a];

             x[a] -= delta;

          }

          else

             xassert(a != a);

       }

       goto loop;

 done: /* free working arrays */

       xfree(ptr);

       xfree(arc);

       xfree(link);

       xfree(list);

       return;

 }

 /* eof */