lemon-project-template-glpk
comparison deps/glpk/src/glpnet06.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 21:43:29 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:66fafe38455a |
---|---|
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, 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 * 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 */ |