1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnet06.c Mon Dec 06 13:09:21 2010 +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 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 */