|
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 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 |
|
77 void 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; |
|
123 loop: /* 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; |
|
172 brkt: /* 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; |
|
214 done: /* free working arrays */ |
|
215 xfree(ptr); |
|
216 xfree(arc); |
|
217 xfree(link); |
|
218 xfree(list); |
|
219 return; |
|
220 } |
|
221 |
|
222 /* eof */ |