src/glpnet05.c
changeset 1 c445c931472f
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/glpnet05.c	Mon Dec 06 13:09:21 2010 +0100
     1.3 @@ -0,0 +1,367 @@
     1.4 +/* glpnet05.c (Goldfarb's maximum flow 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 RMFGEN, a maxflow
    1.10 +*  problem generator developed by D.Goldfarb and M.Grigoriadis, and
    1.11 +*  originally implemented by Tamas Badics <badics@rutcor.rutgers.edu>.
    1.12 +*  The original code is publically available on the DIMACS ftp site at:
    1.13 +*  <ftp://dimacs.rutgers.edu/pub/netflow/generators/network/genrmf>.
    1.14 +*
    1.15 +*  All changes concern only the program interface, so this modified
    1.16 +*  version produces exactly the same instances as the original version.
    1.17 +*
    1.18 +*  Changes were made by Andrew Makhorin <mao@gnu.org>.
    1.19 +*
    1.20 +*  GLPK is free software: you can redistribute it and/or modify it
    1.21 +*  under the terms of the GNU General Public License as published by
    1.22 +*  the Free Software Foundation, either version 3 of the License, or
    1.23 +*  (at your option) any later version.
    1.24 +*
    1.25 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.26 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.27 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.28 +*  License for more details.
    1.29 +*
    1.30 +*  You should have received a copy of the GNU General Public License
    1.31 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.32 +***********************************************************************/
    1.33 +
    1.34 +#include "glpapi.h"
    1.35 +#include "glprng.h"
    1.36 +
    1.37 +/***********************************************************************
    1.38 +*  NAME
    1.39 +*
    1.40 +*  glp_rmfgen - Goldfarb's maximum flow problem generator
    1.41 +*
    1.42 +*  SYNOPSIS
    1.43 +*
    1.44 +*  int glp_rmfgen(glp_graph *G, int *s, int *t, int a_cap,
    1.45 +*     const int parm[1+5]);
    1.46 +*
    1.47 +*  DESCRIPTION
    1.48 +*
    1.49 +*  The routine glp_rmfgen is a maximum flow problem generator developed
    1.50 +*  by D.Goldfarb and M.Grigoriadis.
    1.51 +*
    1.52 +*  The parameter G specifies the graph object, to which the generated
    1.53 +*  problem data have to be stored. Note that on entry the graph object
    1.54 +*  is erased with the routine glp_erase_graph.
    1.55 +*
    1.56 +*  The pointer s specifies a location, to which the routine stores the
    1.57 +*  source node number. If s is NULL, the node number is not stored.
    1.58 +*
    1.59 +*  The pointer t specifies a location, to which the routine stores the
    1.60 +*  sink node number. If t is NULL, the node number is not stored.
    1.61 +*
    1.62 +*  The parameter a_cap specifies an offset of the field of type double
    1.63 +*  in the arc data block, to which the routine stores the arc capacity.
    1.64 +*  If a_cap < 0, the capacity 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]  (seed)   random number seed (a positive integer)
    1.70 +*  parm[2]  (a)      frame size
    1.71 +*  parm[3]  (b)      depth
    1.72 +*  parm[4]  (c1)     minimal arc capacity
    1.73 +*  parm[5]  (c2)     maximal arc capacity
    1.74 +*
    1.75 +*  RETURNS
    1.76 +*
    1.77 +*  If the instance was successfully generated, the routine glp_netgen
    1.78 +*  returns zero; otherwise, if specified parameters are inconsistent,
    1.79 +*  the routine returns a non-zero error code.
    1.80 +*
    1.81 +*  COMMENTS
    1.82 +*
    1.83 +*  The generated network is as follows. It has b pieces of frames of
    1.84 +*  size a * a. (So alltogether the number of vertices is a * a * b)
    1.85 +*
    1.86 +*  In each frame all the vertices are connected with their neighbours
    1.87 +*  (forth and back). In addition the vertices of a frame are connected
    1.88 +*  one to one with the vertices of next frame using a random permutation
    1.89 +*  of those vertices.
    1.90 +*
    1.91 +*  The source is the lower left vertex of the first frame, the sink is
    1.92 +*  the upper right vertex of the b'th frame.
    1.93 +*
    1.94 +*                    t
    1.95 +*           +-------+
    1.96 +*           |      .|
    1.97 +*           |     . |
    1.98 +*        /  |    /  |
    1.99 +*       +-------+/ -+ b
   1.100 +*       |    |  |/.
   1.101 +*     a |   -v- |/
   1.102 +*       |    |  |/
   1.103 +*       +-------+ 1
   1.104 +*      s    a
   1.105 +*
   1.106 +*  The capacities are randomly chosen integers from the range of [c1,c2]
   1.107 +*  in the case of interconnecting edges, and c2 * a * a for the in-frame
   1.108 +*  edges.
   1.109 +*
   1.110 +*  REFERENCES
   1.111 +*
   1.112 +*  D.Goldfarb and M.D.Grigoriadis, "A computational comparison of the
   1.113 +*  Dinic and network simplex methods for maximum flow." Annals of Op.
   1.114 +*  Res. 13 (1988), pp. 83-123.
   1.115 +*
   1.116 +*  U.Derigs and W.Meier, "Implementing Goldberg's max-flow algorithm:
   1.117 +*  A computational investigation." Zeitschrift fuer Operations Research
   1.118 +*  33 (1989), pp. 383-403. */
   1.119 +
   1.120 +typedef struct VERTEX
   1.121 +{     struct EDGE **edgelist;
   1.122 +      /* Pointer to the list of pointers to the adjacent edges.
   1.123 +         (No matter that to or from edges) */
   1.124 +      struct EDGE **current;
   1.125 +      /* Pointer to the current edge */
   1.126 +      int degree;
   1.127 +      /* Number of adjacent edges (both direction) */
   1.128 +      int index;
   1.129 +} vertex;
   1.130 +
   1.131 +typedef struct EDGE
   1.132 +{     int from;
   1.133 +      int to;
   1.134 +      int cap;
   1.135 +      /* Capacity */
   1.136 +} edge;
   1.137 +
   1.138 +typedef struct NETWORK
   1.139 +{     struct NETWORK *next, *prev;
   1.140 +      int vertnum;
   1.141 +      int edgenum;
   1.142 +      vertex *verts;
   1.143 +      /* Vertex array[1..vertnum] */
   1.144 +      edge *edges;
   1.145 +      /* Edge array[1..edgenum] */
   1.146 +      int source;
   1.147 +      /* Pointer to the source */
   1.148 +      int sink;
   1.149 +      /* Pointer to the sink */
   1.150 +} network;
   1.151 +
   1.152 +struct csa
   1.153 +{     /* common storage area */
   1.154 +      glp_graph *G;
   1.155 +      int *s, *t, a_cap;
   1.156 +      RNG *rand;
   1.157 +      network *N;
   1.158 +      int *Parr;
   1.159 +      int A, AA, C2AA, Ec;
   1.160 +};
   1.161 +
   1.162 +#define G      (csa->G)
   1.163 +#define s      (csa->s)
   1.164 +#define t      (csa->t)
   1.165 +#define a_cap  (csa->a_cap)
   1.166 +#define N      (csa->N)
   1.167 +#define Parr   (csa->Parr)
   1.168 +#define A      (csa->A)
   1.169 +#define AA     (csa->AA)
   1.170 +#define C2AA   (csa->C2AA)
   1.171 +#define Ec     (csa->Ec)
   1.172 +
   1.173 +#undef random
   1.174 +#define random(A) (int)(rng_unif_01(csa->rand) * (double)(A))
   1.175 +#define RANDOM(A, B) (int)(random((B) - (A) + 1) + (A))
   1.176 +#define sgn(A) (((A) > 0) ? 1 : ((A) == 0) ? 0 : -1)
   1.177 +
   1.178 +static void make_edge(struct csa *csa, int from, int to, int c1, int c2)
   1.179 +{     Ec++;
   1.180 +      N->edges[Ec].from = from;
   1.181 +      N->edges[Ec].to = to;
   1.182 +      N->edges[Ec].cap = RANDOM(c1, c2);
   1.183 +      return;
   1.184 +}
   1.185 +
   1.186 +static void permute(struct csa *csa)
   1.187 +{     int i, j, tmp;
   1.188 +      for (i = 1; i < AA; i++)
   1.189 +      {  j = RANDOM(i, AA);
   1.190 +         tmp = Parr[i];
   1.191 +         Parr[i] = Parr[j];
   1.192 +         Parr[j] = tmp;
   1.193 +      }
   1.194 +      return;
   1.195 +}
   1.196 +
   1.197 +static void connect(struct csa *csa, int offset, int cv, int x1, int y1)
   1.198 +{     int cv1;
   1.199 +      cv1 = offset + (x1 - 1) * A + y1;
   1.200 +      Ec++;
   1.201 +      N->edges[Ec].from = cv;
   1.202 +      N->edges[Ec].to = cv1;
   1.203 +      N->edges[Ec].cap = C2AA;
   1.204 +      return;
   1.205 +}
   1.206 +
   1.207 +static network *gen_rmf(struct csa *csa, int a, int b, int c1, int c2)
   1.208 +{     /* generates a network with a*a*b nodes and 6a*a*b-4ab-2a*a edges
   1.209 +         random_frame network:
   1.210 +         Derigs & Meier, Methods & Models of OR (1989), 33:383-403 */
   1.211 +      int x, y, z, offset, cv;
   1.212 +      A = a;
   1.213 +      AA = a * a;
   1.214 +      C2AA = c2 * AA;
   1.215 +      Ec = 0;
   1.216 +      N = (network *)xmalloc(sizeof(network));
   1.217 +      N->vertnum = AA * b;
   1.218 +      N->edgenum = 5 * AA * b - 4 * A * b - AA;
   1.219 +      N->edges = (edge *)xcalloc(N->edgenum + 1, sizeof(edge));
   1.220 +      N->source = 1;
   1.221 +      N->sink = N->vertnum;
   1.222 +      Parr = (int *)xcalloc(AA + 1, sizeof(int));
   1.223 +      for (x = 1; x <= AA; x++)
   1.224 +         Parr[x] = x;
   1.225 +      for (z = 1; z <= b; z++)
   1.226 +      {  offset = AA * (z - 1);
   1.227 +         if (z != b)
   1.228 +            permute(csa);
   1.229 +         for (x = 1; x <= A; x++)
   1.230 +         {  for (y = 1; y <= A; y++)
   1.231 +            {  cv = offset + (x - 1) * A + y;
   1.232 +               if (z != b)
   1.233 +                  make_edge(csa, cv, offset + AA + Parr[cv - offset],
   1.234 +                     c1, c2); /* the intermediate edges */
   1.235 +               if (y < A)
   1.236 +                  connect(csa, offset, cv, x, y + 1);
   1.237 +               if (y > 1)
   1.238 +                  connect(csa, offset, cv, x, y - 1);
   1.239 +               if (x < A)
   1.240 +                  connect(csa, offset, cv, x + 1, y);
   1.241 +               if (x > 1)
   1.242 +                  connect(csa, offset, cv, x - 1, y);
   1.243 +            }
   1.244 +         }
   1.245 +      }
   1.246 +      xfree(Parr);
   1.247 +      return N;
   1.248 +}
   1.249 +
   1.250 +static void print_max_format(struct csa *csa, network *n, char *comm[],
   1.251 +      int dim)
   1.252 +{     /* prints a network heading with dim lines of comments (no \n
   1.253 +         needs at the ends) */
   1.254 +      int i, vnum, e_num;
   1.255 +      edge *e;
   1.256 +      vnum = n->vertnum;
   1.257 +      e_num = n->edgenum;
   1.258 +      if (G == NULL)
   1.259 +      {  for (i = 0; i < dim; i++)
   1.260 +            xprintf("c %s\n", comm[i]);
   1.261 +         xprintf("p max %7d %10d\n", vnum, e_num);
   1.262 +         xprintf("n %7d s\n", n->source);
   1.263 +         xprintf("n %7d t\n", n->sink);
   1.264 +      }
   1.265 +      else
   1.266 +      {  glp_add_vertices(G, vnum);
   1.267 +         if (s != NULL) *s = n->source;
   1.268 +         if (t != NULL) *t = n->sink;
   1.269 +      }
   1.270 +      for (i = 1; i <= e_num; i++)
   1.271 +      {  e = &n->edges[i];
   1.272 +         if (G == NULL)
   1.273 +            xprintf("a %7d %7d %10d\n", e->from, e->to, (int)e->cap);
   1.274 +         else
   1.275 +         {  glp_arc *a = glp_add_arc(G, e->from, e->to);
   1.276 +            if (a_cap >= 0)
   1.277 +            {  double temp = (double)e->cap;
   1.278 +               memcpy((char *)a->data + a_cap, &temp, sizeof(double));
   1.279 +            }
   1.280 +         }
   1.281 +      }
   1.282 +      return;
   1.283 +}
   1.284 +
   1.285 +static void gen_free_net(network *n)
   1.286 +{     xfree(n->edges);
   1.287 +      xfree(n);
   1.288 +      return;
   1.289 +}
   1.290 +
   1.291 +int glp_rmfgen(glp_graph *G_, int *_s, int *_t, int _a_cap,
   1.292 +      const int parm[1+5])
   1.293 +{     struct csa _csa, *csa = &_csa;
   1.294 +      network *n;
   1.295 +      char comm[10][80], *com1[10];
   1.296 +      int seed, a, b, c1, c2, ret;
   1.297 +      G = G_;
   1.298 +      s = _s;
   1.299 +      t = _t;
   1.300 +      a_cap = _a_cap;
   1.301 +      if (G != NULL)
   1.302 +      {  if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
   1.303 +           xerror("glp_rmfgen: a_cap = %d; invalid offset\n", a_cap);
   1.304 +      }
   1.305 +      seed = parm[1];
   1.306 +      a = parm[2];
   1.307 +      b = parm[3];
   1.308 +      c1 = parm[4];
   1.309 +      c2 = parm[5];
   1.310 +      if (!(seed > 0 && 1 <= a && a <= 1000 && 1 <= b && b <= 1000 &&
   1.311 +            0 <= c1 && c1 <= c2 && c2 <= 1000))
   1.312 +      {  ret = 1;
   1.313 +         goto done;
   1.314 +      }
   1.315 +      if (G != NULL)
   1.316 +      {  glp_erase_graph(G, G->v_size, G->a_size);
   1.317 +         glp_set_graph_name(G, "RMFGEN");
   1.318 +      }
   1.319 +      csa->rand = rng_create_rand();
   1.320 +      rng_init_rand(csa->rand, seed);
   1.321 +      n = gen_rmf(csa, a, b, c1, c2);
   1.322 +      sprintf(comm[0], "This file was generated by genrmf.");
   1.323 +      sprintf(comm[1], "The parameters are: a: %d b: %d c1: %d c2: %d",
   1.324 +         a, b, c1, c2);
   1.325 +      com1[0] = comm[0];
   1.326 +      com1[1] = comm[1];
   1.327 +      print_max_format(csa, n, com1, 2);
   1.328 +      gen_free_net(n);
   1.329 +      rng_delete_rand(csa->rand);
   1.330 +      ret = 0;
   1.331 +done: return ret;
   1.332 +}
   1.333 +
   1.334 +/**********************************************************************/
   1.335 +
   1.336 +#if 0
   1.337 +int main(int argc, char *argv[])
   1.338 +{     int seed, a, b, c1, c2, i, parm[1+5];
   1.339 +      seed = 123;
   1.340 +      a = b = c1 = c2 = -1;
   1.341 +      for (i = 1; i < argc; i++)
   1.342 +      {  if (strcmp(argv[i], "-seed") == 0)
   1.343 +            seed = atoi(argv[++i]);
   1.344 +         else if (strcmp(argv[i], "-a") == 0)
   1.345 +            a = atoi(argv[++i]);
   1.346 +         else if (strcmp(argv[i], "-b") == 0)
   1.347 +            b = atoi(argv[++i]);
   1.348 +         else if (strcmp(argv[i], "-c1") == 0)
   1.349 +            c1 = atoi(argv[++i]);
   1.350 +         else if (strcmp(argv[i], "-c2") == 0)
   1.351 +            c2 = atoi(argv[++i]);
   1.352 +      }
   1.353 +      if (a < 0 || b < 0 || c1 < 0 || c2 < 0)
   1.354 +      {  xprintf("Usage:\n");
   1.355 +         xprintf("genrmf [-seed seed] -a frame_size -b depth\n");
   1.356 +         xprintf("        -c1 cap_range1 -c2 cap_range2\n");
   1.357 +      }
   1.358 +      else
   1.359 +      {  parm[1] = seed;
   1.360 +         parm[2] = a;
   1.361 +         parm[3] = b;
   1.362 +         parm[4] = c1;
   1.363 +         parm[5] = c2;
   1.364 +         glp_rmfgen(NULL, NULL, NULL, 0, parm);
   1.365 +      }
   1.366 +      return 0;
   1.367 +}
   1.368 +#endif
   1.369 +
   1.370 +/* eof */