[9] | 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 */ |
---|