1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpapi17.c	Sun Nov 06 20:59:10 2011 +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, 2011 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 */