lemon-project-template-glpk

annotate deps/glpk/src/glpnet06.c @ 11:4fc6ad2fb8a6

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