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