alpar@1: /* glpnet03.c (Klingman's network problem generator) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  This code is the result of translation of the Fortran program NETGEN
alpar@1: *  developed by Dr. Darwin Klingman, which is publically available from
alpar@1: *  NETLIB at <http://www.netlib.org/lp/generators>.
alpar@1: *
alpar@1: *  The translation was made by Andrew Makhorin <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpapi.h"
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_netgen - Klingman's network problem generator
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_netgen(glp_graph *G, int v_rhs, int a_cap, int a_cost,
alpar@1: *     const int parm[1+15]);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_netgen is a network problem generator developed by
alpar@1: *  Dr. Darwin Klingman. It can create capacitated and uncapacitated
alpar@1: *  minimum cost flow (or transshipment), transportation, and assignment
alpar@1: *  problems.
alpar@1: *
alpar@1: *  The parameter G specifies the graph object, to which the generated
alpar@1: *  problem data have to be stored. Note that on entry the graph object
alpar@1: *  is erased with the routine glp_erase_graph.
alpar@1: *
alpar@1: *  The parameter v_rhs specifies an offset of the field of type double
alpar@1: *  in the vertex data block, to which the routine stores the supply or
alpar@1: *  demand value. If v_rhs < 0, the value is not stored.
alpar@1: *
alpar@1: *  The parameter a_cap specifies an offset of the field of type double
alpar@1: *  in the arc data block, to which the routine stores the arc capacity.
alpar@1: *  If a_cap < 0, the capacity is not stored.
alpar@1: *
alpar@1: *  The parameter a_cost specifies an offset of the field of type double
alpar@1: *  in the arc data block, to which the routine stores the per-unit cost
alpar@1: *  if the arc flow. If a_cost < 0, the cost is not stored.
alpar@1: *
alpar@1: *  The array parm contains description of the network to be generated:
alpar@1: *
alpar@1: *  parm[0]           not used
alpar@1: *  parm[1]  (iseed)  8-digit positive random number seed
alpar@1: *  parm[2]  (nprob)  8-digit problem id number
alpar@1: *  parm[3]  (nodes)  total number of nodes
alpar@1: *  parm[4]  (nsorc)  total number of source nodes (including
alpar@1: *                    transshipment nodes)
alpar@1: *  parm[5]  (nsink)  total number of sink nodes (including
alpar@1: *                    transshipment nodes)
alpar@1: *  parm[6]  (iarcs)  number of arcs
alpar@1: *  parm[7]  (mincst) minimum cost for arcs
alpar@1: *  parm[8]  (maxcst) maximum cost for arcs
alpar@1: *  parm[9]  (itsup)  total supply
alpar@1: *  parm[10] (ntsorc) number of transshipment source nodes
alpar@1: *  parm[11] (ntsink) number of transshipment sink nodes
alpar@1: *  parm[12] (iphic)  percentage of skeleton arcs to be given
alpar@1: *                    the maximum cost
alpar@1: *  parm[13] (ipcap)  percentage of arcs to be capacitated
alpar@1: *  parm[14] (mincap) minimum upper bound for capacitated arcs
alpar@1: *  parm[15] (maxcap) maximum upper bound for capacitated arcs
alpar@1: *
alpar@1: *  The routine generates a transportation problem if:
alpar@1: *
alpar@1: *     nsorc + nsink = nodes, ntsorc = 0, and ntsink = 0.
alpar@1: *
alpar@1: *  The routine generates an assignment problem if the requirements for
alpar@1: *  a transportation problem are met and:
alpar@1: *
alpar@1: *     nsorc = nsink and itsup = nsorc.
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  If the instance was successfully generated, the routine glp_netgen
alpar@1: *  returns zero; otherwise, if specified parameters are inconsistent,
alpar@1: *  the routine returns a non-zero error code.
alpar@1: *
alpar@1: *  REFERENCES
alpar@1: *
alpar@1: *  D.Klingman, A.Napier, and J.Stutz. NETGEN: A program for generating
alpar@1: *  large scale capacitated assignment, transportation, and minimum cost
alpar@1: *  flow networks. Management Science 20 (1974), 814-20. */
alpar@1: 
alpar@1: struct csa
alpar@1: {     /* common storage area */
alpar@1:       glp_graph *G;
alpar@1:       int v_rhs, a_cap, a_cost;
alpar@1:       int nodes, iarcs, mincst, maxcst, itsup, nsorc, nsink, nonsor,
alpar@1:          nfsink, narcs, nsort, nftsor, ipcap, mincap, maxcap, ktl,
alpar@1:          nodlft, *ipred, *ihead, *itail, *iflag, *isup, *lsinks, mult,
alpar@1:          modul, i15, i16, jran;
alpar@1: };
alpar@1: 
alpar@1: #define G      (csa->G)
alpar@1: #define v_rhs  (csa->v_rhs)
alpar@1: #define a_cap  (csa->a_cap)
alpar@1: #define a_cost (csa->a_cost)
alpar@1: #define nodes  (csa->nodes)
alpar@1: #define iarcs  (csa->iarcs)
alpar@1: #define mincst (csa->mincst)
alpar@1: #define maxcst (csa->maxcst)
alpar@1: #define itsup  (csa->itsup)
alpar@1: #define nsorc  (csa->nsorc)
alpar@1: #define nsink  (csa->nsink)
alpar@1: #define nonsor (csa->nonsor)
alpar@1: #define nfsink (csa->nfsink)
alpar@1: #define narcs  (csa->narcs)
alpar@1: #define nsort  (csa->nsort)
alpar@1: #define nftsor (csa->nftsor)
alpar@1: #define ipcap  (csa->ipcap)
alpar@1: #define mincap (csa->mincap)
alpar@1: #define maxcap (csa->maxcap)
alpar@1: #define ktl    (csa->ktl)
alpar@1: #define nodlft (csa->nodlft)
alpar@1: #if 0
alpar@1: /* spent a day to find out this bug */
alpar@1: #define ist    (csa->ist)
alpar@1: #else
alpar@1: #define ist    (ipred[0])
alpar@1: #endif
alpar@1: #define ipred  (csa->ipred)
alpar@1: #define ihead  (csa->ihead)
alpar@1: #define itail  (csa->itail)
alpar@1: #define iflag  (csa->iflag)
alpar@1: #define isup   (csa->isup)
alpar@1: #define lsinks (csa->lsinks)
alpar@1: #define mult   (csa->mult)
alpar@1: #define modul  (csa->modul)
alpar@1: #define i15    (csa->i15)
alpar@1: #define i16    (csa->i16)
alpar@1: #define jran   (csa->jran)
alpar@1: 
alpar@1: static void cresup(struct csa *csa);
alpar@1: static void chain(struct csa *csa, int lpick, int lsorc);
alpar@1: static void chnarc(struct csa *csa, int lsorc);
alpar@1: static void sort(struct csa *csa);
alpar@1: static void pickj(struct csa *csa, int it);
alpar@1: static void assign(struct csa *csa);
alpar@1: static void setran(struct csa *csa, int iseed);
alpar@1: static int iran(struct csa *csa, int ilow, int ihigh);
alpar@1: 
alpar@1: int glp_netgen(glp_graph *G_, int _v_rhs, int _a_cap, int _a_cost,
alpar@1:       const int parm[1+15])
alpar@1: {     struct csa _csa, *csa = &_csa;
alpar@1:       int iseed, nprob, ntsorc, ntsink, iphic, i, nskel, nltr, ltsink,
alpar@1:          ntrans, npsink, nftr, npsorc, ntravl, ntrrem, lsorc, lpick,
alpar@1:          nsksr, nsrchn, j, item, l, ks, k, ksp, li, n, ii, it, ih, icap,
alpar@1:          jcap, icost, jcost, ret;
alpar@1:       G = G_;
alpar@1:       v_rhs = _v_rhs;
alpar@1:       a_cap = _a_cap;
alpar@1:       a_cost = _a_cost;
alpar@1:       if (G != NULL)
alpar@1:       {  if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
alpar@1:             xerror("glp_netgen: v_rhs = %d; invalid offset\n", v_rhs);
alpar@1:          if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
alpar@1:             xerror("glp_netgen: a_cap = %d; invalid offset\n", a_cap);
alpar@1:          if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
alpar@1:             xerror("glp_netgen: a_cost = %d; invalid offset\n", a_cost);
alpar@1:       }
alpar@1:       /* Input the user's random number seed and fix it if
alpar@1:          non-positive. */
alpar@1:       iseed = parm[1];
alpar@1:       nprob = parm[2];
alpar@1:       if (iseed <= 0) iseed = 13502460;
alpar@1:       setran(csa, iseed);
alpar@1:       /* Input the user's problem characteristics. */
alpar@1:       nodes = parm[3];
alpar@1:       nsorc = parm[4];
alpar@1:       nsink = parm[5];
alpar@1:       iarcs = parm[6];
alpar@1:       mincst = parm[7];
alpar@1:       maxcst = parm[8];
alpar@1:       itsup = parm[9];
alpar@1:       ntsorc = parm[10];
alpar@1:       ntsink = parm[11];
alpar@1:       iphic = parm[12];
alpar@1:       ipcap = parm[13];
alpar@1:       mincap = parm[14];
alpar@1:       maxcap = parm[15];
alpar@1:       /* Check the size of the problem. */
alpar@1:       if (!(10 <= nodes && nodes <= 100000))
alpar@1:       {  ret = 1;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* Check user supplied parameters for consistency. */
alpar@1:       if (!(nsorc >= 0 && nsink >= 0 && nsorc + nsink <= nodes))
alpar@1:       {  ret = 2;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (iarcs < 0)
alpar@1:       {  ret = 3;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (mincst > maxcst)
alpar@1:       {  ret = 4;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (itsup < 0)
alpar@1:       {  ret = 5;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (!(0 <= ntsorc && ntsorc <= nsorc))
alpar@1:       {  ret = 6;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (!(0 <= ntsink && ntsink <= nsink))
alpar@1:       {  ret = 7;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (!(0 <= iphic && iphic <= 100))
alpar@1:       {  ret = 8;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (!(0 <= ipcap && ipcap <= 100))
alpar@1:       {  ret = 9;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       if (mincap > maxcap)
alpar@1:       {  ret = 10;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* Initailize the graph object. */
alpar@1:       if (G != NULL)
alpar@1:       {  glp_erase_graph(G, G->v_size, G->a_size);
alpar@1:          glp_add_vertices(G, nodes);
alpar@1:          if (v_rhs >= 0)
alpar@1:          {  double zero = 0.0;
alpar@1:             for (i = 1; i <= nodes; i++)
alpar@1:             {  glp_vertex *v = G->v[i];
alpar@1:                memcpy((char *)v->data + v_rhs, &zero, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* Allocate working arrays. */
alpar@1:       ipred = xcalloc(1+nodes, sizeof(int));
alpar@1:       ihead = xcalloc(1+nodes, sizeof(int));
alpar@1:       itail = xcalloc(1+nodes, sizeof(int));
alpar@1:       iflag = xcalloc(1+nodes, sizeof(int));
alpar@1:       isup = xcalloc(1+nodes, sizeof(int));
alpar@1:       lsinks = xcalloc(1+nodes, sizeof(int));
alpar@1:       /* Print the problem documentation records. */
alpar@1:       if (G == NULL)
alpar@1:       {  xprintf("BEGIN\n");
alpar@1:          xprintf("NETGEN PROBLEM%8d%10s%10d NODES AND%10d ARCS\n",
alpar@1:             nprob, "", nodes, iarcs);
alpar@1:          xprintf("USER:%11d%11d%11d%11d%11d%11d\nDATA:%11d%11d%11d%11d%"
alpar@1:             "11d%11d\n", iseed, nsorc, nsink, mincst,
alpar@1:             maxcst, itsup, ntsorc, ntsink, iphic, ipcap,
alpar@1:             mincap, maxcap);
alpar@1:       }
alpar@1:       else
alpar@1:          glp_set_graph_name(G, "NETGEN");
alpar@1:       /* Set various constants used in the program. */
alpar@1:       narcs = 0;
alpar@1:       nskel = 0;
alpar@1:       nltr = nodes - nsink;
alpar@1:       ltsink = nltr + ntsink;
alpar@1:       ntrans = nltr - nsorc;
alpar@1:       nfsink = nltr + 1;
alpar@1:       nonsor = nodes - nsorc + ntsorc;
alpar@1:       npsink = nsink - ntsink;
alpar@1:       nodlft = nodes - nsink + ntsink;
alpar@1:       nftr = nsorc + 1;
alpar@1:       nftsor = nsorc - ntsorc + 1;
alpar@1:       npsorc = nsorc - ntsorc;
alpar@1:       /* Randomly distribute the supply among the source nodes. */
alpar@1:       if (npsorc + npsink == nodes && npsorc == npsink &&
alpar@1:           itsup == nsorc)
alpar@1:       {  assign(csa);
alpar@1:          nskel = nsorc;
alpar@1:          goto L390;
alpar@1:       }
alpar@1:       cresup(csa);
alpar@1:       /* Print the supply records. */
alpar@1:       if (G == NULL)
alpar@1:       {  xprintf("SUPPLY\n");
alpar@1:          for (i = 1; i <= nsorc; i++)
alpar@1:             xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
alpar@1:          xprintf("ARCS\n");
alpar@1:       }
alpar@1:       else
alpar@1:       {  if (v_rhs >= 0)
alpar@1:          {  for (i = 1; i <= nsorc; i++)
alpar@1:             {  double temp = (double)isup[i];
alpar@1:                glp_vertex *v = G->v[i];
alpar@1:                memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* Make the sources point to themselves in ipred array. */
alpar@1:       for (i = 1; i <= nsorc; i++)
alpar@1:          ipred[i] = i;
alpar@1:       if (ntrans == 0) goto L170;
alpar@1:       /* Chain the transshipment nodes together in the ipred array. */
alpar@1:       ist = nftr;
alpar@1:       ipred[nltr] = 0;
alpar@1:       for (i = nftr; i < nltr; i++)
alpar@1:          ipred[i] = i+1;
alpar@1:       /* Form even length chains for 60 percent of the transshipments.*/
alpar@1:       ntravl = 6 * ntrans / 10;
alpar@1:       ntrrem = ntrans - ntravl;
alpar@1: L140: lsorc = 1;
alpar@1:       while (ntravl != 0)
alpar@1:       {  lpick = iran(csa, 1, ntravl + ntrrem);
alpar@1:          ntravl--;
alpar@1:          chain(csa, lpick, lsorc);
alpar@1:          if (lsorc == nsorc) goto L140;
alpar@1:          lsorc++;
alpar@1:       }
alpar@1:       /* Add the remaining transshipments to the chains. */
alpar@1:       while (ntrrem != 0)
alpar@1:       {
alpar@1:          lpick = iran(csa, 1, ntrrem);
alpar@1:          ntrrem--;
alpar@1:          lsorc = iran(csa, 1, nsorc);
alpar@1:          chain(csa, lpick, lsorc);
alpar@1:       }
alpar@1: L170: /* Set all demands equal to zero. */
alpar@1:       for (i = nfsink; i <= nodes; i++)
alpar@1:          ipred[i] = 0;
alpar@1:       /* The following loop takes one chain at a time (through the use
alpar@1:          of logic contained in the loop and calls to other routines) and
alpar@1:          creates the remaining network arcs. */
alpar@1:       for (lsorc = 1; lsorc <= nsorc; lsorc++)
alpar@1:       {  chnarc(csa, lsorc);
alpar@1:          for (i = nfsink; i <= nodes; i++)
alpar@1:             iflag[i] = 0;
alpar@1:          /* Choose the number of sinks to be hooked up to the current
alpar@1:             chain. */
alpar@1:          if (ntrans != 0)
alpar@1:             nsksr = (nsort * 2 * nsink) / ntrans;
alpar@1:          else
alpar@1:             nsksr = nsink / nsorc + 1;
alpar@1:          if (nsksr < 2) nsksr = 2;
alpar@1:          if (nsksr > nsink) nsksr = nsink;
alpar@1:          nsrchn = nsort;
alpar@1:          /* Randomly pick nsksr sinks and put their names in lsinks. */
alpar@1:          ktl = nsink;
alpar@1:          for (j = 1; j <= nsksr; j++)
alpar@1:          {  item = iran(csa, 1, ktl);
alpar@1:             ktl--;
alpar@1:             for (l = nfsink; l <= nodes; l++)
alpar@1:             {  if (iflag[l] != 1)
alpar@1:                {  item--;
alpar@1:                   if (item == 0) goto L230;
alpar@1:                }
alpar@1:             }
alpar@1:             break;
alpar@1: L230:       lsinks[j] = l;
alpar@1:             iflag[l] = 1;
alpar@1:          }
alpar@1:          /* If last source chain, add all sinks with zero demand to
alpar@1:             lsinks list. */
alpar@1:          if (lsorc == nsorc)
alpar@1:          {  for (j = nfsink; j <= nodes; j++)
alpar@1:             {  if (ipred[j] == 0 && iflag[j] != 1)
alpar@1:                {  nsksr++;
alpar@1:                   lsinks[nsksr] = j;
alpar@1:                   iflag[j] = 1;
alpar@1:                }
alpar@1:             }
alpar@1:          }
alpar@1:          /* Create demands for group of sinks in lsinks. */
alpar@1:          ks = isup[lsorc] / nsksr;
alpar@1:          k = ipred[lsorc];
alpar@1:          for (i = 1; i <= nsksr; i++)
alpar@1:          {  nsort++;
alpar@1:             ksp = iran(csa, 1, ks);
alpar@1:             j = iran(csa, 1, nsksr);
alpar@1:             itail[nsort] = k;
alpar@1:             li = lsinks[i];
alpar@1:             ihead[nsort] = li;
alpar@1:             ipred[li] += ksp;
alpar@1:             li = lsinks[j];
alpar@1:             ipred[li] += ks - ksp;
alpar@1:             n = iran(csa, 1, nsrchn);
alpar@1:             k = lsorc;
alpar@1:             for (ii = 1; ii <= n; ii++)
alpar@1:                k = ipred[k];
alpar@1:          }
alpar@1:          li = lsinks[1];
alpar@1:          ipred[li] += isup[lsorc] - ks * nsksr;
alpar@1:          nskel += nsort;
alpar@1:          /* Sort the arcs in the chain from source lsorc using itail as
alpar@1:             sort key. */
alpar@1:          sort(csa);
alpar@1:          /* Print this part of skeleton and create the arcs for these
alpar@1:             nodes. */
alpar@1:          i = 1;
alpar@1:          itail[nsort+1] = 0;
alpar@1: L300:    for (j = nftsor; j <= nodes; j++)
alpar@1:             iflag[j] = 0;
alpar@1:          ktl = nonsor - 1;
alpar@1:          it = itail[i];
alpar@1:          iflag[it] = 1;
alpar@1: L320:    ih = ihead[i];
alpar@1:          iflag[ih] = 1;
alpar@1:          narcs++;
alpar@1:          ktl--;
alpar@1:          /* Determine if this skeleton arc should be capacitated. */
alpar@1:          icap = itsup;
alpar@1:          jcap = iran(csa, 1, 100);
alpar@1:          if (jcap <= ipcap)
alpar@1:          {  icap = isup[lsorc];
alpar@1:             if (mincap > icap) icap = mincap;
alpar@1:          }
alpar@1:          /* Determine if this skeleton arc should have the maximum
alpar@1:             cost. */
alpar@1:          icost = maxcst;
alpar@1:          jcost = iran(csa, 1, 100);
alpar@1:          if (jcost > iphic)
alpar@1:             icost = iran(csa, mincst, maxcst);
alpar@1:          if (G == NULL)
alpar@1:             xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ih, "", icost,
alpar@1:                icap);
alpar@1:          else
alpar@1:          {  glp_arc *a = glp_add_arc(G, it, ih);
alpar@1:             if (a_cap >= 0)
alpar@1:             {  double temp = (double)icap;
alpar@1:                memcpy((char *)a->data + a_cap, &temp, sizeof(double));
alpar@1:             }
alpar@1:             if (a_cost >= 0)
alpar@1:             {  double temp = (double)icost;
alpar@1:                memcpy((char *)a->data + a_cost, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:          i++;
alpar@1:          if (itail[i] == it) goto L320;
alpar@1:          pickj(csa, it);
alpar@1:          if (i <= nsort) goto L300;
alpar@1:       }
alpar@1:       /* Create arcs from the transshipment sinks. */
alpar@1:       if (ntsink != 0)
alpar@1:       {  for (i = nfsink; i <= ltsink; i++)
alpar@1:          {  for (j = nftsor; j <= nodes; j++)
alpar@1:                iflag[j] = 0;
alpar@1:             ktl = nonsor - 1;
alpar@1:             iflag[i] = 1;
alpar@1:             pickj(csa, i);
alpar@1:          }
alpar@1:       }
alpar@1: L390: /* Print the demand records and end record. */
alpar@1:       if (G == NULL)
alpar@1:       {  xprintf("DEMAND\n");
alpar@1:          for (i = nfsink; i <= nodes; i++)
alpar@1:             xprintf("%6s%6d%18s%10d\n", "", i, "", ipred[i]);
alpar@1:          xprintf("END\n");
alpar@1:       }
alpar@1:       else
alpar@1:       {  if (v_rhs >= 0)
alpar@1:          {  for (i = nfsink; i <= nodes; i++)
alpar@1:             {  double temp = - (double)ipred[i];
alpar@1:                glp_vertex *v = G->v[i];
alpar@1:                memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* Free working arrays. */
alpar@1:       xfree(ipred);
alpar@1:       xfree(ihead);
alpar@1:       xfree(itail);
alpar@1:       xfree(iflag);
alpar@1:       xfree(isup);
alpar@1:       xfree(lsinks);
alpar@1:       /* The instance has been successfully generated. */
alpar@1:       ret = 0;
alpar@1: done: return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The routine cresup randomly distributes the total supply among the
alpar@1: *  source nodes. */
alpar@1: 
alpar@1: static void cresup(struct csa *csa)
alpar@1: {     int i, j, ks, ksp;
alpar@1:       xassert(itsup > nsorc);
alpar@1:       ks = itsup / nsorc;
alpar@1:       for (i = 1; i <= nsorc; i++)
alpar@1:          isup[i] = 0;
alpar@1:       for (i = 1; i <= nsorc; i++)
alpar@1:       {  ksp = iran(csa, 1, ks);
alpar@1:          j = iran(csa, 1, nsorc);
alpar@1:          isup[i] += ksp;
alpar@1:          isup[j] += ks - ksp;
alpar@1:       }
alpar@1:       j = iran(csa, 1, nsorc);
alpar@1:       isup[j] += itsup - ks * nsorc;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The routine chain adds node lpick to the end of the chain with source
alpar@1: *  node lsorc. */
alpar@1: 
alpar@1: static void chain(struct csa *csa, int lpick, int lsorc)
alpar@1: {     int i, j, k, l, m;
alpar@1:       k = 0;
alpar@1:       m = ist;
alpar@1:       for (i = 1; i <= lpick; i++)
alpar@1:       {  l = k;
alpar@1:          k = m;
alpar@1:          m = ipred[k];
alpar@1:       }
alpar@1:       ipred[l] = m;
alpar@1:       j = ipred[lsorc];
alpar@1:       ipred[k] = j;
alpar@1:       ipred[lsorc] = k;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The routine chnarc puts the arcs in the chain from source lsorc into
alpar@1: *  the ihead and itail arrays for sorting. */
alpar@1: 
alpar@1: static void chnarc(struct csa *csa, int lsorc)
alpar@1: {     int ito, ifrom;
alpar@1:       nsort = 0;
alpar@1:       ito = ipred[lsorc];
alpar@1: L10:  if (ito == lsorc) return;
alpar@1:       nsort++;
alpar@1:       ifrom = ipred[ito];
alpar@1:       ihead[nsort] = ito;
alpar@1:       itail[nsort] = ifrom;
alpar@1:       ito = ifrom;
alpar@1:       goto L10;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The routine sort sorts the nsort arcs in the ihead and itail arrays.
alpar@1: *  ihead is used as the sort key (i.e. forward star sort order). */
alpar@1: 
alpar@1: static void sort(struct csa *csa)
alpar@1: {     int i, j, k, l, m, n, it;
alpar@1:       n = nsort;
alpar@1:       m = n;
alpar@1: L10:  m /= 2;
alpar@1:       if (m == 0) return;
alpar@1:       k = n - m;
alpar@1:       j = 1;
alpar@1: L20:  i = j;
alpar@1: L30:  l = i + m;
alpar@1:       if (itail[i] <= itail[l]) goto L40;
alpar@1:       it = itail[i];
alpar@1:       itail[i] = itail[l];
alpar@1:       itail[l] = it;
alpar@1:       it = ihead[i];
alpar@1:       ihead[i] = ihead[l];
alpar@1:       ihead[l] = it;
alpar@1:       i -= m;
alpar@1:       if (i >= 1) goto L30;
alpar@1: L40:  j++;
alpar@1:       if (j <= k) goto L20;
alpar@1:       goto L10;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The routine pickj creates a random number of arcs out of node 'it'.
alpar@1: *  Various parameters are dynamically adjusted in an attempt to ensure
alpar@1: *  that the generated network has the correct number of arcs. */
alpar@1: 
alpar@1: static void pickj(struct csa *csa, int it)
alpar@1: {     int j, k, l, nn, nupbnd, icap, jcap, icost;
alpar@1:       if ((nodlft - 1) * 2 > iarcs - narcs - 1)
alpar@1:       {  nodlft--;
alpar@1:          return;
alpar@1:       }
alpar@1:       if ((iarcs - narcs + nonsor - ktl - 1) / nodlft - nonsor + 1 >= 0)
alpar@1:          k = nonsor;
alpar@1:       else
alpar@1:       {  nupbnd = (iarcs - narcs - nodlft) / nodlft * 2;
alpar@1: L40:     k = iran(csa, 1, nupbnd);
alpar@1:          if (nodlft == 1) k = iarcs - narcs;
alpar@1:          if ((nodlft - 1) * (nonsor - 1) < iarcs - narcs - k) goto L40;
alpar@1:       }
alpar@1:       nodlft--;
alpar@1:       for (j = 1; j <= k; j++)
alpar@1:       {  nn = iran(csa, 1, ktl);
alpar@1:          ktl--;
alpar@1:          for (l = nftsor; l <= nodes; l++)
alpar@1:          {  if (iflag[l] != 1)
alpar@1:             {  nn--;
alpar@1:                if (nn == 0) goto L70;
alpar@1:             }
alpar@1:          }
alpar@1:          return;
alpar@1: L70:     iflag[l] = 1;
alpar@1:          icap = itsup;
alpar@1:          jcap = iran(csa, 1, 100);
alpar@1:          if (jcap <= ipcap)
alpar@1:             icap = iran(csa, mincap, maxcap);
alpar@1:          icost = iran(csa, mincst, maxcst);
alpar@1:          if (G == NULL)
alpar@1:             xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, l, "", icost,
alpar@1:                icap);
alpar@1:          else
alpar@1:          {  glp_arc *a = glp_add_arc(G, it, l);
alpar@1:             if (a_cap >= 0)
alpar@1:             {  double temp = (double)icap;
alpar@1:                memcpy((char *)a->data + a_cap, &temp, sizeof(double));
alpar@1:             }
alpar@1:             if (a_cost >= 0)
alpar@1:             {  double temp = (double)icost;
alpar@1:                memcpy((char *)a->data + a_cost, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:          narcs++;
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The routine assign generate assignment problems. It defines the unit
alpar@1: *  supplies, builds a skeleton, then calls pickj to create the arcs. */
alpar@1: 
alpar@1: static void assign(struct csa *csa)
alpar@1: {     int i, it, nn, l, ll, icost;
alpar@1:       if (G == NULL)
alpar@1:          xprintf("SUPPLY\n");
alpar@1:       for (i = 1; i <= nsorc; i++)
alpar@1:       {  isup[i] = 1;
alpar@1:          iflag[i] = 0;
alpar@1:          if (G == NULL)
alpar@1:             xprintf("%6s%6d%18s%10d\n", "", i, "", isup[i]);
alpar@1:          else
alpar@1:          {  if (v_rhs >= 0)
alpar@1:             {  double temp = (double)isup[i];
alpar@1:                glp_vertex *v = G->v[i];
alpar@1:                memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       if (G == NULL)
alpar@1:          xprintf("ARCS\n");
alpar@1:       for (i = nfsink; i <= nodes; i++)
alpar@1:          ipred[i] = 1;
alpar@1:       for (it = 1; it <= nsorc; it++)
alpar@1:       {  for (i = nfsink; i <= nodes; i++)
alpar@1:             iflag[i] = 0;
alpar@1:          ktl = nsink - 1;
alpar@1:          nn = iran(csa, 1, nsink - it + 1);
alpar@1:          for (l = 1; l <= nsorc; l++)
alpar@1:          {  if (iflag[l] != 1)
alpar@1:             {  nn--;
alpar@1:                if (nn == 0) break;
alpar@1:             }
alpar@1:          }
alpar@1:          narcs++;
alpar@1:          ll = nsorc + l;
alpar@1:          icost = iran(csa, mincst, maxcst);
alpar@1:          if (G == NULL)
alpar@1:             xprintf("%6s%6d%6d%2s%10d%10d\n", "", it, ll, "", icost,
alpar@1:                isup[1]);
alpar@1:          else
alpar@1:          {  glp_arc *a = glp_add_arc(G, it, ll);
alpar@1:             if (a_cap >= 0)
alpar@1:             {  double temp = (double)isup[1];
alpar@1:                memcpy((char *)a->data + a_cap, &temp, sizeof(double));
alpar@1:             }
alpar@1:             if (a_cost >= 0)
alpar@1:             {  double temp = (double)icost;
alpar@1:                memcpy((char *)a->data + a_cost, &temp, sizeof(double));
alpar@1:             }
alpar@1:          }
alpar@1:          iflag[l] = 1;
alpar@1:          iflag[ll] = 1;
alpar@1:          pickj(csa, it);
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  Portable congruential (uniform) random number generator:
alpar@1: *
alpar@1: *     next_value = ((7**5) * previous_value) modulo ((2**31)-1)
alpar@1: *
alpar@1: *  This generator consists of three routines:
alpar@1: *
alpar@1: *  (1) setran - initializes constants and seed
alpar@1: *  (2) iran   - generates an integer random number
alpar@1: *  (3) rran   - generates a real random number
alpar@1: *
alpar@1: *  The generator requires a machine with at least 32 bits of precision.
alpar@1: *  The seed (iseed) must be in the range [1,(2**31)-1]. */
alpar@1: 
alpar@1: static void setran(struct csa *csa, int iseed)
alpar@1: {     xassert(iseed >= 1);
alpar@1:       mult = 16807;
alpar@1:       modul = 2147483647;
alpar@1:       i15 = 1 << 15;
alpar@1:       i16 = 1 << 16;
alpar@1:       jran = iseed;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The routine iran generates an integer random number between ilow and
alpar@1: *  ihigh. If ilow > ihigh then iran returns ihigh. */
alpar@1: 
alpar@1: static int iran(struct csa *csa, int ilow, int ihigh)
alpar@1: {     int ixhi, ixlo, ixalo, leftlo, ixahi, ifulhi, irtlo, iover,
alpar@1:          irthi, j;
alpar@1:       ixhi = jran / i16;
alpar@1:       ixlo = jran - ixhi * i16;
alpar@1:       ixalo = ixlo * mult;
alpar@1:       leftlo = ixalo / i16;
alpar@1:       ixahi = ixhi * mult;
alpar@1:       ifulhi = ixahi + leftlo;
alpar@1:       irtlo = ixalo - leftlo * i16;
alpar@1:       iover = ifulhi / i15;
alpar@1:       irthi = ifulhi - iover * i15;
alpar@1:       jran = ((irtlo - modul) + irthi * i16) + iover;
alpar@1:       if (jran < 0) jran += modul;
alpar@1:       j = ihigh - ilow + 1;
alpar@1:       if (j > 0)
alpar@1:          return jran % j + ilow;
alpar@1:       else
alpar@1:          return ihigh;
alpar@1: }
alpar@1: 
alpar@1: /**********************************************************************/
alpar@1: 
alpar@1: #if 0
alpar@1: static int scan(char card[80+1], int pos, int len)
alpar@1: {     char buf[10+1];
alpar@1:       memcpy(buf, &card[pos-1], len);
alpar@1:       buf[len] = '\0';
alpar@1:       return atoi(buf);
alpar@1: }
alpar@1: 
alpar@1: int main(void)
alpar@1: {     int parm[1+15];
alpar@1:       char card[80+1];
alpar@1:       xassert(fgets(card, sizeof(card), stdin) == card);
alpar@1:       parm[1] = scan(card, 1, 8);
alpar@1:       parm[2] = scan(card, 9, 8);
alpar@1:       xassert(fgets(card, sizeof(card), stdin) == card);
alpar@1:       parm[3] = scan(card, 1, 5);
alpar@1:       parm[4] = scan(card, 6, 5);
alpar@1:       parm[5] = scan(card, 11, 5);
alpar@1:       parm[6] = scan(card, 16, 5);
alpar@1:       parm[7] = scan(card, 21, 5);
alpar@1:       parm[8] = scan(card, 26, 5);
alpar@1:       parm[9] = scan(card, 31, 10);
alpar@1:       parm[10] = scan(card, 41, 5);
alpar@1:       parm[11] = scan(card, 46, 5);
alpar@1:       parm[12] = scan(card, 51, 5);
alpar@1:       parm[13] = scan(card, 56, 5);
alpar@1:       parm[14] = scan(card, 61, 10);
alpar@1:       parm[15] = scan(card, 71, 10);
alpar@1:       glp_netgen(NULL, 0, 0, 0, parm);
alpar@1:       return 0;
alpar@1: }
alpar@1: #endif
alpar@1: 
alpar@1: /* eof */