/* 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 */