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
overflow(int u,int v)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
okalg(int nv,int na,const int tail[],const int head[],const int low[],const int cap[],const int cost[],int x[],int pi[])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 */
384