COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpnet06.c @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 12.4 KB
Line 
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
93static 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
100int 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;
146loop: /* 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;
309brkt: /* 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;
375done: /* free working arrays */
376      xfree(ptr);
377      xfree(arc);
378      xfree(link);
379      xfree(list);
380      return ret;
381}
382
383/* eof */
Note: See TracBrowser for help on using the repository browser.