alpar@9: /* glpnet07.c (Ford-Fulkerson algorithm) */
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 "glpenv.h"
alpar@9: #include "glpnet.h"
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  NAME
alpar@9: *
alpar@9: *  ffalg - Ford-Fulkerson algorithm
alpar@9: *
alpar@9: *  SYNOPSIS
alpar@9: *
alpar@9: *  #include "glpnet.h"
alpar@9: *  void ffalg(int nv, int na, const int tail[], const int head[],
alpar@9: *     int s, int t, const int cap[], int x[], char cut[]);
alpar@9: *
alpar@9: *  DESCRIPTION
alpar@9: *
alpar@9: *  The routine ffalg implements the Ford-Fulkerson algorithm to find a
alpar@9: *  maximal flow in the specified flow network.
alpar@9: *
alpar@9: *  INPUT PARAMETERS
alpar@9: *
alpar@9: *  nv is the number of nodes, nv >= 2.
alpar@9: *
alpar@9: *  na is the number of arcs, na >= 0.
alpar@9: *
alpar@9: *  tail[a], a = 1,...,na, is the index of tail node of arc a.
alpar@9: *
alpar@9: *  head[a], a = 1,...,na, is the index of head node of arc a.
alpar@9: *
alpar@9: *  s is the source node index, 1 <= s <= nv.
alpar@9: *
alpar@9: *  t is the sink node index, 1 <= t <= nv, t != s.
alpar@9: *
alpar@9: *  cap[a], a = 1,...,na, is the capacity of arc a, cap[a] >= 0.
alpar@9: *
alpar@9: *  NOTE: Multiple arcs are allowed, but self-loops are not allowed.
alpar@9: *
alpar@9: *  OUTPUT PARAMETERS
alpar@9: *
alpar@9: *  x[a], a = 1,...,na, is optimal value of the flow through arc a.
alpar@9: *
alpar@9: *  cut[i], i = 1,...,nv, is 1 if node i is labelled, and 0 otherwise.
alpar@9: *  The set of arcs, whose one endpoint is labelled and other is not,
alpar@9: *  defines the minimal cut corresponding to the maximal flow found.
alpar@9: *  If the parameter cut is NULL, the cut information are not stored.
alpar@9: *
alpar@9: *  REFERENCES
alpar@9: *
alpar@9: *  L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND
alpar@9: *  Corp., Report R-375-PR (August 1962), Chap. I "Static Maximal Flow,"
alpar@9: *  pp.30-33. */
alpar@9: 
alpar@9: void ffalg(int nv, int na, const int tail[], const int head[],
alpar@9:       int s, int t, const int cap[], int x[], char cut[])
alpar@9: {     int a, delta, i, j, k, pos1, pos2, temp,
alpar@9:          *ptr, *arc, *link, *list;
alpar@9:       /* sanity checks */
alpar@9:       xassert(nv >= 2);
alpar@9:       xassert(na >= 0);
alpar@9:       xassert(1 <= s && s <= nv);
alpar@9:       xassert(1 <= t && t <= nv);
alpar@9:       xassert(s != t);
alpar@9:       for (a = 1; a <= na; a++)
alpar@9:       {  i = tail[a], j = head[a];
alpar@9:          xassert(1 <= i && i <= nv);
alpar@9:          xassert(1 <= j && j <= nv);
alpar@9:          xassert(i != j);
alpar@9:          xassert(cap[a] >= 0);
alpar@9:       }
alpar@9:       /* allocate working arrays */
alpar@9:       ptr = xcalloc(1+nv+1, sizeof(int));
alpar@9:       arc = xcalloc(1+na+na, sizeof(int));
alpar@9:       link = xcalloc(1+nv, sizeof(int));
alpar@9:       list = xcalloc(1+nv, sizeof(int));
alpar@9:       /* ptr[i] := (degree of node i) */
alpar@9:       for (i = 1; i <= nv; i++)
alpar@9:          ptr[i] = 0;
alpar@9:       for (a = 1; a <= na; a++)
alpar@9:       {  ptr[tail[a]]++;
alpar@9:          ptr[head[a]]++;
alpar@9:       }
alpar@9:       /* initialize arc pointers */
alpar@9:       ptr[1]++;
alpar@9:       for (i = 1; i < nv; i++)
alpar@9:          ptr[i+1] += ptr[i];
alpar@9:       ptr[nv+1] = ptr[nv];
alpar@9:       /* build arc lists */
alpar@9:       for (a = 1; a <= na; a++)
alpar@9:       {  arc[--ptr[tail[a]]] = a;
alpar@9:          arc[--ptr[head[a]]] = a;
alpar@9:       }
alpar@9:       xassert(ptr[1] == 1);
alpar@9:       xassert(ptr[nv+1] == na+na+1);
alpar@9:       /* now the indices of arcs incident to node i are stored in
alpar@9:          locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */
alpar@9:       /* initialize arc flows */
alpar@9:       for (a = 1; a <= na; a++)
alpar@9:          x[a] = 0;
alpar@9: loop: /* main loop starts here */
alpar@9:       /* build augmenting tree rooted at s */
alpar@9:       /* link[i] = 0 means that node i is not labelled yet;
alpar@9:          link[i] = a means that arc a immediately precedes node i */
alpar@9:       /* initially node s is labelled as the root */
alpar@9:       for (i = 1; i <= nv; i++)
alpar@9:          link[i] = 0;
alpar@9:       link[s] = -1, list[1] = s, pos1 = pos2 = 1;
alpar@9:       /* breadth first search */
alpar@9:       while (pos1 <= pos2)
alpar@9:       {  /* dequeue node i */
alpar@9:          i = list[pos1++];
alpar@9:          /* consider all arcs incident to node i */
alpar@9:          for (k = ptr[i]; k < ptr[i+1]; k++)
alpar@9:          {  a = arc[k];
alpar@9:             if (tail[a] == i)
alpar@9:             {  /* a = i->j is a forward arc from s to t */
alpar@9:                j = head[a];
alpar@9:                /* if node j has been labelled, skip the arc */
alpar@9:                if (link[j] != 0) continue;
alpar@9:                /* if the arc does not allow increasing the flow through
alpar@9:                   it, skip the arc */
alpar@9:                if (x[a] == cap[a]) continue;
alpar@9:             }
alpar@9:             else if (head[a] == i)
alpar@9:             {  /* a = i<-j is a backward arc from s to t */
alpar@9:                j = tail[a];
alpar@9:                /* if node j has been labelled, skip the arc */
alpar@9:                if (link[j] != 0) continue;
alpar@9:                /* if the arc does not allow decreasing the flow through
alpar@9:                   it, skip the arc */
alpar@9:                if (x[a] == 0) continue;
alpar@9:             }
alpar@9:             else
alpar@9:                xassert(a != a);
alpar@9:             /* label node j and enqueue it */
alpar@9:             link[j] = a, list[++pos2] = j;
alpar@9:             /* check for breakthrough */
alpar@9:             if (j == t) goto brkt;
alpar@9:          }
alpar@9:       }
alpar@9:       /* NONBREAKTHROUGH */
alpar@9:       /* no augmenting path exists; current flow is maximal */
alpar@9:       /* store minimal cut information, if necessary */
alpar@9:       if (cut != NULL)
alpar@9:       {  for (i = 1; i <= nv; i++)
alpar@9:             cut[i] = (char)(link[i] != 0);
alpar@9:       }
alpar@9:       goto done;
alpar@9: brkt: /* BREAKTHROUGH */
alpar@9:       /* walk through arcs of the augmenting path (s, ..., t) found in
alpar@9:          the reverse order and determine maximal change of the flow */
alpar@9:       delta = 0;
alpar@9:       for (j = t; j != s; j = i)
alpar@9:       {  /* arc a immediately precedes node j in the path */
alpar@9:          a = link[j];
alpar@9:          if (head[a] == j)
alpar@9:          {  /* a = i->j is a forward arc of the cycle */
alpar@9:             i = tail[a];
alpar@9:             /* x[a] may be increased until its upper bound */
alpar@9:             temp = cap[a] - x[a];
alpar@9:          }
alpar@9:          else if (tail[a] == j)
alpar@9:          {  /* a = i<-j is a backward arc of the cycle */
alpar@9:             i = head[a];
alpar@9:             /* x[a] may be decreased until its lower bound */
alpar@9:             temp = x[a];
alpar@9:          }
alpar@9:          else
alpar@9:             xassert(a != a);
alpar@9:          if (delta == 0 || delta > temp) delta = temp;
alpar@9:       }
alpar@9:       xassert(delta > 0);
alpar@9:       /* increase the flow along the path */
alpar@9:       for (j = t; j != s; j = i)
alpar@9:       {  /* arc a immediately precedes node j in the path */
alpar@9:          a = link[j];
alpar@9:          if (head[a] == j)
alpar@9:          {  /* a = i->j is a forward arc of the cycle */
alpar@9:             i = tail[a];
alpar@9:             x[a] += delta;
alpar@9:          }
alpar@9:          else if (tail[a] == j)
alpar@9:          {  /* a = i<-j is a backward arc of the cycle */
alpar@9:             i = head[a];
alpar@9:             x[a] -= delta;
alpar@9:          }
alpar@9:          else
alpar@9:             xassert(a != a);
alpar@9:       }
alpar@9:       goto loop;
alpar@9: done: /* free working arrays */
alpar@9:       xfree(ptr);
alpar@9:       xfree(arc);
alpar@9:       xfree(link);
alpar@9:       xfree(list);
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: /* eof */