1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpapi17.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,1048 @@
1.4 +/* glpapi17.c (flow network problems) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpapi.h"
1.29 +#include "glpnet.h"
1.30 +
1.31 +/***********************************************************************
1.32 +* NAME
1.33 +*
1.34 +* glp_mincost_lp - convert minimum cost flow problem to LP
1.35 +*
1.36 +* SYNOPSIS
1.37 +*
1.38 +* void glp_mincost_lp(glp_prob *lp, glp_graph *G, int names,
1.39 +* int v_rhs, int a_low, int a_cap, int a_cost);
1.40 +*
1.41 +* DESCRIPTION
1.42 +*
1.43 +* The routine glp_mincost_lp builds an LP problem, which corresponds
1.44 +* to the minimum cost flow problem on the specified network G. */
1.45 +
1.46 +void glp_mincost_lp(glp_prob *lp, glp_graph *G, int names, int v_rhs,
1.47 + int a_low, int a_cap, int a_cost)
1.48 +{ glp_vertex *v;
1.49 + glp_arc *a;
1.50 + int i, j, type, ind[1+2];
1.51 + double rhs, low, cap, cost, val[1+2];
1.52 + if (!(names == GLP_ON || names == GLP_OFF))
1.53 + xerror("glp_mincost_lp: names = %d; invalid parameter\n",
1.54 + names);
1.55 + if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
1.56 + xerror("glp_mincost_lp: v_rhs = %d; invalid offset\n", v_rhs);
1.57 + if (a_low >= 0 && a_low > G->a_size - (int)sizeof(double))
1.58 + xerror("glp_mincost_lp: a_low = %d; invalid offset\n", a_low);
1.59 + if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
1.60 + xerror("glp_mincost_lp: a_cap = %d; invalid offset\n", a_cap);
1.61 + if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
1.62 + xerror("glp_mincost_lp: a_cost = %d; invalid offset\n", a_cost)
1.63 + ;
1.64 + glp_erase_prob(lp);
1.65 + if (names) glp_set_prob_name(lp, G->name);
1.66 + if (G->nv > 0) glp_add_rows(lp, G->nv);
1.67 + for (i = 1; i <= G->nv; i++)
1.68 + { v = G->v[i];
1.69 + if (names) glp_set_row_name(lp, i, v->name);
1.70 + if (v_rhs >= 0)
1.71 + memcpy(&rhs, (char *)v->data + v_rhs, sizeof(double));
1.72 + else
1.73 + rhs = 0.0;
1.74 + glp_set_row_bnds(lp, i, GLP_FX, rhs, rhs);
1.75 + }
1.76 + if (G->na > 0) glp_add_cols(lp, G->na);
1.77 + for (i = 1, j = 0; i <= G->nv; i++)
1.78 + { v = G->v[i];
1.79 + for (a = v->out; a != NULL; a = a->t_next)
1.80 + { j++;
1.81 + if (names)
1.82 + { char name[50+1];
1.83 + sprintf(name, "x[%d,%d]", a->tail->i, a->head->i);
1.84 + xassert(strlen(name) < sizeof(name));
1.85 + glp_set_col_name(lp, j, name);
1.86 + }
1.87 + if (a->tail->i != a->head->i)
1.88 + { ind[1] = a->tail->i, val[1] = +1.0;
1.89 + ind[2] = a->head->i, val[2] = -1.0;
1.90 + glp_set_mat_col(lp, j, 2, ind, val);
1.91 + }
1.92 + if (a_low >= 0)
1.93 + memcpy(&low, (char *)a->data + a_low, sizeof(double));
1.94 + else
1.95 + low = 0.0;
1.96 + if (a_cap >= 0)
1.97 + memcpy(&cap, (char *)a->data + a_cap, sizeof(double));
1.98 + else
1.99 + cap = 1.0;
1.100 + if (cap == DBL_MAX)
1.101 + type = GLP_LO;
1.102 + else if (low != cap)
1.103 + type = GLP_DB;
1.104 + else
1.105 + type = GLP_FX;
1.106 + glp_set_col_bnds(lp, j, type, low, cap);
1.107 + if (a_cost >= 0)
1.108 + memcpy(&cost, (char *)a->data + a_cost, sizeof(double));
1.109 + else
1.110 + cost = 0.0;
1.111 + glp_set_obj_coef(lp, j, cost);
1.112 + }
1.113 + }
1.114 + xassert(j == G->na);
1.115 + return;
1.116 +}
1.117 +
1.118 +/**********************************************************************/
1.119 +
1.120 +int glp_mincost_okalg(glp_graph *G, int v_rhs, int a_low, int a_cap,
1.121 + int a_cost, double *sol, int a_x, int v_pi)
1.122 +{ /* find minimum-cost flow with out-of-kilter algorithm */
1.123 + glp_vertex *v;
1.124 + glp_arc *a;
1.125 + int nv, na, i, k, s, t, *tail, *head, *low, *cap, *cost, *x, *pi,
1.126 + ret;
1.127 + double sum, temp;
1.128 + if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
1.129 + xerror("glp_mincost_okalg: v_rhs = %d; invalid offset\n",
1.130 + v_rhs);
1.131 + if (a_low >= 0 && a_low > G->a_size - (int)sizeof(double))
1.132 + xerror("glp_mincost_okalg: a_low = %d; invalid offset\n",
1.133 + a_low);
1.134 + if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
1.135 + xerror("glp_mincost_okalg: a_cap = %d; invalid offset\n",
1.136 + a_cap);
1.137 + if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
1.138 + xerror("glp_mincost_okalg: a_cost = %d; invalid offset\n",
1.139 + a_cost);
1.140 + if (a_x >= 0 && a_x > G->a_size - (int)sizeof(double))
1.141 + xerror("glp_mincost_okalg: a_x = %d; invalid offset\n", a_x);
1.142 + if (v_pi >= 0 && v_pi > G->v_size - (int)sizeof(double))
1.143 + xerror("glp_mincost_okalg: v_pi = %d; invalid offset\n", v_pi);
1.144 + /* s is artificial source node */
1.145 + s = G->nv + 1;
1.146 + /* t is artificial sink node */
1.147 + t = s + 1;
1.148 + /* nv is the total number of nodes in the resulting network */
1.149 + nv = t;
1.150 + /* na is the total number of arcs in the resulting network */
1.151 + na = G->na + 1;
1.152 + for (i = 1; i <= G->nv; i++)
1.153 + { v = G->v[i];
1.154 + if (v_rhs >= 0)
1.155 + memcpy(&temp, (char *)v->data + v_rhs, sizeof(double));
1.156 + else
1.157 + temp = 0.0;
1.158 + if (temp != 0.0) na++;
1.159 + }
1.160 + /* allocate working arrays */
1.161 + tail = xcalloc(1+na, sizeof(int));
1.162 + head = xcalloc(1+na, sizeof(int));
1.163 + low = xcalloc(1+na, sizeof(int));
1.164 + cap = xcalloc(1+na, sizeof(int));
1.165 + cost = xcalloc(1+na, sizeof(int));
1.166 + x = xcalloc(1+na, sizeof(int));
1.167 + pi = xcalloc(1+nv, sizeof(int));
1.168 + /* construct the resulting network */
1.169 + k = 0;
1.170 + /* (original arcs) */
1.171 + for (i = 1; i <= G->nv; i++)
1.172 + { v = G->v[i];
1.173 + for (a = v->out; a != NULL; a = a->t_next)
1.174 + { k++;
1.175 + tail[k] = a->tail->i;
1.176 + head[k] = a->head->i;
1.177 + if (tail[k] == head[k])
1.178 + { ret = GLP_EDATA;
1.179 + goto done;
1.180 + }
1.181 + if (a_low >= 0)
1.182 + memcpy(&temp, (char *)a->data + a_low, sizeof(double));
1.183 + else
1.184 + temp = 0.0;
1.185 + if (!(0.0 <= temp && temp <= (double)INT_MAX &&
1.186 + temp == floor(temp)))
1.187 + { ret = GLP_EDATA;
1.188 + goto done;
1.189 + }
1.190 + low[k] = (int)temp;
1.191 + if (a_cap >= 0)
1.192 + memcpy(&temp, (char *)a->data + a_cap, sizeof(double));
1.193 + else
1.194 + temp = 1.0;
1.195 + if (!((double)low[k] <= temp && temp <= (double)INT_MAX &&
1.196 + temp == floor(temp)))
1.197 + { ret = GLP_EDATA;
1.198 + goto done;
1.199 + }
1.200 + cap[k] = (int)temp;
1.201 + if (a_cost >= 0)
1.202 + memcpy(&temp, (char *)a->data + a_cost, sizeof(double));
1.203 + else
1.204 + temp = 0.0;
1.205 + if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp)))
1.206 + { ret = GLP_EDATA;
1.207 + goto done;
1.208 + }
1.209 + cost[k] = (int)temp;
1.210 + }
1.211 + }
1.212 + /* (artificial arcs) */
1.213 + sum = 0.0;
1.214 + for (i = 1; i <= G->nv; i++)
1.215 + { v = G->v[i];
1.216 + if (v_rhs >= 0)
1.217 + memcpy(&temp, (char *)v->data + v_rhs, sizeof(double));
1.218 + else
1.219 + temp = 0.0;
1.220 + if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp)))
1.221 + { ret = GLP_EDATA;
1.222 + goto done;
1.223 + }
1.224 + if (temp > 0.0)
1.225 + { /* artificial arc from s to original source i */
1.226 + k++;
1.227 + tail[k] = s;
1.228 + head[k] = i;
1.229 + low[k] = cap[k] = (int)(+temp); /* supply */
1.230 + cost[k] = 0;
1.231 + sum += (double)temp;
1.232 + }
1.233 + else if (temp < 0.0)
1.234 + { /* artificial arc from original sink i to t */
1.235 + k++;
1.236 + tail[k] = i;
1.237 + head[k] = t;
1.238 + low[k] = cap[k] = (int)(-temp); /* demand */
1.239 + cost[k] = 0;
1.240 + }
1.241 + }
1.242 + /* (feedback arc from t to s) */
1.243 + k++;
1.244 + xassert(k == na);
1.245 + tail[k] = t;
1.246 + head[k] = s;
1.247 + if (sum > (double)INT_MAX)
1.248 + { ret = GLP_EDATA;
1.249 + goto done;
1.250 + }
1.251 + low[k] = cap[k] = (int)sum; /* total supply/demand */
1.252 + cost[k] = 0;
1.253 + /* find minimal-cost circulation in the resulting network */
1.254 + ret = okalg(nv, na, tail, head, low, cap, cost, x, pi);
1.255 + switch (ret)
1.256 + { case 0:
1.257 + /* optimal circulation found */
1.258 + ret = 0;
1.259 + break;
1.260 + case 1:
1.261 + /* no feasible circulation exists */
1.262 + ret = GLP_ENOPFS;
1.263 + break;
1.264 + case 2:
1.265 + /* integer overflow occured */
1.266 + ret = GLP_ERANGE;
1.267 + goto done;
1.268 + case 3:
1.269 + /* optimality test failed (logic error) */
1.270 + ret = GLP_EFAIL;
1.271 + goto done;
1.272 + default:
1.273 + xassert(ret != ret);
1.274 + }
1.275 + /* store solution components */
1.276 + /* (objective function = the total cost) */
1.277 + if (sol != NULL)
1.278 + { temp = 0.0;
1.279 + for (k = 1; k <= na; k++)
1.280 + temp += (double)cost[k] * (double)x[k];
1.281 + *sol = temp;
1.282 + }
1.283 + /* (arc flows) */
1.284 + if (a_x >= 0)
1.285 + { k = 0;
1.286 + for (i = 1; i <= G->nv; i++)
1.287 + { v = G->v[i];
1.288 + for (a = v->out; a != NULL; a = a->t_next)
1.289 + { temp = (double)x[++k];
1.290 + memcpy((char *)a->data + a_x, &temp, sizeof(double));
1.291 + }
1.292 + }
1.293 + }
1.294 + /* (node potentials = Lagrange multipliers) */
1.295 + if (v_pi >= 0)
1.296 + { for (i = 1; i <= G->nv; i++)
1.297 + { v = G->v[i];
1.298 + temp = - (double)pi[i];
1.299 + memcpy((char *)v->data + v_pi, &temp, sizeof(double));
1.300 + }
1.301 + }
1.302 +done: /* free working arrays */
1.303 + xfree(tail);
1.304 + xfree(head);
1.305 + xfree(low);
1.306 + xfree(cap);
1.307 + xfree(cost);
1.308 + xfree(x);
1.309 + xfree(pi);
1.310 + return ret;
1.311 +}
1.312 +
1.313 +/***********************************************************************
1.314 +* NAME
1.315 +*
1.316 +* glp_maxflow_lp - convert maximum flow problem to LP
1.317 +*
1.318 +* SYNOPSIS
1.319 +*
1.320 +* void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s,
1.321 +* int t, int a_cap);
1.322 +*
1.323 +* DESCRIPTION
1.324 +*
1.325 +* The routine glp_maxflow_lp builds an LP problem, which corresponds
1.326 +* to the maximum flow problem on the specified network G. */
1.327 +
1.328 +void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s,
1.329 + int t, int a_cap)
1.330 +{ glp_vertex *v;
1.331 + glp_arc *a;
1.332 + int i, j, type, ind[1+2];
1.333 + double cap, val[1+2];
1.334 + if (!(names == GLP_ON || names == GLP_OFF))
1.335 + xerror("glp_maxflow_lp: names = %d; invalid parameter\n",
1.336 + names);
1.337 + if (!(1 <= s && s <= G->nv))
1.338 + xerror("glp_maxflow_lp: s = %d; source node number out of rang"
1.339 + "e\n", s);
1.340 + if (!(1 <= t && t <= G->nv))
1.341 + xerror("glp_maxflow_lp: t = %d: sink node number out of range "
1.342 + "\n", t);
1.343 + if (s == t)
1.344 + xerror("glp_maxflow_lp: s = t = %d; source and sink nodes must"
1.345 + " be distinct\n", s);
1.346 + if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
1.347 + xerror("glp_maxflow_lp: a_cap = %d; invalid offset\n", a_cap);
1.348 + glp_erase_prob(lp);
1.349 + if (names) glp_set_prob_name(lp, G->name);
1.350 + glp_set_obj_dir(lp, GLP_MAX);
1.351 + glp_add_rows(lp, G->nv);
1.352 + for (i = 1; i <= G->nv; i++)
1.353 + { v = G->v[i];
1.354 + if (names) glp_set_row_name(lp, i, v->name);
1.355 + if (i == s)
1.356 + type = GLP_LO;
1.357 + else if (i == t)
1.358 + type = GLP_UP;
1.359 + else
1.360 + type = GLP_FX;
1.361 + glp_set_row_bnds(lp, i, type, 0.0, 0.0);
1.362 + }
1.363 + if (G->na > 0) glp_add_cols(lp, G->na);
1.364 + for (i = 1, j = 0; i <= G->nv; i++)
1.365 + { v = G->v[i];
1.366 + for (a = v->out; a != NULL; a = a->t_next)
1.367 + { j++;
1.368 + if (names)
1.369 + { char name[50+1];
1.370 + sprintf(name, "x[%d,%d]", a->tail->i, a->head->i);
1.371 + xassert(strlen(name) < sizeof(name));
1.372 + glp_set_col_name(lp, j, name);
1.373 + }
1.374 + if (a->tail->i != a->head->i)
1.375 + { ind[1] = a->tail->i, val[1] = +1.0;
1.376 + ind[2] = a->head->i, val[2] = -1.0;
1.377 + glp_set_mat_col(lp, j, 2, ind, val);
1.378 + }
1.379 + if (a_cap >= 0)
1.380 + memcpy(&cap, (char *)a->data + a_cap, sizeof(double));
1.381 + else
1.382 + cap = 1.0;
1.383 + if (cap == DBL_MAX)
1.384 + type = GLP_LO;
1.385 + else if (cap != 0.0)
1.386 + type = GLP_DB;
1.387 + else
1.388 + type = GLP_FX;
1.389 + glp_set_col_bnds(lp, j, type, 0.0, cap);
1.390 + if (a->tail->i == s)
1.391 + glp_set_obj_coef(lp, j, +1.0);
1.392 + else if (a->head->i == s)
1.393 + glp_set_obj_coef(lp, j, -1.0);
1.394 + }
1.395 + }
1.396 + xassert(j == G->na);
1.397 + return;
1.398 +}
1.399 +
1.400 +int glp_maxflow_ffalg(glp_graph *G, int s, int t, int a_cap,
1.401 + double *sol, int a_x, int v_cut)
1.402 +{ /* find maximal flow with Ford-Fulkerson algorithm */
1.403 + glp_vertex *v;
1.404 + glp_arc *a;
1.405 + int nv, na, i, k, flag, *tail, *head, *cap, *x, ret;
1.406 + char *cut;
1.407 + double temp;
1.408 + if (!(1 <= s && s <= G->nv))
1.409 + xerror("glp_maxflow_ffalg: s = %d; source node number out of r"
1.410 + "ange\n", s);
1.411 + if (!(1 <= t && t <= G->nv))
1.412 + xerror("glp_maxflow_ffalg: t = %d: sink node number out of ran"
1.413 + "ge\n", t);
1.414 + if (s == t)
1.415 + xerror("glp_maxflow_ffalg: s = t = %d; source and sink nodes m"
1.416 + "ust be distinct\n", s);
1.417 + if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
1.418 + xerror("glp_maxflow_ffalg: a_cap = %d; invalid offset\n",
1.419 + a_cap);
1.420 + if (v_cut >= 0 && v_cut > G->v_size - (int)sizeof(int))
1.421 + xerror("glp_maxflow_ffalg: v_cut = %d; invalid offset\n",
1.422 + v_cut);
1.423 + /* allocate working arrays */
1.424 + nv = G->nv;
1.425 + na = G->na;
1.426 + tail = xcalloc(1+na, sizeof(int));
1.427 + head = xcalloc(1+na, sizeof(int));
1.428 + cap = xcalloc(1+na, sizeof(int));
1.429 + x = xcalloc(1+na, sizeof(int));
1.430 + if (v_cut < 0)
1.431 + cut = NULL;
1.432 + else
1.433 + cut = xcalloc(1+nv, sizeof(char));
1.434 + /* copy the flow network */
1.435 + k = 0;
1.436 + for (i = 1; i <= G->nv; i++)
1.437 + { v = G->v[i];
1.438 + for (a = v->out; a != NULL; a = a->t_next)
1.439 + { k++;
1.440 + tail[k] = a->tail->i;
1.441 + head[k] = a->head->i;
1.442 + if (tail[k] == head[k])
1.443 + { ret = GLP_EDATA;
1.444 + goto done;
1.445 + }
1.446 + if (a_cap >= 0)
1.447 + memcpy(&temp, (char *)a->data + a_cap, sizeof(double));
1.448 + else
1.449 + temp = 1.0;
1.450 + if (!(0.0 <= temp && temp <= (double)INT_MAX &&
1.451 + temp == floor(temp)))
1.452 + { ret = GLP_EDATA;
1.453 + goto done;
1.454 + }
1.455 + cap[k] = (int)temp;
1.456 + }
1.457 + }
1.458 + xassert(k == na);
1.459 + /* find maximal flow in the flow network */
1.460 + ffalg(nv, na, tail, head, s, t, cap, x, cut);
1.461 + ret = 0;
1.462 + /* store solution components */
1.463 + /* (objective function = total flow through the network) */
1.464 + if (sol != NULL)
1.465 + { temp = 0.0;
1.466 + for (k = 1; k <= na; k++)
1.467 + { if (tail[k] == s)
1.468 + temp += (double)x[k];
1.469 + else if (head[k] == s)
1.470 + temp -= (double)x[k];
1.471 + }
1.472 + *sol = temp;
1.473 + }
1.474 + /* (arc flows) */
1.475 + if (a_x >= 0)
1.476 + { k = 0;
1.477 + for (i = 1; i <= G->nv; i++)
1.478 + { v = G->v[i];
1.479 + for (a = v->out; a != NULL; a = a->t_next)
1.480 + { temp = (double)x[++k];
1.481 + memcpy((char *)a->data + a_x, &temp, sizeof(double));
1.482 + }
1.483 + }
1.484 + }
1.485 + /* (node flags) */
1.486 + if (v_cut >= 0)
1.487 + { for (i = 1; i <= G->nv; i++)
1.488 + { v = G->v[i];
1.489 + flag = cut[i];
1.490 + memcpy((char *)v->data + v_cut, &flag, sizeof(int));
1.491 + }
1.492 + }
1.493 +done: /* free working arrays */
1.494 + xfree(tail);
1.495 + xfree(head);
1.496 + xfree(cap);
1.497 + xfree(x);
1.498 + if (cut != NULL) xfree(cut);
1.499 + return ret;
1.500 +}
1.501 +
1.502 +/***********************************************************************
1.503 +* NAME
1.504 +*
1.505 +* glp_check_asnprob - check correctness of assignment problem data
1.506 +*
1.507 +* SYNOPSIS
1.508 +*
1.509 +* int glp_check_asnprob(glp_graph *G, int v_set);
1.510 +*
1.511 +* RETURNS
1.512 +*
1.513 +* If the specified assignment problem data are correct, the routine
1.514 +* glp_check_asnprob returns zero, otherwise, non-zero. */
1.515 +
1.516 +int glp_check_asnprob(glp_graph *G, int v_set)
1.517 +{ glp_vertex *v;
1.518 + int i, k, ret = 0;
1.519 + if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
1.520 + xerror("glp_check_asnprob: v_set = %d; invalid offset\n",
1.521 + v_set);
1.522 + for (i = 1; i <= G->nv; i++)
1.523 + { v = G->v[i];
1.524 + if (v_set >= 0)
1.525 + { memcpy(&k, (char *)v->data + v_set, sizeof(int));
1.526 + if (k == 0)
1.527 + { if (v->in != NULL)
1.528 + { ret = 1;
1.529 + break;
1.530 + }
1.531 + }
1.532 + else if (k == 1)
1.533 + { if (v->out != NULL)
1.534 + { ret = 2;
1.535 + break;
1.536 + }
1.537 + }
1.538 + else
1.539 + { ret = 3;
1.540 + break;
1.541 + }
1.542 + }
1.543 + else
1.544 + { if (v->in != NULL && v->out != NULL)
1.545 + { ret = 4;
1.546 + break;
1.547 + }
1.548 + }
1.549 + }
1.550 + return ret;
1.551 +}
1.552 +
1.553 +/***********************************************************************
1.554 +* NAME
1.555 +*
1.556 +* glp_asnprob_lp - convert assignment problem to LP
1.557 +*
1.558 +* SYNOPSIS
1.559 +*
1.560 +* int glp_asnprob_lp(glp_prob *P, int form, glp_graph *G, int names,
1.561 +* int v_set, int a_cost);
1.562 +*
1.563 +* DESCRIPTION
1.564 +*
1.565 +* The routine glp_asnprob_lp builds an LP problem, which corresponds
1.566 +* to the assignment problem on the specified graph G.
1.567 +*
1.568 +* RETURNS
1.569 +*
1.570 +* If the LP problem has been successfully built, the routine returns
1.571 +* zero, otherwise, non-zero. */
1.572 +
1.573 +int glp_asnprob_lp(glp_prob *P, int form, glp_graph *G, int names,
1.574 + int v_set, int a_cost)
1.575 +{ glp_vertex *v;
1.576 + glp_arc *a;
1.577 + int i, j, ret, ind[1+2];
1.578 + double cost, val[1+2];
1.579 + if (!(form == GLP_ASN_MIN || form == GLP_ASN_MAX ||
1.580 + form == GLP_ASN_MMP))
1.581 + xerror("glp_asnprob_lp: form = %d; invalid parameter\n",
1.582 + form);
1.583 + if (!(names == GLP_ON || names == GLP_OFF))
1.584 + xerror("glp_asnprob_lp: names = %d; invalid parameter\n",
1.585 + names);
1.586 + if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
1.587 + xerror("glp_asnprob_lp: v_set = %d; invalid offset\n",
1.588 + v_set);
1.589 + if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
1.590 + xerror("glp_asnprob_lp: a_cost = %d; invalid offset\n",
1.591 + a_cost);
1.592 + ret = glp_check_asnprob(G, v_set);
1.593 + if (ret != 0) goto done;
1.594 + glp_erase_prob(P);
1.595 + if (names) glp_set_prob_name(P, G->name);
1.596 + glp_set_obj_dir(P, form == GLP_ASN_MIN ? GLP_MIN : GLP_MAX);
1.597 + if (G->nv > 0) glp_add_rows(P, G->nv);
1.598 + for (i = 1; i <= G->nv; i++)
1.599 + { v = G->v[i];
1.600 + if (names) glp_set_row_name(P, i, v->name);
1.601 + glp_set_row_bnds(P, i, form == GLP_ASN_MMP ? GLP_UP : GLP_FX,
1.602 + 1.0, 1.0);
1.603 + }
1.604 + if (G->na > 0) glp_add_cols(P, G->na);
1.605 + for (i = 1, j = 0; i <= G->nv; i++)
1.606 + { v = G->v[i];
1.607 + for (a = v->out; a != NULL; a = a->t_next)
1.608 + { j++;
1.609 + if (names)
1.610 + { char name[50+1];
1.611 + sprintf(name, "x[%d,%d]", a->tail->i, a->head->i);
1.612 + xassert(strlen(name) < sizeof(name));
1.613 + glp_set_col_name(P, j, name);
1.614 + }
1.615 + ind[1] = a->tail->i, val[1] = +1.0;
1.616 + ind[2] = a->head->i, val[2] = +1.0;
1.617 + glp_set_mat_col(P, j, 2, ind, val);
1.618 + glp_set_col_bnds(P, j, GLP_DB, 0.0, 1.0);
1.619 + if (a_cost >= 0)
1.620 + memcpy(&cost, (char *)a->data + a_cost, sizeof(double));
1.621 + else
1.622 + cost = 1.0;
1.623 + glp_set_obj_coef(P, j, cost);
1.624 + }
1.625 + }
1.626 + xassert(j == G->na);
1.627 +done: return ret;
1.628 +}
1.629 +
1.630 +/**********************************************************************/
1.631 +
1.632 +int glp_asnprob_okalg(int form, glp_graph *G, int v_set, int a_cost,
1.633 + double *sol, int a_x)
1.634 +{ /* solve assignment problem with out-of-kilter algorithm */
1.635 + glp_vertex *v;
1.636 + glp_arc *a;
1.637 + int nv, na, i, k, *tail, *head, *low, *cap, *cost, *x, *pi, ret;
1.638 + double temp;
1.639 + if (!(form == GLP_ASN_MIN || form == GLP_ASN_MAX ||
1.640 + form == GLP_ASN_MMP))
1.641 + xerror("glp_asnprob_okalg: form = %d; invalid parameter\n",
1.642 + form);
1.643 + if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
1.644 + xerror("glp_asnprob_okalg: v_set = %d; invalid offset\n",
1.645 + v_set);
1.646 + if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
1.647 + xerror("glp_asnprob_okalg: a_cost = %d; invalid offset\n",
1.648 + a_cost);
1.649 + if (a_x >= 0 && a_x > G->a_size - (int)sizeof(int))
1.650 + xerror("glp_asnprob_okalg: a_x = %d; invalid offset\n", a_x);
1.651 + if (glp_check_asnprob(G, v_set))
1.652 + return GLP_EDATA;
1.653 + /* nv is the total number of nodes in the resulting network */
1.654 + nv = G->nv + 1;
1.655 + /* na is the total number of arcs in the resulting network */
1.656 + na = G->na + G->nv;
1.657 + /* allocate working arrays */
1.658 + tail = xcalloc(1+na, sizeof(int));
1.659 + head = xcalloc(1+na, sizeof(int));
1.660 + low = xcalloc(1+na, sizeof(int));
1.661 + cap = xcalloc(1+na, sizeof(int));
1.662 + cost = xcalloc(1+na, sizeof(int));
1.663 + x = xcalloc(1+na, sizeof(int));
1.664 + pi = xcalloc(1+nv, sizeof(int));
1.665 + /* construct the resulting network */
1.666 + k = 0;
1.667 + /* (original arcs) */
1.668 + for (i = 1; i <= G->nv; i++)
1.669 + { v = G->v[i];
1.670 + for (a = v->out; a != NULL; a = a->t_next)
1.671 + { k++;
1.672 + tail[k] = a->tail->i;
1.673 + head[k] = a->head->i;
1.674 + low[k] = 0;
1.675 + cap[k] = 1;
1.676 + if (a_cost >= 0)
1.677 + memcpy(&temp, (char *)a->data + a_cost, sizeof(double));
1.678 + else
1.679 + temp = 1.0;
1.680 + if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp)))
1.681 + { ret = GLP_EDATA;
1.682 + goto done;
1.683 + }
1.684 + cost[k] = (int)temp;
1.685 + if (form != GLP_ASN_MIN) cost[k] = - cost[k];
1.686 + }
1.687 + }
1.688 + /* (artificial arcs) */
1.689 + for (i = 1; i <= G->nv; i++)
1.690 + { v = G->v[i];
1.691 + k++;
1.692 + if (v->out == NULL)
1.693 + tail[k] = i, head[k] = nv;
1.694 + else if (v->in == NULL)
1.695 + tail[k] = nv, head[k] = i;
1.696 + else
1.697 + xassert(v != v);
1.698 + low[k] = (form == GLP_ASN_MMP ? 0 : 1);
1.699 + cap[k] = 1;
1.700 + cost[k] = 0;
1.701 + }
1.702 + xassert(k == na);
1.703 + /* find minimal-cost circulation in the resulting network */
1.704 + ret = okalg(nv, na, tail, head, low, cap, cost, x, pi);
1.705 + switch (ret)
1.706 + { case 0:
1.707 + /* optimal circulation found */
1.708 + ret = 0;
1.709 + break;
1.710 + case 1:
1.711 + /* no feasible circulation exists */
1.712 + ret = GLP_ENOPFS;
1.713 + break;
1.714 + case 2:
1.715 + /* integer overflow occured */
1.716 + ret = GLP_ERANGE;
1.717 + goto done;
1.718 + case 3:
1.719 + /* optimality test failed (logic error) */
1.720 + ret = GLP_EFAIL;
1.721 + goto done;
1.722 + default:
1.723 + xassert(ret != ret);
1.724 + }
1.725 + /* store solution components */
1.726 + /* (objective function = the total cost) */
1.727 + if (sol != NULL)
1.728 + { temp = 0.0;
1.729 + for (k = 1; k <= na; k++)
1.730 + temp += (double)cost[k] * (double)x[k];
1.731 + if (form != GLP_ASN_MIN) temp = - temp;
1.732 + *sol = temp;
1.733 + }
1.734 + /* (arc flows) */
1.735 + if (a_x >= 0)
1.736 + { k = 0;
1.737 + for (i = 1; i <= G->nv; i++)
1.738 + { v = G->v[i];
1.739 + for (a = v->out; a != NULL; a = a->t_next)
1.740 + { k++;
1.741 + if (ret == 0)
1.742 + xassert(x[k] == 0 || x[k] == 1);
1.743 + memcpy((char *)a->data + a_x, &x[k], sizeof(int));
1.744 + }
1.745 + }
1.746 + }
1.747 +done: /* free working arrays */
1.748 + xfree(tail);
1.749 + xfree(head);
1.750 + xfree(low);
1.751 + xfree(cap);
1.752 + xfree(cost);
1.753 + xfree(x);
1.754 + xfree(pi);
1.755 + return ret;
1.756 +}
1.757 +
1.758 +/***********************************************************************
1.759 +* NAME
1.760 +*
1.761 +* glp_asnprob_hall - find bipartite matching of maximum cardinality
1.762 +*
1.763 +* SYNOPSIS
1.764 +*
1.765 +* int glp_asnprob_hall(glp_graph *G, int v_set, int a_x);
1.766 +*
1.767 +* DESCRIPTION
1.768 +*
1.769 +* The routine glp_asnprob_hall finds a matching of maximal cardinality
1.770 +* in the specified bipartite graph G. It uses a version of the Fortran
1.771 +* routine MC21A developed by I.S.Duff [1], which implements Hall's
1.772 +* algorithm [2].
1.773 +*
1.774 +* RETURNS
1.775 +*
1.776 +* The routine glp_asnprob_hall returns the cardinality of the matching
1.777 +* found. However, if the specified graph is incorrect (as detected by
1.778 +* the routine glp_check_asnprob), the routine returns negative value.
1.779 +*
1.780 +* REFERENCES
1.781 +*
1.782 +* 1. I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM
1.783 +* Trans. on Math. Softw. 7 (1981), 387-390.
1.784 +*
1.785 +* 2. M.Hall, "An Algorithm for distinct representatives," Amer. Math.
1.786 +* Monthly 63 (1956), 716-717. */
1.787 +
1.788 +int glp_asnprob_hall(glp_graph *G, int v_set, int a_x)
1.789 +{ glp_vertex *v;
1.790 + glp_arc *a;
1.791 + int card, i, k, loc, n, n1, n2, xij;
1.792 + int *num, *icn, *ip, *lenr, *iperm, *pr, *arp, *cv, *out;
1.793 + if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
1.794 + xerror("glp_asnprob_hall: v_set = %d; invalid offset\n",
1.795 + v_set);
1.796 + if (a_x >= 0 && a_x > G->a_size - (int)sizeof(int))
1.797 + xerror("glp_asnprob_hall: a_x = %d; invalid offset\n", a_x);
1.798 + if (glp_check_asnprob(G, v_set))
1.799 + return -1;
1.800 + /* determine the number of vertices in sets R and S and renumber
1.801 + vertices in S which correspond to columns of the matrix; skip
1.802 + all isolated vertices */
1.803 + num = xcalloc(1+G->nv, sizeof(int));
1.804 + n1 = n2 = 0;
1.805 + for (i = 1; i <= G->nv; i++)
1.806 + { v = G->v[i];
1.807 + if (v->in == NULL && v->out != NULL)
1.808 + n1++, num[i] = 0; /* vertex in R */
1.809 + else if (v->in != NULL && v->out == NULL)
1.810 + n2++, num[i] = n2; /* vertex in S */
1.811 + else
1.812 + { xassert(v->in == NULL && v->out == NULL);
1.813 + num[i] = -1; /* isolated vertex */
1.814 + }
1.815 + }
1.816 + /* the matrix must be square, thus, if it has more columns than
1.817 + rows, extra rows will be just empty, and vice versa */
1.818 + n = (n1 >= n2 ? n1 : n2);
1.819 + /* allocate working arrays */
1.820 + icn = xcalloc(1+G->na, sizeof(int));
1.821 + ip = xcalloc(1+n, sizeof(int));
1.822 + lenr = xcalloc(1+n, sizeof(int));
1.823 + iperm = xcalloc(1+n, sizeof(int));
1.824 + pr = xcalloc(1+n, sizeof(int));
1.825 + arp = xcalloc(1+n, sizeof(int));
1.826 + cv = xcalloc(1+n, sizeof(int));
1.827 + out = xcalloc(1+n, sizeof(int));
1.828 + /* build the adjacency matrix of the bipartite graph in row-wise
1.829 + format (rows are vertices in R, columns are vertices in S) */
1.830 + k = 0, loc = 1;
1.831 + for (i = 1; i <= G->nv; i++)
1.832 + { if (num[i] != 0) continue;
1.833 + /* vertex i in R */
1.834 + ip[++k] = loc;
1.835 + v = G->v[i];
1.836 + for (a = v->out; a != NULL; a = a->t_next)
1.837 + { xassert(num[a->head->i] != 0);
1.838 + icn[loc++] = num[a->head->i];
1.839 + }
1.840 + lenr[k] = loc - ip[k];
1.841 + }
1.842 + xassert(loc-1 == G->na);
1.843 + /* make all extra rows empty (all extra columns are empty due to
1.844 + the row-wise format used) */
1.845 + for (k++; k <= n; k++)
1.846 + ip[k] = loc, lenr[k] = 0;
1.847 + /* find a row permutation that maximizes the number of non-zeros
1.848 + on the main diagonal */
1.849 + card = mc21a(n, icn, ip, lenr, iperm, pr, arp, cv, out);
1.850 +#if 1 /* 18/II-2010 */
1.851 + /* FIXED: if card = n, arp remains clobbered on exit */
1.852 + for (i = 1; i <= n; i++)
1.853 + arp[i] = 0;
1.854 + for (i = 1; i <= card; i++)
1.855 + { k = iperm[i];
1.856 + xassert(1 <= k && k <= n);
1.857 + xassert(arp[k] == 0);
1.858 + arp[k] = i;
1.859 + }
1.860 +#endif
1.861 + /* store solution, if necessary */
1.862 + if (a_x < 0) goto skip;
1.863 + k = 0;
1.864 + for (i = 1; i <= G->nv; i++)
1.865 + { if (num[i] != 0) continue;
1.866 + /* vertex i in R */
1.867 + k++;
1.868 + v = G->v[i];
1.869 + for (a = v->out; a != NULL; a = a->t_next)
1.870 + { /* arp[k] is the number of matched column or zero */
1.871 + if (arp[k] == num[a->head->i])
1.872 + { xassert(arp[k] != 0);
1.873 + xij = 1;
1.874 + }
1.875 + else
1.876 + xij = 0;
1.877 + memcpy((char *)a->data + a_x, &xij, sizeof(int));
1.878 + }
1.879 + }
1.880 +skip: /* free working arrays */
1.881 + xfree(num);
1.882 + xfree(icn);
1.883 + xfree(ip);
1.884 + xfree(lenr);
1.885 + xfree(iperm);
1.886 + xfree(pr);
1.887 + xfree(arp);
1.888 + xfree(cv);
1.889 + xfree(out);
1.890 + return card;
1.891 +}
1.892 +
1.893 +/***********************************************************************
1.894 +* NAME
1.895 +*
1.896 +* glp_cpp - solve critical path problem
1.897 +*
1.898 +* SYNOPSIS
1.899 +*
1.900 +* double glp_cpp(glp_graph *G, int v_t, int v_es, int v_ls);
1.901 +*
1.902 +* DESCRIPTION
1.903 +*
1.904 +* The routine glp_cpp solves the critical path problem represented in
1.905 +* the form of the project network.
1.906 +*
1.907 +* The parameter G is a pointer to the graph object, which specifies
1.908 +* the project network. This graph must be acyclic. Multiple arcs are
1.909 +* allowed being considered as single arcs.
1.910 +*
1.911 +* The parameter v_t specifies an offset of the field of type double
1.912 +* in the vertex data block, which contains time t[i] >= 0 needed to
1.913 +* perform corresponding job j. If v_t < 0, it is assumed that t[i] = 1
1.914 +* for all jobs.
1.915 +*
1.916 +* The parameter v_es specifies an offset of the field of type double
1.917 +* in the vertex data block, to which the routine stores earliest start
1.918 +* time for corresponding job. If v_es < 0, this time is not stored.
1.919 +*
1.920 +* The parameter v_ls specifies an offset of the field of type double
1.921 +* in the vertex data block, to which the routine stores latest start
1.922 +* time for corresponding job. If v_ls < 0, this time is not stored.
1.923 +*
1.924 +* RETURNS
1.925 +*
1.926 +* The routine glp_cpp returns the minimal project duration, that is,
1.927 +* minimal time needed to perform all jobs in the project. */
1.928 +
1.929 +static void sorting(glp_graph *G, int list[]);
1.930 +
1.931 +double glp_cpp(glp_graph *G, int v_t, int v_es, int v_ls)
1.932 +{ glp_vertex *v;
1.933 + glp_arc *a;
1.934 + int i, j, k, nv, *list;
1.935 + double temp, total, *t, *es, *ls;
1.936 + if (v_t >= 0 && v_t > G->v_size - (int)sizeof(double))
1.937 + xerror("glp_cpp: v_t = %d; invalid offset\n", v_t);
1.938 + if (v_es >= 0 && v_es > G->v_size - (int)sizeof(double))
1.939 + xerror("glp_cpp: v_es = %d; invalid offset\n", v_es);
1.940 + if (v_ls >= 0 && v_ls > G->v_size - (int)sizeof(double))
1.941 + xerror("glp_cpp: v_ls = %d; invalid offset\n", v_ls);
1.942 + nv = G->nv;
1.943 + if (nv == 0)
1.944 + { total = 0.0;
1.945 + goto done;
1.946 + }
1.947 + /* allocate working arrays */
1.948 + t = xcalloc(1+nv, sizeof(double));
1.949 + es = xcalloc(1+nv, sizeof(double));
1.950 + ls = xcalloc(1+nv, sizeof(double));
1.951 + list = xcalloc(1+nv, sizeof(int));
1.952 + /* retrieve job times */
1.953 + for (i = 1; i <= nv; i++)
1.954 + { v = G->v[i];
1.955 + if (v_t >= 0)
1.956 + { memcpy(&t[i], (char *)v->data + v_t, sizeof(double));
1.957 + if (t[i] < 0.0)
1.958 + xerror("glp_cpp: t[%d] = %g; invalid time\n", i, t[i]);
1.959 + }
1.960 + else
1.961 + t[i] = 1.0;
1.962 + }
1.963 + /* perform topological sorting to determine the list of nodes
1.964 + (jobs) such that if list[k] = i and list[kk] = j and there
1.965 + exists arc (i->j), then k < kk */
1.966 + sorting(G, list);
1.967 + /* FORWARD PASS */
1.968 + /* determine earliest start times */
1.969 + for (k = 1; k <= nv; k++)
1.970 + { j = list[k];
1.971 + es[j] = 0.0;
1.972 + for (a = G->v[j]->in; a != NULL; a = a->h_next)
1.973 + { i = a->tail->i;
1.974 + /* there exists arc (i->j) in the project network */
1.975 + temp = es[i] + t[i];
1.976 + if (es[j] < temp) es[j] = temp;
1.977 + }
1.978 + }
1.979 + /* determine the minimal project duration */
1.980 + total = 0.0;
1.981 + for (i = 1; i <= nv; i++)
1.982 + { temp = es[i] + t[i];
1.983 + if (total < temp) total = temp;
1.984 + }
1.985 + /* BACKWARD PASS */
1.986 + /* determine latest start times */
1.987 + for (k = nv; k >= 1; k--)
1.988 + { i = list[k];
1.989 + ls[i] = total - t[i];
1.990 + for (a = G->v[i]->out; a != NULL; a = a->t_next)
1.991 + { j = a->head->i;
1.992 + /* there exists arc (i->j) in the project network */
1.993 + temp = ls[j] - t[i];
1.994 + if (ls[i] > temp) ls[i] = temp;
1.995 + }
1.996 + /* avoid possible round-off errors */
1.997 + if (ls[i] < es[i]) ls[i] = es[i];
1.998 + }
1.999 + /* store results, if necessary */
1.1000 + if (v_es >= 0)
1.1001 + { for (i = 1; i <= nv; i++)
1.1002 + { v = G->v[i];
1.1003 + memcpy((char *)v->data + v_es, &es[i], sizeof(double));
1.1004 + }
1.1005 + }
1.1006 + if (v_ls >= 0)
1.1007 + { for (i = 1; i <= nv; i++)
1.1008 + { v = G->v[i];
1.1009 + memcpy((char *)v->data + v_ls, &ls[i], sizeof(double));
1.1010 + }
1.1011 + }
1.1012 + /* free working arrays */
1.1013 + xfree(t);
1.1014 + xfree(es);
1.1015 + xfree(ls);
1.1016 + xfree(list);
1.1017 +done: return total;
1.1018 +}
1.1019 +
1.1020 +static void sorting(glp_graph *G, int list[])
1.1021 +{ /* perform topological sorting to determine the list of nodes
1.1022 + (jobs) such that if list[k] = i and list[kk] = j and there
1.1023 + exists arc (i->j), then k < kk */
1.1024 + int i, k, nv, v_size, *num;
1.1025 + void **save;
1.1026 + nv = G->nv;
1.1027 + v_size = G->v_size;
1.1028 + save = xcalloc(1+nv, sizeof(void *));
1.1029 + num = xcalloc(1+nv, sizeof(int));
1.1030 + G->v_size = sizeof(int);
1.1031 + for (i = 1; i <= nv; i++)
1.1032 + { save[i] = G->v[i]->data;
1.1033 + G->v[i]->data = &num[i];
1.1034 + list[i] = 0;
1.1035 + }
1.1036 + if (glp_top_sort(G, 0) != 0)
1.1037 + xerror("glp_cpp: project network is not acyclic\n");
1.1038 + G->v_size = v_size;
1.1039 + for (i = 1; i <= nv; i++)
1.1040 + { G->v[i]->data = save[i];
1.1041 + k = num[i];
1.1042 + xassert(1 <= k && k <= nv);
1.1043 + xassert(list[k] == 0);
1.1044 + list[k] = i;
1.1045 + }
1.1046 + xfree(save);
1.1047 + xfree(num);
1.1048 + return;
1.1049 +}
1.1050 +
1.1051 +/* eof */