COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpnet03.c @ 9:33de93886c88

subpack-glpk
Last change on this file since 9:33de93886c88 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 12 years ago

Import GLPK 4.47

File size: 24.3 KB
Line 
1/* glpnet03.c (Klingman's network problem generator) */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
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>.
9*
10*  The translation was made by Andrew Makhorin <mao@gnu.org>.
11*
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.
16*
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.
21*
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***********************************************************************/
25
26#include "glpapi.h"
27
28/***********************************************************************
29*  NAME
30*
31*  glp_netgen - Klingman's network problem generator
32*
33*  SYNOPSIS
34*
35*  int glp_netgen(glp_graph *G, int v_rhs, int a_cap, int a_cost,
36*     const int parm[1+15]);
37*
38*  DESCRIPTION
39*
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
43*  problems.
44*
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.
48*
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.
52*
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.
56*
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.
60*
61*  The array parm contains description of the network to be generated:
62*
63*  parm[0]           not used
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
78*                    the maximum cost
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
82*
83*  The routine generates a transportation problem if:
84*
85*     nsorc + nsink = nodes, ntsorc = 0, and ntsink = 0.
86*
87*  The routine generates an assignment problem if the requirements for
88*  a transportation problem are met and:
89*
90*     nsorc = nsink and itsup = nsorc.
91*
92*  RETURNS
93*
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.
97*
98*  REFERENCES
99*
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. */
103
104struct csa
105{     /* common storage area */
106      glp_graph *G;
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;
112};
113
114#define G      (csa->G)
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)
135#if 0
136/* spent a day to find out this bug */
137#define ist    (csa->ist)
138#else
139#define ist    (ipred[0])
140#endif
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)
152
153static void cresup(struct csa *csa);
154static void chain(struct csa *csa, int lpick, int lsorc);
155static void chnarc(struct csa *csa, int lsorc);
156static void sort(struct csa *csa);
157static void pickj(struct csa *csa, int it);
158static void assign(struct csa *csa);
159static void setran(struct csa *csa, int iseed);
160static int iran(struct csa *csa, int ilow, int ihigh);
161
162int 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;
169      G = G_;
170      v_rhs = _v_rhs;
171      a_cap = _a_cap;
172      a_cost = _a_cost;
173      if (G != NULL)
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);
180      }
181      /* Input the user's random number seed and fix it if
182         non-positive. */
183      iseed = parm[1];
184      nprob = parm[2];
185      if (iseed <= 0) iseed = 13502460;
186      setran(csa, iseed);
187      /* Input the user's problem characteristics. */
188      nodes = parm[3];
189      nsorc = parm[4];
190      nsink = parm[5];
191      iarcs = parm[6];
192      mincst = parm[7];
193      maxcst = parm[8];
194      itsup = parm[9];
195      ntsorc = parm[10];
196      ntsink = parm[11];
197      iphic = parm[12];
198      ipcap = parm[13];
199      mincap = parm[14];
200      maxcap = parm[15];
201      /* Check the size of the problem. */
202      if (!(10 <= nodes && nodes <= 100000))
203      {  ret = 1;
204         goto done;
205      }
206      /* Check user supplied parameters for consistency. */
207      if (!(nsorc >= 0 && nsink >= 0 && nsorc + nsink <= nodes))
208      {  ret = 2;
209         goto done;
210      }
211      if (iarcs < 0)
212      {  ret = 3;
213         goto done;
214      }
215      if (mincst > maxcst)
216      {  ret = 4;
217         goto done;
218      }
219      if (itsup < 0)
220      {  ret = 5;
221         goto done;
222      }
223      if (!(0 <= ntsorc && ntsorc <= nsorc))
224      {  ret = 6;
225         goto done;
226      }
227      if (!(0 <= ntsink && ntsink <= nsink))
228      {  ret = 7;
229         goto done;
230      }
231      if (!(0 <= iphic && iphic <= 100))
232      {  ret = 8;
233         goto done;
234      }
235      if (!(0 <= ipcap && ipcap <= 100))
236      {  ret = 9;
237         goto done;
238      }
239      if (mincap > maxcap)
240      {  ret = 10;
241         goto done;
242      }
243      /* Initailize the graph object. */
244      if (G != NULL)
245      {  glp_erase_graph(G, G->v_size, G->a_size);
246         glp_add_vertices(G, nodes);
247         if (v_rhs >= 0)
248         {  double zero = 0.0;
249            for (i = 1; i <= nodes; i++)
250            {  glp_vertex *v = G->v[i];
251               memcpy((char *)v->data + v_rhs, &zero, sizeof(double));
252            }
253         }
254      }
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. */
263      if (G == NULL)
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,
270            mincap, maxcap);
271      }
272      else
273         glp_set_graph_name(G, "NETGEN");
274      /* Set various constants used in the program. */
275      narcs = 0;
276      nskel = 0;
277      nltr = nodes - nsink;
278      ltsink = nltr + ntsink;
279      ntrans = nltr - nsorc;
280      nfsink = nltr + 1;
281      nonsor = nodes - nsorc + ntsorc;
282      npsink = nsink - ntsink;
283      nodlft = nodes - nsink + ntsink;
284      nftr = nsorc + 1;
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 &&
289          itsup == nsorc)
290      {  assign(csa);
291         nskel = nsorc;
292         goto L390;
293      }
294      cresup(csa);
295      /* Print the supply records. */
296      if (G == NULL)
297      {  xprintf("SUPPLY\n");
298         for (i = 1; i <= nsorc; i++)
299            xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
300         xprintf("ARCS\n");
301      }
302      else
303      {  if (v_rhs >= 0)
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));
308            }
309         }
310      }
311      /* Make the sources point to themselves in ipred array. */
312      for (i = 1; i <= nsorc; i++)
313         ipred[i] = i;
314      if (ntrans == 0) goto L170;
315      /* Chain the transshipment nodes together in the ipred array. */
316      ist = nftr;
317      ipred[nltr] = 0;
318      for (i = nftr; i < nltr; i++)
319         ipred[i] = i+1;
320      /* Form even length chains for 60 percent of the transshipments.*/
321      ntravl = 6 * ntrans / 10;
322      ntrrem = ntrans - ntravl;
323L140: lsorc = 1;
324      while (ntravl != 0)
325      {  lpick = iran(csa, 1, ntravl + ntrrem);
326         ntravl--;
327         chain(csa, lpick, lsorc);
328         if (lsorc == nsorc) goto L140;
329         lsorc++;
330      }
331      /* Add the remaining transshipments to the chains. */
332      while (ntrrem != 0)
333      {
334         lpick = iran(csa, 1, ntrrem);
335         ntrrem--;
336         lsorc = iran(csa, 1, nsorc);
337         chain(csa, lpick, lsorc);
338      }
339L170: /* Set all demands equal to zero. */
340      for (i = nfsink; i <= nodes; i++)
341         ipred[i] = 0;
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++)
348            iflag[i] = 0;
349         /* Choose the number of sinks to be hooked up to the current
350            chain. */
351         if (ntrans != 0)
352            nsksr = (nsort * 2 * nsink) / ntrans;
353         else
354            nsksr = nsink / nsorc + 1;
355         if (nsksr < 2) nsksr = 2;
356         if (nsksr > nsink) nsksr = nsink;
357         nsrchn = nsort;
358         /* Randomly pick nsksr sinks and put their names in lsinks. */
359         ktl = nsink;
360         for (j = 1; j <= nsksr; j++)
361         {  item = iran(csa, 1, ktl);
362            ktl--;
363            for (l = nfsink; l <= nodes; l++)
364            {  if (iflag[l] != 1)
365               {  item--;
366                  if (item == 0) goto L230;
367               }
368            }
369            break;
370L230:       lsinks[j] = l;
371            iflag[l] = 1;
372         }
373         /* If last source chain, add all sinks with zero demand to
374            lsinks list. */
375         if (lsorc == nsorc)
376         {  for (j = nfsink; j <= nodes; j++)
377            {  if (ipred[j] == 0 && iflag[j] != 1)
378               {  nsksr++;
379                  lsinks[nsksr] = j;
380                  iflag[j] = 1;
381               }
382            }
383         }
384         /* Create demands for group of sinks in lsinks. */
385         ks = isup[lsorc] / nsksr;
386         k = ipred[lsorc];
387         for (i = 1; i <= nsksr; i++)
388         {  nsort++;
389            ksp = iran(csa, 1, ks);
390            j = iran(csa, 1, nsksr);
391            itail[nsort] = k;
392            li = lsinks[i];
393            ihead[nsort] = li;
394            ipred[li] += ksp;
395            li = lsinks[j];
396            ipred[li] += ks - ksp;
397            n = iran(csa, 1, nsrchn);
398            k = lsorc;
399            for (ii = 1; ii <= n; ii++)
400               k = ipred[k];
401         }
402         li = lsinks[1];
403         ipred[li] += isup[lsorc] - ks * nsksr;
404         nskel += nsort;
405         /* Sort the arcs in the chain from source lsorc using itail as
406            sort key. */
407         sort(csa);
408         /* Print this part of skeleton and create the arcs for these
409            nodes. */
410         i = 1;
411         itail[nsort+1] = 0;
412L300:    for (j = nftsor; j <= nodes; j++)
413            iflag[j] = 0;
414         ktl = nonsor - 1;
415         it = itail[i];
416         iflag[it] = 1;
417L320:    ih = ihead[i];
418         iflag[ih] = 1;
419         narcs++;
420         ktl--;
421         /* Determine if this skeleton arc should be capacitated. */
422         icap = itsup;
423         jcap = iran(csa, 1, 100);
424         if (jcap <= ipcap)
425         {  icap = isup[lsorc];
426            if (mincap > icap) icap = mincap;
427         }
428         /* Determine if this skeleton arc should have the maximum
429            cost. */
430         icost = maxcst;
431         jcost = iran(csa, 1, 100);
432         if (jcost > iphic)
433            icost = iran(csa, mincst, maxcst);
434         if (G == NULL)
435            xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ih, "", icost,
436               icap);
437         else
438         {  glp_arc *a = glp_add_arc(G, it, ih);
439            if (a_cap >= 0)
440            {  double temp = (double)icap;
441               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
442            }
443            if (a_cost >= 0)
444            {  double temp = (double)icost;
445               memcpy((char *)a->data + a_cost, &temp, sizeof(double));
446            }
447         }
448         i++;
449         if (itail[i] == it) goto L320;
450         pickj(csa, it);
451         if (i <= nsort) goto L300;
452      }
453      /* Create arcs from the transshipment sinks. */
454      if (ntsink != 0)
455      {  for (i = nfsink; i <= ltsink; i++)
456         {  for (j = nftsor; j <= nodes; j++)
457               iflag[j] = 0;
458            ktl = nonsor - 1;
459            iflag[i] = 1;
460            pickj(csa, i);
461         }
462      }
463L390: /* Print the demand records and end record. */
464      if (G == NULL)
465      {  xprintf("DEMAND\n");
466         for (i = nfsink; i <= nodes; i++)
467            xprintf("%6s%6d%18s%10d\n", "", i, "", ipred[i]);
468         xprintf("END\n");
469      }
470      else
471      {  if (v_rhs >= 0)
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));
476            }
477         }
478      }
479      /* Free working arrays. */
480      xfree(ipred);
481      xfree(ihead);
482      xfree(itail);
483      xfree(iflag);
484      xfree(isup);
485      xfree(lsinks);
486      /* The instance has been successfully generated. */
487      ret = 0;
488done: return ret;
489}
490
491/***********************************************************************
492*  The routine cresup randomly distributes the total supply among the
493*  source nodes. */
494
495static void cresup(struct csa *csa)
496{     int i, j, ks, ksp;
497      xassert(itsup > nsorc);
498      ks = itsup / nsorc;
499      for (i = 1; i <= nsorc; i++)
500         isup[i] = 0;
501      for (i = 1; i <= nsorc; i++)
502      {  ksp = iran(csa, 1, ks);
503         j = iran(csa, 1, nsorc);
504         isup[i] += ksp;
505         isup[j] += ks - ksp;
506      }
507      j = iran(csa, 1, nsorc);
508      isup[j] += itsup - ks * nsorc;
509      return;
510}
511
512/***********************************************************************
513*  The routine chain adds node lpick to the end of the chain with source
514*  node lsorc. */
515
516static void chain(struct csa *csa, int lpick, int lsorc)
517{     int i, j, k, l, m;
518      k = 0;
519      m = ist;
520      for (i = 1; i <= lpick; i++)
521      {  l = k;
522         k = m;
523         m = ipred[k];
524      }
525      ipred[l] = m;
526      j = ipred[lsorc];
527      ipred[k] = j;
528      ipred[lsorc] = k;
529      return;
530}
531
532/***********************************************************************
533*  The routine chnarc puts the arcs in the chain from source lsorc into
534*  the ihead and itail arrays for sorting. */
535
536static void chnarc(struct csa *csa, int lsorc)
537{     int ito, ifrom;
538      nsort = 0;
539      ito = ipred[lsorc];
540L10:  if (ito == lsorc) return;
541      nsort++;
542      ifrom = ipred[ito];
543      ihead[nsort] = ito;
544      itail[nsort] = ifrom;
545      ito = ifrom;
546      goto L10;
547}
548
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). */
552
553static void sort(struct csa *csa)
554{     int i, j, k, l, m, n, it;
555      n = nsort;
556      m = n;
557L10:  m /= 2;
558      if (m == 0) return;
559      k = n - m;
560      j = 1;
561L20:  i = j;
562L30:  l = i + m;
563      if (itail[i] <= itail[l]) goto L40;
564      it = itail[i];
565      itail[i] = itail[l];
566      itail[l] = it;
567      it = ihead[i];
568      ihead[i] = ihead[l];
569      ihead[l] = it;
570      i -= m;
571      if (i >= 1) goto L30;
572L40:  j++;
573      if (j <= k) goto L20;
574      goto L10;
575}
576
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. */
581
582static 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)
585      {  nodlft--;
586         return;
587      }
588      if ((iarcs - narcs + nonsor - ktl - 1) / nodlft - nonsor + 1 >= 0)
589         k = nonsor;
590      else
591      {  nupbnd = (iarcs - narcs - nodlft) / nodlft * 2;
592L40:     k = iran(csa, 1, nupbnd);
593         if (nodlft == 1) k = iarcs - narcs;
594         if ((nodlft - 1) * (nonsor - 1) < iarcs - narcs - k) goto L40;
595      }
596      nodlft--;
597      for (j = 1; j <= k; j++)
598      {  nn = iran(csa, 1, ktl);
599         ktl--;
600         for (l = nftsor; l <= nodes; l++)
601         {  if (iflag[l] != 1)
602            {  nn--;
603               if (nn == 0) goto L70;
604            }
605         }
606         return;
607L70:     iflag[l] = 1;
608         icap = itsup;
609         jcap = iran(csa, 1, 100);
610         if (jcap <= ipcap)
611            icap = iran(csa, mincap, maxcap);
612         icost = iran(csa, mincst, maxcst);
613         if (G == NULL)
614            xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, l, "", icost,
615               icap);
616         else
617         {  glp_arc *a = glp_add_arc(G, it, l);
618            if (a_cap >= 0)
619            {  double temp = (double)icap;
620               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
621            }
622            if (a_cost >= 0)
623            {  double temp = (double)icost;
624               memcpy((char *)a->data + a_cost, &temp, sizeof(double));
625            }
626         }
627         narcs++;
628      }
629      return;
630}
631
632/***********************************************************************
633*  The routine assign generate assignment problems. It defines the unit
634*  supplies, builds a skeleton, then calls pickj to create the arcs. */
635
636static void assign(struct csa *csa)
637{     int i, it, nn, l, ll, icost;
638      if (G == NULL)
639         xprintf("SUPPLY\n");
640      for (i = 1; i <= nsorc; i++)
641      {  isup[i] = 1;
642         iflag[i] = 0;
643         if (G == NULL)
644            xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
645         else
646         {  if (v_rhs >= 0)
647            {  double temp = (double)isup[i];
648               glp_vertex *v = G->v[i];
649               memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
650            }
651         }
652      }
653      if (G == NULL)
654         xprintf("ARCS\n");
655      for (i = nfsink; i <= nodes; i++)
656         ipred[i] = 1;
657      for (it = 1; it <= nsorc; it++)
658      {  for (i = nfsink; i <= nodes; i++)
659            iflag[i] = 0;
660         ktl = nsink - 1;
661         nn = iran(csa, 1, nsink - it + 1);
662         for (l = 1; l <= nsorc; l++)
663         {  if (iflag[l] != 1)
664            {  nn--;
665               if (nn == 0) break;
666            }
667         }
668         narcs++;
669         ll = nsorc + l;
670         icost = iran(csa, mincst, maxcst);
671         if (G == NULL)
672            xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ll, "", icost,
673               isup[1]);
674         else
675         {  glp_arc *a = glp_add_arc(G, it, ll);
676            if (a_cap >= 0)
677            {  double temp = (double)isup[1];
678               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
679            }
680            if (a_cost >= 0)
681            {  double temp = (double)icost;
682               memcpy((char *)a->data + a_cost, &temp, sizeof(double));
683            }
684         }
685         iflag[l] = 1;
686         iflag[ll] = 1;
687         pickj(csa, it);
688      }
689      return;
690}
691
692/***********************************************************************
693*  Portable congruential (uniform) random number generator:
694*
695*     next_value = ((7**5) * previous_value) modulo ((2**31)-1)
696*
697*  This generator consists of three routines:
698*
699*  (1) setran - initializes constants and seed
700*  (2) iran   - generates an integer random number
701*  (3) rran   - generates a real random number
702*
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]. */
705
706static void setran(struct csa *csa, int iseed)
707{     xassert(iseed >= 1);
708      mult = 16807;
709      modul = 2147483647;
710      i15 = 1 << 15;
711      i16 = 1 << 16;
712      jran = iseed;
713      return;
714}
715
716/***********************************************************************
717*  The routine iran generates an integer random number between ilow and
718*  ihigh. If ilow > ihigh then iran returns ihigh. */
719
720static int iran(struct csa *csa, int ilow, int ihigh)
721{     int ixhi, ixlo, ixalo, leftlo, ixahi, ifulhi, irtlo, iover,
722         irthi, j;
723      ixhi = jran / i16;
724      ixlo = jran - ixhi * i16;
725      ixalo = ixlo * mult;
726      leftlo = ixalo / i16;
727      ixahi = ixhi * mult;
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;
735      if (j > 0)
736         return jran % j + ilow;
737      else
738         return ihigh;
739}
740
741/**********************************************************************/
742
743#if 0
744static int scan(char card[80+1], int pos, int len)
745{     char buf[10+1];
746      memcpy(buf, &card[pos-1], len);
747      buf[len] = '\0';
748      return atoi(buf);
749}
750
751int main(void)
752{     int parm[1+15];
753      char card[80+1];
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);
772      return 0;
773}
774#endif
775
776/* eof */
Note: See TracBrowser for help on using the repository browser.