|
1 /* glpapi17.c (flow network problems) */ |
|
2 |
|
3 /*********************************************************************** |
|
4 * This code is part of GLPK (GNU Linear Programming Kit). |
|
5 * |
|
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
|
7 * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
|
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
|
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 ***********************************************************************/ |
|
24 |
|
25 #include "glpapi.h" |
|
26 #include "glpnet.h" |
|
27 |
|
28 /*********************************************************************** |
|
29 * NAME |
|
30 * |
|
31 * glp_mincost_lp - convert minimum cost flow problem to LP |
|
32 * |
|
33 * SYNOPSIS |
|
34 * |
|
35 * void glp_mincost_lp(glp_prob *lp, glp_graph *G, int names, |
|
36 * int v_rhs, int a_low, int a_cap, int a_cost); |
|
37 * |
|
38 * DESCRIPTION |
|
39 * |
|
40 * The routine glp_mincost_lp builds an LP problem, which corresponds |
|
41 * to the minimum cost flow problem on the specified network G. */ |
|
42 |
|
43 void glp_mincost_lp(glp_prob *lp, glp_graph *G, int names, int v_rhs, |
|
44 int a_low, int a_cap, int a_cost) |
|
45 { glp_vertex *v; |
|
46 glp_arc *a; |
|
47 int i, j, type, ind[1+2]; |
|
48 double rhs, low, cap, cost, val[1+2]; |
|
49 if (!(names == GLP_ON || names == GLP_OFF)) |
|
50 xerror("glp_mincost_lp: names = %d; invalid parameter\n", |
|
51 names); |
|
52 if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double)) |
|
53 xerror("glp_mincost_lp: v_rhs = %d; invalid offset\n", v_rhs); |
|
54 if (a_low >= 0 && a_low > G->a_size - (int)sizeof(double)) |
|
55 xerror("glp_mincost_lp: a_low = %d; invalid offset\n", a_low); |
|
56 if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double)) |
|
57 xerror("glp_mincost_lp: a_cap = %d; invalid offset\n", a_cap); |
|
58 if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double)) |
|
59 xerror("glp_mincost_lp: a_cost = %d; invalid offset\n", a_cost) |
|
60 ; |
|
61 glp_erase_prob(lp); |
|
62 if (names) glp_set_prob_name(lp, G->name); |
|
63 if (G->nv > 0) glp_add_rows(lp, G->nv); |
|
64 for (i = 1; i <= G->nv; i++) |
|
65 { v = G->v[i]; |
|
66 if (names) glp_set_row_name(lp, i, v->name); |
|
67 if (v_rhs >= 0) |
|
68 memcpy(&rhs, (char *)v->data + v_rhs, sizeof(double)); |
|
69 else |
|
70 rhs = 0.0; |
|
71 glp_set_row_bnds(lp, i, GLP_FX, rhs, rhs); |
|
72 } |
|
73 if (G->na > 0) glp_add_cols(lp, G->na); |
|
74 for (i = 1, j = 0; i <= G->nv; i++) |
|
75 { v = G->v[i]; |
|
76 for (a = v->out; a != NULL; a = a->t_next) |
|
77 { j++; |
|
78 if (names) |
|
79 { char name[50+1]; |
|
80 sprintf(name, "x[%d,%d]", a->tail->i, a->head->i); |
|
81 xassert(strlen(name) < sizeof(name)); |
|
82 glp_set_col_name(lp, j, name); |
|
83 } |
|
84 if (a->tail->i != a->head->i) |
|
85 { ind[1] = a->tail->i, val[1] = +1.0; |
|
86 ind[2] = a->head->i, val[2] = -1.0; |
|
87 glp_set_mat_col(lp, j, 2, ind, val); |
|
88 } |
|
89 if (a_low >= 0) |
|
90 memcpy(&low, (char *)a->data + a_low, sizeof(double)); |
|
91 else |
|
92 low = 0.0; |
|
93 if (a_cap >= 0) |
|
94 memcpy(&cap, (char *)a->data + a_cap, sizeof(double)); |
|
95 else |
|
96 cap = 1.0; |
|
97 if (cap == DBL_MAX) |
|
98 type = GLP_LO; |
|
99 else if (low != cap) |
|
100 type = GLP_DB; |
|
101 else |
|
102 type = GLP_FX; |
|
103 glp_set_col_bnds(lp, j, type, low, cap); |
|
104 if (a_cost >= 0) |
|
105 memcpy(&cost, (char *)a->data + a_cost, sizeof(double)); |
|
106 else |
|
107 cost = 0.0; |
|
108 glp_set_obj_coef(lp, j, cost); |
|
109 } |
|
110 } |
|
111 xassert(j == G->na); |
|
112 return; |
|
113 } |
|
114 |
|
115 /**********************************************************************/ |
|
116 |
|
117 int glp_mincost_okalg(glp_graph *G, int v_rhs, int a_low, int a_cap, |
|
118 int a_cost, double *sol, int a_x, int v_pi) |
|
119 { /* find minimum-cost flow with out-of-kilter algorithm */ |
|
120 glp_vertex *v; |
|
121 glp_arc *a; |
|
122 int nv, na, i, k, s, t, *tail, *head, *low, *cap, *cost, *x, *pi, |
|
123 ret; |
|
124 double sum, temp; |
|
125 if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double)) |
|
126 xerror("glp_mincost_okalg: v_rhs = %d; invalid offset\n", |
|
127 v_rhs); |
|
128 if (a_low >= 0 && a_low > G->a_size - (int)sizeof(double)) |
|
129 xerror("glp_mincost_okalg: a_low = %d; invalid offset\n", |
|
130 a_low); |
|
131 if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double)) |
|
132 xerror("glp_mincost_okalg: a_cap = %d; invalid offset\n", |
|
133 a_cap); |
|
134 if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double)) |
|
135 xerror("glp_mincost_okalg: a_cost = %d; invalid offset\n", |
|
136 a_cost); |
|
137 if (a_x >= 0 && a_x > G->a_size - (int)sizeof(double)) |
|
138 xerror("glp_mincost_okalg: a_x = %d; invalid offset\n", a_x); |
|
139 if (v_pi >= 0 && v_pi > G->v_size - (int)sizeof(double)) |
|
140 xerror("glp_mincost_okalg: v_pi = %d; invalid offset\n", v_pi); |
|
141 /* s is artificial source node */ |
|
142 s = G->nv + 1; |
|
143 /* t is artificial sink node */ |
|
144 t = s + 1; |
|
145 /* nv is the total number of nodes in the resulting network */ |
|
146 nv = t; |
|
147 /* na is the total number of arcs in the resulting network */ |
|
148 na = G->na + 1; |
|
149 for (i = 1; i <= G->nv; i++) |
|
150 { v = G->v[i]; |
|
151 if (v_rhs >= 0) |
|
152 memcpy(&temp, (char *)v->data + v_rhs, sizeof(double)); |
|
153 else |
|
154 temp = 0.0; |
|
155 if (temp != 0.0) na++; |
|
156 } |
|
157 /* allocate working arrays */ |
|
158 tail = xcalloc(1+na, sizeof(int)); |
|
159 head = xcalloc(1+na, sizeof(int)); |
|
160 low = xcalloc(1+na, sizeof(int)); |
|
161 cap = xcalloc(1+na, sizeof(int)); |
|
162 cost = xcalloc(1+na, sizeof(int)); |
|
163 x = xcalloc(1+na, sizeof(int)); |
|
164 pi = xcalloc(1+nv, sizeof(int)); |
|
165 /* construct the resulting network */ |
|
166 k = 0; |
|
167 /* (original arcs) */ |
|
168 for (i = 1; i <= G->nv; i++) |
|
169 { v = G->v[i]; |
|
170 for (a = v->out; a != NULL; a = a->t_next) |
|
171 { k++; |
|
172 tail[k] = a->tail->i; |
|
173 head[k] = a->head->i; |
|
174 if (tail[k] == head[k]) |
|
175 { ret = GLP_EDATA; |
|
176 goto done; |
|
177 } |
|
178 if (a_low >= 0) |
|
179 memcpy(&temp, (char *)a->data + a_low, sizeof(double)); |
|
180 else |
|
181 temp = 0.0; |
|
182 if (!(0.0 <= temp && temp <= (double)INT_MAX && |
|
183 temp == floor(temp))) |
|
184 { ret = GLP_EDATA; |
|
185 goto done; |
|
186 } |
|
187 low[k] = (int)temp; |
|
188 if (a_cap >= 0) |
|
189 memcpy(&temp, (char *)a->data + a_cap, sizeof(double)); |
|
190 else |
|
191 temp = 1.0; |
|
192 if (!((double)low[k] <= temp && temp <= (double)INT_MAX && |
|
193 temp == floor(temp))) |
|
194 { ret = GLP_EDATA; |
|
195 goto done; |
|
196 } |
|
197 cap[k] = (int)temp; |
|
198 if (a_cost >= 0) |
|
199 memcpy(&temp, (char *)a->data + a_cost, sizeof(double)); |
|
200 else |
|
201 temp = 0.0; |
|
202 if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp))) |
|
203 { ret = GLP_EDATA; |
|
204 goto done; |
|
205 } |
|
206 cost[k] = (int)temp; |
|
207 } |
|
208 } |
|
209 /* (artificial arcs) */ |
|
210 sum = 0.0; |
|
211 for (i = 1; i <= G->nv; i++) |
|
212 { v = G->v[i]; |
|
213 if (v_rhs >= 0) |
|
214 memcpy(&temp, (char *)v->data + v_rhs, sizeof(double)); |
|
215 else |
|
216 temp = 0.0; |
|
217 if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp))) |
|
218 { ret = GLP_EDATA; |
|
219 goto done; |
|
220 } |
|
221 if (temp > 0.0) |
|
222 { /* artificial arc from s to original source i */ |
|
223 k++; |
|
224 tail[k] = s; |
|
225 head[k] = i; |
|
226 low[k] = cap[k] = (int)(+temp); /* supply */ |
|
227 cost[k] = 0; |
|
228 sum += (double)temp; |
|
229 } |
|
230 else if (temp < 0.0) |
|
231 { /* artificial arc from original sink i to t */ |
|
232 k++; |
|
233 tail[k] = i; |
|
234 head[k] = t; |
|
235 low[k] = cap[k] = (int)(-temp); /* demand */ |
|
236 cost[k] = 0; |
|
237 } |
|
238 } |
|
239 /* (feedback arc from t to s) */ |
|
240 k++; |
|
241 xassert(k == na); |
|
242 tail[k] = t; |
|
243 head[k] = s; |
|
244 if (sum > (double)INT_MAX) |
|
245 { ret = GLP_EDATA; |
|
246 goto done; |
|
247 } |
|
248 low[k] = cap[k] = (int)sum; /* total supply/demand */ |
|
249 cost[k] = 0; |
|
250 /* find minimal-cost circulation in the resulting network */ |
|
251 ret = okalg(nv, na, tail, head, low, cap, cost, x, pi); |
|
252 switch (ret) |
|
253 { case 0: |
|
254 /* optimal circulation found */ |
|
255 ret = 0; |
|
256 break; |
|
257 case 1: |
|
258 /* no feasible circulation exists */ |
|
259 ret = GLP_ENOPFS; |
|
260 break; |
|
261 case 2: |
|
262 /* integer overflow occured */ |
|
263 ret = GLP_ERANGE; |
|
264 goto done; |
|
265 case 3: |
|
266 /* optimality test failed (logic error) */ |
|
267 ret = GLP_EFAIL; |
|
268 goto done; |
|
269 default: |
|
270 xassert(ret != ret); |
|
271 } |
|
272 /* store solution components */ |
|
273 /* (objective function = the total cost) */ |
|
274 if (sol != NULL) |
|
275 { temp = 0.0; |
|
276 for (k = 1; k <= na; k++) |
|
277 temp += (double)cost[k] * (double)x[k]; |
|
278 *sol = temp; |
|
279 } |
|
280 /* (arc flows) */ |
|
281 if (a_x >= 0) |
|
282 { k = 0; |
|
283 for (i = 1; i <= G->nv; i++) |
|
284 { v = G->v[i]; |
|
285 for (a = v->out; a != NULL; a = a->t_next) |
|
286 { temp = (double)x[++k]; |
|
287 memcpy((char *)a->data + a_x, &temp, sizeof(double)); |
|
288 } |
|
289 } |
|
290 } |
|
291 /* (node potentials = Lagrange multipliers) */ |
|
292 if (v_pi >= 0) |
|
293 { for (i = 1; i <= G->nv; i++) |
|
294 { v = G->v[i]; |
|
295 temp = - (double)pi[i]; |
|
296 memcpy((char *)v->data + v_pi, &temp, sizeof(double)); |
|
297 } |
|
298 } |
|
299 done: /* free working arrays */ |
|
300 xfree(tail); |
|
301 xfree(head); |
|
302 xfree(low); |
|
303 xfree(cap); |
|
304 xfree(cost); |
|
305 xfree(x); |
|
306 xfree(pi); |
|
307 return ret; |
|
308 } |
|
309 |
|
310 /*********************************************************************** |
|
311 * NAME |
|
312 * |
|
313 * glp_maxflow_lp - convert maximum flow problem to LP |
|
314 * |
|
315 * SYNOPSIS |
|
316 * |
|
317 * void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s, |
|
318 * int t, int a_cap); |
|
319 * |
|
320 * DESCRIPTION |
|
321 * |
|
322 * The routine glp_maxflow_lp builds an LP problem, which corresponds |
|
323 * to the maximum flow problem on the specified network G. */ |
|
324 |
|
325 void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s, |
|
326 int t, int a_cap) |
|
327 { glp_vertex *v; |
|
328 glp_arc *a; |
|
329 int i, j, type, ind[1+2]; |
|
330 double cap, val[1+2]; |
|
331 if (!(names == GLP_ON || names == GLP_OFF)) |
|
332 xerror("glp_maxflow_lp: names = %d; invalid parameter\n", |
|
333 names); |
|
334 if (!(1 <= s && s <= G->nv)) |
|
335 xerror("glp_maxflow_lp: s = %d; source node number out of rang" |
|
336 "e\n", s); |
|
337 if (!(1 <= t && t <= G->nv)) |
|
338 xerror("glp_maxflow_lp: t = %d: sink node number out of range " |
|
339 "\n", t); |
|
340 if (s == t) |
|
341 xerror("glp_maxflow_lp: s = t = %d; source and sink nodes must" |
|
342 " be distinct\n", s); |
|
343 if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double)) |
|
344 xerror("glp_maxflow_lp: a_cap = %d; invalid offset\n", a_cap); |
|
345 glp_erase_prob(lp); |
|
346 if (names) glp_set_prob_name(lp, G->name); |
|
347 glp_set_obj_dir(lp, GLP_MAX); |
|
348 glp_add_rows(lp, G->nv); |
|
349 for (i = 1; i <= G->nv; i++) |
|
350 { v = G->v[i]; |
|
351 if (names) glp_set_row_name(lp, i, v->name); |
|
352 if (i == s) |
|
353 type = GLP_LO; |
|
354 else if (i == t) |
|
355 type = GLP_UP; |
|
356 else |
|
357 type = GLP_FX; |
|
358 glp_set_row_bnds(lp, i, type, 0.0, 0.0); |
|
359 } |
|
360 if (G->na > 0) glp_add_cols(lp, G->na); |
|
361 for (i = 1, j = 0; i <= G->nv; i++) |
|
362 { v = G->v[i]; |
|
363 for (a = v->out; a != NULL; a = a->t_next) |
|
364 { j++; |
|
365 if (names) |
|
366 { char name[50+1]; |
|
367 sprintf(name, "x[%d,%d]", a->tail->i, a->head->i); |
|
368 xassert(strlen(name) < sizeof(name)); |
|
369 glp_set_col_name(lp, j, name); |
|
370 } |
|
371 if (a->tail->i != a->head->i) |
|
372 { ind[1] = a->tail->i, val[1] = +1.0; |
|
373 ind[2] = a->head->i, val[2] = -1.0; |
|
374 glp_set_mat_col(lp, j, 2, ind, val); |
|
375 } |
|
376 if (a_cap >= 0) |
|
377 memcpy(&cap, (char *)a->data + a_cap, sizeof(double)); |
|
378 else |
|
379 cap = 1.0; |
|
380 if (cap == DBL_MAX) |
|
381 type = GLP_LO; |
|
382 else if (cap != 0.0) |
|
383 type = GLP_DB; |
|
384 else |
|
385 type = GLP_FX; |
|
386 glp_set_col_bnds(lp, j, type, 0.0, cap); |
|
387 if (a->tail->i == s) |
|
388 glp_set_obj_coef(lp, j, +1.0); |
|
389 else if (a->head->i == s) |
|
390 glp_set_obj_coef(lp, j, -1.0); |
|
391 } |
|
392 } |
|
393 xassert(j == G->na); |
|
394 return; |
|
395 } |
|
396 |
|
397 int glp_maxflow_ffalg(glp_graph *G, int s, int t, int a_cap, |
|
398 double *sol, int a_x, int v_cut) |
|
399 { /* find maximal flow with Ford-Fulkerson algorithm */ |
|
400 glp_vertex *v; |
|
401 glp_arc *a; |
|
402 int nv, na, i, k, flag, *tail, *head, *cap, *x, ret; |
|
403 char *cut; |
|
404 double temp; |
|
405 if (!(1 <= s && s <= G->nv)) |
|
406 xerror("glp_maxflow_ffalg: s = %d; source node number out of r" |
|
407 "ange\n", s); |
|
408 if (!(1 <= t && t <= G->nv)) |
|
409 xerror("glp_maxflow_ffalg: t = %d: sink node number out of ran" |
|
410 "ge\n", t); |
|
411 if (s == t) |
|
412 xerror("glp_maxflow_ffalg: s = t = %d; source and sink nodes m" |
|
413 "ust be distinct\n", s); |
|
414 if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double)) |
|
415 xerror("glp_maxflow_ffalg: a_cap = %d; invalid offset\n", |
|
416 a_cap); |
|
417 if (v_cut >= 0 && v_cut > G->v_size - (int)sizeof(int)) |
|
418 xerror("glp_maxflow_ffalg: v_cut = %d; invalid offset\n", |
|
419 v_cut); |
|
420 /* allocate working arrays */ |
|
421 nv = G->nv; |
|
422 na = G->na; |
|
423 tail = xcalloc(1+na, sizeof(int)); |
|
424 head = xcalloc(1+na, sizeof(int)); |
|
425 cap = xcalloc(1+na, sizeof(int)); |
|
426 x = xcalloc(1+na, sizeof(int)); |
|
427 if (v_cut < 0) |
|
428 cut = NULL; |
|
429 else |
|
430 cut = xcalloc(1+nv, sizeof(char)); |
|
431 /* copy the flow network */ |
|
432 k = 0; |
|
433 for (i = 1; i <= G->nv; i++) |
|
434 { v = G->v[i]; |
|
435 for (a = v->out; a != NULL; a = a->t_next) |
|
436 { k++; |
|
437 tail[k] = a->tail->i; |
|
438 head[k] = a->head->i; |
|
439 if (tail[k] == head[k]) |
|
440 { ret = GLP_EDATA; |
|
441 goto done; |
|
442 } |
|
443 if (a_cap >= 0) |
|
444 memcpy(&temp, (char *)a->data + a_cap, sizeof(double)); |
|
445 else |
|
446 temp = 1.0; |
|
447 if (!(0.0 <= temp && temp <= (double)INT_MAX && |
|
448 temp == floor(temp))) |
|
449 { ret = GLP_EDATA; |
|
450 goto done; |
|
451 } |
|
452 cap[k] = (int)temp; |
|
453 } |
|
454 } |
|
455 xassert(k == na); |
|
456 /* find maximal flow in the flow network */ |
|
457 ffalg(nv, na, tail, head, s, t, cap, x, cut); |
|
458 ret = 0; |
|
459 /* store solution components */ |
|
460 /* (objective function = total flow through the network) */ |
|
461 if (sol != NULL) |
|
462 { temp = 0.0; |
|
463 for (k = 1; k <= na; k++) |
|
464 { if (tail[k] == s) |
|
465 temp += (double)x[k]; |
|
466 else if (head[k] == s) |
|
467 temp -= (double)x[k]; |
|
468 } |
|
469 *sol = temp; |
|
470 } |
|
471 /* (arc flows) */ |
|
472 if (a_x >= 0) |
|
473 { k = 0; |
|
474 for (i = 1; i <= G->nv; i++) |
|
475 { v = G->v[i]; |
|
476 for (a = v->out; a != NULL; a = a->t_next) |
|
477 { temp = (double)x[++k]; |
|
478 memcpy((char *)a->data + a_x, &temp, sizeof(double)); |
|
479 } |
|
480 } |
|
481 } |
|
482 /* (node flags) */ |
|
483 if (v_cut >= 0) |
|
484 { for (i = 1; i <= G->nv; i++) |
|
485 { v = G->v[i]; |
|
486 flag = cut[i]; |
|
487 memcpy((char *)v->data + v_cut, &flag, sizeof(int)); |
|
488 } |
|
489 } |
|
490 done: /* free working arrays */ |
|
491 xfree(tail); |
|
492 xfree(head); |
|
493 xfree(cap); |
|
494 xfree(x); |
|
495 if (cut != NULL) xfree(cut); |
|
496 return ret; |
|
497 } |
|
498 |
|
499 /*********************************************************************** |
|
500 * NAME |
|
501 * |
|
502 * glp_check_asnprob - check correctness of assignment problem data |
|
503 * |
|
504 * SYNOPSIS |
|
505 * |
|
506 * int glp_check_asnprob(glp_graph *G, int v_set); |
|
507 * |
|
508 * RETURNS |
|
509 * |
|
510 * If the specified assignment problem data are correct, the routine |
|
511 * glp_check_asnprob returns zero, otherwise, non-zero. */ |
|
512 |
|
513 int glp_check_asnprob(glp_graph *G, int v_set) |
|
514 { glp_vertex *v; |
|
515 int i, k, ret = 0; |
|
516 if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int)) |
|
517 xerror("glp_check_asnprob: v_set = %d; invalid offset\n", |
|
518 v_set); |
|
519 for (i = 1; i <= G->nv; i++) |
|
520 { v = G->v[i]; |
|
521 if (v_set >= 0) |
|
522 { memcpy(&k, (char *)v->data + v_set, sizeof(int)); |
|
523 if (k == 0) |
|
524 { if (v->in != NULL) |
|
525 { ret = 1; |
|
526 break; |
|
527 } |
|
528 } |
|
529 else if (k == 1) |
|
530 { if (v->out != NULL) |
|
531 { ret = 2; |
|
532 break; |
|
533 } |
|
534 } |
|
535 else |
|
536 { ret = 3; |
|
537 break; |
|
538 } |
|
539 } |
|
540 else |
|
541 { if (v->in != NULL && v->out != NULL) |
|
542 { ret = 4; |
|
543 break; |
|
544 } |
|
545 } |
|
546 } |
|
547 return ret; |
|
548 } |
|
549 |
|
550 /*********************************************************************** |
|
551 * NAME |
|
552 * |
|
553 * glp_asnprob_lp - convert assignment problem to LP |
|
554 * |
|
555 * SYNOPSIS |
|
556 * |
|
557 * int glp_asnprob_lp(glp_prob *P, int form, glp_graph *G, int names, |
|
558 * int v_set, int a_cost); |
|
559 * |
|
560 * DESCRIPTION |
|
561 * |
|
562 * The routine glp_asnprob_lp builds an LP problem, which corresponds |
|
563 * to the assignment problem on the specified graph G. |
|
564 * |
|
565 * RETURNS |
|
566 * |
|
567 * If the LP problem has been successfully built, the routine returns |
|
568 * zero, otherwise, non-zero. */ |
|
569 |
|
570 int glp_asnprob_lp(glp_prob *P, int form, glp_graph *G, int names, |
|
571 int v_set, int a_cost) |
|
572 { glp_vertex *v; |
|
573 glp_arc *a; |
|
574 int i, j, ret, ind[1+2]; |
|
575 double cost, val[1+2]; |
|
576 if (!(form == GLP_ASN_MIN || form == GLP_ASN_MAX || |
|
577 form == GLP_ASN_MMP)) |
|
578 xerror("glp_asnprob_lp: form = %d; invalid parameter\n", |
|
579 form); |
|
580 if (!(names == GLP_ON || names == GLP_OFF)) |
|
581 xerror("glp_asnprob_lp: names = %d; invalid parameter\n", |
|
582 names); |
|
583 if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int)) |
|
584 xerror("glp_asnprob_lp: v_set = %d; invalid offset\n", |
|
585 v_set); |
|
586 if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double)) |
|
587 xerror("glp_asnprob_lp: a_cost = %d; invalid offset\n", |
|
588 a_cost); |
|
589 ret = glp_check_asnprob(G, v_set); |
|
590 if (ret != 0) goto done; |
|
591 glp_erase_prob(P); |
|
592 if (names) glp_set_prob_name(P, G->name); |
|
593 glp_set_obj_dir(P, form == GLP_ASN_MIN ? GLP_MIN : GLP_MAX); |
|
594 if (G->nv > 0) glp_add_rows(P, G->nv); |
|
595 for (i = 1; i <= G->nv; i++) |
|
596 { v = G->v[i]; |
|
597 if (names) glp_set_row_name(P, i, v->name); |
|
598 glp_set_row_bnds(P, i, form == GLP_ASN_MMP ? GLP_UP : GLP_FX, |
|
599 1.0, 1.0); |
|
600 } |
|
601 if (G->na > 0) glp_add_cols(P, G->na); |
|
602 for (i = 1, j = 0; i <= G->nv; i++) |
|
603 { v = G->v[i]; |
|
604 for (a = v->out; a != NULL; a = a->t_next) |
|
605 { j++; |
|
606 if (names) |
|
607 { char name[50+1]; |
|
608 sprintf(name, "x[%d,%d]", a->tail->i, a->head->i); |
|
609 xassert(strlen(name) < sizeof(name)); |
|
610 glp_set_col_name(P, j, name); |
|
611 } |
|
612 ind[1] = a->tail->i, val[1] = +1.0; |
|
613 ind[2] = a->head->i, val[2] = +1.0; |
|
614 glp_set_mat_col(P, j, 2, ind, val); |
|
615 glp_set_col_bnds(P, j, GLP_DB, 0.0, 1.0); |
|
616 if (a_cost >= 0) |
|
617 memcpy(&cost, (char *)a->data + a_cost, sizeof(double)); |
|
618 else |
|
619 cost = 1.0; |
|
620 glp_set_obj_coef(P, j, cost); |
|
621 } |
|
622 } |
|
623 xassert(j == G->na); |
|
624 done: return ret; |
|
625 } |
|
626 |
|
627 /**********************************************************************/ |
|
628 |
|
629 int glp_asnprob_okalg(int form, glp_graph *G, int v_set, int a_cost, |
|
630 double *sol, int a_x) |
|
631 { /* solve assignment problem with out-of-kilter algorithm */ |
|
632 glp_vertex *v; |
|
633 glp_arc *a; |
|
634 int nv, na, i, k, *tail, *head, *low, *cap, *cost, *x, *pi, ret; |
|
635 double temp; |
|
636 if (!(form == GLP_ASN_MIN || form == GLP_ASN_MAX || |
|
637 form == GLP_ASN_MMP)) |
|
638 xerror("glp_asnprob_okalg: form = %d; invalid parameter\n", |
|
639 form); |
|
640 if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int)) |
|
641 xerror("glp_asnprob_okalg: v_set = %d; invalid offset\n", |
|
642 v_set); |
|
643 if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double)) |
|
644 xerror("glp_asnprob_okalg: a_cost = %d; invalid offset\n", |
|
645 a_cost); |
|
646 if (a_x >= 0 && a_x > G->a_size - (int)sizeof(int)) |
|
647 xerror("glp_asnprob_okalg: a_x = %d; invalid offset\n", a_x); |
|
648 if (glp_check_asnprob(G, v_set)) |
|
649 return GLP_EDATA; |
|
650 /* nv is the total number of nodes in the resulting network */ |
|
651 nv = G->nv + 1; |
|
652 /* na is the total number of arcs in the resulting network */ |
|
653 na = G->na + G->nv; |
|
654 /* allocate working arrays */ |
|
655 tail = xcalloc(1+na, sizeof(int)); |
|
656 head = xcalloc(1+na, sizeof(int)); |
|
657 low = xcalloc(1+na, sizeof(int)); |
|
658 cap = xcalloc(1+na, sizeof(int)); |
|
659 cost = xcalloc(1+na, sizeof(int)); |
|
660 x = xcalloc(1+na, sizeof(int)); |
|
661 pi = xcalloc(1+nv, sizeof(int)); |
|
662 /* construct the resulting network */ |
|
663 k = 0; |
|
664 /* (original arcs) */ |
|
665 for (i = 1; i <= G->nv; i++) |
|
666 { v = G->v[i]; |
|
667 for (a = v->out; a != NULL; a = a->t_next) |
|
668 { k++; |
|
669 tail[k] = a->tail->i; |
|
670 head[k] = a->head->i; |
|
671 low[k] = 0; |
|
672 cap[k] = 1; |
|
673 if (a_cost >= 0) |
|
674 memcpy(&temp, (char *)a->data + a_cost, sizeof(double)); |
|
675 else |
|
676 temp = 1.0; |
|
677 if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp))) |
|
678 { ret = GLP_EDATA; |
|
679 goto done; |
|
680 } |
|
681 cost[k] = (int)temp; |
|
682 if (form != GLP_ASN_MIN) cost[k] = - cost[k]; |
|
683 } |
|
684 } |
|
685 /* (artificial arcs) */ |
|
686 for (i = 1; i <= G->nv; i++) |
|
687 { v = G->v[i]; |
|
688 k++; |
|
689 if (v->out == NULL) |
|
690 tail[k] = i, head[k] = nv; |
|
691 else if (v->in == NULL) |
|
692 tail[k] = nv, head[k] = i; |
|
693 else |
|
694 xassert(v != v); |
|
695 low[k] = (form == GLP_ASN_MMP ? 0 : 1); |
|
696 cap[k] = 1; |
|
697 cost[k] = 0; |
|
698 } |
|
699 xassert(k == na); |
|
700 /* find minimal-cost circulation in the resulting network */ |
|
701 ret = okalg(nv, na, tail, head, low, cap, cost, x, pi); |
|
702 switch (ret) |
|
703 { case 0: |
|
704 /* optimal circulation found */ |
|
705 ret = 0; |
|
706 break; |
|
707 case 1: |
|
708 /* no feasible circulation exists */ |
|
709 ret = GLP_ENOPFS; |
|
710 break; |
|
711 case 2: |
|
712 /* integer overflow occured */ |
|
713 ret = GLP_ERANGE; |
|
714 goto done; |
|
715 case 3: |
|
716 /* optimality test failed (logic error) */ |
|
717 ret = GLP_EFAIL; |
|
718 goto done; |
|
719 default: |
|
720 xassert(ret != ret); |
|
721 } |
|
722 /* store solution components */ |
|
723 /* (objective function = the total cost) */ |
|
724 if (sol != NULL) |
|
725 { temp = 0.0; |
|
726 for (k = 1; k <= na; k++) |
|
727 temp += (double)cost[k] * (double)x[k]; |
|
728 if (form != GLP_ASN_MIN) temp = - temp; |
|
729 *sol = temp; |
|
730 } |
|
731 /* (arc flows) */ |
|
732 if (a_x >= 0) |
|
733 { k = 0; |
|
734 for (i = 1; i <= G->nv; i++) |
|
735 { v = G->v[i]; |
|
736 for (a = v->out; a != NULL; a = a->t_next) |
|
737 { k++; |
|
738 if (ret == 0) |
|
739 xassert(x[k] == 0 || x[k] == 1); |
|
740 memcpy((char *)a->data + a_x, &x[k], sizeof(int)); |
|
741 } |
|
742 } |
|
743 } |
|
744 done: /* free working arrays */ |
|
745 xfree(tail); |
|
746 xfree(head); |
|
747 xfree(low); |
|
748 xfree(cap); |
|
749 xfree(cost); |
|
750 xfree(x); |
|
751 xfree(pi); |
|
752 return ret; |
|
753 } |
|
754 |
|
755 /*********************************************************************** |
|
756 * NAME |
|
757 * |
|
758 * glp_asnprob_hall - find bipartite matching of maximum cardinality |
|
759 * |
|
760 * SYNOPSIS |
|
761 * |
|
762 * int glp_asnprob_hall(glp_graph *G, int v_set, int a_x); |
|
763 * |
|
764 * DESCRIPTION |
|
765 * |
|
766 * The routine glp_asnprob_hall finds a matching of maximal cardinality |
|
767 * in the specified bipartite graph G. It uses a version of the Fortran |
|
768 * routine MC21A developed by I.S.Duff [1], which implements Hall's |
|
769 * algorithm [2]. |
|
770 * |
|
771 * RETURNS |
|
772 * |
|
773 * The routine glp_asnprob_hall returns the cardinality of the matching |
|
774 * found. However, if the specified graph is incorrect (as detected by |
|
775 * the routine glp_check_asnprob), the routine returns negative value. |
|
776 * |
|
777 * REFERENCES |
|
778 * |
|
779 * 1. I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM |
|
780 * Trans. on Math. Softw. 7 (1981), 387-390. |
|
781 * |
|
782 * 2. M.Hall, "An Algorithm for distinct representatives," Amer. Math. |
|
783 * Monthly 63 (1956), 716-717. */ |
|
784 |
|
785 int glp_asnprob_hall(glp_graph *G, int v_set, int a_x) |
|
786 { glp_vertex *v; |
|
787 glp_arc *a; |
|
788 int card, i, k, loc, n, n1, n2, xij; |
|
789 int *num, *icn, *ip, *lenr, *iperm, *pr, *arp, *cv, *out; |
|
790 if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int)) |
|
791 xerror("glp_asnprob_hall: v_set = %d; invalid offset\n", |
|
792 v_set); |
|
793 if (a_x >= 0 && a_x > G->a_size - (int)sizeof(int)) |
|
794 xerror("glp_asnprob_hall: a_x = %d; invalid offset\n", a_x); |
|
795 if (glp_check_asnprob(G, v_set)) |
|
796 return -1; |
|
797 /* determine the number of vertices in sets R and S and renumber |
|
798 vertices in S which correspond to columns of the matrix; skip |
|
799 all isolated vertices */ |
|
800 num = xcalloc(1+G->nv, sizeof(int)); |
|
801 n1 = n2 = 0; |
|
802 for (i = 1; i <= G->nv; i++) |
|
803 { v = G->v[i]; |
|
804 if (v->in == NULL && v->out != NULL) |
|
805 n1++, num[i] = 0; /* vertex in R */ |
|
806 else if (v->in != NULL && v->out == NULL) |
|
807 n2++, num[i] = n2; /* vertex in S */ |
|
808 else |
|
809 { xassert(v->in == NULL && v->out == NULL); |
|
810 num[i] = -1; /* isolated vertex */ |
|
811 } |
|
812 } |
|
813 /* the matrix must be square, thus, if it has more columns than |
|
814 rows, extra rows will be just empty, and vice versa */ |
|
815 n = (n1 >= n2 ? n1 : n2); |
|
816 /* allocate working arrays */ |
|
817 icn = xcalloc(1+G->na, sizeof(int)); |
|
818 ip = xcalloc(1+n, sizeof(int)); |
|
819 lenr = xcalloc(1+n, sizeof(int)); |
|
820 iperm = xcalloc(1+n, sizeof(int)); |
|
821 pr = xcalloc(1+n, sizeof(int)); |
|
822 arp = xcalloc(1+n, sizeof(int)); |
|
823 cv = xcalloc(1+n, sizeof(int)); |
|
824 out = xcalloc(1+n, sizeof(int)); |
|
825 /* build the adjacency matrix of the bipartite graph in row-wise |
|
826 format (rows are vertices in R, columns are vertices in S) */ |
|
827 k = 0, loc = 1; |
|
828 for (i = 1; i <= G->nv; i++) |
|
829 { if (num[i] != 0) continue; |
|
830 /* vertex i in R */ |
|
831 ip[++k] = loc; |
|
832 v = G->v[i]; |
|
833 for (a = v->out; a != NULL; a = a->t_next) |
|
834 { xassert(num[a->head->i] != 0); |
|
835 icn[loc++] = num[a->head->i]; |
|
836 } |
|
837 lenr[k] = loc - ip[k]; |
|
838 } |
|
839 xassert(loc-1 == G->na); |
|
840 /* make all extra rows empty (all extra columns are empty due to |
|
841 the row-wise format used) */ |
|
842 for (k++; k <= n; k++) |
|
843 ip[k] = loc, lenr[k] = 0; |
|
844 /* find a row permutation that maximizes the number of non-zeros |
|
845 on the main diagonal */ |
|
846 card = mc21a(n, icn, ip, lenr, iperm, pr, arp, cv, out); |
|
847 #if 1 /* 18/II-2010 */ |
|
848 /* FIXED: if card = n, arp remains clobbered on exit */ |
|
849 for (i = 1; i <= n; i++) |
|
850 arp[i] = 0; |
|
851 for (i = 1; i <= card; i++) |
|
852 { k = iperm[i]; |
|
853 xassert(1 <= k && k <= n); |
|
854 xassert(arp[k] == 0); |
|
855 arp[k] = i; |
|
856 } |
|
857 #endif |
|
858 /* store solution, if necessary */ |
|
859 if (a_x < 0) goto skip; |
|
860 k = 0; |
|
861 for (i = 1; i <= G->nv; i++) |
|
862 { if (num[i] != 0) continue; |
|
863 /* vertex i in R */ |
|
864 k++; |
|
865 v = G->v[i]; |
|
866 for (a = v->out; a != NULL; a = a->t_next) |
|
867 { /* arp[k] is the number of matched column or zero */ |
|
868 if (arp[k] == num[a->head->i]) |
|
869 { xassert(arp[k] != 0); |
|
870 xij = 1; |
|
871 } |
|
872 else |
|
873 xij = 0; |
|
874 memcpy((char *)a->data + a_x, &xij, sizeof(int)); |
|
875 } |
|
876 } |
|
877 skip: /* free working arrays */ |
|
878 xfree(num); |
|
879 xfree(icn); |
|
880 xfree(ip); |
|
881 xfree(lenr); |
|
882 xfree(iperm); |
|
883 xfree(pr); |
|
884 xfree(arp); |
|
885 xfree(cv); |
|
886 xfree(out); |
|
887 return card; |
|
888 } |
|
889 |
|
890 /*********************************************************************** |
|
891 * NAME |
|
892 * |
|
893 * glp_cpp - solve critical path problem |
|
894 * |
|
895 * SYNOPSIS |
|
896 * |
|
897 * double glp_cpp(glp_graph *G, int v_t, int v_es, int v_ls); |
|
898 * |
|
899 * DESCRIPTION |
|
900 * |
|
901 * The routine glp_cpp solves the critical path problem represented in |
|
902 * the form of the project network. |
|
903 * |
|
904 * The parameter G is a pointer to the graph object, which specifies |
|
905 * the project network. This graph must be acyclic. Multiple arcs are |
|
906 * allowed being considered as single arcs. |
|
907 * |
|
908 * The parameter v_t specifies an offset of the field of type double |
|
909 * in the vertex data block, which contains time t[i] >= 0 needed to |
|
910 * perform corresponding job j. If v_t < 0, it is assumed that t[i] = 1 |
|
911 * for all jobs. |
|
912 * |
|
913 * The parameter v_es specifies an offset of the field of type double |
|
914 * in the vertex data block, to which the routine stores earliest start |
|
915 * time for corresponding job. If v_es < 0, this time is not stored. |
|
916 * |
|
917 * The parameter v_ls specifies an offset of the field of type double |
|
918 * in the vertex data block, to which the routine stores latest start |
|
919 * time for corresponding job. If v_ls < 0, this time is not stored. |
|
920 * |
|
921 * RETURNS |
|
922 * |
|
923 * The routine glp_cpp returns the minimal project duration, that is, |
|
924 * minimal time needed to perform all jobs in the project. */ |
|
925 |
|
926 static void sorting(glp_graph *G, int list[]); |
|
927 |
|
928 double glp_cpp(glp_graph *G, int v_t, int v_es, int v_ls) |
|
929 { glp_vertex *v; |
|
930 glp_arc *a; |
|
931 int i, j, k, nv, *list; |
|
932 double temp, total, *t, *es, *ls; |
|
933 if (v_t >= 0 && v_t > G->v_size - (int)sizeof(double)) |
|
934 xerror("glp_cpp: v_t = %d; invalid offset\n", v_t); |
|
935 if (v_es >= 0 && v_es > G->v_size - (int)sizeof(double)) |
|
936 xerror("glp_cpp: v_es = %d; invalid offset\n", v_es); |
|
937 if (v_ls >= 0 && v_ls > G->v_size - (int)sizeof(double)) |
|
938 xerror("glp_cpp: v_ls = %d; invalid offset\n", v_ls); |
|
939 nv = G->nv; |
|
940 if (nv == 0) |
|
941 { total = 0.0; |
|
942 goto done; |
|
943 } |
|
944 /* allocate working arrays */ |
|
945 t = xcalloc(1+nv, sizeof(double)); |
|
946 es = xcalloc(1+nv, sizeof(double)); |
|
947 ls = xcalloc(1+nv, sizeof(double)); |
|
948 list = xcalloc(1+nv, sizeof(int)); |
|
949 /* retrieve job times */ |
|
950 for (i = 1; i <= nv; i++) |
|
951 { v = G->v[i]; |
|
952 if (v_t >= 0) |
|
953 { memcpy(&t[i], (char *)v->data + v_t, sizeof(double)); |
|
954 if (t[i] < 0.0) |
|
955 xerror("glp_cpp: t[%d] = %g; invalid time\n", i, t[i]); |
|
956 } |
|
957 else |
|
958 t[i] = 1.0; |
|
959 } |
|
960 /* perform topological sorting to determine the list of nodes |
|
961 (jobs) such that if list[k] = i and list[kk] = j and there |
|
962 exists arc (i->j), then k < kk */ |
|
963 sorting(G, list); |
|
964 /* FORWARD PASS */ |
|
965 /* determine earliest start times */ |
|
966 for (k = 1; k <= nv; k++) |
|
967 { j = list[k]; |
|
968 es[j] = 0.0; |
|
969 for (a = G->v[j]->in; a != NULL; a = a->h_next) |
|
970 { i = a->tail->i; |
|
971 /* there exists arc (i->j) in the project network */ |
|
972 temp = es[i] + t[i]; |
|
973 if (es[j] < temp) es[j] = temp; |
|
974 } |
|
975 } |
|
976 /* determine the minimal project duration */ |
|
977 total = 0.0; |
|
978 for (i = 1; i <= nv; i++) |
|
979 { temp = es[i] + t[i]; |
|
980 if (total < temp) total = temp; |
|
981 } |
|
982 /* BACKWARD PASS */ |
|
983 /* determine latest start times */ |
|
984 for (k = nv; k >= 1; k--) |
|
985 { i = list[k]; |
|
986 ls[i] = total - t[i]; |
|
987 for (a = G->v[i]->out; a != NULL; a = a->t_next) |
|
988 { j = a->head->i; |
|
989 /* there exists arc (i->j) in the project network */ |
|
990 temp = ls[j] - t[i]; |
|
991 if (ls[i] > temp) ls[i] = temp; |
|
992 } |
|
993 /* avoid possible round-off errors */ |
|
994 if (ls[i] < es[i]) ls[i] = es[i]; |
|
995 } |
|
996 /* store results, if necessary */ |
|
997 if (v_es >= 0) |
|
998 { for (i = 1; i <= nv; i++) |
|
999 { v = G->v[i]; |
|
1000 memcpy((char *)v->data + v_es, &es[i], sizeof(double)); |
|
1001 } |
|
1002 } |
|
1003 if (v_ls >= 0) |
|
1004 { for (i = 1; i <= nv; i++) |
|
1005 { v = G->v[i]; |
|
1006 memcpy((char *)v->data + v_ls, &ls[i], sizeof(double)); |
|
1007 } |
|
1008 } |
|
1009 /* free working arrays */ |
|
1010 xfree(t); |
|
1011 xfree(es); |
|
1012 xfree(ls); |
|
1013 xfree(list); |
|
1014 done: return total; |
|
1015 } |
|
1016 |
|
1017 static void sorting(glp_graph *G, int list[]) |
|
1018 { /* perform topological sorting to determine the list of nodes |
|
1019 (jobs) such that if list[k] = i and list[kk] = j and there |
|
1020 exists arc (i->j), then k < kk */ |
|
1021 int i, k, nv, v_size, *num; |
|
1022 void **save; |
|
1023 nv = G->nv; |
|
1024 v_size = G->v_size; |
|
1025 save = xcalloc(1+nv, sizeof(void *)); |
|
1026 num = xcalloc(1+nv, sizeof(int)); |
|
1027 G->v_size = sizeof(int); |
|
1028 for (i = 1; i <= nv; i++) |
|
1029 { save[i] = G->v[i]->data; |
|
1030 G->v[i]->data = &num[i]; |
|
1031 list[i] = 0; |
|
1032 } |
|
1033 if (glp_top_sort(G, 0) != 0) |
|
1034 xerror("glp_cpp: project network is not acyclic\n"); |
|
1035 G->v_size = v_size; |
|
1036 for (i = 1; i <= nv; i++) |
|
1037 { G->v[i]->data = save[i]; |
|
1038 k = num[i]; |
|
1039 xassert(1 <= k && k <= nv); |
|
1040 xassert(list[k] == 0); |
|
1041 list[k] = i; |
|
1042 } |
|
1043 xfree(save); |
|
1044 xfree(num); |
|
1045 return; |
|
1046 } |
|
1047 |
|
1048 /* eof */ |