1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnet06.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,383 @@
     1.4 +/* glpnet06.c (out-of-kilter algorithm) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpenv.h"
    1.29 +#include "glpnet.h"
    1.30 +
    1.31 +/***********************************************************************
    1.32 +*  NAME
    1.33 +*
    1.34 +*  okalg - out-of-kilter algorithm
    1.35 +*
    1.36 +*  SYNOPSIS
    1.37 +*
    1.38 +*  #include "glpnet.h"
    1.39 +*  int okalg(int nv, int na, const int tail[], const int head[],
    1.40 +*     const int low[], const int cap[], const int cost[], int x[],
    1.41 +*     int pi[]);
    1.42 +*
    1.43 +*  DESCRIPTION
    1.44 +*
    1.45 +*  The routine okalg implements the out-of-kilter algorithm to find a
    1.46 +*  minimal-cost circulation in the specified flow network.
    1.47 +*
    1.48 +*  INPUT PARAMETERS
    1.49 +*
    1.50 +*  nv is the number of nodes, nv >= 0.
    1.51 +*
    1.52 +*  na is the number of arcs, na >= 0.
    1.53 +*
    1.54 +*  tail[a], a = 1,...,na, is the index of tail node of arc a.
    1.55 +*
    1.56 +*  head[a], a = 1,...,na, is the index of head node of arc a.
    1.57 +*
    1.58 +*  low[a], a = 1,...,na, is an lower bound to the flow through arc a.
    1.59 +*
    1.60 +*  cap[a], a = 1,...,na, is an upper bound to the flow through arc a,
    1.61 +*  which is the capacity of the arc.
    1.62 +*
    1.63 +*  cost[a], a = 1,...,na, is a per-unit cost of the flow through arc a.
    1.64 +*
    1.65 +*  NOTES
    1.66 +*
    1.67 +*  1. Multiple arcs are allowed, but self-loops are not allowed.
    1.68 +*
    1.69 +*  2. It is required that 0 <= low[a] <= cap[a] for all arcs.
    1.70 +*
    1.71 +*  3. Arc costs may have any sign.
    1.72 +*
    1.73 +*  OUTPUT PARAMETERS
    1.74 +*
    1.75 +*  x[a], a = 1,...,na, is optimal value of the flow through arc a.
    1.76 +*
    1.77 +*  pi[i], i = 1,...,nv, is Lagrange multiplier for flow conservation
    1.78 +*  equality constraint corresponding to node i (the node potential).
    1.79 +*
    1.80 +*  RETURNS
    1.81 +*
    1.82 +*  0  optimal circulation found;
    1.83 +*
    1.84 +*  1  there is no feasible circulation;
    1.85 +*
    1.86 +*  2  integer overflow occured;
    1.87 +*
    1.88 +*  3  optimality test failed (logic error).
    1.89 +*
    1.90 +*  REFERENCES
    1.91 +*
    1.92 +*  L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND
    1.93 +*  Corp., Report R-375-PR (August 1962), Chap. III "Minimal Cost Flow
    1.94 +*  Problems," pp.113-26. */
    1.95 +
    1.96 +static int overflow(int u, int v)
    1.97 +{     /* check for integer overflow on computing u + v */
    1.98 +      if (u > 0 && v > 0 && u + v < 0) return 1;
    1.99 +      if (u < 0 && v < 0 && u + v > 0) return 1;
   1.100 +      return 0;
   1.101 +}
   1.102 +
   1.103 +int okalg(int nv, int na, const int tail[], const int head[],
   1.104 +      const int low[], const int cap[], const int cost[], int x[],
   1.105 +      int pi[])
   1.106 +{     int a, aok, delta, i, j, k, lambda, pos1, pos2, s, t, temp, ret,
   1.107 +         *ptr, *arc, *link, *list;
   1.108 +      /* sanity checks */
   1.109 +      xassert(nv >= 0);
   1.110 +      xassert(na >= 0);
   1.111 +      for (a = 1; a <= na; a++)
   1.112 +      {  i = tail[a], j = head[a];
   1.113 +         xassert(1 <= i && i <= nv);
   1.114 +         xassert(1 <= j && j <= nv);
   1.115 +         xassert(i != j);
   1.116 +         xassert(0 <= low[a] && low[a] <= cap[a]);
   1.117 +      }
   1.118 +      /* allocate working arrays */
   1.119 +      ptr = xcalloc(1+nv+1, sizeof(int));
   1.120 +      arc = xcalloc(1+na+na, sizeof(int));
   1.121 +      link = xcalloc(1+nv, sizeof(int));
   1.122 +      list = xcalloc(1+nv, sizeof(int));
   1.123 +      /* ptr[i] := (degree of node i) */
   1.124 +      for (i = 1; i <= nv; i++)
   1.125 +         ptr[i] = 0;
   1.126 +      for (a = 1; a <= na; a++)
   1.127 +      {  ptr[tail[a]]++;
   1.128 +         ptr[head[a]]++;
   1.129 +      }
   1.130 +      /* initialize arc pointers */
   1.131 +      ptr[1]++;
   1.132 +      for (i = 1; i < nv; i++)
   1.133 +         ptr[i+1] += ptr[i];
   1.134 +      ptr[nv+1] = ptr[nv];
   1.135 +      /* build arc lists */
   1.136 +      for (a = 1; a <= na; a++)
   1.137 +      {  arc[--ptr[tail[a]]] = a;
   1.138 +         arc[--ptr[head[a]]] = a;
   1.139 +      }
   1.140 +      xassert(ptr[1] == 1);
   1.141 +      xassert(ptr[nv+1] == na+na+1);
   1.142 +      /* now the indices of arcs incident to node i are stored in
   1.143 +         locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */
   1.144 +      /* initialize arc flows and node potentials */
   1.145 +      for (a = 1; a <= na; a++)
   1.146 +         x[a] = 0;
   1.147 +      for (i = 1; i <= nv; i++)
   1.148 +         pi[i] = 0;
   1.149 +loop: /* main loop starts here */
   1.150 +      /* find out-of-kilter arc */
   1.151 +      aok = 0;
   1.152 +      for (a = 1; a <= na; a++)
   1.153 +      {  i = tail[a], j = head[a];
   1.154 +         if (overflow(cost[a], pi[i] - pi[j]))
   1.155 +         {  ret = 2;
   1.156 +            goto done;
   1.157 +         }
   1.158 +         lambda = cost[a] + (pi[i] - pi[j]);
   1.159 +         if (x[a] < low[a] || lambda < 0 && x[a] < cap[a])
   1.160 +         {  /* arc a = i->j is out of kilter, and we need to increase
   1.161 +               the flow through this arc */
   1.162 +            aok = a, s = j, t = i;
   1.163 +            break;
   1.164 +         }
   1.165 +         if (x[a] > cap[a] || lambda > 0 && x[a] > low[a])
   1.166 +         {  /* arc a = i->j is out of kilter, and we need to decrease
   1.167 +               the flow through this arc */
   1.168 +            aok = a, s = i, t = j;
   1.169 +            break;
   1.170 +         }
   1.171 +      }
   1.172 +      if (aok == 0)
   1.173 +      {  /* all arcs are in kilter */
   1.174 +         /* check for feasibility */
   1.175 +         for (a = 1; a <= na; a++)
   1.176 +         {  if (!(low[a] <= x[a] && x[a] <= cap[a]))
   1.177 +            {  ret = 3;
   1.178 +               goto done;
   1.179 +            }
   1.180 +         }
   1.181 +         for (i = 1; i <= nv; i++)
   1.182 +         {  temp = 0;
   1.183 +            for (k = ptr[i]; k < ptr[i+1]; k++)
   1.184 +            {  a = arc[k];
   1.185 +               if (tail[a] == i)
   1.186 +               {  /* a is outgoing arc */
   1.187 +                  temp += x[a];
   1.188 +               }
   1.189 +               else if (head[a] == i)
   1.190 +               {  /* a is incoming arc */
   1.191 +                  temp -= x[a];
   1.192 +               }
   1.193 +               else
   1.194 +                  xassert(a != a);
   1.195 +            }
   1.196 +            if (temp != 0)
   1.197 +            {  ret = 3;
   1.198 +               goto done;
   1.199 +            }
   1.200 +         }
   1.201 +         /* check for optimality */
   1.202 +         for (a = 1; a <= na; a++)
   1.203 +         {  i = tail[a], j = head[a];
   1.204 +            lambda = cost[a] + (pi[i] - pi[j]);
   1.205 +            if (lambda > 0 && x[a] != low[a] ||
   1.206 +                lambda < 0 && x[a] != cap[a])
   1.207 +            {  ret = 3;
   1.208 +               goto done;
   1.209 +            }
   1.210 +         }
   1.211 +         /* current circulation is optimal */
   1.212 +         ret = 0;
   1.213 +         goto done;
   1.214 +      }
   1.215 +      /* now we need to find a cycle (t, a, s, ..., t), which allows
   1.216 +         increasing the flow along it, where a is the out-of-kilter arc
   1.217 +         just found */
   1.218 +      /* link[i] = 0 means that node i is not labelled yet;
   1.219 +         link[i] = a means that arc a immediately precedes node i */
   1.220 +      /* initially only node s is labelled */
   1.221 +      for (i = 1; i <= nv; i++)
   1.222 +         link[i] = 0;
   1.223 +      link[s] = aok, list[1] = s, pos1 = pos2 = 1;
   1.224 +      /* breadth first search */
   1.225 +      while (pos1 <= pos2)
   1.226 +      {  /* dequeue node i */
   1.227 +         i = list[pos1++];
   1.228 +         /* consider all arcs incident to node i */
   1.229 +         for (k = ptr[i]; k < ptr[i+1]; k++)
   1.230 +         {  a = arc[k];
   1.231 +            if (tail[a] == i)
   1.232 +            {  /* a = i->j is a forward arc from s to t */
   1.233 +               j = head[a];
   1.234 +               /* if node j has been labelled, skip the arc */
   1.235 +               if (link[j] != 0) continue;
   1.236 +               /* if the arc does not allow increasing the flow through
   1.237 +                  it, skip the arc */
   1.238 +               if (x[a] >= cap[a]) continue;
   1.239 +               if (overflow(cost[a], pi[i] - pi[j]))
   1.240 +               {  ret = 2;
   1.241 +                  goto done;
   1.242 +               }
   1.243 +               lambda = cost[a] + (pi[i] - pi[j]);
   1.244 +               if (lambda > 0 && x[a] >= low[a]) continue;
   1.245 +            }
   1.246 +            else if (head[a] == i)
   1.247 +            {  /* a = i<-j is a backward arc from s to t */
   1.248 +               j = tail[a];
   1.249 +               /* if node j has been labelled, skip the arc */
   1.250 +               if (link[j] != 0) continue;
   1.251 +               /* if the arc does not allow decreasing the flow through
   1.252 +                  it, skip the arc */
   1.253 +               if (x[a] <= low[a]) continue;
   1.254 +               if (overflow(cost[a], pi[j] - pi[i]))
   1.255 +               {  ret = 2;
   1.256 +                  goto done;
   1.257 +               }
   1.258 +               lambda = cost[a] + (pi[j] - pi[i]);
   1.259 +               if (lambda < 0 && x[a] <= cap[a]) continue;
   1.260 +            }
   1.261 +            else
   1.262 +               xassert(a != a);
   1.263 +            /* label node j and enqueue it */
   1.264 +            link[j] = a, list[++pos2] = j;
   1.265 +            /* check for breakthrough */
   1.266 +            if (j == t) goto brkt;
   1.267 +         }
   1.268 +      }
   1.269 +      /* NONBREAKTHROUGH */
   1.270 +      /* consider all arcs, whose one endpoint is labelled and other is
   1.271 +         not, and determine maximal change of node potentials */
   1.272 +      delta = 0;
   1.273 +      for (a = 1; a <= na; a++)
   1.274 +      {  i = tail[a], j = head[a];
   1.275 +         if (link[i] != 0 && link[j] == 0)
   1.276 +         {  /* a = i->j, where node i is labelled, node j is not */
   1.277 +            if (overflow(cost[a], pi[i] - pi[j]))
   1.278 +            {  ret = 2;
   1.279 +               goto done;
   1.280 +            }
   1.281 +            lambda = cost[a] + (pi[i] - pi[j]);
   1.282 +            if (x[a] <= cap[a] && lambda > 0)
   1.283 +               if (delta == 0 || delta > + lambda) delta = + lambda;
   1.284 +         }
   1.285 +         else if (link[i] == 0 && link[j] != 0)
   1.286 +         {  /* a = j<-i, where node j is labelled, node i is not */
   1.287 +            if (overflow(cost[a], pi[i] - pi[j]))
   1.288 +            {  ret = 2;
   1.289 +               goto done;
   1.290 +            }
   1.291 +            lambda = cost[a] + (pi[i] - pi[j]);
   1.292 +            if (x[a] >= low[a] && lambda < 0)
   1.293 +               if (delta == 0 || delta > - lambda) delta = - lambda;
   1.294 +         }
   1.295 +      }
   1.296 +      if (delta == 0)
   1.297 +      {  /* there is no feasible circulation */
   1.298 +         ret = 1;
   1.299 +         goto done;
   1.300 +      }
   1.301 +      /* increase potentials of all unlabelled nodes */
   1.302 +      for (i = 1; i <= nv; i++)
   1.303 +      {  if (link[i] == 0)
   1.304 +         {  if (overflow(pi[i], delta))
   1.305 +            {  ret = 2;
   1.306 +               goto done;
   1.307 +            }
   1.308 +            pi[i] += delta;
   1.309 +         }
   1.310 +      }
   1.311 +      goto loop;
   1.312 +brkt: /* BREAKTHROUGH */
   1.313 +      /* walk through arcs of the cycle (t, a, s, ..., t) found in the
   1.314 +         reverse order and determine maximal change of the flow */
   1.315 +      delta = 0;
   1.316 +      for (j = t;; j = i)
   1.317 +      {  /* arc a immediately precedes node j in the cycle */
   1.318 +         a = link[j];
   1.319 +         if (head[a] == j)
   1.320 +         {  /* a = i->j is a forward arc of the cycle */
   1.321 +            i = tail[a];
   1.322 +            lambda = cost[a] + (pi[i] - pi[j]);
   1.323 +            if (lambda > 0 && x[a] < low[a])
   1.324 +            {  /* x[a] may be increased until its lower bound */
   1.325 +               temp = low[a] - x[a];
   1.326 +            }
   1.327 +            else if (lambda <= 0 && x[a] < cap[a])
   1.328 +            {  /* x[a] may be increased until its upper bound */
   1.329 +               temp = cap[a] - x[a];
   1.330 +            }
   1.331 +            else
   1.332 +               xassert(a != a);
   1.333 +         }
   1.334 +         else if (tail[a] == j)
   1.335 +         {  /* a = i<-j is a backward arc of the cycle */
   1.336 +            i = head[a];
   1.337 +            lambda = cost[a] + (pi[j] - pi[i]);
   1.338 +            if (lambda < 0 && x[a] > cap[a])
   1.339 +            {  /* x[a] may be decreased until its upper bound */
   1.340 +               temp = x[a] - cap[a];
   1.341 +            }
   1.342 +            else if (lambda >= 0 && x[a] > low[a])
   1.343 +            {  /* x[a] may be decreased until its lower bound */
   1.344 +               temp = x[a] - low[a];
   1.345 +            }
   1.346 +            else
   1.347 +               xassert(a != a);
   1.348 +         }
   1.349 +         else
   1.350 +            xassert(a != a);
   1.351 +         if (delta == 0 || delta > temp) delta = temp;
   1.352 +         /* check for end of the cycle */
   1.353 +         if (i == t) break;
   1.354 +      }
   1.355 +      xassert(delta > 0);
   1.356 +      /* increase the flow along the cycle */
   1.357 +      for (j = t;; j = i)
   1.358 +      {  /* arc a immediately precedes node j in the cycle */
   1.359 +         a = link[j];
   1.360 +         if (head[a] == j)
   1.361 +         {  /* a = i->j is a forward arc of the cycle */
   1.362 +            i = tail[a];
   1.363 +            /* overflow cannot occur */
   1.364 +            x[a] += delta;
   1.365 +         }
   1.366 +         else if (tail[a] == j)
   1.367 +         {  /* a = i<-j is a backward arc of the cycle */
   1.368 +            i = head[a];
   1.369 +            /* overflow cannot occur */
   1.370 +            x[a] -= delta;
   1.371 +         }
   1.372 +         else
   1.373 +            xassert(a != a);
   1.374 +         /* check for end of the cycle */
   1.375 +         if (i == t) break;
   1.376 +      }
   1.377 +      goto loop;
   1.378 +done: /* free working arrays */
   1.379 +      xfree(ptr);
   1.380 +      xfree(arc);
   1.381 +      xfree(link);
   1.382 +      xfree(list);
   1.383 +      return ret;
   1.384 +}
   1.385 +
   1.386 +/* eof */