1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnet03.c Mon Dec 06 13:09:21 2010 +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 */