COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpnet07.c

subpack-glpk
Last change on this file was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 7.3 KB
Line 
1/* glpnet07.c (Ford-Fulkerson algorithm) */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
6*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9*  E-mail: <mao@gnu.org>.
10*
11*  GLPK is free software: you can redistribute it and/or modify it
12*  under the terms of the GNU General Public License as published by
13*  the Free Software Foundation, either version 3 of the License, or
14*  (at your option) any later version.
15*
16*  GLPK is distributed in the hope that it will be useful, but WITHOUT
17*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19*  License for more details.
20*
21*  You should have received a copy of the GNU General Public License
22*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23***********************************************************************/
24
25#include "glpenv.h"
26#include "glpnet.h"
27
28/***********************************************************************
29*  NAME
30*
31*  ffalg - Ford-Fulkerson algorithm
32*
33*  SYNOPSIS
34*
35*  #include "glpnet.h"
36*  void ffalg(int nv, int na, const int tail[], const int head[],
37*     int s, int t, const int cap[], int x[], char cut[]);
38*
39*  DESCRIPTION
40*
41*  The routine ffalg implements the Ford-Fulkerson algorithm to find a
42*  maximal flow in the specified flow network.
43*
44*  INPUT PARAMETERS
45*
46*  nv is the number of nodes, nv >= 2.
47*
48*  na is the number of arcs, na >= 0.
49*
50*  tail[a], a = 1,...,na, is the index of tail node of arc a.
51*
52*  head[a], a = 1,...,na, is the index of head node of arc a.
53*
54*  s is the source node index, 1 <= s <= nv.
55*
56*  t is the sink node index, 1 <= t <= nv, t != s.
57*
58*  cap[a], a = 1,...,na, is the capacity of arc a, cap[a] >= 0.
59*
60*  NOTE: Multiple arcs are allowed, but self-loops are not allowed.
61*
62*  OUTPUT PARAMETERS
63*
64*  x[a], a = 1,...,na, is optimal value of the flow through arc a.
65*
66*  cut[i], i = 1,...,nv, is 1 if node i is labelled, and 0 otherwise.
67*  The set of arcs, whose one endpoint is labelled and other is not,
68*  defines the minimal cut corresponding to the maximal flow found.
69*  If the parameter cut is NULL, the cut information are not stored.
70*
71*  REFERENCES
72*
73*  L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND
74*  Corp., Report R-375-PR (August 1962), Chap. I "Static Maximal Flow,"
75*  pp.30-33. */
76
77void ffalg(int nv, int na, const int tail[], const int head[],
78      int s, int t, const int cap[], int x[], char cut[])
79{     int a, delta, i, j, k, pos1, pos2, temp,
80         *ptr, *arc, *link, *list;
81      /* sanity checks */
82      xassert(nv >= 2);
83      xassert(na >= 0);
84      xassert(1 <= s && s <= nv);
85      xassert(1 <= t && t <= nv);
86      xassert(s != t);
87      for (a = 1; a <= na; a++)
88      {  i = tail[a], j = head[a];
89         xassert(1 <= i && i <= nv);
90         xassert(1 <= j && j <= nv);
91         xassert(i != j);
92         xassert(cap[a] >= 0);
93      }
94      /* allocate working arrays */
95      ptr = xcalloc(1+nv+1, sizeof(int));
96      arc = xcalloc(1+na+na, sizeof(int));
97      link = xcalloc(1+nv, sizeof(int));
98      list = xcalloc(1+nv, sizeof(int));
99      /* ptr[i] := (degree of node i) */
100      for (i = 1; i <= nv; i++)
101         ptr[i] = 0;
102      for (a = 1; a <= na; a++)
103      {  ptr[tail[a]]++;
104         ptr[head[a]]++;
105      }
106      /* initialize arc pointers */
107      ptr[1]++;
108      for (i = 1; i < nv; i++)
109         ptr[i+1] += ptr[i];
110      ptr[nv+1] = ptr[nv];
111      /* build arc lists */
112      for (a = 1; a <= na; a++)
113      {  arc[--ptr[tail[a]]] = a;
114         arc[--ptr[head[a]]] = a;
115      }
116      xassert(ptr[1] == 1);
117      xassert(ptr[nv+1] == na+na+1);
118      /* now the indices of arcs incident to node i are stored in
119         locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */
120      /* initialize arc flows */
121      for (a = 1; a <= na; a++)
122         x[a] = 0;
123loop: /* main loop starts here */
124      /* build augmenting tree rooted at s */
125      /* link[i] = 0 means that node i is not labelled yet;
126         link[i] = a means that arc a immediately precedes node i */
127      /* initially node s is labelled as the root */
128      for (i = 1; i <= nv; i++)
129         link[i] = 0;
130      link[s] = -1, list[1] = s, pos1 = pos2 = 1;
131      /* breadth first search */
132      while (pos1 <= pos2)
133      {  /* dequeue node i */
134         i = list[pos1++];
135         /* consider all arcs incident to node i */
136         for (k = ptr[i]; k < ptr[i+1]; k++)
137         {  a = arc[k];
138            if (tail[a] == i)
139            {  /* a = i->j is a forward arc from s to t */
140               j = head[a];
141               /* if node j has been labelled, skip the arc */
142               if (link[j] != 0) continue;
143               /* if the arc does not allow increasing the flow through
144                  it, skip the arc */
145               if (x[a] == cap[a]) continue;
146            }
147            else if (head[a] == i)
148            {  /* a = i<-j is a backward arc from s to t */
149               j = tail[a];
150               /* if node j has been labelled, skip the arc */
151               if (link[j] != 0) continue;
152               /* if the arc does not allow decreasing the flow through
153                  it, skip the arc */
154               if (x[a] == 0) continue;
155            }
156            else
157               xassert(a != a);
158            /* label node j and enqueue it */
159            link[j] = a, list[++pos2] = j;
160            /* check for breakthrough */
161            if (j == t) goto brkt;
162         }
163      }
164      /* NONBREAKTHROUGH */
165      /* no augmenting path exists; current flow is maximal */
166      /* store minimal cut information, if necessary */
167      if (cut != NULL)
168      {  for (i = 1; i <= nv; i++)
169            cut[i] = (char)(link[i] != 0);
170      }
171      goto done;
172brkt: /* BREAKTHROUGH */
173      /* walk through arcs of the augmenting path (s, ..., t) found in
174         the reverse order and determine maximal change of the flow */
175      delta = 0;
176      for (j = t; j != s; j = i)
177      {  /* arc a immediately precedes node j in the path */
178         a = link[j];
179         if (head[a] == j)
180         {  /* a = i->j is a forward arc of the cycle */
181            i = tail[a];
182            /* x[a] may be increased until its upper bound */
183            temp = cap[a] - x[a];
184         }
185         else if (tail[a] == j)
186         {  /* a = i<-j is a backward arc of the cycle */
187            i = head[a];
188            /* x[a] may be decreased until its lower bound */
189            temp = x[a];
190         }
191         else
192            xassert(a != a);
193         if (delta == 0 || delta > temp) delta = temp;
194      }
195      xassert(delta > 0);
196      /* increase the flow along the path */
197      for (j = t; j != s; j = i)
198      {  /* arc a immediately precedes node j in the path */
199         a = link[j];
200         if (head[a] == j)
201         {  /* a = i->j is a forward arc of the cycle */
202            i = tail[a];
203            x[a] += delta;
204         }
205         else if (tail[a] == j)
206         {  /* a = i<-j is a backward arc of the cycle */
207            i = head[a];
208            x[a] -= delta;
209         }
210         else
211            xassert(a != a);
212      }
213      goto loop;
214done: /* free working arrays */
215      xfree(ptr);
216      xfree(arc);
217      xfree(link);
218      xfree(list);
219      return;
220}
221
222/* eof */
Note: See TracBrowser for help on using the repository browser.