1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnet04.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,768 @@
     1.4 +/* glpnet04.c (grid-like 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 a modified version of the program GRIDGEN, a grid-like
    1.10 +*  network problem generator developed by Yusin Lee and Jim Orlin.
    1.11 +*  The original code is publically available on the DIMACS ftp site at:
    1.12 +*  <ftp://dimacs.rutgers.edu/pub/netflow/generators/network/gridgen>.
    1.13 +*
    1.14 +*  All changes concern only the program interface, so this modified
    1.15 +*  version produces exactly the same instances as the original version.
    1.16 +*
    1.17 +*  Changes were made by Andrew Makhorin <mao@gnu.org>.
    1.18 +*
    1.19 +*  GLPK is free software: you can redistribute it and/or modify it
    1.20 +*  under the terms of the GNU General Public License as published by
    1.21 +*  the Free Software Foundation, either version 3 of the License, or
    1.22 +*  (at your option) any later version.
    1.23 +*
    1.24 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.25 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.26 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.27 +*  License for more details.
    1.28 +*
    1.29 +*  You should have received a copy of the GNU General Public License
    1.30 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.31 +***********************************************************************/
    1.32 +
    1.33 +#include "glpapi.h"
    1.34 +
    1.35 +/***********************************************************************
    1.36 +*  NAME
    1.37 +*
    1.38 +*  glp_gridgen - grid-like network problem generator
    1.39 +*
    1.40 +*  SYNOPSIS
    1.41 +*
    1.42 +*  int glp_gridgen(glp_graph *G, int v_rhs, int a_cap, int a_cost,
    1.43 +*     const int parm[1+14]);
    1.44 +*
    1.45 +*  DESCRIPTION
    1.46 +*
    1.47 +*  The routine glp_gridgen is a grid-like network problem generator
    1.48 +*  developed by Yusin Lee and Jim Orlin.
    1.49 +*
    1.50 +*  The parameter G specifies the graph object, to which the generated
    1.51 +*  problem data have to be stored. Note that on entry the graph object
    1.52 +*  is erased with the routine glp_erase_graph.
    1.53 +*
    1.54 +*  The parameter v_rhs specifies an offset of the field of type double
    1.55 +*  in the vertex data block, to which the routine stores the supply or
    1.56 +*  demand value. If v_rhs < 0, the value is not stored.
    1.57 +*
    1.58 +*  The parameter a_cap specifies an offset of the field of type double
    1.59 +*  in the arc data block, to which the routine stores the arc capacity.
    1.60 +*  If a_cap < 0, the capacity is not stored.
    1.61 +*
    1.62 +*  The parameter a_cost specifies an offset of the field of type double
    1.63 +*  in the arc data block, to which the routine stores the per-unit cost
    1.64 +*  if the arc flow. If a_cost < 0, the cost is not stored.
    1.65 +*
    1.66 +*  The array parm contains description of the network to be generated:
    1.67 +*
    1.68 +*  parm[0]  not used
    1.69 +*  parm[1]  two-ways arcs indicator:
    1.70 +*           1 - if links in both direction should be generated
    1.71 +*           0 - otherwise
    1.72 +*  parm[2]  random number seed (a positive integer)
    1.73 +*  parm[3]  number of nodes (the number of nodes generated might be
    1.74 +*           slightly different to make the network a grid)
    1.75 +*  parm[4]  grid width
    1.76 +*  parm[5]  number of sources
    1.77 +*  parm[6]  number of sinks
    1.78 +*  parm[7]  average degree
    1.79 +*  parm[8]  total flow
    1.80 +*  parm[9]  distribution of arc costs:
    1.81 +*           1 - uniform
    1.82 +*           2 - exponential
    1.83 +*  parm[10] lower bound for arc cost (uniform)
    1.84 +*           100 * lambda (exponential)
    1.85 +*  parm[11] upper bound for arc cost (uniform)
    1.86 +*           not used (exponential)
    1.87 +*  parm[12] distribution of arc capacities:
    1.88 +*           1 - uniform
    1.89 +*           2 - exponential
    1.90 +*  parm[13] lower bound for arc capacity (uniform)
    1.91 +*           100 * lambda (exponential)
    1.92 +*  parm[14] upper bound for arc capacity (uniform)
    1.93 +*           not used (exponential)
    1.94 +*
    1.95 +*  RETURNS
    1.96 +*
    1.97 +*  If the instance was successfully generated, the routine glp_gridgen
    1.98 +*  returns zero; otherwise, if specified parameters are inconsistent,
    1.99 +*  the routine returns a non-zero error code.
   1.100 +*
   1.101 +*  COMMENTS
   1.102 +*
   1.103 +*  This network generator generates a grid-like network plus a super
   1.104 +*  node. In additional to the arcs connecting the nodes in the grid,
   1.105 +*  there is an arc from each supply node to the super node and from the
   1.106 +*  super node to each demand node to guarantee feasiblity. These arcs
   1.107 +*  have very high costs and very big capacities.
   1.108 +*
   1.109 +*  The idea of this network generator is as follows: First, a grid of
   1.110 +*  n1 * n2 is generated. For example, 5 * 3. The nodes are numbered as
   1.111 +*  1 to 15, and the supernode is numbered as n1*n2+1. Then arcs between
   1.112 +*  adjacent nodes are generated. For these arcs, the user is allowed to
   1.113 +*  specify either to generate two-way arcs or one-way arcs. If two-way
   1.114 +*  arcs are to be generated, two arcs, one in each direction, will be
   1.115 +*  generated between each adjacent node pairs. Otherwise, only one arc
   1.116 +*  will be generated. If this is the case, the arcs will be generated
   1.117 +*  in alterntive directions as shown below.
   1.118 +*
   1.119 +*      1 ---> 2 ---> 3 ---> 4 ---> 5
   1.120 +*      |      ^      |      ^      |
   1.121 +*      |      |      |      |      |
   1.122 +*      V      |      V      |      V
   1.123 +*      6 <--- 7 <--- 8 <--- 9 <--- 10
   1.124 +*      |      ^      |      ^      |
   1.125 +*      |      |      |      |      |
   1.126 +*      V      |      V      |      V
   1.127 +*     11 --->12 --->13 --->14 ---> 15
   1.128 +*
   1.129 +*  Then the arcs between the super node and the source/sink nodes are
   1.130 +*  added as mentioned before. If the number of arcs still doesn't reach
   1.131 +*  the requirement, additional arcs will be added by uniformly picking
   1.132 +*  random node pairs. There is no checking to prevent multiple arcs
   1.133 +*  between any pair of nodes. However, there will be no self-arcs (arcs
   1.134 +*  that poins back to its tail node) in the network.
   1.135 +*
   1.136 +*  The source and sink nodes are selected uniformly in the network, and
   1.137 +*  the imbalances of each source/sink node are also assigned by uniform
   1.138 +*  distribution. */
   1.139 +
   1.140 +struct stat_para
   1.141 +{     /* structure for statistical distributions */
   1.142 +      int distribution;
   1.143 +      /* the distribution: */
   1.144 +#define UNIFORM      1  /* uniform distribution */
   1.145 +#define EXPONENTIAL  2  /* exponential distribution */
   1.146 +      double parameter[5];
   1.147 +      /* the parameters of the distribution */
   1.148 +};
   1.149 +
   1.150 +struct arcs
   1.151 +{     int from;
   1.152 +      /* the FROM node of that arc */
   1.153 +      int to;
   1.154 +      /* the TO node of that arc */
   1.155 +      int cost;
   1.156 +      /* original cost of that arc */
   1.157 +      int u;
   1.158 +      /* capacity of the arc */
   1.159 +};
   1.160 +
   1.161 +struct imbalance
   1.162 +{     int node;
   1.163 +      /* Node ID */
   1.164 +      int supply;
   1.165 +      /* Supply of that node */
   1.166 +};
   1.167 +
   1.168 +struct csa
   1.169 +{     /* common storage area */
   1.170 +      glp_graph *G;
   1.171 +      int v_rhs, a_cap, a_cost;
   1.172 +      int seed;
   1.173 +      /* random number seed */
   1.174 +      int seed_original;
   1.175 +      /* the original seed from input */
   1.176 +      int two_way;
   1.177 +      /* 0: generate arcs in both direction for the basic grid, except
   1.178 +         for the arcs to/from the super node.  1: o/w */
   1.179 +      int n_node;
   1.180 +      /* total number of nodes in the network, numbered 1 to n_node,
   1.181 +         including the super node, which is the last one */
   1.182 +      int n_arc;
   1.183 +      /* total number of arcs in the network, counting EVERY arc. */
   1.184 +      int n_grid_arc;
   1.185 +      /* number of arcs in the basic grid, including the arcs to/from
   1.186 +         the super node */
   1.187 +      int n_source, n_sink;
   1.188 +      /* number of source and sink nodes */
   1.189 +      int avg_degree;
   1.190 +      /* average degree, arcs to and from the super node are counted */
   1.191 +      int t_supply;
   1.192 +      /* total supply in the network */
   1.193 +      int n1, n2;
   1.194 +      /* the two edges of the network grid.  n1 >= n2 */
   1.195 +      struct imbalance *source_list, *sink_list;
   1.196 +      /* head of the array of source/sink nodes */
   1.197 +      struct stat_para arc_costs;
   1.198 +      /* the distribution of arc costs */
   1.199 +      struct stat_para capacities;
   1.200 +      /* distribution of the capacities of the arcs */
   1.201 +      struct arcs *arc_list;
   1.202 +      /* head of the arc list array.  Arcs in this array are in the
   1.203 +         order of grid_arcs, arcs to/from super node, and other arcs */
   1.204 +};
   1.205 +
   1.206 +#define G (csa->G)
   1.207 +#define v_rhs (csa->v_rhs)
   1.208 +#define a_cap (csa->a_cap)
   1.209 +#define a_cost (csa->a_cost)
   1.210 +#define seed (csa->seed)
   1.211 +#define seed_original (csa->seed_original)
   1.212 +#define two_way (csa->two_way)
   1.213 +#define n_node (csa->n_node)
   1.214 +#define n_arc (csa->n_arc)
   1.215 +#define n_grid_arc (csa->n_grid_arc)
   1.216 +#define n_source (csa->n_source)
   1.217 +#define n_sink (csa->n_sink)
   1.218 +#define avg_degree (csa->avg_degree)
   1.219 +#define t_supply (csa->t_supply)
   1.220 +#define n1 (csa->n1)
   1.221 +#define n2 (csa->n2)
   1.222 +#define source_list (csa->source_list)
   1.223 +#define sink_list (csa->sink_list)
   1.224 +#define arc_costs (csa->arc_costs)
   1.225 +#define capacities (csa->capacities)
   1.226 +#define arc_list (csa->arc_list)
   1.227 +
   1.228 +static void assign_capacities(struct csa *csa);
   1.229 +static void assign_costs(struct csa *csa);
   1.230 +static void assign_imbalance(struct csa *csa);
   1.231 +static int exponential(struct csa *csa, double lambda[1]);
   1.232 +static struct arcs *gen_additional_arcs(struct csa *csa, struct arcs
   1.233 +      *arc_ptr);
   1.234 +static struct arcs *gen_basic_grid(struct csa *csa, struct arcs
   1.235 +      *arc_ptr);
   1.236 +static void gen_more_arcs(struct csa *csa, struct arcs *arc_ptr);
   1.237 +static void generate(struct csa *csa);
   1.238 +static void output(struct csa *csa);
   1.239 +static double randy(struct csa *csa);
   1.240 +static void select_source_sinks(struct csa *csa);
   1.241 +static int uniform(struct csa *csa, double a[2]);
   1.242 +
   1.243 +int glp_gridgen(glp_graph *G_, int _v_rhs, int _a_cap, int _a_cost,
   1.244 +      const int parm[1+14])
   1.245 +{     struct csa _csa, *csa = &_csa;
   1.246 +      int n, ret;
   1.247 +      G = G_;
   1.248 +      v_rhs = _v_rhs;
   1.249 +      a_cap = _a_cap;
   1.250 +      a_cost = _a_cost;
   1.251 +      if (G != NULL)
   1.252 +      {  if (v_rhs >= 0 && v_rhs > G->v_size - (int)sizeof(double))
   1.253 +            xerror("glp_gridgen: v_rhs = %d; invalid offset\n", v_rhs);
   1.254 +         if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
   1.255 +            xerror("glp_gridgen: a_cap = %d; invalid offset\n", a_cap);
   1.256 +         if (a_cost >= 0 && a_cost > G->a_size - (int)sizeof(double))
   1.257 +            xerror("glp_gridgen: a_cost = %d; invalid offset\n", a_cost)
   1.258 +               ;
   1.259 +      }
   1.260 +      /* Check the parameters for consistency. */
   1.261 +      if (!(parm[1] == 0 || parm[1] == 1))
   1.262 +      {  ret = 1;
   1.263 +         goto done;
   1.264 +      }
   1.265 +      if (parm[2] < 1)
   1.266 +      {  ret = 2;
   1.267 +         goto done;
   1.268 +      }
   1.269 +      if (!(10 <= parm[3] && parm[3] <= 40000))
   1.270 +      {  ret = 3;
   1.271 +         goto done;
   1.272 +      }
   1.273 +      if (!(1 <= parm[4] && parm[4] <= 40000))
   1.274 +      {  ret = 4;
   1.275 +         goto done;
   1.276 +      }
   1.277 +      if (!(parm[5] >= 0 && parm[6] >= 0 && parm[5] + parm[6] <=
   1.278 +         parm[3]))
   1.279 +      {  ret = 5;
   1.280 +         goto done;
   1.281 +      }
   1.282 +      if (!(1 <= parm[7] && parm[7] <= parm[3]))
   1.283 +      {  ret = 6;
   1.284 +         goto done;
   1.285 +      }
   1.286 +      if (parm[8] < 0)
   1.287 +      {  ret = 7;
   1.288 +         goto done;
   1.289 +      }
   1.290 +      if (!(parm[9] == 1 || parm[9] == 2))
   1.291 +      {  ret = 8;
   1.292 +         goto done;
   1.293 +      }
   1.294 +      if (parm[9] == 1 && parm[10] > parm[11] ||
   1.295 +          parm[9] == 2 && parm[10] < 1)
   1.296 +      {  ret = 9;
   1.297 +         goto done;
   1.298 +      }
   1.299 +      if (!(parm[12] == 1 || parm[12] == 2))
   1.300 +      {  ret = 10;
   1.301 +         goto done;
   1.302 +      }
   1.303 +      if (parm[12] == 1 && !(0 <= parm[13] && parm[13] <= parm[14]) ||
   1.304 +          parm[12] == 2 && parm[13] < 1)
   1.305 +      {  ret = 11;
   1.306 +         goto done;
   1.307 +      }
   1.308 +      /* Initialize the graph object. */
   1.309 +      if (G != NULL)
   1.310 +      {  glp_erase_graph(G, G->v_size, G->a_size);
   1.311 +         glp_set_graph_name(G, "GRIDGEN");
   1.312 +      }
   1.313 +      /* Copy the generator parameters. */
   1.314 +      two_way = parm[1];
   1.315 +      seed_original = seed = parm[2];
   1.316 +      n_node = parm[3];
   1.317 +      n = parm[4];
   1.318 +      n_source = parm[5];
   1.319 +      n_sink = parm[6];
   1.320 +      avg_degree = parm[7];
   1.321 +      t_supply = parm[8];
   1.322 +      arc_costs.distribution = parm[9];
   1.323 +      if (parm[9] == 1)
   1.324 +      {  arc_costs.parameter[0] = parm[10];
   1.325 +         arc_costs.parameter[1] = parm[11];
   1.326 +      }
   1.327 +      else
   1.328 +      {  arc_costs.parameter[0] = (double)parm[10] / 100.0;
   1.329 +         arc_costs.parameter[1] = 0.0;
   1.330 +      }
   1.331 +      capacities.distribution = parm[12];
   1.332 +      if (parm[12] == 1)
   1.333 +      {  capacities.parameter[0] = parm[13];
   1.334 +         capacities.parameter[1] = parm[14];
   1.335 +      }
   1.336 +      else
   1.337 +      {  capacities.parameter[0] = (double)parm[13] / 100.0;
   1.338 +         capacities.parameter[1] = 0.0;
   1.339 +      }
   1.340 +      /* Calculate the edge lengths of the grid according to the
   1.341 +         input. */
   1.342 +      if (n * n >= n_node)
   1.343 +      {  n1 = n;
   1.344 +         n2 = (int)((double)n_node / (double)n + 0.5);
   1.345 +      }
   1.346 +      else
   1.347 +      {  n2 = n;
   1.348 +         n1 = (int)((double)n_node / (double)n + 0.5);
   1.349 +      }
   1.350 +      /* Recalculate the total number of nodes and plus 1 for the super
   1.351 +         node. */
   1.352 +      n_node = n1 * n2 + 1;
   1.353 +      n_arc = n_node * avg_degree;
   1.354 +      n_grid_arc = (two_way + 1) * ((n1 - 1) * n2 + (n2 - 1) * n1) +
   1.355 +         n_source + n_sink;
   1.356 +      if (n_grid_arc > n_arc) n_arc = n_grid_arc;
   1.357 +      arc_list = xcalloc(n_arc, sizeof(struct arcs));
   1.358 +      source_list = xcalloc(n_source, sizeof(struct imbalance));
   1.359 +      sink_list = xcalloc(n_sink, sizeof(struct imbalance));
   1.360 +      /* Generate a random network. */
   1.361 +      generate(csa);
   1.362 +      /* Output the network. */
   1.363 +      output(csa);
   1.364 +      /* Free all allocated memory. */
   1.365 +      xfree(arc_list);
   1.366 +      xfree(source_list);
   1.367 +      xfree(sink_list);
   1.368 +      /* The instance has been successfully generated. */
   1.369 +      ret = 0;
   1.370 +done: return ret;
   1.371 +}
   1.372 +
   1.373 +#undef random
   1.374 +
   1.375 +static void assign_capacities(struct csa *csa)
   1.376 +{     /* Assign a capacity to each arc. */
   1.377 +      struct arcs *arc_ptr = arc_list;
   1.378 +      int (*random)(struct csa *csa, double *);
   1.379 +      int i;
   1.380 +      /* Determine the random number generator to use. */
   1.381 +      switch (arc_costs.distribution)
   1.382 +      {  case UNIFORM:
   1.383 +            random = uniform;
   1.384 +            break;
   1.385 +         case EXPONENTIAL:
   1.386 +            random = exponential;
   1.387 +            break;
   1.388 +         default:
   1.389 +            xassert(csa != csa);
   1.390 +      }
   1.391 +      /* Assign capacities to grid arcs. */
   1.392 +      for (i = n_source + n_sink; i < n_grid_arc; i++, arc_ptr++)
   1.393 +         arc_ptr->u = random(csa, capacities.parameter);
   1.394 +      i = i - n_source - n_sink;
   1.395 +      /* Assign capacities to arcs to/from supernode. */
   1.396 +      for (; i < n_grid_arc; i++, arc_ptr++)
   1.397 +         arc_ptr->u = t_supply;
   1.398 +      /* Assign capacities to all other arcs. */
   1.399 +      for (; i < n_arc; i++, arc_ptr++)
   1.400 +         arc_ptr->u = random(csa, capacities.parameter);
   1.401 +      return;
   1.402 +}
   1.403 +
   1.404 +static void assign_costs(struct csa *csa)
   1.405 +{     /* Assign a cost to each arc. */
   1.406 +      struct arcs *arc_ptr = arc_list;
   1.407 +      int (*random)(struct csa *csa, double *);
   1.408 +      int i;
   1.409 +      /* A high cost assigned to arcs to/from the supernode. */
   1.410 +      int high_cost;
   1.411 +      /* The maximum cost assigned to arcs in the base grid. */
   1.412 +      int max_cost = 0;
   1.413 +      /* Determine the random number generator to use. */
   1.414 +      switch (arc_costs.distribution)
   1.415 +      {  case UNIFORM:
   1.416 +            random = uniform;
   1.417 +            break;
   1.418 +         case EXPONENTIAL:
   1.419 +            random = exponential;
   1.420 +            break;
   1.421 +         default:
   1.422 +            xassert(csa != csa);
   1.423 +      }
   1.424 +      /* Assign costs to arcs in the base grid. */
   1.425 +      for (i = n_source + n_sink; i < n_grid_arc; i++, arc_ptr++)
   1.426 +      {  arc_ptr->cost = random(csa, arc_costs.parameter);
   1.427 +         if (max_cost < arc_ptr->cost) max_cost = arc_ptr->cost;
   1.428 +      }
   1.429 +      i = i - n_source - n_sink;
   1.430 +      /* Assign costs to arcs to/from the super node. */
   1.431 +      high_cost = max_cost * 2;
   1.432 +      for (; i < n_grid_arc; i++, arc_ptr++)
   1.433 +         arc_ptr->cost = high_cost;
   1.434 +      /* Assign costs to all other arcs. */
   1.435 +      for (; i < n_arc; i++, arc_ptr++)
   1.436 +         arc_ptr->cost = random(csa, arc_costs.parameter);
   1.437 +      return;
   1.438 +}
   1.439 +
   1.440 +static void assign_imbalance(struct csa *csa)
   1.441 +{     /* Assign an imbalance to each node. */
   1.442 +      int total, i;
   1.443 +      double avg;
   1.444 +      struct imbalance *ptr;
   1.445 +      /* assign the supply nodes */
   1.446 +      avg = 2.0 * t_supply / n_source;
   1.447 +      do
   1.448 +      {  for (i = 1, total = t_supply, ptr = source_list + 1;
   1.449 +            i < n_source; i++, ptr++)
   1.450 +         {  ptr->supply = (int)(randy(csa) * avg + 0.5);
   1.451 +            total -= ptr->supply;
   1.452 +         }
   1.453 +         source_list->supply = total;
   1.454 +      }
   1.455 +      /* redo all if the assignment "overshooted" */
   1.456 +      while (total <= 0);
   1.457 +      /* assign the demand nodes */
   1.458 +      avg = -2.0 * t_supply / n_sink;
   1.459 +      do
   1.460 +      {  for (i = 1, total = t_supply, ptr = sink_list + 1;
   1.461 +            i < n_sink; i++, ptr++)
   1.462 +         {  ptr->supply = (int)(randy(csa) * avg - 0.5);
   1.463 +            total += ptr->supply;
   1.464 +         }
   1.465 +         sink_list->supply = - total;
   1.466 +      }
   1.467 +      while (total <= 0);
   1.468 +      return;
   1.469 +}
   1.470 +
   1.471 +static int exponential(struct csa *csa, double lambda[1])
   1.472 +{     /* Returns an "exponentially distributed" integer with parameter
   1.473 +         lambda. */
   1.474 +      return ((int)(- lambda[0] * log((double)randy(csa)) + 0.5));
   1.475 +}
   1.476 +
   1.477 +static struct arcs *gen_additional_arcs(struct csa *csa, struct arcs
   1.478 +      *arc_ptr)
   1.479 +{     /* Generate an arc from each source to the supernode and from
   1.480 +         supernode to each sink. */
   1.481 +      int i;
   1.482 +      for (i = 0; i < n_source; i++, arc_ptr++)
   1.483 +      {  arc_ptr->from = source_list[i].node;
   1.484 +         arc_ptr->to = n_node;
   1.485 +      }
   1.486 +      for (i = 0; i < n_sink; i++, arc_ptr++)
   1.487 +      {  arc_ptr->to = sink_list[i].node;
   1.488 +         arc_ptr->from = n_node;
   1.489 +      }
   1.490 +      return arc_ptr;
   1.491 +}
   1.492 +
   1.493 +static struct arcs *gen_basic_grid(struct csa *csa, struct arcs
   1.494 +      *arc_ptr)
   1.495 +{     /* Generate the basic grid. */
   1.496 +      int direction = 1, i, j, k;
   1.497 +      if (two_way)
   1.498 +      {  /* Generate an arc in each direction. */
   1.499 +         for (i = 1; i < n_node; i += n1)
   1.500 +         {  for (j = i, k = j + n1 - 1; j < k; j++)
   1.501 +            {  arc_ptr->from = j;
   1.502 +               arc_ptr->to = j + 1;
   1.503 +               arc_ptr++;
   1.504 +               arc_ptr->from = j + 1;
   1.505 +               arc_ptr->to = j;
   1.506 +               arc_ptr++;
   1.507 +            }
   1.508 +         }
   1.509 +         for (i = 1; i <= n1; i++)
   1.510 +         {  for (j = i + n1; j < n_node; j += n1)
   1.511 +            {  arc_ptr->from = j;
   1.512 +               arc_ptr->to = j - n1;
   1.513 +               arc_ptr++;
   1.514 +               arc_ptr->from = j - n1;
   1.515 +               arc_ptr->to = j;
   1.516 +               arc_ptr++;
   1.517 +            }
   1.518 +         }
   1.519 +      }
   1.520 +      else
   1.521 +      {  /* Generate one arc in each direction. */
   1.522 +         for (i = 1; i < n_node; i += n1)
   1.523 +         {  if (direction == 1)
   1.524 +               j = i;
   1.525 +            else
   1.526 +               j = i + 1;
   1.527 +            for (k = j + n1 - 1; j < k; j++)
   1.528 +            {  arc_ptr->from = j;
   1.529 +               arc_ptr->to = j + direction;
   1.530 +               arc_ptr++;
   1.531 +            }
   1.532 +            direction = - direction;
   1.533 +         }
   1.534 +         for (i = 1; i <= n1; i++)
   1.535 +         {  j = i + n1;
   1.536 +            if (direction == 1)
   1.537 +            {  for (; j < n_node; j += n1)
   1.538 +               {  arc_ptr->from = j - n1;
   1.539 +                  arc_ptr->to = j;
   1.540 +                  arc_ptr++;
   1.541 +               }
   1.542 +            }
   1.543 +            else
   1.544 +            {  for (; j < n_node; j += n1)
   1.545 +               {  arc_ptr->from = j - n1;
   1.546 +                  arc_ptr->to = j;
   1.547 +                  arc_ptr++;
   1.548 +               }
   1.549 +            }
   1.550 +            direction = - direction;
   1.551 +         }
   1.552 +      }
   1.553 +      return arc_ptr;
   1.554 +}
   1.555 +
   1.556 +static void gen_more_arcs(struct csa *csa, struct arcs *arc_ptr)
   1.557 +{     /* Generate random arcs to meet the specified density. */
   1.558 +      int i;
   1.559 +      double ab[2];
   1.560 +      ab[0] = 0.9;
   1.561 +      ab[1] = n_node - 0.99;  /* upper limit is n_node-1 because the
   1.562 +                                 supernode cannot be selected */
   1.563 +      for (i = n_grid_arc; i < n_arc; i++, arc_ptr++)
   1.564 +      {  arc_ptr->from = uniform(csa, ab);
   1.565 +         arc_ptr->to = uniform(csa, ab);
   1.566 +         if (arc_ptr->from == arc_ptr->to)
   1.567 +         {  arc_ptr--;
   1.568 +            i--;
   1.569 +         }
   1.570 +      }
   1.571 +      return;
   1.572 +}
   1.573 +
   1.574 +static void generate(struct csa *csa)
   1.575 +{     /* Generate a random network. */
   1.576 +      struct arcs *arc_ptr = arc_list;
   1.577 +      arc_ptr = gen_basic_grid(csa, arc_ptr);
   1.578 +      select_source_sinks(csa);
   1.579 +      arc_ptr = gen_additional_arcs(csa, arc_ptr);
   1.580 +      gen_more_arcs(csa, arc_ptr);
   1.581 +      assign_costs(csa);
   1.582 +      assign_capacities(csa);
   1.583 +      assign_imbalance(csa);
   1.584 +      return;
   1.585 +}
   1.586 +
   1.587 +static void output(struct csa *csa)
   1.588 +{     /* Output the network in DIMACS format. */
   1.589 +      struct arcs *arc_ptr;
   1.590 +      struct imbalance *imb_ptr;
   1.591 +      int i;
   1.592 +      if (G != NULL) goto skip;
   1.593 +      /* Output "c", "p" records. */
   1.594 +      xprintf("c generated by GRIDGEN\n");
   1.595 +      xprintf("c seed %d\n", seed_original);
   1.596 +      xprintf("c nodes %d\n", n_node);
   1.597 +      xprintf("c grid size %d X %d\n", n1, n2);
   1.598 +      xprintf("c sources %d sinks %d\n", n_source, n_sink);
   1.599 +      xprintf("c avg. degree %d\n", avg_degree);
   1.600 +      xprintf("c supply %d\n", t_supply);
   1.601 +      switch (arc_costs.distribution)
   1.602 +      {  case UNIFORM:
   1.603 +            xprintf("c arc costs: UNIFORM distr. min %d max %d\n",
   1.604 +               (int)arc_costs.parameter[0],
   1.605 +               (int)arc_costs.parameter[1]);
   1.606 +            break;
   1.607 +         case EXPONENTIAL:
   1.608 +            xprintf("c arc costs: EXPONENTIAL distr. lambda %d\n",
   1.609 +               (int)arc_costs.parameter[0]);
   1.610 +            break;
   1.611 +         default:
   1.612 +            xassert(csa != csa);
   1.613 +      }
   1.614 +      switch (capacities.distribution)
   1.615 +      {  case UNIFORM:
   1.616 +            xprintf("c arc caps :  UNIFORM distr. min %d max %d\n",
   1.617 +               (int)capacities.parameter[0],
   1.618 +               (int)capacities.parameter[1]);
   1.619 +            break;
   1.620 +         case EXPONENTIAL:
   1.621 +            xprintf("c arc caps :  EXPONENTIAL distr. %d lambda %d\n",
   1.622 +               (int)capacities.parameter[0]);
   1.623 +            break;
   1.624 +         default:
   1.625 +            xassert(csa != csa);
   1.626 +      }
   1.627 +skip: if (G == NULL)
   1.628 +         xprintf("p min %d %d\n", n_node, n_arc);
   1.629 +      else
   1.630 +      {  glp_add_vertices(G, n_node);
   1.631 +         if (v_rhs >= 0)
   1.632 +         {  double zero = 0.0;
   1.633 +            for (i = 1; i <= n_node; i++)
   1.634 +            {  glp_vertex *v = G->v[i];
   1.635 +               memcpy((char *)v->data + v_rhs, &zero, sizeof(double));
   1.636 +            }
   1.637 +         }
   1.638 +      }
   1.639 +      /* Output "n node supply". */
   1.640 +      for (i = 0, imb_ptr = source_list; i < n_source; i++, imb_ptr++)
   1.641 +      {  if (G == NULL)
   1.642 +            xprintf("n %d %d\n", imb_ptr->node, imb_ptr->supply);
   1.643 +         else
   1.644 +         {  if (v_rhs >= 0)
   1.645 +            {  double temp = (double)imb_ptr->supply;
   1.646 +               glp_vertex *v = G->v[imb_ptr->node];
   1.647 +               memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
   1.648 +            }
   1.649 +         }
   1.650 +      }
   1.651 +      for (i = 0, imb_ptr = sink_list; i < n_sink; i++, imb_ptr++)
   1.652 +      {  if (G == NULL)
   1.653 +            xprintf("n %d %d\n", imb_ptr->node, imb_ptr->supply);
   1.654 +         else
   1.655 +         {  if (v_rhs >= 0)
   1.656 +            {  double temp = (double)imb_ptr->supply;
   1.657 +               glp_vertex *v = G->v[imb_ptr->node];
   1.658 +               memcpy((char *)v->data + v_rhs, &temp, sizeof(double));
   1.659 +            }
   1.660 +         }
   1.661 +      }
   1.662 +      /* Output "a from to lowcap=0 hicap cost". */
   1.663 +      for (i = 0, arc_ptr = arc_list; i < n_arc; i++, arc_ptr++)
   1.664 +      {  if (G == NULL)
   1.665 +            xprintf("a %d %d 0 %d %d\n", arc_ptr->from, arc_ptr->to,
   1.666 +               arc_ptr->u, arc_ptr->cost);
   1.667 +         else
   1.668 +         {  glp_arc *a = glp_add_arc(G, arc_ptr->from, arc_ptr->to);
   1.669 +            if (a_cap >= 0)
   1.670 +            {  double temp = (double)arc_ptr->u;
   1.671 +               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
   1.672 +            }
   1.673 +            if (a_cost >= 0)
   1.674 +            {  double temp = (double)arc_ptr->cost;
   1.675 +               memcpy((char *)a->data + a_cost, &temp, sizeof(double));
   1.676 +            }
   1.677 +         }
   1.678 +      }
   1.679 +      return;
   1.680 +}
   1.681 +
   1.682 +static double randy(struct csa *csa)
   1.683 +{     /* Returns a random number between 0.0 and 1.0.
   1.684 +         See Ward Cheney & David Kincaid, "Numerical Mathematics and
   1.685 +         Computing," 2Ed, pp. 335. */
   1.686 +      seed = 16807 * seed % 2147483647;
   1.687 +      if (seed < 0) seed = - seed;
   1.688 +      return seed * 4.6566128752459e-10;
   1.689 +}
   1.690 +
   1.691 +static void select_source_sinks(struct csa *csa)
   1.692 +{     /* Randomly select the source nodes and sink nodes. */
   1.693 +      int i, *int_ptr;
   1.694 +      int *temp_list;   /* a temporary list of nodes */
   1.695 +      struct imbalance *ptr;
   1.696 +      double ab[2];     /* parameter for random number generator */
   1.697 +      ab[0] = 0.9;
   1.698 +      ab[1] = n_node - 0.99;  /* upper limit is n_node-1 because the
   1.699 +                                 supernode cannot be selected */
   1.700 +      temp_list = xcalloc(n_node, sizeof(int));
   1.701 +      for (i = 0, int_ptr = temp_list; i < n_node; i++, int_ptr++)
   1.702 +         *int_ptr = 0;
   1.703 +      /* Select the source nodes. */
   1.704 +      for (i = 0, ptr = source_list; i < n_source; i++, ptr++)
   1.705 +      {  ptr->node = uniform(csa, ab);
   1.706 +         if (temp_list[ptr->node] == 1) /* check for duplicates */
   1.707 +         {  ptr--;
   1.708 +            i--;
   1.709 +         }
   1.710 +         else
   1.711 +            temp_list[ptr->node] = 1;
   1.712 +      }
   1.713 +      /* Select the sink nodes. */
   1.714 +      for (i = 0, ptr = sink_list; i < n_sink; i++, ptr++)
   1.715 +      {  ptr->node = uniform(csa, ab);
   1.716 +         if (temp_list[ptr->node] == 1)
   1.717 +         {  ptr--;
   1.718 +            i--;
   1.719 +         }
   1.720 +         else
   1.721 +            temp_list[ptr->node] = 1;
   1.722 +      }
   1.723 +      xfree(temp_list);
   1.724 +      return;
   1.725 +}
   1.726 +
   1.727 +int uniform(struct csa *csa, double a[2])
   1.728 +{     /* Generates an integer uniformly selected from [a[0],a[1]]. */
   1.729 +      return (int)((a[1] - a[0]) * randy(csa) + a[0] + 0.5);
   1.730 +}
   1.731 +
   1.732 +/**********************************************************************/
   1.733 +
   1.734 +#if 0
   1.735 +int main(void)
   1.736 +{     int parm[1+14];
   1.737 +      double temp;
   1.738 +      scanf("%d", &parm[1]);
   1.739 +      scanf("%d", &parm[2]);
   1.740 +      scanf("%d", &parm[3]);
   1.741 +      scanf("%d", &parm[4]);
   1.742 +      scanf("%d", &parm[5]);
   1.743 +      scanf("%d", &parm[6]);
   1.744 +      scanf("%d", &parm[7]);
   1.745 +      scanf("%d", &parm[8]);
   1.746 +      scanf("%d", &parm[9]);
   1.747 +      if (parm[9] == 1)
   1.748 +      {  scanf("%d", &parm[10]);
   1.749 +         scanf("%d", &parm[11]);
   1.750 +      }
   1.751 +      else
   1.752 +      {  scanf("%le", &temp);
   1.753 +         parm[10] = (int)(100.0 * temp + .5);
   1.754 +         parm[11] = 0;
   1.755 +      }
   1.756 +      scanf("%d", &parm[12]);
   1.757 +      if (parm[12] == 1)
   1.758 +      {  scanf("%d", &parm[13]);
   1.759 +         scanf("%d", &parm[14]);
   1.760 +      }
   1.761 +      else
   1.762 +      {  scanf("%le", &temp);
   1.763 +         parm[13] = (int)(100.0 * temp + .5);
   1.764 +         parm[14] = 0;
   1.765 +      }
   1.766 +      glp_gridgen(NULL, 0, 0, 0, parm);
   1.767 +      return 0;
   1.768 +}
   1.769 +#endif
   1.770 +
   1.771 +/* eof */