1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnet03.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,776 @@
     1.4 +/* glpnet03.c (Klingman's network problem generator) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  This code is the result of translation of the Fortran program NETGEN
    1.10 +*  developed by Dr. Darwin Klingman, which is publically available from
    1.11 +*  NETLIB at <http://www.netlib.org/lp/generators>.
    1.12 +*
    1.13 +*  The translation was made by Andrew Makhorin <mao@gnu.org>.
    1.14 +*
    1.15 +*  GLPK is free software: you can redistribute it and/or modify it
    1.16 +*  under the terms of the GNU General Public License as published by
    1.17 +*  the Free Software Foundation, either version 3 of the License, or
    1.18 +*  (at your option) any later version.
    1.19 +*
    1.20 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.21 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.22 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.23 +*  License for more details.
    1.24 +*
    1.25 +*  You should have received a copy of the GNU General Public License
    1.26 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.27 +***********************************************************************/
    1.28 +
    1.29 +#include "glpapi.h"
    1.30 +
    1.31 +/***********************************************************************
    1.32 +*  NAME
    1.33 +*
    1.34 +*  glp_netgen - Klingman's network problem generator
    1.35 +*
    1.36 +*  SYNOPSIS
    1.37 +*
    1.38 +*  int glp_netgen(glp_graph *G, int v_rhs, int a_cap, int a_cost,
    1.39 +*     const int parm[1+15]);
    1.40 +*
    1.41 +*  DESCRIPTION
    1.42 +*
    1.43 +*  The routine glp_netgen is a network problem generator developed by
    1.44 +*  Dr. Darwin Klingman. It can create capacitated and uncapacitated
    1.45 +*  minimum cost flow (or transshipment), transportation, and assignment
    1.46 +*  problems.
    1.47 +*
    1.48 +*  The parameter G specifies the graph object, to which the generated
    1.49 +*  problem data have to be stored. Note that on entry the graph object
    1.50 +*  is erased with the routine glp_erase_graph.
    1.51 +*
    1.52 +*  The parameter v_rhs specifies an offset of the field of type double
    1.53 +*  in the vertex data block, to which the routine stores the supply or
    1.54 +*  demand value. If v_rhs < 0, the value is not stored.
    1.55 +*
    1.56 +*  The parameter a_cap specifies an offset of the field of type double
    1.57 +*  in the arc data block, to which the routine stores the arc capacity.
    1.58 +*  If a_cap < 0, the capacity is not stored.
    1.59 +*
    1.60 +*  The parameter a_cost specifies an offset of the field of type double
    1.61 +*  in the arc data block, to which the routine stores the per-unit cost
    1.62 +*  if the arc flow. If a_cost < 0, the cost is not stored.
    1.63 +*
    1.64 +*  The array parm contains description of the network to be generated:
    1.65 +*
    1.66 +*  parm[0]           not used
    1.67 +*  parm[1]  (iseed)  8-digit positive random number seed
    1.68 +*  parm[2]  (nprob)  8-digit problem id number
    1.69 +*  parm[3]  (nodes)  total number of nodes
    1.70 +*  parm[4]  (nsorc)  total number of source nodes (including
    1.71 +*                    transshipment nodes)
    1.72 +*  parm[5]  (nsink)  total number of sink nodes (including
    1.73 +*                    transshipment nodes)
    1.74 +*  parm[6]  (iarcs)  number of arcs
    1.75 +*  parm[7]  (mincst) minimum cost for arcs
    1.76 +*  parm[8]  (maxcst) maximum cost for arcs
    1.77 +*  parm[9]  (itsup)  total supply
    1.78 +*  parm[10] (ntsorc) number of transshipment source nodes
    1.79 +*  parm[11] (ntsink) number of transshipment sink nodes
    1.80 +*  parm[12] (iphic)  percentage of skeleton arcs to be given
    1.81 +*                    the maximum cost
    1.82 +*  parm[13] (ipcap)  percentage of arcs to be capacitated
    1.83 +*  parm[14] (mincap) minimum upper bound for capacitated arcs
    1.84 +*  parm[15] (maxcap) maximum upper bound for capacitated arcs
    1.85 +*
    1.86 +*  The routine generates a transportation problem if:
    1.87 +*
    1.88 +*     nsorc + nsink = nodes, ntsorc = 0, and ntsink = 0.
    1.89 +*
    1.90 +*  The routine generates an assignment problem if the requirements for
    1.91 +*  a transportation problem are met and:
    1.92 +*
    1.93 +*     nsorc = nsink and itsup = nsorc.
    1.94 +*
    1.95 +*  RETURNS
    1.96 +*
    1.97 +*  If the instance was successfully generated, the routine glp_netgen
    1.98 +*  returns zero; otherwise, if specified parameters are inconsistent,
    1.99 +*  the routine returns a non-zero error code.
   1.100 +*
   1.101 +*  REFERENCES
   1.102 +*
   1.103 +*  D.Klingman, A.Napier, and J.Stutz. NETGEN: A program for generating
   1.104 +*  large scale capacitated assignment, transportation, and minimum cost
   1.105 +*  flow networks. Management Science 20 (1974), 814-20. */
   1.106 +
   1.107 +struct csa
   1.108 +{     /* common storage area */
   1.109 +      glp_graph *G;
   1.110 +      int v_rhs, a_cap, a_cost;
   1.111 +      int nodes, iarcs, mincst, maxcst, itsup, nsorc, nsink, nonsor,
   1.112 +         nfsink, narcs, nsort, nftsor, ipcap, mincap, maxcap, ktl,
   1.113 +         nodlft, *ipred, *ihead, *itail, *iflag, *isup, *lsinks, mult,
   1.114 +         modul, i15, i16, jran;
   1.115 +};
   1.116 +
   1.117 +#define G      (csa->G)
   1.118 +#define v_rhs  (csa->v_rhs)
   1.119 +#define a_cap  (csa->a_cap)
   1.120 +#define a_cost (csa->a_cost)
   1.121 +#define nodes  (csa->nodes)
   1.122 +#define iarcs  (csa->iarcs)
   1.123 +#define mincst (csa->mincst)
   1.124 +#define maxcst (csa->maxcst)
   1.125 +#define itsup  (csa->itsup)
   1.126 +#define nsorc  (csa->nsorc)
   1.127 +#define nsink  (csa->nsink)
   1.128 +#define nonsor (csa->nonsor)
   1.129 +#define nfsink (csa->nfsink)
   1.130 +#define narcs  (csa->narcs)
   1.131 +#define nsort  (csa->nsort)
   1.132 +#define nftsor (csa->nftsor)
   1.133 +#define ipcap  (csa->ipcap)
   1.134 +#define mincap (csa->mincap)
   1.135 +#define maxcap (csa->maxcap)
   1.136 +#define ktl    (csa->ktl)
   1.137 +#define nodlft (csa->nodlft)
   1.138 +#if 0
   1.139 +/* spent a day to find out this bug */
   1.140 +#define ist    (csa->ist)
   1.141 +#else
   1.142 +#define ist    (ipred[0])
   1.143 +#endif
   1.144 +#define ipred  (csa->ipred)
   1.145 +#define ihead  (csa->ihead)
   1.146 +#define itail  (csa->itail)
   1.147 +#define iflag  (csa->iflag)
   1.148 +#define isup   (csa->isup)
   1.149 +#define lsinks (csa->lsinks)
   1.150 +#define mult   (csa->mult)
   1.151 +#define modul  (csa->modul)
   1.152 +#define i15    (csa->i15)
   1.153 +#define i16    (csa->i16)
   1.154 +#define jran   (csa->jran)
   1.155 +
   1.156 +static void cresup(struct csa *csa);
   1.157 +static void chain(struct csa *csa, int lpick, int lsorc);
   1.158 +static void chnarc(struct csa *csa, int lsorc);
   1.159 +static void sort(struct csa *csa);
   1.160 +static void pickj(struct csa *csa, int it);
   1.161 +static void assign(struct csa *csa);
   1.162 +static void setran(struct csa *csa, int iseed);
   1.163 +static int iran(struct csa *csa, int ilow, int ihigh);
   1.164 +
   1.165 +int glp_netgen(glp_graph *G_, int _v_rhs, int _a_cap, int _a_cost,
   1.166 +      const int parm[1+15])
   1.167 +{     struct csa _csa, *csa = &_csa;
   1.168 +      int iseed, nprob, ntsorc, ntsink, iphic, i, nskel, nltr, ltsink,
   1.169 +         ntrans, npsink, nftr, npsorc, ntravl, ntrrem, lsorc, lpick,
   1.170 +         nsksr, nsrchn, j, item, l, ks, k, ksp, li, n, ii, it, ih, icap,
   1.171 +         jcap, icost, jcost, ret;
   1.172 +      G = G_;
   1.173 +      v_rhs = _v_rhs;
   1.174 +      a_cap = _a_cap;
   1.175 +      a_cost = _a_cost;
   1.176 +      if (G != NULL)
   1.177 +      {  if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
   1.178 +            xerror("glp_netgen: v_rhs = %d; invalid offset\n", v_rhs);
   1.179 +         if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
   1.180 +            xerror("glp_netgen: a_cap = %d; invalid offset\n", a_cap);
   1.181 +         if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
   1.182 +            xerror("glp_netgen: a_cost = %d; invalid offset\n", a_cost);
   1.183 +      }
   1.184 +      /* Input the user's random number seed and fix it if
   1.185 +         non-positive. */
   1.186 +      iseed = parm[1];
   1.187 +      nprob = parm[2];
   1.188 +      if (iseed <= 0) iseed = 13502460;
   1.189 +      setran(csa, iseed);
   1.190 +      /* Input the user's problem characteristics. */
   1.191 +      nodes = parm[3];
   1.192 +      nsorc = parm[4];
   1.193 +      nsink = parm[5];
   1.194 +      iarcs = parm[6];
   1.195 +      mincst = parm[7];
   1.196 +      maxcst = parm[8];
   1.197 +      itsup = parm[9];
   1.198 +      ntsorc = parm[10];
   1.199 +      ntsink = parm[11];
   1.200 +      iphic = parm[12];
   1.201 +      ipcap = parm[13];
   1.202 +      mincap = parm[14];
   1.203 +      maxcap = parm[15];
   1.204 +      /* Check the size of the problem. */
   1.205 +      if (!(10 <= nodes && nodes <= 100000))
   1.206 +      {  ret = 1;
   1.207 +         goto done;
   1.208 +      }
   1.209 +      /* Check user supplied parameters for consistency. */
   1.210 +      if (!(nsorc >= 0 && nsink >= 0 && nsorc + nsink <= nodes))
   1.211 +      {  ret = 2;
   1.212 +         goto done;
   1.213 +      }
   1.214 +      if (iarcs < 0)
   1.215 +      {  ret = 3;
   1.216 +         goto done;
   1.217 +      }
   1.218 +      if (mincst > maxcst)
   1.219 +      {  ret = 4;
   1.220 +         goto done;
   1.221 +      }
   1.222 +      if (itsup < 0)
   1.223 +      {  ret = 5;
   1.224 +         goto done;
   1.225 +      }
   1.226 +      if (!(0 <= ntsorc && ntsorc <= nsorc))
   1.227 +      {  ret = 6;
   1.228 +         goto done;
   1.229 +      }
   1.230 +      if (!(0 <= ntsink && ntsink <= nsink))
   1.231 +      {  ret = 7;
   1.232 +         goto done;
   1.233 +      }
   1.234 +      if (!(0 <= iphic && iphic <= 100))
   1.235 +      {  ret = 8;
   1.236 +         goto done;
   1.237 +      }
   1.238 +      if (!(0 <= ipcap && ipcap <= 100))
   1.239 +      {  ret = 9;
   1.240 +         goto done;
   1.241 +      }
   1.242 +      if (mincap > maxcap)
   1.243 +      {  ret = 10;
   1.244 +         goto done;
   1.245 +      }
   1.246 +      /* Initailize the graph object. */
   1.247 +      if (G != NULL)
   1.248 +      {  glp_erase_graph(G, G->v_size, G->a_size);
   1.249 +         glp_add_vertices(G, nodes);
   1.250 +         if (v_rhs >= 0)
   1.251 +         {  double zero = 0.0;
   1.252 +            for (i = 1; i <= nodes; i++)
   1.253 +            {  glp_vertex *v = G->v[i];
   1.254 +               memcpy((char *)v->data + v_rhs, &zero, sizeof(double));
   1.255 +            }
   1.256 +         }
   1.257 +      }
   1.258 +      /* Allocate working arrays. */
   1.259 +      ipred = xcalloc(1+nodes, sizeof(int));
   1.260 +      ihead = xcalloc(1+nodes, sizeof(int));
   1.261 +      itail = xcalloc(1+nodes, sizeof(int));
   1.262 +      iflag = xcalloc(1+nodes, sizeof(int));
   1.263 +      isup = xcalloc(1+nodes, sizeof(int));
   1.264 +      lsinks = xcalloc(1+nodes, sizeof(int));
   1.265 +      /* Print the problem documentation records. */
   1.266 +      if (G == NULL)
   1.267 +      {  xprintf("BEGIN\n");
   1.268 +         xprintf("NETGEN PROBLEM%8d%10s%10d NODES AND%10d ARCS\n",
   1.269 +            nprob, "", nodes, iarcs);
   1.270 +         xprintf("USER:%11d%11d%11d%11d%11d%11d\nDATA:%11d%11d%11d%11d%"
   1.271 +            "11d%11d\n", iseed, nsorc, nsink, mincst,
   1.272 +            maxcst, itsup, ntsorc, ntsink, iphic, ipcap,
   1.273 +            mincap, maxcap);
   1.274 +      }
   1.275 +      else
   1.276 +         glp_set_graph_name(G, "NETGEN");
   1.277 +      /* Set various constants used in the program. */
   1.278 +      narcs = 0;
   1.279 +      nskel = 0;
   1.280 +      nltr = nodes - nsink;
   1.281 +      ltsink = nltr + ntsink;
   1.282 +      ntrans = nltr - nsorc;
   1.283 +      nfsink = nltr + 1;
   1.284 +      nonsor = nodes - nsorc + ntsorc;
   1.285 +      npsink = nsink - ntsink;
   1.286 +      nodlft = nodes - nsink + ntsink;
   1.287 +      nftr = nsorc + 1;
   1.288 +      nftsor = nsorc - ntsorc + 1;
   1.289 +      npsorc = nsorc - ntsorc;
   1.290 +      /* Randomly distribute the supply among the source nodes. */
   1.291 +      if (npsorc + npsink == nodes && npsorc == npsink &&
   1.292 +          itsup == nsorc)
   1.293 +      {  assign(csa);
   1.294 +         nskel = nsorc;
   1.295 +         goto L390;
   1.296 +      }
   1.297 +      cresup(csa);
   1.298 +      /* Print the supply records. */
   1.299 +      if (G == NULL)
   1.300 +      {  xprintf("SUPPLY\n");
   1.301 +         for (i = 1; i <= nsorc; i++)
   1.302 +            xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
   1.303 +         xprintf("ARCS\n");
   1.304 +      }
   1.305 +      else
   1.306 +      {  if (v_rhs >= 0)
   1.307 +         {  for (i = 1; i <= nsorc; i++)
   1.308 +            {  double temp = (double)isup[i];
   1.309 +               glp_vertex *v = G->v[i];
   1.310 +               memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
   1.311 +            }
   1.312 +         }
   1.313 +      }
   1.314 +      /* Make the sources point to themselves in ipred array. */
   1.315 +      for (i = 1; i <= nsorc; i++)
   1.316 +         ipred[i] = i;
   1.317 +      if (ntrans == 0) goto L170;
   1.318 +      /* Chain the transshipment nodes together in the ipred array. */
   1.319 +      ist = nftr;
   1.320 +      ipred[nltr] = 0;
   1.321 +      for (i = nftr; i < nltr; i++)
   1.322 +         ipred[i] = i+1;
   1.323 +      /* Form even length chains for 60 percent of the transshipments.*/
   1.324 +      ntravl = 6 * ntrans / 10;
   1.325 +      ntrrem = ntrans - ntravl;
   1.326 +L140: lsorc = 1;
   1.327 +      while (ntravl != 0)
   1.328 +      {  lpick = iran(csa, 1, ntravl + ntrrem);
   1.329 +         ntravl--;
   1.330 +         chain(csa, lpick, lsorc);
   1.331 +         if (lsorc == nsorc) goto L140;
   1.332 +         lsorc++;
   1.333 +      }
   1.334 +      /* Add the remaining transshipments to the chains. */
   1.335 +      while (ntrrem != 0)
   1.336 +      {
   1.337 +         lpick = iran(csa, 1, ntrrem);
   1.338 +         ntrrem--;
   1.339 +         lsorc = iran(csa, 1, nsorc);
   1.340 +         chain(csa, lpick, lsorc);
   1.341 +      }
   1.342 +L170: /* Set all demands equal to zero. */
   1.343 +      for (i = nfsink; i <= nodes; i++)
   1.344 +         ipred[i] = 0;
   1.345 +      /* The following loop takes one chain at a time (through the use
   1.346 +         of logic contained in the loop and calls to other routines) and
   1.347 +         creates the remaining network arcs. */
   1.348 +      for (lsorc = 1; lsorc <= nsorc; lsorc++)
   1.349 +      {  chnarc(csa, lsorc);
   1.350 +         for (i = nfsink; i <= nodes; i++)
   1.351 +            iflag[i] = 0;
   1.352 +         /* Choose the number of sinks to be hooked up to the current
   1.353 +            chain. */
   1.354 +         if (ntrans != 0)
   1.355 +            nsksr = (nsort * 2 * nsink) / ntrans;
   1.356 +         else
   1.357 +            nsksr = nsink / nsorc + 1;
   1.358 +         if (nsksr < 2) nsksr = 2;
   1.359 +         if (nsksr > nsink) nsksr = nsink;
   1.360 +         nsrchn = nsort;
   1.361 +         /* Randomly pick nsksr sinks and put their names in lsinks. */
   1.362 +         ktl = nsink;
   1.363 +         for (j = 1; j <= nsksr; j++)
   1.364 +         {  item = iran(csa, 1, ktl);
   1.365 +            ktl--;
   1.366 +            for (l = nfsink; l <= nodes; l++)
   1.367 +            {  if (iflag[l] != 1)
   1.368 +               {  item--;
   1.369 +                  if (item == 0) goto L230;
   1.370 +               }
   1.371 +            }
   1.372 +            break;
   1.373 +L230:       lsinks[j] = l;
   1.374 +            iflag[l] = 1;
   1.375 +         }
   1.376 +         /* If last source chain, add all sinks with zero demand to
   1.377 +            lsinks list. */
   1.378 +         if (lsorc == nsorc)
   1.379 +         {  for (j = nfsink; j <= nodes; j++)
   1.380 +            {  if (ipred[j] == 0 && iflag[j] != 1)
   1.381 +               {  nsksr++;
   1.382 +                  lsinks[nsksr] = j;
   1.383 +                  iflag[j] = 1;
   1.384 +               }
   1.385 +            }
   1.386 +         }
   1.387 +         /* Create demands for group of sinks in lsinks. */
   1.388 +         ks = isup[lsorc] / nsksr;
   1.389 +         k = ipred[lsorc];
   1.390 +         for (i = 1; i <= nsksr; i++)
   1.391 +         {  nsort++;
   1.392 +            ksp = iran(csa, 1, ks);
   1.393 +            j = iran(csa, 1, nsksr);
   1.394 +            itail[nsort] = k;
   1.395 +            li = lsinks[i];
   1.396 +            ihead[nsort] = li;
   1.397 +            ipred[li] += ksp;
   1.398 +            li = lsinks[j];
   1.399 +            ipred[li] += ks - ksp;
   1.400 +            n = iran(csa, 1, nsrchn);
   1.401 +            k = lsorc;
   1.402 +            for (ii = 1; ii <= n; ii++)
   1.403 +               k = ipred[k];
   1.404 +         }
   1.405 +         li = lsinks[1];
   1.406 +         ipred[li] += isup[lsorc] - ks * nsksr;
   1.407 +         nskel += nsort;
   1.408 +         /* Sort the arcs in the chain from source lsorc using itail as
   1.409 +            sort key. */
   1.410 +         sort(csa);
   1.411 +         /* Print this part of skeleton and create the arcs for these
   1.412 +            nodes. */
   1.413 +         i = 1;
   1.414 +         itail[nsort+1] = 0;
   1.415 +L300:    for (j = nftsor; j <= nodes; j++)
   1.416 +            iflag[j] = 0;
   1.417 +         ktl = nonsor - 1;
   1.418 +         it = itail[i];
   1.419 +         iflag[it] = 1;
   1.420 +L320:    ih = ihead[i];
   1.421 +         iflag[ih] = 1;
   1.422 +         narcs++;
   1.423 +         ktl--;
   1.424 +         /* Determine if this skeleton arc should be capacitated. */
   1.425 +         icap = itsup;
   1.426 +         jcap = iran(csa, 1, 100);
   1.427 +         if (jcap <= ipcap)
   1.428 +         {  icap = isup[lsorc];
   1.429 +            if (mincap > icap) icap = mincap;
   1.430 +         }
   1.431 +         /* Determine if this skeleton arc should have the maximum
   1.432 +            cost. */
   1.433 +         icost = maxcst;
   1.434 +         jcost = iran(csa, 1, 100);
   1.435 +         if (jcost > iphic)
   1.436 +            icost = iran(csa, mincst, maxcst);
   1.437 +         if (G == NULL)
   1.438 +            xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ih, "", icost,
   1.439 +               icap);
   1.440 +         else
   1.441 +         {  glp_arc *a = glp_add_arc(G, it, ih);
   1.442 +            if (a_cap >= 0)
   1.443 +            {  double temp = (double)icap;
   1.444 +               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
   1.445 +            }
   1.446 +            if (a_cost >= 0)
   1.447 +            {  double temp = (double)icost;
   1.448 +               memcpy((char *)a->data + a_cost, &temp, sizeof(double));
   1.449 +            }
   1.450 +         }
   1.451 +         i++;
   1.452 +         if (itail[i] == it) goto L320;
   1.453 +         pickj(csa, it);
   1.454 +         if (i <= nsort) goto L300;
   1.455 +      }
   1.456 +      /* Create arcs from the transshipment sinks. */
   1.457 +      if (ntsink != 0)
   1.458 +      {  for (i = nfsink; i <= ltsink; i++)
   1.459 +         {  for (j = nftsor; j <= nodes; j++)
   1.460 +               iflag[j] = 0;
   1.461 +            ktl = nonsor - 1;
   1.462 +            iflag[i] = 1;
   1.463 +            pickj(csa, i);
   1.464 +         }
   1.465 +      }
   1.466 +L390: /* Print the demand records and end record. */
   1.467 +      if (G == NULL)
   1.468 +      {  xprintf("DEMAND\n");
   1.469 +         for (i = nfsink; i <= nodes; i++)
   1.470 +            xprintf("%6s%6d%18s%10d\n", "", i, "", ipred[i]);
   1.471 +         xprintf("END\n");
   1.472 +      }
   1.473 +      else
   1.474 +      {  if (v_rhs >= 0)
   1.475 +         {  for (i = nfsink; i <= nodes; i++)
   1.476 +            {  double temp = - (double)ipred[i];
   1.477 +               glp_vertex *v = G->v[i];
   1.478 +               memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
   1.479 +            }
   1.480 +         }
   1.481 +      }
   1.482 +      /* Free working arrays. */
   1.483 +      xfree(ipred);
   1.484 +      xfree(ihead);
   1.485 +      xfree(itail);
   1.486 +      xfree(iflag);
   1.487 +      xfree(isup);
   1.488 +      xfree(lsinks);
   1.489 +      /* The instance has been successfully generated. */
   1.490 +      ret = 0;
   1.491 +done: return ret;
   1.492 +}
   1.493 +
   1.494 +/***********************************************************************
   1.495 +*  The routine cresup randomly distributes the total supply among the
   1.496 +*  source nodes. */
   1.497 +
   1.498 +static void cresup(struct csa *csa)
   1.499 +{     int i, j, ks, ksp;
   1.500 +      xassert(itsup > nsorc);
   1.501 +      ks = itsup / nsorc;
   1.502 +      for (i = 1; i <= nsorc; i++)
   1.503 +         isup[i] = 0;
   1.504 +      for (i = 1; i <= nsorc; i++)
   1.505 +      {  ksp = iran(csa, 1, ks);
   1.506 +         j = iran(csa, 1, nsorc);
   1.507 +         isup[i] += ksp;
   1.508 +         isup[j] += ks - ksp;
   1.509 +      }
   1.510 +      j = iran(csa, 1, nsorc);
   1.511 +      isup[j] += itsup - ks * nsorc;
   1.512 +      return;
   1.513 +}
   1.514 +
   1.515 +/***********************************************************************
   1.516 +*  The routine chain adds node lpick to the end of the chain with source
   1.517 +*  node lsorc. */
   1.518 +
   1.519 +static void chain(struct csa *csa, int lpick, int lsorc)
   1.520 +{     int i, j, k, l, m;
   1.521 +      k = 0;
   1.522 +      m = ist;
   1.523 +      for (i = 1; i <= lpick; i++)
   1.524 +      {  l = k;
   1.525 +         k = m;
   1.526 +         m = ipred[k];
   1.527 +      }
   1.528 +      ipred[l] = m;
   1.529 +      j = ipred[lsorc];
   1.530 +      ipred[k] = j;
   1.531 +      ipred[lsorc] = k;
   1.532 +      return;
   1.533 +}
   1.534 +
   1.535 +/***********************************************************************
   1.536 +*  The routine chnarc puts the arcs in the chain from source lsorc into
   1.537 +*  the ihead and itail arrays for sorting. */
   1.538 +
   1.539 +static void chnarc(struct csa *csa, int lsorc)
   1.540 +{     int ito, ifrom;
   1.541 +      nsort = 0;
   1.542 +      ito = ipred[lsorc];
   1.543 +L10:  if (ito == lsorc) return;
   1.544 +      nsort++;
   1.545 +      ifrom = ipred[ito];
   1.546 +      ihead[nsort] = ito;
   1.547 +      itail[nsort] = ifrom;
   1.548 +      ito = ifrom;
   1.549 +      goto L10;
   1.550 +}
   1.551 +
   1.552 +/***********************************************************************
   1.553 +*  The routine sort sorts the nsort arcs in the ihead and itail arrays.
   1.554 +*  ihead is used as the sort key (i.e. forward star sort order). */
   1.555 +
   1.556 +static void sort(struct csa *csa)
   1.557 +{     int i, j, k, l, m, n, it;
   1.558 +      n = nsort;
   1.559 +      m = n;
   1.560 +L10:  m /= 2;
   1.561 +      if (m == 0) return;
   1.562 +      k = n - m;
   1.563 +      j = 1;
   1.564 +L20:  i = j;
   1.565 +L30:  l = i + m;
   1.566 +      if (itail[i] <= itail[l]) goto L40;
   1.567 +      it = itail[i];
   1.568 +      itail[i] = itail[l];
   1.569 +      itail[l] = it;
   1.570 +      it = ihead[i];
   1.571 +      ihead[i] = ihead[l];
   1.572 +      ihead[l] = it;
   1.573 +      i -= m;
   1.574 +      if (i >= 1) goto L30;
   1.575 +L40:  j++;
   1.576 +      if (j <= k) goto L20;
   1.577 +      goto L10;
   1.578 +}
   1.579 +
   1.580 +/***********************************************************************
   1.581 +*  The routine pickj creates a random number of arcs out of node 'it'.
   1.582 +*  Various parameters are dynamically adjusted in an attempt to ensure
   1.583 +*  that the generated network has the correct number of arcs. */
   1.584 +
   1.585 +static void pickj(struct csa *csa, int it)
   1.586 +{     int j, k, l, nn, nupbnd, icap, jcap, icost;
   1.587 +      if ((nodlft - 1) * 2 > iarcs - narcs - 1)
   1.588 +      {  nodlft--;
   1.589 +         return;
   1.590 +      }
   1.591 +      if ((iarcs - narcs + nonsor - ktl - 1) / nodlft - nonsor + 1 >= 0)
   1.592 +         k = nonsor;
   1.593 +      else
   1.594 +      {  nupbnd = (iarcs - narcs - nodlft) / nodlft * 2;
   1.595 +L40:     k = iran(csa, 1, nupbnd);
   1.596 +         if (nodlft == 1) k = iarcs - narcs;
   1.597 +         if ((nodlft - 1) * (nonsor - 1) < iarcs - narcs - k) goto L40;
   1.598 +      }
   1.599 +      nodlft--;
   1.600 +      for (j = 1; j <= k; j++)
   1.601 +      {  nn = iran(csa, 1, ktl);
   1.602 +         ktl--;
   1.603 +         for (l = nftsor; l <= nodes; l++)
   1.604 +         {  if (iflag[l] != 1)
   1.605 +            {  nn--;
   1.606 +               if (nn == 0) goto L70;
   1.607 +            }
   1.608 +         }
   1.609 +         return;
   1.610 +L70:     iflag[l] = 1;
   1.611 +         icap = itsup;
   1.612 +         jcap = iran(csa, 1, 100);
   1.613 +         if (jcap <= ipcap)
   1.614 +            icap = iran(csa, mincap, maxcap);
   1.615 +         icost = iran(csa, mincst, maxcst);
   1.616 +         if (G == NULL)
   1.617 +            xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, l, "", icost,
   1.618 +               icap);
   1.619 +         else
   1.620 +         {  glp_arc *a = glp_add_arc(G, it, l);
   1.621 +            if (a_cap >= 0)
   1.622 +            {  double temp = (double)icap;
   1.623 +               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
   1.624 +            }
   1.625 +            if (a_cost >= 0)
   1.626 +            {  double temp = (double)icost;
   1.627 +               memcpy((char *)a->data + a_cost, &temp, sizeof(double));
   1.628 +            }
   1.629 +         }
   1.630 +         narcs++;
   1.631 +      }
   1.632 +      return;
   1.633 +}
   1.634 +
   1.635 +/***********************************************************************
   1.636 +*  The routine assign generate assignment problems. It defines the unit
   1.637 +*  supplies, builds a skeleton, then calls pickj to create the arcs. */
   1.638 +
   1.639 +static void assign(struct csa *csa)
   1.640 +{     int i, it, nn, l, ll, icost;
   1.641 +      if (G == NULL)
   1.642 +         xprintf("SUPPLY\n");
   1.643 +      for (i = 1; i <= nsorc; i++)
   1.644 +      {  isup[i] = 1;
   1.645 +         iflag[i] = 0;
   1.646 +         if (G == NULL)
   1.647 +            xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
   1.648 +         else
   1.649 +         {  if (v_rhs >= 0)
   1.650 +            {  double temp = (double)isup[i];
   1.651 +               glp_vertex *v = G->v[i];
   1.652 +               memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
   1.653 +            }
   1.654 +         }
   1.655 +      }
   1.656 +      if (G == NULL)
   1.657 +         xprintf("ARCS\n");
   1.658 +      for (i = nfsink; i <= nodes; i++)
   1.659 +         ipred[i] = 1;
   1.660 +      for (it = 1; it <= nsorc; it++)
   1.661 +      {  for (i = nfsink; i <= nodes; i++)
   1.662 +            iflag[i] = 0;
   1.663 +         ktl = nsink - 1;
   1.664 +         nn = iran(csa, 1, nsink - it + 1);
   1.665 +         for (l = 1; l <= nsorc; l++)
   1.666 +         {  if (iflag[l] != 1)
   1.667 +            {  nn--;
   1.668 +               if (nn == 0) break;
   1.669 +            }
   1.670 +         }
   1.671 +         narcs++;
   1.672 +         ll = nsorc + l;
   1.673 +         icost = iran(csa, mincst, maxcst);
   1.674 +         if (G == NULL)
   1.675 +            xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ll, "", icost,
   1.676 +               isup[1]);
   1.677 +         else
   1.678 +         {  glp_arc *a = glp_add_arc(G, it, ll);
   1.679 +            if (a_cap >= 0)
   1.680 +            {  double temp = (double)isup[1];
   1.681 +               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
   1.682 +            }
   1.683 +            if (a_cost >= 0)
   1.684 +            {  double temp = (double)icost;
   1.685 +               memcpy((char *)a->data + a_cost, &temp, sizeof(double));
   1.686 +            }
   1.687 +         }
   1.688 +         iflag[l] = 1;
   1.689 +         iflag[ll] = 1;
   1.690 +         pickj(csa, it);
   1.691 +      }
   1.692 +      return;
   1.693 +}
   1.694 +
   1.695 +/***********************************************************************
   1.696 +*  Portable congruential (uniform) random number generator:
   1.697 +*
   1.698 +*     next_value = ((7**5) * previous_value) modulo ((2**31)-1)
   1.699 +*
   1.700 +*  This generator consists of three routines:
   1.701 +*
   1.702 +*  (1) setran - initializes constants and seed
   1.703 +*  (2) iran   - generates an integer random number
   1.704 +*  (3) rran   - generates a real random number
   1.705 +*
   1.706 +*  The generator requires a machine with at least 32 bits of precision.
   1.707 +*  The seed (iseed) must be in the range [1,(2**31)-1]. */
   1.708 +
   1.709 +static void setran(struct csa *csa, int iseed)
   1.710 +{     xassert(iseed >= 1);
   1.711 +      mult = 16807;
   1.712 +      modul = 2147483647;
   1.713 +      i15 = 1 << 15;
   1.714 +      i16 = 1 << 16;
   1.715 +      jran = iseed;
   1.716 +      return;
   1.717 +}
   1.718 +
   1.719 +/***********************************************************************
   1.720 +*  The routine iran generates an integer random number between ilow and
   1.721 +*  ihigh. If ilow > ihigh then iran returns ihigh. */
   1.722 +
   1.723 +static int iran(struct csa *csa, int ilow, int ihigh)
   1.724 +{     int ixhi, ixlo, ixalo, leftlo, ixahi, ifulhi, irtlo, iover,
   1.725 +         irthi, j;
   1.726 +      ixhi = jran / i16;
   1.727 +      ixlo = jran - ixhi * i16;
   1.728 +      ixalo = ixlo * mult;
   1.729 +      leftlo = ixalo / i16;
   1.730 +      ixahi = ixhi * mult;
   1.731 +      ifulhi = ixahi + leftlo;
   1.732 +      irtlo = ixalo - leftlo * i16;
   1.733 +      iover = ifulhi / i15;
   1.734 +      irthi = ifulhi - iover * i15;
   1.735 +      jran = ((irtlo - modul) + irthi * i16) + iover;
   1.736 +      if (jran < 0) jran += modul;
   1.737 +      j = ihigh - ilow + 1;
   1.738 +      if (j > 0)
   1.739 +         return jran % j + ilow;
   1.740 +      else
   1.741 +         return ihigh;
   1.742 +}
   1.743 +
   1.744 +/**********************************************************************/
   1.745 +
   1.746 +#if 0
   1.747 +static int scan(char card[80+1], int pos, int len)
   1.748 +{     char buf[10+1];
   1.749 +      memcpy(buf, &card[pos-1], len);
   1.750 +      buf[len] = '\0';
   1.751 +      return atoi(buf);
   1.752 +}
   1.753 +
   1.754 +int main(void)
   1.755 +{     int parm[1+15];
   1.756 +      char card[80+1];
   1.757 +      xassert(fgets(card, sizeof(card), stdin) == card);
   1.758 +      parm[1] = scan(card, 1, 8);
   1.759 +      parm[2] = scan(card, 9, 8);
   1.760 +      xassert(fgets(card, sizeof(card), stdin) == card);
   1.761 +      parm[3] = scan(card, 1, 5);
   1.762 +      parm[4] = scan(card, 6, 5);
   1.763 +      parm[5] = scan(card, 11, 5);
   1.764 +      parm[6] = scan(card, 16, 5);
   1.765 +      parm[7] = scan(card, 21, 5);
   1.766 +      parm[8] = scan(card, 26, 5);
   1.767 +      parm[9] = scan(card, 31, 10);
   1.768 +      parm[10] = scan(card, 41, 5);
   1.769 +      parm[11] = scan(card, 46, 5);
   1.770 +      parm[12] = scan(card, 51, 5);
   1.771 +      parm[13] = scan(card, 56, 5);
   1.772 +      parm[14] = scan(card, 61, 10);
   1.773 +      parm[15] = scan(card, 71, 10);
   1.774 +      glp_netgen(NULL, 0, 0, 0, parm);
   1.775 +      return 0;
   1.776 +}
   1.777 +#endif
   1.778 +
   1.779 +/* eof */