|
1 /* glpnet06.c (out-of-kilter 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 * okalg - out-of-kilter algorithm |
|
32 * |
|
33 * SYNOPSIS |
|
34 * |
|
35 * #include "glpnet.h" |
|
36 * int okalg(int nv, int na, const int tail[], const int head[], |
|
37 * const int low[], const int cap[], const int cost[], int x[], |
|
38 * int pi[]); |
|
39 * |
|
40 * DESCRIPTION |
|
41 * |
|
42 * The routine okalg implements the out-of-kilter algorithm to find a |
|
43 * minimal-cost circulation in the specified flow network. |
|
44 * |
|
45 * INPUT PARAMETERS |
|
46 * |
|
47 * nv is the number of nodes, nv >= 0. |
|
48 * |
|
49 * na is the number of arcs, na >= 0. |
|
50 * |
|
51 * tail[a], a = 1,...,na, is the index of tail node of arc a. |
|
52 * |
|
53 * head[a], a = 1,...,na, is the index of head node of arc a. |
|
54 * |
|
55 * low[a], a = 1,...,na, is an lower bound to the flow through arc a. |
|
56 * |
|
57 * cap[a], a = 1,...,na, is an upper bound to the flow through arc a, |
|
58 * which is the capacity of the arc. |
|
59 * |
|
60 * cost[a], a = 1,...,na, is a per-unit cost of the flow through arc a. |
|
61 * |
|
62 * NOTES |
|
63 * |
|
64 * 1. Multiple arcs are allowed, but self-loops are not allowed. |
|
65 * |
|
66 * 2. It is required that 0 <= low[a] <= cap[a] for all arcs. |
|
67 * |
|
68 * 3. Arc costs may have any sign. |
|
69 * |
|
70 * OUTPUT PARAMETERS |
|
71 * |
|
72 * x[a], a = 1,...,na, is optimal value of the flow through arc a. |
|
73 * |
|
74 * pi[i], i = 1,...,nv, is Lagrange multiplier for flow conservation |
|
75 * equality constraint corresponding to node i (the node potential). |
|
76 * |
|
77 * RETURNS |
|
78 * |
|
79 * 0 optimal circulation found; |
|
80 * |
|
81 * 1 there is no feasible circulation; |
|
82 * |
|
83 * 2 integer overflow occured; |
|
84 * |
|
85 * 3 optimality test failed (logic error). |
|
86 * |
|
87 * REFERENCES |
|
88 * |
|
89 * L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND |
|
90 * Corp., Report R-375-PR (August 1962), Chap. III "Minimal Cost Flow |
|
91 * Problems," pp.113-26. */ |
|
92 |
|
93 static int overflow(int u, int v) |
|
94 { /* check for integer overflow on computing u + v */ |
|
95 if (u > 0 && v > 0 && u + v < 0) return 1; |
|
96 if (u < 0 && v < 0 && u + v > 0) return 1; |
|
97 return 0; |
|
98 } |
|
99 |
|
100 int okalg(int nv, int na, const int tail[], const int head[], |
|
101 const int low[], const int cap[], const int cost[], int x[], |
|
102 int pi[]) |
|
103 { int a, aok, delta, i, j, k, lambda, pos1, pos2, s, t, temp, ret, |
|
104 *ptr, *arc, *link, *list; |
|
105 /* sanity checks */ |
|
106 xassert(nv >= 0); |
|
107 xassert(na >= 0); |
|
108 for (a = 1; a <= na; a++) |
|
109 { i = tail[a], j = head[a]; |
|
110 xassert(1 <= i && i <= nv); |
|
111 xassert(1 <= j && j <= nv); |
|
112 xassert(i != j); |
|
113 xassert(0 <= low[a] && low[a] <= cap[a]); |
|
114 } |
|
115 /* allocate working arrays */ |
|
116 ptr = xcalloc(1+nv+1, sizeof(int)); |
|
117 arc = xcalloc(1+na+na, sizeof(int)); |
|
118 link = xcalloc(1+nv, sizeof(int)); |
|
119 list = xcalloc(1+nv, sizeof(int)); |
|
120 /* ptr[i] := (degree of node i) */ |
|
121 for (i = 1; i <= nv; i++) |
|
122 ptr[i] = 0; |
|
123 for (a = 1; a <= na; a++) |
|
124 { ptr[tail[a]]++; |
|
125 ptr[head[a]]++; |
|
126 } |
|
127 /* initialize arc pointers */ |
|
128 ptr[1]++; |
|
129 for (i = 1; i < nv; i++) |
|
130 ptr[i+1] += ptr[i]; |
|
131 ptr[nv+1] = ptr[nv]; |
|
132 /* build arc lists */ |
|
133 for (a = 1; a <= na; a++) |
|
134 { arc[--ptr[tail[a]]] = a; |
|
135 arc[--ptr[head[a]]] = a; |
|
136 } |
|
137 xassert(ptr[1] == 1); |
|
138 xassert(ptr[nv+1] == na+na+1); |
|
139 /* now the indices of arcs incident to node i are stored in |
|
140 locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */ |
|
141 /* initialize arc flows and node potentials */ |
|
142 for (a = 1; a <= na; a++) |
|
143 x[a] = 0; |
|
144 for (i = 1; i <= nv; i++) |
|
145 pi[i] = 0; |
|
146 loop: /* main loop starts here */ |
|
147 /* find out-of-kilter arc */ |
|
148 aok = 0; |
|
149 for (a = 1; a <= na; a++) |
|
150 { i = tail[a], j = head[a]; |
|
151 if (overflow(cost[a], pi[i] - pi[j])) |
|
152 { ret = 2; |
|
153 goto done; |
|
154 } |
|
155 lambda = cost[a] + (pi[i] - pi[j]); |
|
156 if (x[a] < low[a] || lambda < 0 && x[a] < cap[a]) |
|
157 { /* arc a = i->j is out of kilter, and we need to increase |
|
158 the flow through this arc */ |
|
159 aok = a, s = j, t = i; |
|
160 break; |
|
161 } |
|
162 if (x[a] > cap[a] || lambda > 0 && x[a] > low[a]) |
|
163 { /* arc a = i->j is out of kilter, and we need to decrease |
|
164 the flow through this arc */ |
|
165 aok = a, s = i, t = j; |
|
166 break; |
|
167 } |
|
168 } |
|
169 if (aok == 0) |
|
170 { /* all arcs are in kilter */ |
|
171 /* check for feasibility */ |
|
172 for (a = 1; a <= na; a++) |
|
173 { if (!(low[a] <= x[a] && x[a] <= cap[a])) |
|
174 { ret = 3; |
|
175 goto done; |
|
176 } |
|
177 } |
|
178 for (i = 1; i <= nv; i++) |
|
179 { temp = 0; |
|
180 for (k = ptr[i]; k < ptr[i+1]; k++) |
|
181 { a = arc[k]; |
|
182 if (tail[a] == i) |
|
183 { /* a is outgoing arc */ |
|
184 temp += x[a]; |
|
185 } |
|
186 else if (head[a] == i) |
|
187 { /* a is incoming arc */ |
|
188 temp -= x[a]; |
|
189 } |
|
190 else |
|
191 xassert(a != a); |
|
192 } |
|
193 if (temp != 0) |
|
194 { ret = 3; |
|
195 goto done; |
|
196 } |
|
197 } |
|
198 /* check for optimality */ |
|
199 for (a = 1; a <= na; a++) |
|
200 { i = tail[a], j = head[a]; |
|
201 lambda = cost[a] + (pi[i] - pi[j]); |
|
202 if (lambda > 0 && x[a] != low[a] || |
|
203 lambda < 0 && x[a] != cap[a]) |
|
204 { ret = 3; |
|
205 goto done; |
|
206 } |
|
207 } |
|
208 /* current circulation is optimal */ |
|
209 ret = 0; |
|
210 goto done; |
|
211 } |
|
212 /* now we need to find a cycle (t, a, s, ..., t), which allows |
|
213 increasing the flow along it, where a is the out-of-kilter arc |
|
214 just found */ |
|
215 /* link[i] = 0 means that node i is not labelled yet; |
|
216 link[i] = a means that arc a immediately precedes node i */ |
|
217 /* initially only node s is labelled */ |
|
218 for (i = 1; i <= nv; i++) |
|
219 link[i] = 0; |
|
220 link[s] = aok, list[1] = s, pos1 = pos2 = 1; |
|
221 /* breadth first search */ |
|
222 while (pos1 <= pos2) |
|
223 { /* dequeue node i */ |
|
224 i = list[pos1++]; |
|
225 /* consider all arcs incident to node i */ |
|
226 for (k = ptr[i]; k < ptr[i+1]; k++) |
|
227 { a = arc[k]; |
|
228 if (tail[a] == i) |
|
229 { /* a = i->j is a forward arc from s to t */ |
|
230 j = head[a]; |
|
231 /* if node j has been labelled, skip the arc */ |
|
232 if (link[j] != 0) continue; |
|
233 /* if the arc does not allow increasing the flow through |
|
234 it, skip the arc */ |
|
235 if (x[a] >= cap[a]) continue; |
|
236 if (overflow(cost[a], pi[i] - pi[j])) |
|
237 { ret = 2; |
|
238 goto done; |
|
239 } |
|
240 lambda = cost[a] + (pi[i] - pi[j]); |
|
241 if (lambda > 0 && x[a] >= low[a]) continue; |
|
242 } |
|
243 else if (head[a] == i) |
|
244 { /* a = i<-j is a backward arc from s to t */ |
|
245 j = tail[a]; |
|
246 /* if node j has been labelled, skip the arc */ |
|
247 if (link[j] != 0) continue; |
|
248 /* if the arc does not allow decreasing the flow through |
|
249 it, skip the arc */ |
|
250 if (x[a] <= low[a]) continue; |
|
251 if (overflow(cost[a], pi[j] - pi[i])) |
|
252 { ret = 2; |
|
253 goto done; |
|
254 } |
|
255 lambda = cost[a] + (pi[j] - pi[i]); |
|
256 if (lambda < 0 && x[a] <= cap[a]) continue; |
|
257 } |
|
258 else |
|
259 xassert(a != a); |
|
260 /* label node j and enqueue it */ |
|
261 link[j] = a, list[++pos2] = j; |
|
262 /* check for breakthrough */ |
|
263 if (j == t) goto brkt; |
|
264 } |
|
265 } |
|
266 /* NONBREAKTHROUGH */ |
|
267 /* consider all arcs, whose one endpoint is labelled and other is |
|
268 not, and determine maximal change of node potentials */ |
|
269 delta = 0; |
|
270 for (a = 1; a <= na; a++) |
|
271 { i = tail[a], j = head[a]; |
|
272 if (link[i] != 0 && link[j] == 0) |
|
273 { /* a = i->j, where node i is labelled, node j is not */ |
|
274 if (overflow(cost[a], pi[i] - pi[j])) |
|
275 { ret = 2; |
|
276 goto done; |
|
277 } |
|
278 lambda = cost[a] + (pi[i] - pi[j]); |
|
279 if (x[a] <= cap[a] && lambda > 0) |
|
280 if (delta == 0 || delta > + lambda) delta = + lambda; |
|
281 } |
|
282 else if (link[i] == 0 && link[j] != 0) |
|
283 { /* a = j<-i, where node j is labelled, node i is not */ |
|
284 if (overflow(cost[a], pi[i] - pi[j])) |
|
285 { ret = 2; |
|
286 goto done; |
|
287 } |
|
288 lambda = cost[a] + (pi[i] - pi[j]); |
|
289 if (x[a] >= low[a] && lambda < 0) |
|
290 if (delta == 0 || delta > - lambda) delta = - lambda; |
|
291 } |
|
292 } |
|
293 if (delta == 0) |
|
294 { /* there is no feasible circulation */ |
|
295 ret = 1; |
|
296 goto done; |
|
297 } |
|
298 /* increase potentials of all unlabelled nodes */ |
|
299 for (i = 1; i <= nv; i++) |
|
300 { if (link[i] == 0) |
|
301 { if (overflow(pi[i], delta)) |
|
302 { ret = 2; |
|
303 goto done; |
|
304 } |
|
305 pi[i] += delta; |
|
306 } |
|
307 } |
|
308 goto loop; |
|
309 brkt: /* BREAKTHROUGH */ |
|
310 /* walk through arcs of the cycle (t, a, s, ..., t) found in the |
|
311 reverse order and determine maximal change of the flow */ |
|
312 delta = 0; |
|
313 for (j = t;; j = i) |
|
314 { /* arc a immediately precedes node j in the cycle */ |
|
315 a = link[j]; |
|
316 if (head[a] == j) |
|
317 { /* a = i->j is a forward arc of the cycle */ |
|
318 i = tail[a]; |
|
319 lambda = cost[a] + (pi[i] - pi[j]); |
|
320 if (lambda > 0 && x[a] < low[a]) |
|
321 { /* x[a] may be increased until its lower bound */ |
|
322 temp = low[a] - x[a]; |
|
323 } |
|
324 else if (lambda <= 0 && x[a] < cap[a]) |
|
325 { /* x[a] may be increased until its upper bound */ |
|
326 temp = cap[a] - x[a]; |
|
327 } |
|
328 else |
|
329 xassert(a != a); |
|
330 } |
|
331 else if (tail[a] == j) |
|
332 { /* a = i<-j is a backward arc of the cycle */ |
|
333 i = head[a]; |
|
334 lambda = cost[a] + (pi[j] - pi[i]); |
|
335 if (lambda < 0 && x[a] > cap[a]) |
|
336 { /* x[a] may be decreased until its upper bound */ |
|
337 temp = x[a] - cap[a]; |
|
338 } |
|
339 else if (lambda >= 0 && x[a] > low[a]) |
|
340 { /* x[a] may be decreased until its lower bound */ |
|
341 temp = x[a] - low[a]; |
|
342 } |
|
343 else |
|
344 xassert(a != a); |
|
345 } |
|
346 else |
|
347 xassert(a != a); |
|
348 if (delta == 0 || delta > temp) delta = temp; |
|
349 /* check for end of the cycle */ |
|
350 if (i == t) break; |
|
351 } |
|
352 xassert(delta > 0); |
|
353 /* increase the flow along the cycle */ |
|
354 for (j = t;; j = i) |
|
355 { /* arc a immediately precedes node j in the cycle */ |
|
356 a = link[j]; |
|
357 if (head[a] == j) |
|
358 { /* a = i->j is a forward arc of the cycle */ |
|
359 i = tail[a]; |
|
360 /* overflow cannot occur */ |
|
361 x[a] += delta; |
|
362 } |
|
363 else if (tail[a] == j) |
|
364 { /* a = i<-j is a backward arc of the cycle */ |
|
365 i = head[a]; |
|
366 /* overflow cannot occur */ |
|
367 x[a] -= delta; |
|
368 } |
|
369 else |
|
370 xassert(a != a); |
|
371 /* check for end of the cycle */ |
|
372 if (i == t) break; |
|
373 } |
|
374 goto loop; |
|
375 done: /* free working arrays */ |
|
376 xfree(ptr); |
|
377 xfree(arc); |
|
378 xfree(link); |
|
379 xfree(list); |
|
380 return ret; |
|
381 } |
|
382 |
|
383 /* eof */ |