1 /****************************************************************************/
2 /* */
3 /* This file is part of CONCORDE */
4 /* */
5 /* (c) Copyright 1995--1999 by David Applegate, Robert Bixby, */
6 /* Vasek Chvatal, and William Cook */
7 /* */
8 /* Permission is granted for academic research use. For other uses, */
9 /* contact the authors for licensing options. */
10 /* */
11 /* Use at your own risk. We make no guarantees about the */
12 /* correctness or usefulness of this code. */
13 /* */
14 /****************************************************************************/
15
16 /****************************************************************************/
17 /* */
18 /* LINEAR SUBTOUR SEPARATION ROUTINES */
19 /* preliminary */
20 /* */
21 /* TSP CODE */
22 /* Written by: Applegate, Bixby, Chvatal, and Cook */
23 /* Date: July 25, 1995 */
24 /* Thanks: Phil Gibbons, for helpful ideas */
25 /* */
26 /* EXPORTED FUNCTIONS: */
27 /* */
28 /* int CCcut_linsub (int ncount, int ecount, int *endmark, int *elist, */
29 /* double *x, double maxval, void *u_data, */
30 /* int (*cut_callback) (double cut_val, int cut_start, int cut_end, */
31 /* void *u_data)) */
32 /* -ncount is the number of nodes */
33 /* -ecount is the number of edges */
34 /* -endmark indicates which segments are of interest, by indicating */
35 /* whether a node can be a right or left end of a segment */
36 /* -elist contains the LP edges in node node format */
37 /* -x is an LP solution */
38 /* -maxval is the maximum cut value desired */
39 /* -u_data is user data to be passed to cut_callback */
40 /* -cut_callback is a function to be called for segments which */
41 /* define a cut of value < cutlim. The cut is cut_start, */
42 /* cut_start+1, ..., cut_end, and has value cut_val. cut_callback */
43 /* will be called for the minimum segment cut starting at each */
44 /* endpoint marked as a right end, provided that cut has value < */
45 /* cutlim. */
46 /* */
47 /* int CCcut_linsub_allcuts (int ncount, int ecount, int *perm, */
48 /* int *endmark, int *elist, double *x, double maxval, */
49 /* void *u_data, int (*cut_callback) (double cut_val, */
50 /* int cut_start, int cut_end, void *u_data)) */
51 /* -ncount is the number of nodes */
52 /* -ecount is the number of edges */
53 /* -perm is a permutation of the nodes (if perm == (int *) NULL, */
54 /* the identity permutation will be used) */
55 /* -elist contains the LP edges in node node format */
56 /* -endmark indicates which segments are of interest, by indicating */
57 /* whether a node can be a right or left end of a segment */
58 /* -x is an LP solution */
59 /* -maxval is the maximum cut value desired */
60 /* -u_data is data to be passed to the callback */
61 /* -cut_callback is a function to be called for every segment which */
62 /* defines a cut of value <= cutlim. The cut is perm[cut_start], */
63 /* perm[cut_start+1], ..., perm[cut_end], and has value cut_val. */
64 /* if cut_callback returns a nonzero value, CCcut_linsub_allcuts */
65 /* will terminate. */
66 /* */
67 /* NOTES: */
68 /* CCcut_linsub runs in time O(m log n). */
69 /* CCcut_linsub_allcuts runs in time O(m log n + |C| log n) where */
70 /* |C| is the number of cuts found. */
71 /* */
72 /****************************************************************************/
73
74 #include "machdefs.h"
75 #include "util.h"
76 #include "cut.h"
77
78 #define LINSUB_INF 1e20
79
80 typedef struct psh {
81 int base;
82 int size;
83 double *sum;
84 double *minpre;
85 } psh;
86
87
88 static void
89 psh_free (psh *p),
90 psh_add (psh *p, int i, double v);
91
92 static int
93 psh_init (psh *p, int k, int *endmark),
94 psh_minloc (psh *p),
95 psh_enum (psh *p, int cut_start, double maxval, void *u_data,
96 int (*cut_callback)(double cut_val, int cut_start, int cut_end,
97 void *u_data)),
98 psh_enum_work (psh *p, int n, int mul, double pre_sum, int cut_start,
99 double maxslack, void *u_data, int (*cut_callback)(double cut_val,
100 int cut_start, int cut_end, void *u_data));
101
102 static double
103 psh_minval (psh *p);
104
105
CCcut_linsub(int ncount,int ecount,int * endmark,int * elist,double * x,double maxval,void * u_data,int (* cut_callback)(double cut_val,int cut_start,int cut_end,void * u_data))106 int CCcut_linsub (int ncount, int ecount, int *endmark, int *elist, double *x,
107 double maxval, void *u_data, int (*cut_callback) (double cut_val,
108 int cut_start, int cut_end, void *u_data))
109 {
110 psh p;
111 int i, j;
112 double v;
113 int rval = 0;
114 int *perm = (int *) NULL;
115 int *eperm = (int *) NULL;
116 int *ends = (int *) NULL;
117 double *xends = (double *) NULL;
118
119 if (psh_init (&p, ncount, endmark)) {
120 return -1;
121 }
122
123 /* arrange elist into the array ends (with xends being the x values */
124 /* so that ends[2*i] is less than ends[2*i+1], and ends is sorted in */
125 /* increasing order of ends[2*i] */
126
127 perm = CC_SAFE_MALLOC (ecount, int);
128 eperm = CC_SAFE_MALLOC (ecount, int);
129 if (!perm || !eperm) {
130 fprintf (stderr, "out of memory in CCcut_linsub\n");
131 rval = 1; goto CLEANUP;
132 }
133 for (i = 0; i < ecount; i++) {
134 eperm[i] = (elist[2*i] < elist[2*i+1] ? elist[2*i] : elist[2*i+1]);
135 perm[i] = i;
136 }
137 CCutil_int_perm_quicksort (perm, eperm, ecount);
138
139 ends = CC_SAFE_MALLOC (2*ecount, int);
140 xends = CC_SAFE_MALLOC (ecount, double);
141 if (!ends || !xends) {
142 fprintf (stderr, "out of memory in CCcut_linsub\n");
143 rval = 1; goto CLEANUP;
144 }
145 for (i = 0; i < ecount; i++) {
146 j = perm[i];
147 if (elist[2*j] < elist[2*j+1]) {
148 ends[2*i] = elist[2*j];
149 ends[2*i+1] = elist[2*j+1];
150 } else {
151 ends[2*i] = elist[2*j+1];
152 ends[2*i+1] = elist[2*j];
153 }
154 xends[i] = x[j];
155 }
156 CC_FREE (perm, int);
157 CC_FREE (eperm, int);
158
159 for (i = ncount - 1, j = ecount - 1; i > 0; i--) {
160 while (j >= 0 && ends[2*j] == i) {
161 psh_add (&p, ends[2*j+1], -xends[j]);
162
163 j--;
164 }
165 if (endmark[i] & CC_LINSUB_LEFT_END) {
166 v = 2.0 + 2.0 * psh_minval (&p);
167 if (v < maxval) {
168 rval = (*cut_callback) (v, i, psh_minloc (&p), u_data);
169 if (rval) {
170 fprintf (stderr, "cut_callback failed\n"); goto CLEANUP;
171 }
172 }
173 }
174 psh_add (&p, i, 1.0);
175 }
176
177 CLEANUP:
178
179 psh_free (&p);
180 CC_IFFREE (ends, int);
181 CC_IFFREE (xends, double);
182 CC_IFFREE (perm, int);
183 CC_IFFREE (eperm, int);
184
185 return rval;
186 }
187
CCcut_linsub_allcuts(int ncount,int ecount,int * perm,int * endmark,int * elist,double * x,double maxval,void * u_data,int (* cut_callback)(double cut_val,int cut_start,int cut_end,void * u_data))188 int CCcut_linsub_allcuts (int ncount, int ecount, int *perm, int *endmark,
189 int *elist, double *x, double maxval, void *u_data,
190 int (*cut_callback) (double cut_val, int cut_start, int cut_end,
191 void *u_data))
192 {
193 psh p;
194 int i, j;
195 int rval = 0;
196 int *perm_inv = (int *) NULL;
197 int *esort = (int *) NULL;
198 int *eperm = (int *) NULL;
199 int *ends = (int *) NULL;
200 int *pendmark = (int *) NULL;
201 double *xends = (double *) NULL;
202
203 perm_inv = CC_SAFE_MALLOC (ncount, int);
204 pendmark = CC_SAFE_MALLOC (ncount, int);
205 if (perm_inv == (int *) NULL ||
206 pendmark == (int *) NULL) {
207 fprintf (stderr, "Out of memory in CCcut_linsub_allcuts\n");
208 CC_IFFREE (perm_inv, int);
209 CC_IFFREE (pendmark, int);
210 return -1;
211 }
212
213 if (perm == (int *) NULL) {
214 for (i=0; i<ncount; i++) {
215 perm_inv[i] = i;
216 pendmark[i] = endmark[i];
217 }
218 } else {
219 for (i=0; i<ncount; i++) {
220 perm_inv[perm[i]] = i;
221 pendmark[i] = endmark[perm[i]];
222 }
223 }
224
225 if (psh_init (&p, ncount, pendmark)) {
226 fprintf (stderr, "psh_init failed\n");
227 CC_IFFREE (perm_inv, int);
228 CC_IFFREE (pendmark, int);
229 return -1;
230 }
231
232 /* arrange elist into the array ends (with xends being the x values */
233 /* so that ends[2*i] is less than ends[2*i+1], and ends is sorted in */
234 /* increasing order of ends[2*i] */
235
236 eperm = CC_SAFE_MALLOC (ecount, int);
237 esort = CC_SAFE_MALLOC (ecount, int);
238 if (eperm == (int *) NULL || esort == (int *) NULL) {
239 fprintf (stderr, "out of memory in CCcut_linsub_allcuts\n");
240 rval = 1; goto CLEANUP;
241 }
242 for (i = 0; i < ecount; i++) {
243 esort[i] = (perm_inv[elist[2*i]] < perm_inv[elist[2*i+1]]
244 ? perm_inv[elist[2*i]]
245 : perm_inv[elist[2*i+1]]);
246 eperm[i] = i;
247 }
248 CCutil_int_perm_quicksort (eperm, esort, ecount);
249
250 ends = CC_SAFE_MALLOC (2*ecount, int);
251 xends = CC_SAFE_MALLOC (ecount, double);
252 if (ends == (int *) NULL || xends == (double *) NULL) {
253 fprintf (stderr, "out of memory in CCcut_linsub_allcuts\n");
254 rval = 1; goto CLEANUP;
255 }
256 for (i = 0; i < ecount; i++) {
257 j = eperm[i];
258 if (perm_inv[elist[2*j]] < perm_inv[elist[2*j+1]]) {
259 ends[2*i] = perm_inv[elist[2*j]];
260 ends[2*i+1] = perm_inv[elist[2*j+1]];
261 } else {
262 ends[2*i] = perm_inv[elist[2*j+1]];
263 ends[2*i+1] = perm_inv[elist[2*j]];
264 }
265 xends[i] = x[j];
266 }
267 CC_FREE (eperm, int);
268 CC_FREE (esort, int);
269
270 for (i = ncount - 1, j = ecount - 1; i > 0; i--) {
271 while (j >= 0 && ends[2*j] == i) {
272 psh_add (&p, ends[2*j+1], -xends[j]);
273 j--;
274 }
275 if (pendmark[i] & CC_LINSUB_LEFT_END) {
276 rval = psh_enum (&p, i, maxval, u_data, cut_callback);
277 if (rval) {
278 fprintf (stderr, "psh_enum failed\n");
279 goto CLEANUP;
280 }
281 }
282 psh_add (&p, i, 1.0);
283 }
284
285 rval = 0;
286
287 CLEANUP:
288
289 psh_free (&p);
290 CC_IFFREE (ends, int);
291 CC_IFFREE (xends, double);
292 CC_IFFREE (eperm, int);
293 CC_IFFREE (esort, int);
294 CC_IFFREE (pendmark, int);
295 CC_IFFREE (perm_inv, int);
296
297 return rval;
298 }
299
300 /****************************************************************************/
301 /* */
302 /* PREFIX SUM HEAP ROUTINES */
303 /* */
304 /* TSP CODE */
305 /* Written by: Applegate, Bixby, Chvatal, and Cook */
306 /* Date: July 25, 1995 */
307 /* Thanks: Phil Gibbons, for the algorithm */
308 /* */
309 /* These routines implement a heap of prefix sums. The routines are */
310 /* self-contained, but unlikely to be useful outside of the linear */
311 /* subtour separation code. */
312 /* */
313 /* PREFIX SUM HEAP FUNCTIONS: */
314 /* int psh_init (psh *p, int k, int *endmark) */
315 /* -p should point to a psh struct. */
316 /* -k is the number of elements in the psh. */
317 /* -endmark marks which nodes can have prefix sums of interest */
318 /* -this initializes the value of all k elements to 0.0. */
319 /* -return value 0 means success, 1 means failure. */
320 /* void psh_free (psh *p) */
321 /* -frees the spaces allocated by psh_init. */
322 /* void psh_add (psh *p, int i, double v) */
323 /* -adds v to the value of element i. */
324 /* double psh_minval (psh *p) */
325 /* -returns min_j sum_{i=0}^j value[i]. */
326 /* over j with endmark[j] & CC_LINSUB_RIGHT_END */
327 /* int psh_minloc (psh *p) */
328 /* -returns the smallest j which achieves psh_minval(p). */
329 /* over j with endmark[j] & CC_LINSUB_RIGHT_END */
330 /* int psh_enum (psh *p, int cut_start, double maxval, void *u_data, */
331 /* int (*cut_callback)(double cut_val, int cut_start, */
332 /* int cut_end, void *u_data)) */
333 /* -calls cut_callback for each interval starting at cut_start */
334 /* defining a cut <= maxval */
335 /* with endmark[cut_start] & CC_LINSUB_LEFT_END, and */
336 /* with endmark[cut_end] & CC_LINSUB_RIGHT_END. */
337 /* */
338 /* NOTES: */
339 /* A k-element heap will malloc 32k bytes of memory. If memory is */
340 /* tight, using integer values (instead of doubles) brings it down to */
341 /* 16k bytes. */
342 /* psh_init takes O(k) time. psh_add and psh_minloc take O(log k) */
343 /* time. psh_free and psh_minval take O(1) time. */
344 /* */
345 /* It is likely that using a ternary tree instead of binary would */
346 /* improve performance. Also, psh_add could take advantage of knowing */
347 /* which child has changed. */
348 /* psh_minloc could be changed to return all elements which achieve */
349 /* the min, in time O(log k) per element returned. */
350 /* */
351 /****************************************************************************/
352
psh_init(psh * p,int k,int * endmark)353 static int psh_init (psh *p, int k, int *endmark)
354 {
355 int i;
356 int space;
357
358 p->size = k;
359 p->base = 1;
360 while (p->base < k) p->base *= 2;
361
362 space = p->base*2;
363
364 p->sum = CC_SAFE_MALLOC (space, double);
365 if (!p->sum)
366 return 1;
367 p->minpre = CC_SAFE_MALLOC (space, double);
368 if (!p->minpre) {
369 CC_FREE (p->sum, double);
370 return 1;
371 }
372
373 for (i=0; i<k; i++) {
374 p->sum[p->base + i] = 0.0;
375 if (endmark[i] & CC_LINSUB_RIGHT_END) {
376 p->minpre[p->base + i] = 0.0;
377 } else {
378 p->minpre[p->base + i] = LINSUB_INF;
379 }
380 }
381 for (i=k; i<p->base; i++) {
382 p->sum[p->base + i] = 0.0;
383 p->minpre[p->base + i] = LINSUB_INF;
384 }
385 for (i=p->base - 1; i>=1; i--) {
386 p->sum[i] = p->sum[2*i] + p->sum[2*i+1];
387 if (p->minpre[2*i] < p->sum[2*i] + p->minpre[2*i+1]) {
388 p->minpre[i] = p->minpre[2*i];
389 } else {
390 p->minpre[i] = p->sum[2*i] + p->minpre[2*i+1];
391 }
392 }
393
394 return 0;
395 }
396
psh_free(psh * p)397 static void psh_free (psh *p)
398 {
399 CC_FREE (p->minpre, double);
400 CC_FREE (p->sum, double);
401 p->size = 0;
402 p->base = 0;
403 }
404
psh_add(psh * p,int i,double v)405 static void psh_add (psh *p, int i, double v)
406 {
407 double *s = p->sum;
408 double *m = p->minpre;
409
410 i += p->base;
411 s[i] += v;
412 m[i] += v;
413 i /= 2;
414 while (i >= 1) {
415 s[i] += v;
416 if (m[2*i] < s[2*i] + m[2*i+1]) {
417 m[i] = m[2*i];
418 } else {
419 m[i] = s[2*i] + m[2*i+1];
420 }
421 i /= 2;
422 }
423 }
424
psh_minloc(psh * p)425 static int psh_minloc (psh *p)
426 {
427 int i = 1;
428 double *m = p->minpre;
429
430 while (i < p->base) {
431 if (m[i] == m[2*i]) {
432 i = 2*i;
433 } else {
434 i = 2*i+1;
435 }
436 }
437 if (i - p->base < p->size) {
438 return i - p->base;
439 } else {
440 /* we're lost */
441 return p->size - 1;
442 }
443 }
444
psh_minval(psh * p)445 static double psh_minval (psh *p)
446 {
447 return p->minpre[1];
448 }
449
psh_enum(psh * p,int cut_start,double maxval,void * u_data,int (* cut_callback)(double cut_val,int cut_start,int cut_end,void * u_data))450 static int psh_enum (psh *p, int cut_start, double maxval, void *u_data,
451 int (*cut_callback)(double cut_val, int cut_start, int cut_end,
452 void *u_data))
453 {
454 int rval;
455
456 rval = psh_enum_work (p, 1, p->base, 0.0, cut_start,
457 (maxval - 2.0) / 2.0, u_data, cut_callback);
458 if (rval) {
459 fprintf (stderr, "psh_enum_work failed\n");
460 return rval;
461 }
462 return 0;
463 }
464
psh_enum_work(psh * p,int n,int mul,double pre_sum,int cut_start,double maxslack,void * u_data,int (* cut_callback)(double cut_val,int cut_start,int cut_end,void * u_data))465 static int psh_enum_work (psh *p, int n, int mul, double pre_sum,
466 int cut_start, double maxslack, void *u_data,
467 int (*cut_callback)(double cut_val, int cut_start, int cut_end,
468 void *u_data))
469 {
470 int rval;
471
472 if (n >= p->base) {
473 rval = (*cut_callback) (2.0 + 2.0*(pre_sum + p->sum[n]),
474 cut_start, n - p->base, u_data);
475 if (rval) {
476 fprintf (stderr, "cut_callback failed\n");
477 return rval;
478 }
479 return 0;
480 }
481
482 mul /= 2;
483 n *= 2;
484 if (n*mul < p->base + p->size &&
485 n*mul + mul > p->base + cut_start &&
486 pre_sum + p->minpre[n] <= maxslack) {
487 rval = psh_enum_work (p, n, mul, pre_sum, cut_start,
488 maxslack, u_data, cut_callback);
489 if (rval) return rval;
490 }
491 if (n*mul + mul < p->base + p->size &&
492 pre_sum + p->sum[n] + p->minpre[n+1] <= maxslack) {
493 rval = psh_enum_work (p, n+1, mul, pre_sum + p->sum[n], cut_start,
494 maxslack, u_data, cut_callback);
495 if (rval) return rval;
496 }
497 return 0;
498 }
499
500