src/glpnet03.c
author Alpar Juttner <alpar@cs.elte.hu>
Sun, 05 Dec 2010 17:35:23 +0100
changeset 2 4c8956a7bdf4
permissions -rw-r--r--
Set up CMAKE build environment
     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 
   104 struct 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 
   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);
   161 
   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;
   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;
   323 L140: 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       }
   339 L170: /* 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;
   370 L230:       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;
   412 L300:    for (j = nftsor; j <= nodes; j++)
   413             iflag[j] = 0;
   414          ktl = nonsor - 1;
   415          it = itail[i];
   416          iflag[it] = 1;
   417 L320:    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       }
   463 L390: /* 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;
   488 done: return ret;
   489 }
   490 
   491 /***********************************************************************
   492 *  The routine cresup randomly distributes the total supply among the
   493 *  source nodes. */
   494 
   495 static 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 
   516 static 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 
   536 static void chnarc(struct csa *csa, int lsorc)
   537 {     int ito, ifrom;
   538       nsort = 0;
   539       ito = ipred[lsorc];
   540 L10:  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 
   553 static void sort(struct csa *csa)
   554 {     int i, j, k, l, m, n, it;
   555       n = nsort;
   556       m = n;
   557 L10:  m /= 2;
   558       if (m == 0) return;
   559       k = n - m;
   560       j = 1;
   561 L20:  i = j;
   562 L30:  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;
   572 L40:  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 
   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)
   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;
   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;
   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;
   607 L70:     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 
   636 static 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 
   706 static 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 
   720 static 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
   744 static 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 
   751 int 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 */