PageRenderTime 47ms CodeModel.GetById 13ms app.highlight 30ms RepoModel.GetById 1ms app.codeStats 0ms

/farmR/src/glpnet06.c

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