1 /* glpnet03.c (Klingman's network problem generator) */
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
6 * This code is the result of translation of the Fortran program NETGEN
7 * developed by Dr. Darwin Klingman, which is publically available from
8 * NETLIB at <http://www.netlib.org/lp/generators>.
10 * The translation was made by Andrew Makhorin <mao@gnu.org>.
12 * GLPK is free software: you can redistribute it and/or modify it
13 * under the terms of the GNU General Public License as published by
14 * the Free Software Foundation, either version 3 of the License, or
15 * (at your option) any later version.
17 * GLPK is distributed in the hope that it will be useful, but WITHOUT
18 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
19 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
20 * License for more details.
22 * You should have received a copy of the GNU General Public License
23 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
24 ***********************************************************************/
28 /***********************************************************************
31 * glp_netgen - Klingman's network problem generator
35 * int glp_netgen(glp_graph *G, int v_rhs, int a_cap, int a_cost,
36 * const int parm[1+15]);
40 * The routine glp_netgen is a network problem generator developed by
41 * Dr. Darwin Klingman. It can create capacitated and uncapacitated
42 * minimum cost flow (or transshipment), transportation, and assignment
45 * The parameter G specifies the graph object, to which the generated
46 * problem data have to be stored. Note that on entry the graph object
47 * is erased with the routine glp_erase_graph.
49 * The parameter v_rhs specifies an offset of the field of type double
50 * in the vertex data block, to which the routine stores the supply or
51 * demand value. If v_rhs < 0, the value is not stored.
53 * The parameter a_cap specifies an offset of the field of type double
54 * in the arc data block, to which the routine stores the arc capacity.
55 * If a_cap < 0, the capacity is not stored.
57 * The parameter a_cost specifies an offset of the field of type double
58 * in the arc data block, to which the routine stores the per-unit cost
59 * if the arc flow. If a_cost < 0, the cost is not stored.
61 * The array parm contains description of the network to be generated:
64 * parm[1] (iseed) 8-digit positive random number seed
65 * parm[2] (nprob) 8-digit problem id number
66 * parm[3] (nodes) total number of nodes
67 * parm[4] (nsorc) total number of source nodes (including
68 * transshipment nodes)
69 * parm[5] (nsink) total number of sink nodes (including
70 * transshipment nodes)
71 * parm[6] (iarcs) number of arcs
72 * parm[7] (mincst) minimum cost for arcs
73 * parm[8] (maxcst) maximum cost for arcs
74 * parm[9] (itsup) total supply
75 * parm[10] (ntsorc) number of transshipment source nodes
76 * parm[11] (ntsink) number of transshipment sink nodes
77 * parm[12] (iphic) percentage of skeleton arcs to be given
79 * parm[13] (ipcap) percentage of arcs to be capacitated
80 * parm[14] (mincap) minimum upper bound for capacitated arcs
81 * parm[15] (maxcap) maximum upper bound for capacitated arcs
83 * The routine generates a transportation problem if:
85 * nsorc + nsink = nodes, ntsorc = 0, and ntsink = 0.
87 * The routine generates an assignment problem if the requirements for
88 * a transportation problem are met and:
90 * nsorc = nsink and itsup = nsorc.
94 * If the instance was successfully generated, the routine glp_netgen
95 * returns zero; otherwise, if specified parameters are inconsistent,
96 * the routine returns a non-zero error code.
100 * D.Klingman, A.Napier, and J.Stutz. NETGEN: A program for generating
101 * large scale capacitated assignment, transportation, and minimum cost
102 * flow networks. Management Science 20 (1974), 814-20. */
105 { /* common storage area */
107 int v_rhs, a_cap, a_cost;
108 int nodes, iarcs, mincst, maxcst, itsup, nsorc, nsink, nonsor,
109 nfsink, narcs, nsort, nftsor, ipcap, mincap, maxcap, ktl,
110 nodlft, *ipred, *ihead, *itail, *iflag, *isup, *lsinks, mult,
111 modul, i15, i16, jran;
115 #define v_rhs (csa->v_rhs)
116 #define a_cap (csa->a_cap)
117 #define a_cost (csa->a_cost)
118 #define nodes (csa->nodes)
119 #define iarcs (csa->iarcs)
120 #define mincst (csa->mincst)
121 #define maxcst (csa->maxcst)
122 #define itsup (csa->itsup)
123 #define nsorc (csa->nsorc)
124 #define nsink (csa->nsink)
125 #define nonsor (csa->nonsor)
126 #define nfsink (csa->nfsink)
127 #define narcs (csa->narcs)
128 #define nsort (csa->nsort)
129 #define nftsor (csa->nftsor)
130 #define ipcap (csa->ipcap)
131 #define mincap (csa->mincap)
132 #define maxcap (csa->maxcap)
133 #define ktl (csa->ktl)
134 #define nodlft (csa->nodlft)
136 /* spent a day to find out this bug */
137 #define ist (csa->ist)
139 #define ist (ipred[0])
141 #define ipred (csa->ipred)
142 #define ihead (csa->ihead)
143 #define itail (csa->itail)
144 #define iflag (csa->iflag)
145 #define isup (csa->isup)
146 #define lsinks (csa->lsinks)
147 #define mult (csa->mult)
148 #define modul (csa->modul)
149 #define i15 (csa->i15)
150 #define i16 (csa->i16)
151 #define jran (csa->jran)
153 static void cresup(struct csa *csa);
154 static void chain(struct csa *csa, int lpick, int lsorc);
155 static void chnarc(struct csa *csa, int lsorc);
156 static void sort(struct csa *csa);
157 static void pickj(struct csa *csa, int it);
158 static void assign(struct csa *csa);
159 static void setran(struct csa *csa, int iseed);
160 static int iran(struct csa *csa, int ilow, int ihigh);
162 int glp_netgen(glp_graph *G_, int _v_rhs, int _a_cap, int _a_cost,
163 const int parm[1+15])
164 { struct csa _csa, *csa = &_csa;
165 int iseed, nprob, ntsorc, ntsink, iphic, i, nskel, nltr, ltsink,
166 ntrans, npsink, nftr, npsorc, ntravl, ntrrem, lsorc, lpick,
167 nsksr, nsrchn, j, item, l, ks, k, ksp, li, n, ii, it, ih, icap,
168 jcap, icost, jcost, ret;
174 { if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
175 xerror("glp_netgen: v_rhs = %d; invalid offset\n", v_rhs);
176 if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
177 xerror("glp_netgen: a_cap = %d; invalid offset\n", a_cap);
178 if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
179 xerror("glp_netgen: a_cost = %d; invalid offset\n", a_cost);
181 /* Input the user's random number seed and fix it if
185 if (iseed <= 0) iseed = 13502460;
187 /* Input the user's problem characteristics. */
201 /* Check the size of the problem. */
202 if (!(10 <= nodes && nodes <= 100000))
206 /* Check user supplied parameters for consistency. */
207 if (!(nsorc >= 0 && nsink >= 0 && nsorc + nsink <= nodes))
223 if (!(0 <= ntsorc && ntsorc <= nsorc))
227 if (!(0 <= ntsink && ntsink <= nsink))
231 if (!(0 <= iphic && iphic <= 100))
235 if (!(0 <= ipcap && ipcap <= 100))
243 /* Initailize the graph object. */
245 { glp_erase_graph(G, G->v_size, G->a_size);
246 glp_add_vertices(G, nodes);
249 for (i = 1; i <= nodes; i++)
250 { glp_vertex *v = G->v[i];
251 memcpy((char *)v->data + v_rhs, &zero, sizeof(double));
255 /* Allocate working arrays. */
256 ipred = xcalloc(1+nodes, sizeof(int));
257 ihead = xcalloc(1+nodes, sizeof(int));
258 itail = xcalloc(1+nodes, sizeof(int));
259 iflag = xcalloc(1+nodes, sizeof(int));
260 isup = xcalloc(1+nodes, sizeof(int));
261 lsinks = xcalloc(1+nodes, sizeof(int));
262 /* Print the problem documentation records. */
264 { xprintf("BEGIN\n");
265 xprintf("NETGEN PROBLEM%8d%10s%10d NODES AND%10d ARCS\n",
266 nprob, "", nodes, iarcs);
267 xprintf("USER:%11d%11d%11d%11d%11d%11d\nDATA:%11d%11d%11d%11d%"
268 "11d%11d\n", iseed, nsorc, nsink, mincst,
269 maxcst, itsup, ntsorc, ntsink, iphic, ipcap,
273 glp_set_graph_name(G, "NETGEN");
274 /* Set various constants used in the program. */
277 nltr = nodes - nsink;
278 ltsink = nltr + ntsink;
279 ntrans = nltr - nsorc;
281 nonsor = nodes - nsorc + ntsorc;
282 npsink = nsink - ntsink;
283 nodlft = nodes - nsink + ntsink;
285 nftsor = nsorc - ntsorc + 1;
286 npsorc = nsorc - ntsorc;
287 /* Randomly distribute the supply among the source nodes. */
288 if (npsorc + npsink == nodes && npsorc == npsink &&
295 /* Print the supply records. */
297 { xprintf("SUPPLY\n");
298 for (i = 1; i <= nsorc; i++)
299 xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
304 { for (i = 1; i <= nsorc; i++)
305 { double temp = (double)isup[i];
306 glp_vertex *v = G->v[i];
307 memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
311 /* Make the sources point to themselves in ipred array. */
312 for (i = 1; i <= nsorc; i++)
314 if (ntrans == 0) goto L170;
315 /* Chain the transshipment nodes together in the ipred array. */
318 for (i = nftr; i < nltr; i++)
320 /* Form even length chains for 60 percent of the transshipments.*/
321 ntravl = 6 * ntrans / 10;
322 ntrrem = ntrans - ntravl;
325 { lpick = iran(csa, 1, ntravl + ntrrem);
327 chain(csa, lpick, lsorc);
328 if (lsorc == nsorc) goto L140;
331 /* Add the remaining transshipments to the chains. */
334 lpick = iran(csa, 1, ntrrem);
336 lsorc = iran(csa, 1, nsorc);
337 chain(csa, lpick, lsorc);
339 L170: /* Set all demands equal to zero. */
340 for (i = nfsink; i <= nodes; i++)
342 /* The following loop takes one chain at a time (through the use
343 of logic contained in the loop and calls to other routines) and
344 creates the remaining network arcs. */
345 for (lsorc = 1; lsorc <= nsorc; lsorc++)
346 { chnarc(csa, lsorc);
347 for (i = nfsink; i <= nodes; i++)
349 /* Choose the number of sinks to be hooked up to the current
352 nsksr = (nsort * 2 * nsink) / ntrans;
354 nsksr = nsink / nsorc + 1;
355 if (nsksr < 2) nsksr = 2;
356 if (nsksr > nsink) nsksr = nsink;
358 /* Randomly pick nsksr sinks and put their names in lsinks. */
360 for (j = 1; j <= nsksr; j++)
361 { item = iran(csa, 1, ktl);
363 for (l = nfsink; l <= nodes; l++)
366 if (item == 0) goto L230;
373 /* If last source chain, add all sinks with zero demand to
376 { for (j = nfsink; j <= nodes; j++)
377 { if (ipred[j] == 0 && iflag[j] != 1)
384 /* Create demands for group of sinks in lsinks. */
385 ks = isup[lsorc] / nsksr;
387 for (i = 1; i <= nsksr; i++)
389 ksp = iran(csa, 1, ks);
390 j = iran(csa, 1, nsksr);
396 ipred[li] += ks - ksp;
397 n = iran(csa, 1, nsrchn);
399 for (ii = 1; ii <= n; ii++)
403 ipred[li] += isup[lsorc] - ks * nsksr;
405 /* Sort the arcs in the chain from source lsorc using itail as
408 /* Print this part of skeleton and create the arcs for these
412 L300: for (j = nftsor; j <= nodes; j++)
421 /* Determine if this skeleton arc should be capacitated. */
423 jcap = iran(csa, 1, 100);
425 { icap = isup[lsorc];
426 if (mincap > icap) icap = mincap;
428 /* Determine if this skeleton arc should have the maximum
431 jcost = iran(csa, 1, 100);
433 icost = iran(csa, mincst, maxcst);
435 xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ih, "", icost,
438 { glp_arc *a = glp_add_arc(G, it, ih);
440 { double temp = (double)icap;
441 memcpy((char *)a->data + a_cap, &temp, sizeof(double));
444 { double temp = (double)icost;
445 memcpy((char *)a->data + a_cost, &temp, sizeof(double));
449 if (itail[i] == it) goto L320;
451 if (i <= nsort) goto L300;
453 /* Create arcs from the transshipment sinks. */
455 { for (i = nfsink; i <= ltsink; i++)
456 { for (j = nftsor; j <= nodes; j++)
463 L390: /* Print the demand records and end record. */
465 { xprintf("DEMAND\n");
466 for (i = nfsink; i <= nodes; i++)
467 xprintf("%6s%6d%18s%10d\n", "", i, "", ipred[i]);
472 { for (i = nfsink; i <= nodes; i++)
473 { double temp = - (double)ipred[i];
474 glp_vertex *v = G->v[i];
475 memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
479 /* Free working arrays. */
486 /* The instance has been successfully generated. */
491 /***********************************************************************
492 * The routine cresup randomly distributes the total supply among the
495 static void cresup(struct csa *csa)
497 xassert(itsup > nsorc);
499 for (i = 1; i <= nsorc; i++)
501 for (i = 1; i <= nsorc; i++)
502 { ksp = iran(csa, 1, ks);
503 j = iran(csa, 1, nsorc);
507 j = iran(csa, 1, nsorc);
508 isup[j] += itsup - ks * nsorc;
512 /***********************************************************************
513 * The routine chain adds node lpick to the end of the chain with source
516 static void chain(struct csa *csa, int lpick, int lsorc)
520 for (i = 1; i <= lpick; i++)
532 /***********************************************************************
533 * The routine chnarc puts the arcs in the chain from source lsorc into
534 * the ihead and itail arrays for sorting. */
536 static void chnarc(struct csa *csa, int lsorc)
540 L10: if (ito == lsorc) return;
544 itail[nsort] = ifrom;
549 /***********************************************************************
550 * The routine sort sorts the nsort arcs in the ihead and itail arrays.
551 * ihead is used as the sort key (i.e. forward star sort order). */
553 static void sort(struct csa *csa)
554 { int i, j, k, l, m, n, it;
563 if (itail[i] <= itail[l]) goto L40;
571 if (i >= 1) goto L30;
573 if (j <= k) goto L20;
577 /***********************************************************************
578 * The routine pickj creates a random number of arcs out of node 'it'.
579 * Various parameters are dynamically adjusted in an attempt to ensure
580 * that the generated network has the correct number of arcs. */
582 static void pickj(struct csa *csa, int it)
583 { int j, k, l, nn, nupbnd, icap, jcap, icost;
584 if ((nodlft - 1) * 2 > iarcs - narcs - 1)
588 if ((iarcs - narcs + nonsor - ktl - 1) / nodlft - nonsor + 1 >= 0)
591 { nupbnd = (iarcs - narcs - nodlft) / nodlft * 2;
592 L40: k = iran(csa, 1, nupbnd);
593 if (nodlft == 1) k = iarcs - narcs;
594 if ((nodlft - 1) * (nonsor - 1) < iarcs - narcs - k) goto L40;
597 for (j = 1; j <= k; j++)
598 { nn = iran(csa, 1, ktl);
600 for (l = nftsor; l <= nodes; l++)
603 if (nn == 0) goto L70;
609 jcap = iran(csa, 1, 100);
611 icap = iran(csa, mincap, maxcap);
612 icost = iran(csa, mincst, maxcst);
614 xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, l, "", icost,
617 { glp_arc *a = glp_add_arc(G, it, l);
619 { double temp = (double)icap;
620 memcpy((char *)a->data + a_cap, &temp, sizeof(double));
623 { double temp = (double)icost;
624 memcpy((char *)a->data + a_cost, &temp, sizeof(double));
632 /***********************************************************************
633 * The routine assign generate assignment problems. It defines the unit
634 * supplies, builds a skeleton, then calls pickj to create the arcs. */
636 static void assign(struct csa *csa)
637 { int i, it, nn, l, ll, icost;
640 for (i = 1; i <= nsorc; i++)
644 xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
647 { double temp = (double)isup[i];
648 glp_vertex *v = G->v[i];
649 memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
655 for (i = nfsink; i <= nodes; i++)
657 for (it = 1; it <= nsorc; it++)
658 { for (i = nfsink; i <= nodes; i++)
661 nn = iran(csa, 1, nsink - it + 1);
662 for (l = 1; l <= nsorc; l++)
670 icost = iran(csa, mincst, maxcst);
672 xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ll, "", icost,
675 { glp_arc *a = glp_add_arc(G, it, ll);
677 { double temp = (double)isup[1];
678 memcpy((char *)a->data + a_cap, &temp, sizeof(double));
681 { double temp = (double)icost;
682 memcpy((char *)a->data + a_cost, &temp, sizeof(double));
692 /***********************************************************************
693 * Portable congruential (uniform) random number generator:
695 * next_value = ((7**5) * previous_value) modulo ((2**31)-1)
697 * This generator consists of three routines:
699 * (1) setran - initializes constants and seed
700 * (2) iran - generates an integer random number
701 * (3) rran - generates a real random number
703 * The generator requires a machine with at least 32 bits of precision.
704 * The seed (iseed) must be in the range [1,(2**31)-1]. */
706 static void setran(struct csa *csa, int iseed)
707 { xassert(iseed >= 1);
716 /***********************************************************************
717 * The routine iran generates an integer random number between ilow and
718 * ihigh. If ilow > ihigh then iran returns ihigh. */
720 static int iran(struct csa *csa, int ilow, int ihigh)
721 { int ixhi, ixlo, ixalo, leftlo, ixahi, ifulhi, irtlo, iover,
724 ixlo = jran - ixhi * i16;
726 leftlo = ixalo / i16;
728 ifulhi = ixahi + leftlo;
729 irtlo = ixalo - leftlo * i16;
730 iover = ifulhi / i15;
731 irthi = ifulhi - iover * i15;
732 jran = ((irtlo - modul) + irthi * i16) + iover;
733 if (jran < 0) jran += modul;
734 j = ihigh - ilow + 1;
736 return jran % j + ilow;
741 /**********************************************************************/
744 static int scan(char card[80+1], int pos, int len)
746 memcpy(buf, &card[pos-1], len);
754 xassert(fgets(card, sizeof(card), stdin) == card);
755 parm[1] = scan(card, 1, 8);
756 parm[2] = scan(card, 9, 8);
757 xassert(fgets(card, sizeof(card), stdin) == card);
758 parm[3] = scan(card, 1, 5);
759 parm[4] = scan(card, 6, 5);
760 parm[5] = scan(card, 11, 5);
761 parm[6] = scan(card, 16, 5);
762 parm[7] = scan(card, 21, 5);
763 parm[8] = scan(card, 26, 5);
764 parm[9] = scan(card, 31, 10);
765 parm[10] = scan(card, 41, 5);
766 parm[11] = scan(card, 46, 5);
767 parm[12] = scan(card, 51, 5);
768 parm[13] = scan(card, 56, 5);
769 parm[14] = scan(card, 61, 10);
770 parm[15] = scan(card, 71, 10);
771 glp_netgen(NULL, 0, 0, 0, parm);