lemon-project-template-glpk

annotate deps/glpk/src/glpapi17.c @ 11:4fc6ad2fb8a6

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