lemon-project-template-glpk
comparison deps/glpk/src/glpapi17.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 21:43:29 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:50328f0754e5 |
---|---|
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, 2011 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 */ |