lemon-project-template-glpk: 33de93886c88 deps/glpk/src/glpnet06.c

lemon-project-template-glpk

view deps/glpk/src/glpnet06.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

line source

1 /* glpnet06.c (out-of-kilter algorithm) */

3 /***********************************************************************

4 * This code is part of GLPK (GNU Linear Programming Kit).

5 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

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 ***********************************************************************/

25 #include "glpenv.h"

26 #include "glpnet.h"

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. */

93 static 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 }

100 int 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;

146 loop: /* 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;

309 brkt: /* 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;

375 done: /* free working arrays */

376 xfree(ptr);

377 xfree(arc);

378 xfree(link);

379 xfree(list);

380 return ret;

381 }

382

383 /* eof */