alpar@1: /* glpapi17.c (flow network problems) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpapi.h"
alpar@1: #include "glpnet.h"
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_mincost_lp - convert minimum cost flow problem to LP
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  void glp_mincost_lp(glp_prob *lp, glp_graph *G, int names,
alpar@1: *     int v_rhs, int a_low, int a_cap, int a_cost);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_mincost_lp builds an LP problem, which corresponds
alpar@1: *  to the minimum cost flow problem on the specified network G. */
alpar@1: 
alpar@1: void glp_mincost_lp(glp_prob *lp, glp_graph *G, int names, int v_rhs,
alpar@1:       int a_low, int a_cap, int a_cost)
alpar@1: {     glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int i, j, type, ind[1+2];
alpar@1:       double rhs, low, cap, cost, val[1+2];
alpar@1:       if (!(names == GLP_ON || names == GLP_OFF))
alpar@1:          xerror("glp_mincost_lp: names = %d; invalid parameter\n",
alpar@1:             names);
alpar@1:       if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_lp: v_rhs = %d; invalid offset\n", v_rhs);
alpar@1:       if (a_low >= 0 && a_low > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_lp: a_low = %d; invalid offset\n", a_low);
alpar@1:       if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_lp: a_cap = %d; invalid offset\n", a_cap);
alpar@1:       if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_lp: a_cost = %d; invalid offset\n", a_cost)
alpar@1:             ;
alpar@1:       glp_erase_prob(lp);
alpar@1:       if (names) glp_set_prob_name(lp, G->name);
alpar@1:       if (G->nv > 0) glp_add_rows(lp, G->nv);
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (names) glp_set_row_name(lp, i, v->name);
alpar@1:          if (v_rhs >= 0)
alpar@1:             memcpy(&rhs, (char *)v->data + v_rhs, sizeof(double));
alpar@1:          else
alpar@1:             rhs = 0.0;
alpar@1:          glp_set_row_bnds(lp, i, GLP_FX, rhs, rhs);
alpar@1:       }
alpar@1:       if (G->na > 0) glp_add_cols(lp, G->na);
alpar@1:       for (i = 1, j = 0; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  j++;
alpar@1:             if (names)
alpar@1:             {  char name[50+1];
alpar@1:                sprintf(name, "x[%d,%d]", a->tail->i, a->head->i);
alpar@1:                xassert(strlen(name) < sizeof(name));
alpar@1:                glp_set_col_name(lp, j, name);
alpar@1:             }
alpar@1:             if (a->tail->i != a->head->i)
alpar@1:             {  ind[1] = a->tail->i, val[1] = +1.0;
alpar@1:                ind[2] = a->head->i, val[2] = -1.0;
alpar@1:                glp_set_mat_col(lp, j, 2, ind, val);
alpar@1:             }
alpar@1:             if (a_low >= 0)
alpar@1:                memcpy(&low, (char *)a->data + a_low, sizeof(double));
alpar@1:             else
alpar@1:                low = 0.0;
alpar@1:             if (a_cap >= 0)
alpar@1:                memcpy(&cap, (char *)a->data + a_cap, sizeof(double));
alpar@1:             else
alpar@1:                cap = 1.0;
alpar@1:             if (cap == DBL_MAX)
alpar@1:                type = GLP_LO;
alpar@1:             else if (low != cap)
alpar@1:                type = GLP_DB;
alpar@1:             else
alpar@1:                type = GLP_FX;
alpar@1:             glp_set_col_bnds(lp, j, type, low, cap);
alpar@1:             if (a_cost >= 0)
alpar@1:                memcpy(&cost, (char *)a->data + a_cost, sizeof(double));
alpar@1:             else
alpar@1:                cost = 0.0;
alpar@1:             glp_set_obj_coef(lp, j, cost);
alpar@1:          }
alpar@1:       }
alpar@1:       xassert(j == G->na);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /**********************************************************************/
alpar@1: 
alpar@1: int glp_mincost_okalg(glp_graph *G, int v_rhs, int a_low, int a_cap,
alpar@1:       int a_cost, double *sol, int a_x, int v_pi)
alpar@1: {     /* find minimum-cost flow with out-of-kilter algorithm */
alpar@1:       glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int nv, na, i, k, s, t, *tail, *head, *low, *cap, *cost, *x, *pi,
alpar@1:          ret;
alpar@1:       double sum, temp;
alpar@1:       if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_okalg: v_rhs = %d; invalid offset\n",
alpar@1:             v_rhs);
alpar@1:       if (a_low >= 0 && a_low > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_okalg: a_low = %d; invalid offset\n",
alpar@1:             a_low);
alpar@1:       if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_okalg: a_cap = %d; invalid offset\n",
alpar@1:             a_cap);
alpar@1:       if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_okalg: a_cost = %d; invalid offset\n",
alpar@1:             a_cost);
alpar@1:       if (a_x >= 0 && a_x > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_okalg: a_x = %d; invalid offset\n", a_x);
alpar@1:       if (v_pi >= 0 && v_pi > G->v_size - (int)sizeof(double))
alpar@1:          xerror("glp_mincost_okalg: v_pi = %d; invalid offset\n", v_pi);
alpar@1:       /* s is artificial source node */
alpar@1:       s = G->nv + 1;
alpar@1:       /* t is artificial sink node */
alpar@1:       t = s + 1;
alpar@1:       /* nv is the total number of nodes in the resulting network */
alpar@1:       nv = t;
alpar@1:       /* na is the total number of arcs in the resulting network */
alpar@1:       na = G->na + 1;
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (v_rhs >= 0)
alpar@1:             memcpy(&temp, (char *)v->data + v_rhs, sizeof(double));
alpar@1:          else
alpar@1:             temp = 0.0;
alpar@1:          if (temp != 0.0) na++;
alpar@1:       }
alpar@1:       /* allocate working arrays */
alpar@1:       tail = xcalloc(1+na, sizeof(int));
alpar@1:       head = xcalloc(1+na, sizeof(int));
alpar@1:       low = xcalloc(1+na, sizeof(int));
alpar@1:       cap = xcalloc(1+na, sizeof(int));
alpar@1:       cost = xcalloc(1+na, sizeof(int));
alpar@1:       x = xcalloc(1+na, sizeof(int));
alpar@1:       pi = xcalloc(1+nv, sizeof(int));
alpar@1:       /* construct the resulting network */
alpar@1:       k = 0;
alpar@1:       /* (original arcs) */
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  k++;
alpar@1:             tail[k] = a->tail->i;
alpar@1:             head[k] = a->head->i;
alpar@1:             if (tail[k] == head[k])
alpar@1:             {  ret = GLP_EDATA;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             if (a_low >= 0)
alpar@1:                memcpy(&temp, (char *)a->data + a_low, sizeof(double));
alpar@1:             else
alpar@1:                temp = 0.0;
alpar@1:             if (!(0.0 <= temp && temp <= (double)INT_MAX &&
alpar@1:                   temp == floor(temp)))
alpar@1:             {  ret = GLP_EDATA;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             low[k] = (int)temp;
alpar@1:             if (a_cap >= 0)
alpar@1:                memcpy(&temp, (char *)a->data + a_cap, sizeof(double));
alpar@1:             else
alpar@1:                temp = 1.0;
alpar@1:             if (!((double)low[k] <= temp && temp <= (double)INT_MAX &&
alpar@1:                   temp == floor(temp)))
alpar@1:             {  ret = GLP_EDATA;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             cap[k] = (int)temp;
alpar@1:             if (a_cost >= 0)
alpar@1:                memcpy(&temp, (char *)a->data + a_cost, sizeof(double));
alpar@1:             else
alpar@1:                temp = 0.0;
alpar@1:             if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp)))
alpar@1:             {  ret = GLP_EDATA;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             cost[k] = (int)temp;
alpar@1:          }
alpar@1:       }
alpar@1:       /* (artificial arcs) */
alpar@1:       sum = 0.0;
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (v_rhs >= 0)
alpar@1:             memcpy(&temp, (char *)v->data + v_rhs, sizeof(double));
alpar@1:          else
alpar@1:             temp = 0.0;
alpar@1:          if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp)))
alpar@1:          {  ret = GLP_EDATA;
alpar@1:             goto done;
alpar@1:          }
alpar@1:          if (temp > 0.0)
alpar@1:          {  /* artificial arc from s to original source i */
alpar@1:             k++;
alpar@1:             tail[k] = s;
alpar@1:             head[k] = i;
alpar@1:             low[k] = cap[k] = (int)(+temp); /* supply */
alpar@1:             cost[k] = 0;
alpar@1:             sum += (double)temp;
alpar@1:          }
alpar@1:          else if (temp < 0.0)
alpar@1:          {  /* artificial arc from original sink i to t */
alpar@1:             k++;
alpar@1:             tail[k] = i;
alpar@1:             head[k] = t;
alpar@1:             low[k] = cap[k] = (int)(-temp); /* demand */
alpar@1:             cost[k] = 0;
alpar@1:          }
alpar@1:       }
alpar@1:       /* (feedback arc from t to s) */
alpar@1:       k++;
alpar@1:       xassert(k == na);
alpar@1:       tail[k] = t;
alpar@1:       head[k] = s;
alpar@1:       if (sum > (double)INT_MAX)
alpar@1:       {  ret = GLP_EDATA;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       low[k] = cap[k] = (int)sum; /* total supply/demand */
alpar@1:       cost[k] = 0;
alpar@1:       /* find minimal-cost circulation in the resulting network */
alpar@1:       ret = okalg(nv, na, tail, head, low, cap, cost, x, pi);
alpar@1:       switch (ret)
alpar@1:       {  case 0:
alpar@1:             /* optimal circulation found */
alpar@1:             ret = 0;
alpar@1:             break;
alpar@1:          case 1:
alpar@1:             /* no feasible circulation exists */
alpar@1:             ret = GLP_ENOPFS;
alpar@1:             break;
alpar@1:          case 2:
alpar@1:             /* integer overflow occured */
alpar@1:             ret = GLP_ERANGE;
alpar@1:             goto done;
alpar@1:          case 3:
alpar@1:             /* optimality test failed (logic error) */
alpar@1:             ret = GLP_EFAIL;
alpar@1:             goto done;
alpar@1:          default:
alpar@1:             xassert(ret != ret);
alpar@1:       }
alpar@1:       /* store solution components */
alpar@1:       /* (objective function = the total cost) */
alpar@1:       if (sol != NULL)
alpar@1:       {  temp = 0.0;
alpar@1:          for (k = 1; k <= na; k++)
alpar@1:             temp += (double)cost[k] * (double)x[k];
alpar@1:          *sol = temp;
alpar@1:       }
alpar@1:       /* (arc flows) */
alpar@1:       if (a_x >= 0)
alpar@1:       {  k = 0;
alpar@1:          for (i = 1; i <= G->nv; i++)
alpar@1:          {  v = G->v[i];
alpar@1:             for (a = v->out; a != NULL; a = a->t_next)
alpar@1:             {  temp = (double)x[++k];
alpar@1:                memcpy((char *)a->data + a_x, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* (node potentials = Lagrange multipliers) */
alpar@1:       if (v_pi >= 0)
alpar@1:       {  for (i = 1; i <= G->nv; i++)
alpar@1:          {  v = G->v[i];
alpar@1:             temp = - (double)pi[i];
alpar@1:             memcpy((char *)v->data + v_pi, &temp, sizeof(double));
alpar@1:          }
alpar@1:       }
alpar@1: done: /* free working arrays */
alpar@1:       xfree(tail);
alpar@1:       xfree(head);
alpar@1:       xfree(low);
alpar@1:       xfree(cap);
alpar@1:       xfree(cost);
alpar@1:       xfree(x);
alpar@1:       xfree(pi);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_maxflow_lp - convert maximum flow problem to LP
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s,
alpar@1: *     int t, int a_cap);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_maxflow_lp builds an LP problem, which corresponds
alpar@1: *  to the maximum flow problem on the specified network G. */
alpar@1: 
alpar@1: void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s,
alpar@1:       int t, int a_cap)
alpar@1: {     glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int i, j, type, ind[1+2];
alpar@1:       double cap, val[1+2];
alpar@1:       if (!(names == GLP_ON || names == GLP_OFF))
alpar@1:          xerror("glp_maxflow_lp: names = %d; invalid parameter\n",
alpar@1:             names);
alpar@1:       if (!(1 <= s && s <= G->nv))
alpar@1:          xerror("glp_maxflow_lp: s = %d; source node number out of rang"
alpar@1:             "e\n", s);
alpar@1:       if (!(1 <= t && t <= G->nv))
alpar@1:          xerror("glp_maxflow_lp: t = %d: sink node number out of range "
alpar@1:             "\n", t);
alpar@1:       if (s == t)
alpar@1:          xerror("glp_maxflow_lp: s = t = %d; source and sink nodes must"
alpar@1:             " be distinct\n", s);
alpar@1:       if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_maxflow_lp: a_cap = %d; invalid offset\n", a_cap);
alpar@1:       glp_erase_prob(lp);
alpar@1:       if (names) glp_set_prob_name(lp, G->name);
alpar@1:       glp_set_obj_dir(lp, GLP_MAX);
alpar@1:       glp_add_rows(lp, G->nv);
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (names) glp_set_row_name(lp, i, v->name);
alpar@1:          if (i == s)
alpar@1:             type = GLP_LO;
alpar@1:          else if (i == t)
alpar@1:             type = GLP_UP;
alpar@1:          else
alpar@1:             type = GLP_FX;
alpar@1:          glp_set_row_bnds(lp, i, type, 0.0, 0.0);
alpar@1:       }
alpar@1:       if (G->na > 0) glp_add_cols(lp, G->na);
alpar@1:       for (i = 1, j = 0; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  j++;
alpar@1:             if (names)
alpar@1:             {  char name[50+1];
alpar@1:                sprintf(name, "x[%d,%d]", a->tail->i, a->head->i);
alpar@1:                xassert(strlen(name) < sizeof(name));
alpar@1:                glp_set_col_name(lp, j, name);
alpar@1:             }
alpar@1:             if (a->tail->i != a->head->i)
alpar@1:             {  ind[1] = a->tail->i, val[1] = +1.0;
alpar@1:                ind[2] = a->head->i, val[2] = -1.0;
alpar@1:                glp_set_mat_col(lp, j, 2, ind, val);
alpar@1:             }
alpar@1:             if (a_cap >= 0)
alpar@1:                memcpy(&cap, (char *)a->data + a_cap, sizeof(double));
alpar@1:             else
alpar@1:                cap = 1.0;
alpar@1:             if (cap == DBL_MAX)
alpar@1:                type = GLP_LO;
alpar@1:             else if (cap != 0.0)
alpar@1:                type = GLP_DB;
alpar@1:             else
alpar@1:                type = GLP_FX;
alpar@1:             glp_set_col_bnds(lp, j, type, 0.0, cap);
alpar@1:             if (a->tail->i == s)
alpar@1:                glp_set_obj_coef(lp, j, +1.0);
alpar@1:             else if (a->head->i == s)
alpar@1:                glp_set_obj_coef(lp, j, -1.0);
alpar@1:          }
alpar@1:       }
alpar@1:       xassert(j == G->na);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: int glp_maxflow_ffalg(glp_graph *G, int s, int t, int a_cap,
alpar@1:       double *sol, int a_x, int v_cut)
alpar@1: {     /* find maximal flow with Ford-Fulkerson algorithm */
alpar@1:       glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int nv, na, i, k, flag, *tail, *head, *cap, *x, ret;
alpar@1:       char *cut;
alpar@1:       double temp;
alpar@1:       if (!(1 <= s && s <= G->nv))
alpar@1:          xerror("glp_maxflow_ffalg: s = %d; source node number out of r"
alpar@1:             "ange\n", s);
alpar@1:       if (!(1 <= t && t <= G->nv))
alpar@1:          xerror("glp_maxflow_ffalg: t = %d: sink node number out of ran"
alpar@1:             "ge\n", t);
alpar@1:       if (s == t)
alpar@1:          xerror("glp_maxflow_ffalg: s = t = %d; source and sink nodes m"
alpar@1:             "ust be distinct\n", s);
alpar@1:       if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_maxflow_ffalg: a_cap = %d; invalid offset\n",
alpar@1:             a_cap);
alpar@1:       if (v_cut >= 0 && v_cut > G->v_size - (int)sizeof(int))
alpar@1:          xerror("glp_maxflow_ffalg: v_cut = %d; invalid offset\n",
alpar@1:             v_cut);
alpar@1:       /* allocate working arrays */
alpar@1:       nv = G->nv;
alpar@1:       na = G->na;
alpar@1:       tail = xcalloc(1+na, sizeof(int));
alpar@1:       head = xcalloc(1+na, sizeof(int));
alpar@1:       cap = xcalloc(1+na, sizeof(int));
alpar@1:       x = xcalloc(1+na, sizeof(int));
alpar@1:       if (v_cut < 0)
alpar@1:          cut = NULL;
alpar@1:       else
alpar@1:          cut = xcalloc(1+nv, sizeof(char));
alpar@1:       /* copy the flow network */
alpar@1:       k = 0;
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  k++;
alpar@1:             tail[k] = a->tail->i;
alpar@1:             head[k] = a->head->i;
alpar@1:             if (tail[k] == head[k])
alpar@1:             {  ret = GLP_EDATA;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             if (a_cap >= 0)
alpar@1:                memcpy(&temp, (char *)a->data + a_cap, sizeof(double));
alpar@1:             else
alpar@1:                temp = 1.0;
alpar@1:             if (!(0.0 <= temp && temp <= (double)INT_MAX &&
alpar@1:                   temp == floor(temp)))
alpar@1:             {  ret = GLP_EDATA;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             cap[k] = (int)temp;
alpar@1:          }
alpar@1:       }
alpar@1:       xassert(k == na);
alpar@1:       /* find maximal flow in the flow network */
alpar@1:       ffalg(nv, na, tail, head, s, t, cap, x, cut);
alpar@1:       ret = 0;
alpar@1:       /* store solution components */
alpar@1:       /* (objective function = total flow through the network) */
alpar@1:       if (sol != NULL)
alpar@1:       {  temp = 0.0;
alpar@1:          for (k = 1; k <= na; k++)
alpar@1:          {  if (tail[k] == s)
alpar@1:                temp += (double)x[k];
alpar@1:             else if (head[k] == s)
alpar@1:                temp -= (double)x[k];
alpar@1:          }
alpar@1:          *sol = temp;
alpar@1:       }
alpar@1:       /* (arc flows) */
alpar@1:       if (a_x >= 0)
alpar@1:       {  k = 0;
alpar@1:          for (i = 1; i <= G->nv; i++)
alpar@1:          {  v = G->v[i];
alpar@1:             for (a = v->out; a != NULL; a = a->t_next)
alpar@1:             {  temp = (double)x[++k];
alpar@1:                memcpy((char *)a->data + a_x, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* (node flags) */
alpar@1:       if (v_cut >= 0)
alpar@1:       {  for (i = 1; i <= G->nv; i++)
alpar@1:          {  v = G->v[i];
alpar@1:             flag = cut[i];
alpar@1:             memcpy((char *)v->data + v_cut, &flag, sizeof(int));
alpar@1:          }
alpar@1:       }
alpar@1: done: /* free working arrays */
alpar@1:       xfree(tail);
alpar@1:       xfree(head);
alpar@1:       xfree(cap);
alpar@1:       xfree(x);
alpar@1:       if (cut != NULL) xfree(cut);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_check_asnprob - check correctness of assignment problem data
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_check_asnprob(glp_graph *G, int v_set);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  If the specified assignment problem data are correct, the routine
alpar@1: *  glp_check_asnprob returns zero, otherwise, non-zero. */
alpar@1: 
alpar@1: int glp_check_asnprob(glp_graph *G, int v_set)
alpar@1: {     glp_vertex *v;
alpar@1:       int i, k, ret = 0;
alpar@1:       if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
alpar@1:          xerror("glp_check_asnprob: v_set = %d; invalid offset\n",
alpar@1:             v_set);
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (v_set >= 0)
alpar@1:          {  memcpy(&k, (char *)v->data + v_set, sizeof(int));
alpar@1:             if (k == 0)
alpar@1:             {  if (v->in != NULL)
alpar@1:                {  ret = 1;
alpar@1:                   break;
alpar@1:                }
alpar@1:             }
alpar@1:             else if (k == 1)
alpar@1:             {  if (v->out != NULL)
alpar@1:                {  ret = 2;
alpar@1:                   break;
alpar@1:                }
alpar@1:             }
alpar@1:             else
alpar@1:             {  ret = 3;
alpar@1:                break;
alpar@1:             }
alpar@1:          }
alpar@1:          else
alpar@1:          {  if (v->in != NULL && v->out != NULL)
alpar@1:             {  ret = 4;
alpar@1:                break;
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_asnprob_lp - convert assignment problem to LP
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_asnprob_lp(glp_prob *P, int form, glp_graph *G, int names,
alpar@1: *     int v_set, int a_cost);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_asnprob_lp builds an LP problem, which corresponds
alpar@1: *  to the assignment problem on the specified graph G.
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  If the LP problem has been successfully built, the routine returns
alpar@1: *  zero, otherwise, non-zero. */
alpar@1: 
alpar@1: int glp_asnprob_lp(glp_prob *P, int form, glp_graph *G, int names,
alpar@1:       int v_set, int a_cost)
alpar@1: {     glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int i, j, ret, ind[1+2];
alpar@1:       double cost, val[1+2];
alpar@1:       if (!(form == GLP_ASN_MIN || form == GLP_ASN_MAX ||
alpar@1:             form == GLP_ASN_MMP))
alpar@1:          xerror("glp_asnprob_lp: form = %d; invalid parameter\n",
alpar@1:             form);
alpar@1:       if (!(names == GLP_ON || names == GLP_OFF))
alpar@1:          xerror("glp_asnprob_lp: names = %d; invalid parameter\n",
alpar@1:             names);
alpar@1:       if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
alpar@1:          xerror("glp_asnprob_lp: v_set = %d; invalid offset\n",
alpar@1:             v_set);
alpar@1:       if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_asnprob_lp: a_cost = %d; invalid offset\n",
alpar@1:             a_cost);
alpar@1:       ret = glp_check_asnprob(G, v_set);
alpar@1:       if (ret != 0) goto done;
alpar@1:       glp_erase_prob(P);
alpar@1:       if (names) glp_set_prob_name(P, G->name);
alpar@1:       glp_set_obj_dir(P, form == GLP_ASN_MIN ? GLP_MIN : GLP_MAX);
alpar@1:       if (G->nv > 0) glp_add_rows(P, G->nv);
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (names) glp_set_row_name(P, i, v->name);
alpar@1:          glp_set_row_bnds(P, i, form == GLP_ASN_MMP ? GLP_UP : GLP_FX,
alpar@1:             1.0, 1.0);
alpar@1:       }
alpar@1:       if (G->na > 0) glp_add_cols(P, G->na);
alpar@1:       for (i = 1, j = 0; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  j++;
alpar@1:             if (names)
alpar@1:             {  char name[50+1];
alpar@1:                sprintf(name, "x[%d,%d]", a->tail->i, a->head->i);
alpar@1:                xassert(strlen(name) < sizeof(name));
alpar@1:                glp_set_col_name(P, j, name);
alpar@1:             }
alpar@1:             ind[1] = a->tail->i, val[1] = +1.0;
alpar@1:             ind[2] = a->head->i, val[2] = +1.0;
alpar@1:             glp_set_mat_col(P, j, 2, ind, val);
alpar@1:             glp_set_col_bnds(P, j, GLP_DB, 0.0, 1.0);
alpar@1:             if (a_cost >= 0)
alpar@1:                memcpy(&cost, (char *)a->data + a_cost, sizeof(double));
alpar@1:             else
alpar@1:                cost = 1.0;
alpar@1:             glp_set_obj_coef(P, j, cost);
alpar@1:          }
alpar@1:       }
alpar@1:       xassert(j == G->na);
alpar@1: done: return ret;
alpar@1: }
alpar@1: 
alpar@1: /**********************************************************************/
alpar@1: 
alpar@1: int glp_asnprob_okalg(int form, glp_graph *G, int v_set, int a_cost,
alpar@1:       double *sol, int a_x)
alpar@1: {     /* solve assignment problem with out-of-kilter algorithm */
alpar@1:       glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int nv, na, i, k, *tail, *head, *low, *cap, *cost, *x, *pi, ret;
alpar@1:       double temp;
alpar@1:       if (!(form == GLP_ASN_MIN || form == GLP_ASN_MAX ||
alpar@1:             form == GLP_ASN_MMP))
alpar@1:          xerror("glp_asnprob_okalg: form = %d; invalid parameter\n",
alpar@1:             form);
alpar@1:       if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
alpar@1:          xerror("glp_asnprob_okalg: v_set = %d; invalid offset\n",
alpar@1:             v_set);
alpar@1:       if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
alpar@1:          xerror("glp_asnprob_okalg: a_cost = %d; invalid offset\n",
alpar@1:             a_cost);
alpar@1:       if (a_x >= 0 && a_x > G->a_size - (int)sizeof(int))
alpar@1:          xerror("glp_asnprob_okalg: a_x = %d; invalid offset\n", a_x);
alpar@1:       if (glp_check_asnprob(G, v_set))
alpar@1:          return GLP_EDATA;
alpar@1:       /* nv is the total number of nodes in the resulting network */
alpar@1:       nv = G->nv + 1;
alpar@1:       /* na is the total number of arcs in the resulting network */
alpar@1:       na = G->na + G->nv;
alpar@1:       /* allocate working arrays */
alpar@1:       tail = xcalloc(1+na, sizeof(int));
alpar@1:       head = xcalloc(1+na, sizeof(int));
alpar@1:       low = xcalloc(1+na, sizeof(int));
alpar@1:       cap = xcalloc(1+na, sizeof(int));
alpar@1:       cost = xcalloc(1+na, sizeof(int));
alpar@1:       x = xcalloc(1+na, sizeof(int));
alpar@1:       pi = xcalloc(1+nv, sizeof(int));
alpar@1:       /* construct the resulting network */
alpar@1:       k = 0;
alpar@1:       /* (original arcs) */
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  k++;
alpar@1:             tail[k] = a->tail->i;
alpar@1:             head[k] = a->head->i;
alpar@1:             low[k] = 0;
alpar@1:             cap[k] = 1;
alpar@1:             if (a_cost >= 0)
alpar@1:                memcpy(&temp, (char *)a->data + a_cost, sizeof(double));
alpar@1:             else
alpar@1:                temp = 1.0;
alpar@1:             if (!(fabs(temp) <= (double)INT_MAX && temp == floor(temp)))
alpar@1:             {  ret = GLP_EDATA;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             cost[k] = (int)temp;
alpar@1:             if (form != GLP_ASN_MIN) cost[k] = - cost[k];
alpar@1:          }
alpar@1:       }
alpar@1:       /* (artificial arcs) */
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          k++;
alpar@1:          if (v->out == NULL)
alpar@1:             tail[k] = i, head[k] = nv;
alpar@1:          else if (v->in == NULL)
alpar@1:             tail[k] = nv, head[k] = i;
alpar@1:          else
alpar@1:             xassert(v != v);
alpar@1:          low[k] = (form == GLP_ASN_MMP ? 0 : 1);
alpar@1:          cap[k] = 1;
alpar@1:          cost[k] = 0;
alpar@1:       }
alpar@1:       xassert(k == na);
alpar@1:       /* find minimal-cost circulation in the resulting network */
alpar@1:       ret = okalg(nv, na, tail, head, low, cap, cost, x, pi);
alpar@1:       switch (ret)
alpar@1:       {  case 0:
alpar@1:             /* optimal circulation found */
alpar@1:             ret = 0;
alpar@1:             break;
alpar@1:          case 1:
alpar@1:             /* no feasible circulation exists */
alpar@1:             ret = GLP_ENOPFS;
alpar@1:             break;
alpar@1:          case 2:
alpar@1:             /* integer overflow occured */
alpar@1:             ret = GLP_ERANGE;
alpar@1:             goto done;
alpar@1:          case 3:
alpar@1:             /* optimality test failed (logic error) */
alpar@1:             ret = GLP_EFAIL;
alpar@1:             goto done;
alpar@1:          default:
alpar@1:             xassert(ret != ret);
alpar@1:       }
alpar@1:       /* store solution components */
alpar@1:       /* (objective function = the total cost) */
alpar@1:       if (sol != NULL)
alpar@1:       {  temp = 0.0;
alpar@1:          for (k = 1; k <= na; k++)
alpar@1:             temp += (double)cost[k] * (double)x[k];
alpar@1:          if (form != GLP_ASN_MIN) temp = - temp;
alpar@1:          *sol = temp;
alpar@1:       }
alpar@1:       /* (arc flows) */
alpar@1:       if (a_x >= 0)
alpar@1:       {  k = 0;
alpar@1:          for (i = 1; i <= G->nv; i++)
alpar@1:          {  v = G->v[i];
alpar@1:             for (a = v->out; a != NULL; a = a->t_next)
alpar@1:             {  k++;
alpar@1:                if (ret == 0)
alpar@1:                   xassert(x[k] == 0 || x[k] == 1);
alpar@1:                memcpy((char *)a->data + a_x, &x[k], sizeof(int));
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1: done: /* free working arrays */
alpar@1:       xfree(tail);
alpar@1:       xfree(head);
alpar@1:       xfree(low);
alpar@1:       xfree(cap);
alpar@1:       xfree(cost);
alpar@1:       xfree(x);
alpar@1:       xfree(pi);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_asnprob_hall - find bipartite matching of maximum cardinality
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_asnprob_hall(glp_graph *G, int v_set, int a_x);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_asnprob_hall finds a matching of maximal cardinality
alpar@1: *  in the specified bipartite graph G. It uses a version of the Fortran
alpar@1: *  routine MC21A developed by I.S.Duff [1], which implements Hall's
alpar@1: *  algorithm [2].
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_asnprob_hall returns the cardinality of the matching
alpar@1: *  found. However, if the specified graph is incorrect (as detected by
alpar@1: *  the routine glp_check_asnprob), the routine returns negative value.
alpar@1: *
alpar@1: *  REFERENCES
alpar@1: *
alpar@1: *  1. I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM
alpar@1: *     Trans. on Math. Softw. 7 (1981), 387-390.
alpar@1: *
alpar@1: *  2. M.Hall, "An Algorithm for distinct representatives," Amer. Math.
alpar@1: *     Monthly 63 (1956), 716-717. */
alpar@1: 
alpar@1: int glp_asnprob_hall(glp_graph *G, int v_set, int a_x)
alpar@1: {     glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int card, i, k, loc, n, n1, n2, xij;
alpar@1:       int *num, *icn, *ip, *lenr, *iperm, *pr, *arp, *cv, *out;
alpar@1:       if (v_set >= 0 && v_set > G->v_size - (int)sizeof(int))
alpar@1:          xerror("glp_asnprob_hall: v_set = %d; invalid offset\n",
alpar@1:             v_set);
alpar@1:       if (a_x >= 0 && a_x > G->a_size - (int)sizeof(int))
alpar@1:          xerror("glp_asnprob_hall: a_x = %d; invalid offset\n", a_x);
alpar@1:       if (glp_check_asnprob(G, v_set))
alpar@1:          return -1;
alpar@1:       /* determine the number of vertices in sets R and S and renumber
alpar@1:          vertices in S which correspond to columns of the matrix; skip
alpar@1:          all isolated vertices */
alpar@1:       num = xcalloc(1+G->nv, sizeof(int));
alpar@1:       n1 = n2 = 0;
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (v->in == NULL && v->out != NULL)
alpar@1:             n1++, num[i] = 0; /* vertex in R */
alpar@1:          else if (v->in != NULL && v->out == NULL)
alpar@1:             n2++, num[i] = n2; /* vertex in S */
alpar@1:          else
alpar@1:          {  xassert(v->in == NULL && v->out == NULL);
alpar@1:             num[i] = -1; /* isolated vertex */
alpar@1:          }
alpar@1:       }
alpar@1:       /* the matrix must be square, thus, if it has more columns than
alpar@1:          rows, extra rows will be just empty, and vice versa */
alpar@1:       n = (n1 >= n2 ? n1 : n2);
alpar@1:       /* allocate working arrays */
alpar@1:       icn = xcalloc(1+G->na, sizeof(int));
alpar@1:       ip = xcalloc(1+n, sizeof(int));
alpar@1:       lenr = xcalloc(1+n, sizeof(int));
alpar@1:       iperm = xcalloc(1+n, sizeof(int));
alpar@1:       pr = xcalloc(1+n, sizeof(int));
alpar@1:       arp = xcalloc(1+n, sizeof(int));
alpar@1:       cv = xcalloc(1+n, sizeof(int));
alpar@1:       out = xcalloc(1+n, sizeof(int));
alpar@1:       /* build the adjacency matrix of the bipartite graph in row-wise
alpar@1:          format (rows are vertices in R, columns are vertices in S) */
alpar@1:       k = 0, loc = 1;
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  if (num[i] != 0) continue;
alpar@1:          /* vertex i in R */
alpar@1:          ip[++k] = loc;
alpar@1:          v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  xassert(num[a->head->i] != 0);
alpar@1:             icn[loc++] = num[a->head->i];
alpar@1:          }
alpar@1:          lenr[k] = loc - ip[k];
alpar@1:       }
alpar@1:       xassert(loc-1 == G->na);
alpar@1:       /* make all extra rows empty (all extra columns are empty due to
alpar@1:          the row-wise format used) */
alpar@1:       for (k++; k <= n; k++)
alpar@1:          ip[k] = loc, lenr[k] = 0;
alpar@1:       /* find a row permutation that maximizes the number of non-zeros
alpar@1:          on the main diagonal */
alpar@1:       card = mc21a(n, icn, ip, lenr, iperm, pr, arp, cv, out);
alpar@1: #if 1 /* 18/II-2010 */
alpar@1:       /* FIXED: if card = n, arp remains clobbered on exit */
alpar@1:       for (i = 1; i <= n; i++)
alpar@1:          arp[i] = 0;
alpar@1:       for (i = 1; i <= card; i++)
alpar@1:       {  k = iperm[i];
alpar@1:          xassert(1 <= k && k <= n);
alpar@1:          xassert(arp[k] == 0);
alpar@1:          arp[k] = i;
alpar@1:       }
alpar@1: #endif
alpar@1:       /* store solution, if necessary */
alpar@1:       if (a_x < 0) goto skip;
alpar@1:       k = 0;
alpar@1:       for (i = 1; i <= G->nv; i++)
alpar@1:       {  if (num[i] != 0) continue;
alpar@1:          /* vertex i in R */
alpar@1:          k++;
alpar@1:          v = G->v[i];
alpar@1:          for (a = v->out; a != NULL; a = a->t_next)
alpar@1:          {  /* arp[k] is the number of matched column or zero */
alpar@1:             if (arp[k] == num[a->head->i])
alpar@1:             {  xassert(arp[k] != 0);
alpar@1:                xij = 1;
alpar@1:             }
alpar@1:             else
alpar@1:                xij = 0;
alpar@1:             memcpy((char *)a->data + a_x, &xij, sizeof(int));
alpar@1:          }
alpar@1:       }
alpar@1: skip: /* free working arrays */
alpar@1:       xfree(num);
alpar@1:       xfree(icn);
alpar@1:       xfree(ip);
alpar@1:       xfree(lenr);
alpar@1:       xfree(iperm);
alpar@1:       xfree(pr);
alpar@1:       xfree(arp);
alpar@1:       xfree(cv);
alpar@1:       xfree(out);
alpar@1:       return card;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_cpp - solve critical path problem
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  double glp_cpp(glp_graph *G, int v_t, int v_es, int v_ls);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_cpp solves the critical path problem represented in
alpar@1: *  the form of the project network.
alpar@1: *
alpar@1: *  The parameter G is a pointer to the graph object, which specifies
alpar@1: *  the project network. This graph must be acyclic. Multiple arcs are
alpar@1: *  allowed being considered as single arcs.
alpar@1: *
alpar@1: *  The parameter v_t specifies an offset of the field of type double
alpar@1: *  in the vertex data block, which contains time t[i] >= 0 needed to
alpar@1: *  perform corresponding job j. If v_t < 0, it is assumed that t[i] = 1
alpar@1: *  for all jobs.
alpar@1: *
alpar@1: *  The parameter v_es specifies an offset of the field of type double
alpar@1: *  in the vertex data block, to which the routine stores earliest start
alpar@1: *  time for corresponding job. If v_es < 0, this time is not stored.
alpar@1: *
alpar@1: *  The parameter v_ls specifies an offset of the field of type double
alpar@1: *  in the vertex data block, to which the routine stores latest start
alpar@1: *  time for corresponding job. If v_ls < 0, this time is not stored.
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_cpp returns the minimal project duration, that is,
alpar@1: *  minimal time needed to perform all jobs in the project. */
alpar@1: 
alpar@1: static void sorting(glp_graph *G, int list[]);
alpar@1: 
alpar@1: double glp_cpp(glp_graph *G, int v_t, int v_es, int v_ls)
alpar@1: {     glp_vertex *v;
alpar@1:       glp_arc *a;
alpar@1:       int i, j, k, nv, *list;
alpar@1:       double temp, total, *t, *es, *ls;
alpar@1:       if (v_t >= 0 && v_t > G->v_size - (int)sizeof(double))
alpar@1:          xerror("glp_cpp: v_t = %d; invalid offset\n", v_t);
alpar@1:       if (v_es >= 0 && v_es > G->v_size - (int)sizeof(double))
alpar@1:          xerror("glp_cpp: v_es = %d; invalid offset\n", v_es);
alpar@1:       if (v_ls >= 0 && v_ls > G->v_size - (int)sizeof(double))
alpar@1:          xerror("glp_cpp: v_ls = %d; invalid offset\n", v_ls);
alpar@1:       nv = G->nv;
alpar@1:       if (nv == 0)
alpar@1:       {  total = 0.0;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* allocate working arrays */
alpar@1:       t = xcalloc(1+nv, sizeof(double));
alpar@1:       es = xcalloc(1+nv, sizeof(double));
alpar@1:       ls = xcalloc(1+nv, sizeof(double));
alpar@1:       list = xcalloc(1+nv, sizeof(int));
alpar@1:       /* retrieve job times */
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:       {  v = G->v[i];
alpar@1:          if (v_t >= 0)
alpar@1:          {  memcpy(&t[i], (char *)v->data + v_t, sizeof(double));
alpar@1:             if (t[i] < 0.0)
alpar@1:                xerror("glp_cpp: t[%d] = %g; invalid time\n", i, t[i]);
alpar@1:          }
alpar@1:          else
alpar@1:             t[i] = 1.0;
alpar@1:       }
alpar@1:       /* perform topological sorting to determine the list of nodes
alpar@1:          (jobs) such that if list[k] = i and list[kk] = j and there
alpar@1:          exists arc (i->j), then k < kk */
alpar@1:       sorting(G, list);
alpar@1:       /* FORWARD PASS */
alpar@1:       /* determine earliest start times */
alpar@1:       for (k = 1; k <= nv; k++)
alpar@1:       {  j = list[k];
alpar@1:          es[j] = 0.0;
alpar@1:          for (a = G->v[j]->in; a != NULL; a = a->h_next)
alpar@1:          {  i = a->tail->i;
alpar@1:             /* there exists arc (i->j) in the project network */
alpar@1:             temp = es[i] + t[i];
alpar@1:             if (es[j] < temp) es[j] = temp;
alpar@1:          }
alpar@1:       }
alpar@1:       /* determine the minimal project duration */
alpar@1:       total = 0.0;
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:       {  temp = es[i] + t[i];
alpar@1:          if (total < temp) total = temp;
alpar@1:       }
alpar@1:       /* BACKWARD PASS */
alpar@1:       /* determine latest start times */
alpar@1:       for (k = nv; k >= 1; k--)
alpar@1:       {  i = list[k];
alpar@1:          ls[i] = total - t[i];
alpar@1:          for (a = G->v[i]->out; a != NULL; a = a->t_next)
alpar@1:          {  j = a->head->i;
alpar@1:             /* there exists arc (i->j) in the project network */
alpar@1:             temp = ls[j] - t[i];
alpar@1:             if (ls[i] > temp) ls[i] = temp;
alpar@1:          }
alpar@1:          /* avoid possible round-off errors */
alpar@1:          if (ls[i] < es[i]) ls[i] = es[i];
alpar@1:       }
alpar@1:       /* store results, if necessary */
alpar@1:       if (v_es >= 0)
alpar@1:       {  for (i = 1; i <= nv; i++)
alpar@1:          {  v = G->v[i];
alpar@1:             memcpy((char *)v->data + v_es, &es[i], sizeof(double));
alpar@1:          }
alpar@1:       }
alpar@1:       if (v_ls >= 0)
alpar@1:       {  for (i = 1; i <= nv; i++)
alpar@1:          {  v = G->v[i];
alpar@1:             memcpy((char *)v->data + v_ls, &ls[i], sizeof(double));
alpar@1:          }
alpar@1:       }
alpar@1:       /* free working arrays */
alpar@1:       xfree(t);
alpar@1:       xfree(es);
alpar@1:       xfree(ls);
alpar@1:       xfree(list);
alpar@1: done: return total;
alpar@1: }
alpar@1: 
alpar@1: static void sorting(glp_graph *G, int list[])
alpar@1: {     /* perform topological sorting to determine the list of nodes
alpar@1:          (jobs) such that if list[k] = i and list[kk] = j and there
alpar@1:          exists arc (i->j), then k < kk */
alpar@1:       int i, k, nv, v_size, *num;
alpar@1:       void **save;
alpar@1:       nv = G->nv;
alpar@1:       v_size = G->v_size;
alpar@1:       save = xcalloc(1+nv, sizeof(void *));
alpar@1:       num = xcalloc(1+nv, sizeof(int));
alpar@1:       G->v_size = sizeof(int);
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:       {  save[i] = G->v[i]->data;
alpar@1:          G->v[i]->data = &num[i];
alpar@1:          list[i] = 0;
alpar@1:       }
alpar@1:       if (glp_top_sort(G, 0) != 0)
alpar@1:          xerror("glp_cpp: project network is not acyclic\n");
alpar@1:       G->v_size = v_size;
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:       {  G->v[i]->data = save[i];
alpar@1:          k = num[i];
alpar@1:          xassert(1 <= k && k <= nv);
alpar@1:          xassert(list[k] == 0);
alpar@1:          list[k] = i;
alpar@1:       }
alpar@1:       xfree(save);
alpar@1:       xfree(num);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /* eof */