lemon-project-template-glpk
view deps/glpk/src/glpapi17.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 /* glpapi17.c (flow network problems) */
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 ***********************************************************************/
25 #include "glpapi.h"
26 #include "glpnet.h"
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. */
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 }
115 /**********************************************************************/
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 }
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. */
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 }
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 }
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. */
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 }
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. */
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 }
627 /**********************************************************************/
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 }
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. */
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 }
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. */
926 static void sorting(glp_graph *G, int list[]);
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 }
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 }
1048 /* eof */