COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpapi17.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

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