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 /* Exact Seperation of Blossoms (Padberg-Rao) */
19 /* */
20 /* TSP CODE */
21 /* */
22 /* Written by: Applegate, Bixby, Chvatal, and Cook */
23 /* Date: Spring 1989 (Bico) */
24 /* Modified January 11, 1999 (Bico) */
25 /* */
26 /* EXPORTED FUNCTIONS: */
27 /* */
28 /* int CCtsp_exactblossom (CCtsp_lpcut_in **cuts, int *cutcount, */
29 /* int ncount, int ecount, int *elist, double *x, */
30 /* CCrandstate *rstate) */
31 /* RUNS Padberg-Rao to seperate over blossom inequalities. */
32 /* -cuts (new cutting plans will be added to front of this list) */
33 /* -cutcount will return the number of new cuts found (can be NULL) */
34 /* -ncount is the number of nodes */
35 /* -ecount is the number of edges */
36 /* -elist is the edge list in node node format */
37 /* -x is an lp solution vector */
38 /* NOTES: */
39 /* The exactblossom code was written very early in our TSP project. */
40 /* In January 1999 it was updated to fit into the current concorde, */
41 /* but the guts of the code are still written in our old sloppy */
42 /* style. This is a good candidate for a rewrite (big speedups are */
43 /* probably possible without too much effort). */
44 /* */
45 /* int CCtsp_fastblossom (CCtsp_lpcut_in **cuts, int *cutcount, */
46 /* int ncount, int ecount, int *elist, double *x) */
47 /* FINDS blossoms by looking at 0 < x < 1 graph for connected comps */
48 /* meeting an odd number of 1 edges. */
49 /* */
50 /* int CCtsp_ghfastblossom (CCtsp_lpcut_in **cuts, int *cutcount, */
51 /* int ncount, int ecount, int *elist, double *x) */
52 /* FINDS blossoms using a heuristic described by Groetschel and */
53 /* Holland. It works with the 0 < x < 1-EPS (with EPS = .3) graph, */
54 /* builds components, and picks a greedy set of teeth. */
55 /* */
56 /****************************************************************************/
57
58 /****************************************************************************/
59 /* */
60 /* The code is set up to allow blossoms with single teeth. To forbid */
61 /* these, undef ALLOW_SINGLE_TEETH. */
62 /* */
63 /****************************************************************************/
64
65
66 #include "machdefs.h"
67 #include "util.h"
68 #include "cut.h"
69 #include "tsp.h"
70
71 #define ALLOW_SINGLE_TEETH
72 #define BLOTOLERANCE .01
73 #define OTHEREND(e,n) ((e)->ends[0] == (n) ? (e)->ends[1] \
74 : (e)->ends[0])
75
76 typedef struct edge {
77 struct node *ends[2];
78 struct node *splitter;
79 double x;
80 int magiclabel;
81 struct edge *next;
82 } edge;
83
84 typedef struct edgeptr {
85 struct edge *this;
86 struct edgeptr *next;
87 } edgeptr;
88
89 typedef struct node {
90 edgeptr *adj;
91 int magiclabel;
92 struct node *next, *prev;
93 struct node *oddnode;
94 edge *pe;
95 int mark;
96 int num;
97 int name;
98 } node;
99
100 typedef struct nodeptr {
101 struct node *this;
102 struct nodeptr *next;
103 } nodeptr;
104
105 typedef struct graph {
106 int ncount;
107 node *nodelist;
108 int ecount;
109 edge *edgelist;
110 node *pseudonodelist;
111 edge *pseudoedgelist;
112 int magicnum;
113 node pseudonodedummy;
114 edge pseudoedgedummy;
115 CCptrworld edge_world;
116 CCptrworld edgeptr_world;
117 CCptrworld node_world;
118 CCptrworld nodeptr_world;
119 } graph;
120
121 typedef struct toothobj {
122 int in;
123 int out;
124 } toothobj;
125
126 CC_PTRWORLD_ROUTINES (edge, edgealloc, edge_bulkalloc, edgefree)
127 CC_PTRWORLD_LEAKS_ROUTINE (edge, edge_check_leaks, x, double)
128
129 CC_PTRWORLD_LIST_ROUTINES (edgeptr, edge *, edgeptralloc, edgeptr_bulkalloc,
130 edgeptrfree, edgeptr_listadd, edgeptr_listfree)
131 CC_PTRWORLD_LEAKS_ROUTINE (edgeptr, edgeptr_check_leaks, this, edge *)
132
133 CC_PTRWORLD_ROUTINES (node, nodealloc, node_bulkalloc, nodefree)
134 CC_PTRWORLD_LEAKS_ROUTINE (node, node_check_leaks, name, int)
135
136 CC_PTRWORLD_LIST_ROUTINES (nodeptr, node *, nodeptralloc, nodeptr_bulkalloc,
137 nodeptrfree, nodeptr_listadd, nodeptr_listfree)
138 CC_PTRWORLD_LEAKS_ROUTINE (nodeptr, nodeptr_check_leaks, this, node *)
139
140
141 static void
142 blolink_init (graph *G),
143 blolink_free (graph *G),
144 removedegreezero (graph *G),
145 splitem (graph *G, node *n),
146 splitedge (graph *G, edge *e, node *n),
147 destroysplitgraph (graph *G),
148 free_adj (graph *G),
149 freesplitedges (graph *G),
150 freesplitters (graph *G),
151 markcuttree_cut (CC_GHnode *n, int v, node **names),
152 initgraph (graph *G),
153 freegraph (graph *G);
154
155 static int
156 buildadj_from_pseudoedgelist (graph *G),
157 searchtree (graph *G, CC_GHnode *n, node **names, CCtsp_lpcut_in **cuts,
158 int *cutcount),
159 loadcuttree_blossom (graph *G, int v, CCtsp_lpcut_in **cuts,
160 int *cutcount),
161 work_blossom (graph *G, nodeptr *handle, int tcount, edgeptr *teeth,
162 CCtsp_lpcut_in **cuts, int *cutcount),
163 add_blossom (graph *G, int hcount, int *handle, int tcount,
164 toothobj *teeth, CCtsp_lpcut_in **cuts, int *cutcount),
165 cuttree_tooth (edge *e, int v),
166 oneend (edge *e, int v),
167 buildgraph (graph *G, int ncount, int ecount, int *elist, double *x),
168 grab_component (graph *G, node *n, int label, nodeptr **comp, double lbd,
169 double ubd),
170 grow_teeth (graph *G, nodeptr *handle, CCtsp_lpcut_in **cuts, int *cutcount),
171 grow_ghteeth (graph *G, nodeptr *handle, CCtsp_lpcut_in **cuts,
172 int *cutcount);
173
174
175 #define ONEMINUS 0.999999
176 #define ZEROPLUS 0.000001
177
178
CCtsp_exactblossom(CCtsp_lpcut_in ** cuts,int * cutcount,int ncount,int ecount,int * elist,double * x,CCrandstate * rstate)179 int CCtsp_exactblossom (CCtsp_lpcut_in **cuts, int *cutcount, int ncount,
180 int ecount, int *elist, double *x, CCrandstate *rstate)
181 {
182 int i, k;
183 node *n;
184 edge *e;
185 edgeptr *ep;
186 CC_GHtree T;
187 int gncount, gecount, markcount;
188 int *marks = (int *) NULL;
189 int *gelist = (int *) NULL;
190 node **names = (node **) NULL;
191 double *gecap = (double *) NULL;
192 graph G;
193 int rval = 0;
194
195 /*
196 printf ("CCtsp_exactblossom (%d, %d)...\n", ncount, ecount);
197 fflush (stdout);
198 */
199
200 *cutcount = 0;
201 CCcut_GHtreeinit (&T);
202 initgraph (&G);
203 blolink_init (&G);
204
205 rval = buildgraph (&G, ncount, ecount, elist, x);
206 if (rval) {
207 fprintf (stderr, "buildgraph failed\n"); goto CLEANUP;
208 }
209
210 for (i = G.ecount, e = G.edgelist; i; i--, e++) {
211 if (e->x > ONEMINUS) {
212 e->splitter = G.nodelist; /* just to kill the edge */
213 e->ends[0]->mark = 1 - e->ends[0]->mark;
214 e->ends[1]->mark = 1 - e->ends[1]->mark;
215 } else if (e->x < ZEROPLUS) {
216 e->splitter = G.nodelist;
217 } else {
218 e->splitter = (node *) NULL;
219 }
220 }
221
222 G.pseudoedgelist = &G.pseudoedgedummy;
223 G.pseudoedgelist->next = (edge *) NULL;
224
225 G.magicnum++;
226 for (n = G.pseudonodelist->next; n; n = n->next) {
227 if (n->magiclabel != G.magicnum) {
228 splitem (&G, n);
229 }
230 }
231
232 free_adj (&G);
233 rval = buildadj_from_pseudoedgelist (&G);
234 if (rval) goto CLEANUP;
235
236 removedegreezero (&G);
237
238 gncount = 0;
239 gecount = 0;
240 markcount = 0;
241 for (n = G.pseudonodelist->next; n; n = n->next) {
242 for (ep = n->adj; ep; ep = ep->next) {
243 ep->this->magiclabel = 0;
244 gecount++;
245 }
246 n->num = gncount++;
247 if (n->mark) markcount++;
248 }
249 gecount /= 2;
250
251 if (gecount == 0) {
252 /* printf ("No edges in blossom graph\n"); fflush (stdout); */
253 rval = 0; goto CLEANUP;
254 }
255
256 names = CC_SAFE_MALLOC (gncount, node *);
257 gelist = CC_SAFE_MALLOC (2*gecount, int);
258 gecap = CC_SAFE_MALLOC (gecount, double);
259 if (!names || !gelist || !gecap) {
260 fprintf (stderr, "out of memory in Xblossom\n");
261 rval = 1; goto CLEANUP;
262 }
263 if (markcount) {
264 marks = CC_SAFE_MALLOC (markcount, int);
265 if (!marks) {
266 fprintf (stderr, "out of memory in Xblossom\n");
267 rval = 1; goto CLEANUP;
268 }
269 }
270
271 k = 0;
272 markcount = 0;
273 for (i = 0, n = G.pseudonodelist->next; n; n = n->next, i++) {
274 names[i] = n;
275 for (ep = n->adj; ep; ep = ep->next) {
276 if (ep->this->magiclabel == 0) {
277 gelist[2*k] = ep->this->ends[0]->num;
278 gelist[2*k + 1] = ep->this->ends[1]->num;
279 gecap[k] = ep->this->x;
280 k++;
281 ep->this->magiclabel = 1;
282 }
283 }
284 if (n->mark) marks[markcount++] = n->num;
285 }
286
287 if (markcount > 0) {
288 rval = CCcut_gomory_hu (&T, gncount, gecount, gelist, gecap,
289 markcount, marks, rstate);
290 if (rval) {
291 fprintf (stderr, "CCcut_gomory_hu failed\n"); goto CLEANUP;
292 }
293 }
294
295 CC_IFFREE (marks, int);
296 CC_IFFREE (gelist, int);
297 CC_IFFREE (gecap, double);
298
299 if (T.root) {
300 rval = searchtree (&G, T.root, names, cuts, cutcount);
301 if (rval) {
302 fprintf (stderr, "searchtree failed\n"); goto CLEANUP;
303 }
304 }
305
306 CLEANUP:
307
308 CCcut_GHtreefree (&T);
309 destroysplitgraph (&G);
310 freegraph (&G);
311 blolink_free (&G);
312
313 CC_IFFREE (marks, int);
314 CC_IFFREE (gelist, int);
315 CC_IFFREE (gecap, double);
316 CC_IFFREE (names, node *);
317
318 return rval;
319 }
320
blolink_init(graph * G)321 static void blolink_init (graph *G)
322 {
323 CCptrworld_init (&G->edge_world);
324 CCptrworld_init (&G->edgeptr_world);
325 CCptrworld_init (&G->node_world);
326 CCptrworld_init (&G->nodeptr_world);
327 }
328
blolink_free(graph * G)329 static void blolink_free (graph *G)
330 {
331 int total, onlist;
332
333 if (edge_check_leaks (&G->edge_world, &total, &onlist)) {
334 fprintf (stderr, "WARNING: %d outstanding BLOSSOM-edges\n",
335 total - onlist);
336 }
337 CCptrworld_delete (&G->edge_world);
338
339 if (edgeptr_check_leaks (&G->edgeptr_world, &total, &onlist)) {
340 fprintf (stderr, "WARNING: %d outstanding BLOSSOM-edgeptrs\n",
341 total - onlist);
342 }
343 CCptrworld_delete (&G->edgeptr_world);
344
345 if (node_check_leaks (&G->node_world, &total, &onlist)) {
346 fprintf (stderr, "WARNING: %d outstanding BLOSSOM-nodes\n",
347 total - onlist);
348 }
349 CCptrworld_delete (&G->node_world);
350
351 if (nodeptr_check_leaks (&G->nodeptr_world, &total, &onlist)) {
352 fprintf (stderr, "WARNING: %d outstanding BLOSSOM-nodeptrs\n",
353 total - onlist);
354 }
355 CCptrworld_delete (&G->nodeptr_world);
356 }
357
buildadj_from_pseudoedgelist(graph * G)358 static int buildadj_from_pseudoedgelist (graph *G)
359 {
360 edge *e;
361 int rval;
362
363 for (e = G->pseudoedgelist->next; e; e = e->next) {
364 rval = edgeptr_listadd (&(e->ends[0]->adj), e, &G->edgeptr_world);
365 if (rval) return rval;
366 rval = edgeptr_listadd (&(e->ends[1]->adj), e, &G->edgeptr_world);
367 if (rval) return rval;
368 }
369 return 0;
370 }
371
removedegreezero(graph * G)372 static void removedegreezero (graph *G)
373 {
374 node *n;
375
376 for (n = G->pseudonodelist->next; n; n = n->next) {
377 if (n->adj == (edgeptr *) NULL) {
378 if (n->next != (node *) NULL) {
379 n->next->prev = n->prev;
380 n->prev->next = n->next;
381 } else {
382 n->prev->next = (node *) NULL;
383 }
384 }
385 }
386 }
387
splitem(graph * G,node * n)388 static void splitem (graph *G, node *n)
389 {
390 edgeptr *ep;
391 edge *e, *last;
392 node *child;
393
394 n->magiclabel = G->magicnum;
395 for (ep = n->adj; ep; ep = ep->next) {
396 e = ep->this;
397 if (!e->splitter) {
398 child = OTHEREND (e, n);
399 if (child->magiclabel != G->magicnum)
400 splitem (G, child);
401 }
402 }
403 for (ep = n->adj, last = (edge *) NULL; ep; ep = ep->next) {
404 e = ep->this;
405 if (!e->splitter) {
406 if (last)
407 splitedge (G, last, n);
408 last = e;
409 }
410 }
411
412 if (last) {
413 if (n->mark)
414 splitedge (G, last, n);
415 else
416 splitedge (G, last, OTHEREND (last, n));
417 }
418 }
419
splitedge(graph * G,edge * e,node * n)420 static void splitedge (graph *G, edge *e, node *n)
421 {
422 node *n1;
423 edge *e1;
424
425 n->mark = 1 - n->mark;
426
427 n1 = nodealloc (&G->node_world);
428 n1->magiclabel = 0;
429 e->splitter = n1;
430 n1->pe = e;
431 n1->oddnode = n;
432 n1->mark = 1;
433 n1->adj = (edgeptr *) NULL;
434
435 n1->next = G->pseudonodelist->next;
436 n1->prev = G->pseudonodelist;
437 G->pseudonodelist->next->prev = n1;
438 G->pseudonodelist->next = n1;
439
440
441 e1 = edgealloc (&G->edge_world);
442 e1->ends[0] = n;
443 e1->ends[1] = n1;
444 e1->x = 1.0 - e->x;
445 e1->next = G->pseudoedgelist->next;
446 G->pseudoedgelist->next = e1;
447
448
449 e1 = edgealloc (&G->edge_world);
450 e1->ends[0] = OTHEREND (e, n);
451 e1->ends[1] = n1;
452 e1->x = e->x;
453 e1->next = G->pseudoedgelist->next;
454 G->pseudoedgelist->next = e1;
455 }
456
destroysplitgraph(graph * G)457 static void destroysplitgraph (graph *G)
458 {
459 free_adj (G);
460 freesplitedges (G);
461 freesplitters (G);
462 }
463
free_adj(graph * G)464 static void free_adj (graph *G)
465 {
466 node *n;
467
468 for (n = G->pseudonodelist->next; n; n = n->next) {
469 edgeptr_listfree (&G->edgeptr_world, n->adj);
470 n->adj = (edgeptr *) NULL;
471 }
472 }
473
freesplitedges(graph * G)474 static void freesplitedges (graph *G)
475 {
476 edge *e, *enext;
477
478 for (e = G->pseudoedgelist->next; e; e = enext) {
479 enext = e->next;
480 edgefree (&G->edge_world, e);
481 }
482 }
483
freesplitters(graph * G)484 static void freesplitters (graph *G)
485 {
486 node *n, *nnext, *prev;
487
488 for (n = G->pseudonodelist->next, prev = G->pseudonodelist; n; n = nnext) {
489 nnext = n->next;
490 if (n->pe) {
491 prev->next = nnext;
492 if (nnext) nnext->prev = prev;
493 nodefree (&G->node_world, n);
494 } else {
495 prev = n;
496 }
497 }
498 }
499
searchtree(graph * G,CC_GHnode * n,node ** names,CCtsp_lpcut_in ** cuts,int * cutcount)500 static int searchtree (graph *G, CC_GHnode *n, node **names,
501 CCtsp_lpcut_in **cuts, int *cutcount)
502 {
503 CC_GHnode *c;
504 int rval = 0;
505
506 if (n->ndescendants % 2 == 1 && n->ndescendants > 1 ) {
507 if (n->cutval < 1.0 - BLOTOLERANCE) {
508 G->magicnum++;
509 markcuttree_cut (n, G->magicnum, names);
510 rval = loadcuttree_blossom (G, G->magicnum, cuts, cutcount);
511 if (rval) {
512 fprintf (stderr, "loadcuttree_blossom failed\n");
513 goto CLEANUP;
514 }
515 }
516 }
517 for (c = n->child; c; c = c->sibling) {
518 rval = searchtree (G, c, names, cuts, cutcount);
519 if (rval) {
520 goto CLEANUP;
521 }
522 }
523
524 CLEANUP:
525
526 return rval;
527 }
528
markcuttree_cut(CC_GHnode * n,int v,node ** names)529 static void markcuttree_cut (CC_GHnode *n, int v, node **names)
530 {
531 CC_GHnode *c;
532 int i;
533
534 for (i = 0; i < n->listcount; i++) {
535 names[n->nlist[i]]->magiclabel = v;
536 }
537 for (c = n->child; c; c = c->sibling) markcuttree_cut (c, v, names);
538 }
539
loadcuttree_blossom(graph * G,int v,CCtsp_lpcut_in ** cuts,int * cutcount)540 static int loadcuttree_blossom (graph *G, int v, CCtsp_lpcut_in **cuts,
541 int *cutcount)
542 {
543 int i, rval = 0;
544 edge *e;
545 node *n;
546 nodeptr *handle = (nodeptr *) NULL;
547 edgeptr *teeth = (edgeptr *) NULL;
548 int tcount = 0;
549 int hcount = 0;
550
551 for (i = G->ecount, e = G->edgelist; i; i--, e++) {
552 if (oneend (e, v) && cuttree_tooth (e, v)) {
553 rval = edgeptr_listadd (&teeth, e, &G->edgeptr_world);
554 if (rval) goto CLEANUP;
555 tcount++;
556 }
557 }
558 for (i = G->ncount, n = G->nodelist; i; i--, n++) {
559 if (n->magiclabel == v) {
560 rval = nodeptr_listadd (&handle, n, &G->nodeptr_world);
561 if (rval) goto CLEANUP;
562 hcount++;
563 }
564 }
565
566 if (hcount >= 3 && (tcount % 2 == 1)) {
567 rval = work_blossom (G, handle, tcount, teeth, cuts,
568 cutcount);
569 if (rval) {
570 fprintf (stderr, "work_blossom failed\n"); goto CLEANUP;
571 }
572 }
573
574 CLEANUP:
575
576 edgeptr_listfree (&G->edgeptr_world, teeth);
577 nodeptr_listfree (&G->nodeptr_world, handle);
578
579 return rval;
580 }
581
work_blossom(graph * G,nodeptr * handle,int tcount,edgeptr * teeth,CCtsp_lpcut_in ** cuts,int * cutcount)582 static int work_blossom (graph *G, nodeptr *handle, int tcount,
583 edgeptr *teeth, CCtsp_lpcut_in **cuts, int *cutcount)
584 {
585 edge *e;
586 edgeptr *ep;
587 nodeptr *np;
588 toothobj *t = (toothobj *) NULL;
589 toothobj *newteeth = (toothobj *) NULL;
590 int i, k, newhcount, newtcount, rval = 0;
591 int *del = (int *) NULL;
592 int *add = (int *) NULL;
593 int *hit = (int *) NULL;
594 int *newhandle = (int *) NULL;
595
596 /* Clean up intersecting teeth */
597
598 t = CC_SAFE_MALLOC (tcount, toothobj);
599 hit = CC_SAFE_MALLOC (G->ncount, int);
600 del = CC_SAFE_MALLOC (G->ncount, int);
601 add = CC_SAFE_MALLOC (G->ncount, int);
602 if (!t || !hit || !del || !add) {
603 fprintf (stderr, "out of memory in work_blossom\n");
604 }
605
606 G->magicnum++;
607 for (np = handle; np; np = np->next) {
608 np->this->magiclabel = G->magicnum;
609 }
610 for (ep = teeth, i = 0; ep; ep = ep->next, i++) {
611 e = ep->this;
612 if (e->ends[0]->magiclabel == e->ends[1]->magiclabel) goto CLEANUP;
613 if (e->ends[0]->magiclabel == G->magicnum) {
614 t[i].in = e->ends[0]->name;
615 t[i].out = e->ends[1]->name;
616 } else {
617 t[i].in = e->ends[1]->name;
618 t[i].out = e->ends[0]->name;
619 }
620 }
621
622 for (i = 0; i < tcount; i++) {
623 hit[t[i].in] = 0;
624 del[t[i].in] = 0;
625 }
626 for (i = 0; i < tcount; i++) {
627 if (hit[t[i].in]) del[t[i].in] = 1;
628 else hit[t[i].in] = 1;
629 }
630
631 for (i = 0; i < tcount; i++) {
632 hit[t[i].out] = 0;
633 add[t[i].out] = 0;
634 }
635 for (i = 0; i < tcount; i++) {
636 if (hit[t[i].out]) add[t[i].out] = 1;
637 else hit[t[i].out] = 1;
638 }
639
640 for (i = 0; i < G->ncount; i++) hit[i] = 0;
641
642 for (np = handle; np; np = np->next) {
643 hit[np->this->name] = 1;
644 }
645 for (i = 0; i < tcount; i++) {
646 if (del[t[i].in]) hit[t[i].in] = 0;
647 }
648 for (i = 0; i < tcount; i++) {
649 if (add[t[i].out]) hit[t[i].out] = 1;
650 }
651
652 for (i = 0, newhcount = 0; i < G->ncount; i++) {
653 if (hit[i]) newhcount++;
654 }
655 for (i = 0, newtcount = 0; i < tcount; i++) {
656 if (hit[t[i].in] != hit[t[i].out]) newtcount++;
657 }
658
659 #ifdef ALLOW_SINGLE_TEETH
660 if (newhcount >= 3 && (newtcount % 2 == 1)) {
661 #else
662 if (newhcount >= 3 && (newtcount % 2 == 1) && newtcount >= 3) {
663 #endif
664 newhandle = CC_SAFE_MALLOC (newhcount, int);
665 newteeth = CC_SAFE_MALLOC (newtcount, toothobj);
666 if (!newhandle || !newteeth) {
667 fprintf (stderr, "out of memory in work_blossom\n");
668 rval = 1; goto CLEANUP;
669 }
670 k = 0;
671 for (i = 0; i < G->ncount; i++) {
672 if (hit[i]) newhandle[k++] = i;
673 }
674 k = 0;
675 for (i = 0; i < tcount; i++) {
676 if (hit[t[i].in] != hit[t[i].out]) {
677 newteeth[k].in = t[i].in;
678 newteeth[k].out = t[i].out;
679 k++;
680 }
681 }
682 rval = add_blossom (G, newhcount, newhandle, newtcount, newteeth, cuts,
683 cutcount);
684 CCcheck_rval (rval, "add_blossom failed");
685 }
686
687 CLEANUP:
688
689 CC_IFFREE (t, toothobj);
690 CC_IFFREE (hit, int);
691 CC_IFFREE (del, int);
692 CC_IFFREE (add, int);
693 CC_IFFREE (newhandle, int);
694 CC_IFFREE (newteeth, toothobj);
695
696 return rval;
697 }
698
699 static int add_blossom (graph *G, int hcount, int *handle, int tcount,
700 toothobj *teeth, CCtsp_lpcut_in **cuts, int *cutcount)
701 {
702 int itooth[2];
703 CCtsp_lpcut_in *lc;
704 int i, rval = 0;
705
706 lc = CC_SAFE_MALLOC (1, CCtsp_lpcut_in);
707 CCcheck_NULL (lc, "out of memory in add_blossom");
708 CCtsp_init_lpcut_in (lc);
709
710 rval = CCtsp_create_lpcliques (lc, tcount+1);
711 CCcheck_rval (rval, "CCtsp_creat_lpcliques_failed");
712
713 rval = CCtsp_array_to_lpclique (handle, hcount, &(lc->cliques[0]));
714 CCcheck_rval (rval, "CCtsp_array_to_lpclique_failed");
715
716 for (i = 0; i < tcount; i++) {
717 itooth[0] = teeth[i].in;
718 itooth[1] = teeth[i].out;
719 rval = CCtsp_array_to_lpclique (itooth, 2, &(lc->cliques[i+1]));
720 CCcheck_rval (rval, "CCtsp_array_to_lpclique_failed");
721 }
722
723 lc->rhs = CCtsp_COMBRHS (lc);
724 lc->sense = 'G';
725 lc->branch = 0;
726
727 rval = CCtsp_construct_skeleton (lc, G->ncount);
728 if (rval) {
729 fprintf (stderr, "CCtsp_construct_skeleton failed\n"); goto CLEANUP;
730 }
731
732 lc->next = *cuts;
733 *cuts = lc;
734 /* CCtsp_print_lpcut_in (lc); */
735
736 (*cutcount)++;
737
738 CLEANUP:
739
740 if (rval) {
741 CCtsp_free_lpcut_in (lc);
742 CC_IFFREE (lc, CCtsp_lpcut_in);
743 }
744
745 return rval;
746 }
747
748 static int cuttree_tooth (edge *e, int v)
749 {
750 if (e->x > ONEMINUS) return 1;
751 if (e->x < ZEROPLUS) return 0;
752
753 if (e->splitter->magiclabel == v) {
754 if (e->splitter->oddnode->magiclabel == v) return 0;
755 else return 1;
756 } else {
757 if (e->splitter->oddnode->magiclabel == v) return 1;
758 else return 0;
759 }
760 }
761
762 static int oneend (edge *e, int v)
763 {
764 if ((e->ends[0]->magiclabel == v && e->ends[1]->magiclabel != v) ||
765 (e->ends[0]->magiclabel != v && e->ends[1]->magiclabel == v)) {
766 return 1;
767 } else {
768 return 0;
769 }
770 }
771
772 static int buildgraph (graph *G, int ncount, int ecount, int *elist, double *x)
773 {
774 node *n;
775 edge *e;
776 int i, k, n1, n2;
777 int rval = 0;
778
779 initgraph (G);
780
781 G->nodelist = CC_SAFE_MALLOC (ncount, node);
782 G->edgelist = CC_SAFE_MALLOC (ecount, edge);
783
784 if (!G->nodelist || !G->edgelist) {
785 fprintf (stderr, "out of memory in buildgraph\n");
786 rval = 1; goto CLEANUP;
787 }
788 G->ncount = ncount;
789 G->ecount = ecount;
790
791 for (i = 0; i < ncount; i++) {
792 G->nodelist[i].adj = (edgeptr *) NULL;
793 G->nodelist[i].magiclabel = 0;
794 G->nodelist[i].name = i;
795 G->nodelist[i].pe = (edge *) NULL;
796 G->nodelist[i].mark = 0;
797 }
798 for (i = 0, e = G->edgelist, k = 0; i < ecount; i++, e++) {
799 n1 = elist[k++];
800 n2 = elist[k++];
801 e->ends[0] = G->nodelist + n1;
802 e->ends[1] = G->nodelist + n2;
803 e->x = x[i];
804 e->magiclabel = 0;
805 }
806
807 for (i = ecount, e = G->edgelist; i; i--, e++) {
808 rval = edgeptr_listadd (&(e->ends[0]->adj), e, &G->edgeptr_world);
809 if (rval) goto CLEANUP;
810 rval = edgeptr_listadd (&(e->ends[1]->adj), e, &G->edgeptr_world);
811 if (rval) goto CLEANUP;
812 }
813
814 G->pseudonodelist = &G->pseudonodedummy;
815 G->pseudonodedummy.prev = (node *) NULL;
816 G->pseudonodedummy.next = G->nodelist;
817 for (i = 0, n = G->nodelist; i < G->ncount; i++, n++) {
818 n->prev = n - 1;
819 n->next = n + 1;
820 }
821 G->nodelist->prev = G->pseudonodelist;
822 G->nodelist[G->ncount - 1].next = (node *) NULL;
823
824 CLEANUP:
825
826 if (rval) freegraph (G);
827 return rval;
828 }
829
830 static void initgraph (graph *G)
831 {
832 if (G) {
833 G->ncount = 0;
834 G->nodelist = (node *) NULL;
835 G->ecount = 0;
836 G->edgelist = (edge *) NULL;
837 G->pseudonodelist = (node *) NULL;
838 G->pseudoedgelist = (edge *) NULL;
839 G->magicnum = 0;
840 }
841 }
842
843 static void freegraph (graph *G)
844 {
845 int i;
846 node *n;
847
848 if (G->nodelist) {
849 for (i = G->ncount, n = G->nodelist; i; i--, n++) {
850 edgeptr_listfree (&G->edgeptr_world, n->adj);
851 n->adj = (edgeptr *) NULL;
852 }
853 CC_FREE (G->nodelist, node);
854 }
855 CC_IFFREE (G->edgelist, edge);
856 }
857
858 int CCtsp_fastblossom (CCtsp_lpcut_in **cuts, int *cutcount, int ncount,
859 int ecount, int *elist, double *x)
860 {
861 graph G;
862 int rval = 0;
863 nodeptr *handle;
864 int i, k;
865
866 *cutcount = 0;
867 initgraph (&G);
868 blolink_init (&G);
869
870 rval = buildgraph (&G, ncount, ecount, elist, x);
871 if (rval) {
872 fprintf (stderr, "buildgraph failed\n"); goto CLEANUP;
873 }
874
875 k = 0;
876 for (i = 0; i < G.ncount; i++) {
877 if (G.nodelist[i].mark == 0) {
878 handle = (nodeptr *) NULL;
879 rval = grab_component (&G, &(G.nodelist[i]), ++k, &handle,
880 ZEROPLUS, ONEMINUS);
881 if (rval) {
882 fprintf (stderr, "grab_component failed\n"); goto CLEANUP;
883 }
884 grow_teeth (&G, handle, cuts, cutcount);
885 nodeptr_listfree (&G.nodeptr_world, handle);
886 }
887 }
888
889 CLEANUP:
890
891 freegraph (&G);
892 blolink_free (&G);
893
894 return rval;
895 }
896
897 static int grab_component (graph *G, node *n, int label, nodeptr **comp,
898 double lbd, double ubd)
899 {
900 node *v;
901 edgeptr *ep;
902 int rval;
903
904 n->mark = label;
905 rval = nodeptr_listadd (comp, n, &G->nodeptr_world);
906 if (rval) return 1;
907
908 for (ep = n->adj; ep; ep = ep->next) {
909 if (ep->this->x > lbd && ep->this->x < ubd) {
910 v = OTHEREND (ep->this, n);
911 if (v->mark == 0) {
912 if (grab_component (G, v, label, comp, lbd, ubd)) {
913 return 1;
914 }
915 }
916 }
917 }
918 return 0;
919 }
920
921 static int grow_teeth (graph *G, nodeptr *handle, CCtsp_lpcut_in **cuts,
922 int *cutcount)
923 {
924 edgeptr *ep, *teeth = (edgeptr *) NULL;
925 nodeptr *np;
926 node *n;
927 int tcount = 0, hcount = 0, rval = 0;
928
929 for (np = handle; np; np = np->next) {
930 hcount++;
931 }
932 if (hcount < 3) goto CLEANUP;
933
934 for (np = handle; np; np = np->next) {
935 n = np->this;
936 for (ep = n->adj; ep; ep = ep->next) {
937 if (ep->this->x >= ONEMINUS) {
938 if (ep->this->ends[0]->mark != ep->this->ends[1]->mark) {
939 rval = edgeptr_listadd (&teeth, ep->this, &G->edgeptr_world);
940 if (rval) goto CLEANUP;
941 tcount++;
942 }
943 }
944 }
945 }
946
947 if (tcount % 2) {
948 rval = work_blossom (G, handle, tcount, teeth, cuts, cutcount);
949 if (rval) {
950 fprintf (stderr, "work_blossom failed\n"); goto CLEANUP;
951 }
952 }
953
954 CLEANUP:
955
956 edgeptr_listfree (&G->edgeptr_world, teeth);
957 return rval;
958 }
959
960 #define GH_EPS 0.3
961
962 int CCtsp_ghfastblossom (CCtsp_lpcut_in **cuts, int *cutcount, int ncount,
963 int ecount, int *elist, double *x)
964 {
965 graph G;
966 int rval = 0;
967 nodeptr *handle;
968 int i, k;
969
970 /* NOTE: Groetchel and Holland use a lowerbound of GH_EPS for the */
971 /* edges allowed in the graph, while we use ZEROPLUS (a much */
972 /* smaller number). */
973
974 *cutcount = 0;
975 initgraph (&G);
976 blolink_init (&G);
977
978 rval = buildgraph (&G, ncount, ecount, elist, x);
979 if (rval) {
980 fprintf (stderr, "buildgraph failed\n"); goto CLEANUP;
981 }
982
983 k = 0;
984 for (i = 0; i < G.ncount; i++) {
985 if (G.nodelist[i].mark == 0) {
986 handle = (nodeptr *) NULL;
987 rval = grab_component (&G, &(G.nodelist[i]), ++k, &handle,
988 ZEROPLUS, 1.0 - GH_EPS);
989 if (rval) {
990 fprintf (stderr, "grab_component failed\n"); goto CLEANUP;
991 }
992 grow_ghteeth (&G, handle, cuts, cutcount);
993 nodeptr_listfree (&G.nodeptr_world, handle);
994 }
995 }
996
997 CLEANUP:
998
999 freegraph (&G);
1000 blolink_free (&G);
1001
1002 return rval;
1003 }
1004
1005 static int grow_ghteeth (graph *G, nodeptr *handle, CCtsp_lpcut_in **cuts,
1006 int *cutcount)
1007 {
1008 edgeptr *ep, *teeth = (edgeptr *) NULL, *pteeth = (edgeptr *) NULL;
1009 edge *emax = (edge *) NULL, **tlist = (edge **) NULL;
1010 nodeptr *np;
1011 node *n;
1012 int i, ptcount = 0, hcount = 0, rval = 0;
1013 int *tperm = (int *) NULL;
1014 double z = 0.0, xemax = 0.0;
1015 double *xtlist = (double *) NULL;
1016
1017 for (np = handle; np; np = np->next) {
1018 hcount++;
1019 }
1020 if (hcount < 3) goto CLEANUP;
1021
1022 for (np = handle; np; np = np->next) {
1023 n = np->this;
1024 for (ep = n->adj; ep; ep = ep->next) {
1025 if (ep->this->x > ZEROPLUS) {
1026 if (ep->this->ends[0]->mark != ep->this->ends[1]->mark) {
1027 if (ep->this->x >= 1.0 - GH_EPS) {
1028 rval = edgeptr_listadd (&pteeth, ep->this, &G->edgeptr_world);
1029 if (rval) goto CLEANUP;
1030 ptcount++;
1031 } else if (ep->this->x <= ZEROPLUS) {
1032 if (ep->this->x > xemax) {
1033 xemax = ep->this->x;
1034 emax = ep->this;
1035 }
1036 }
1037 } else {
1038 z += ep->this->x;
1039 }
1040 }
1041 }
1042 }
1043 z *= 0.5;
1044
1045 if (ptcount == 0) {
1046 goto CLEANUP;
1047 } else if (ptcount % 2 == 0) {
1048 if (emax) {
1049 rval = edgeptr_listadd (&pteeth, emax, &G->edgeptr_world);
1050 if (rval) goto CLEANUP;
1051 ptcount++;
1052 }
1053 }
1054
1055 tlist = CC_SAFE_MALLOC (ptcount, edge *);
1056 xtlist = CC_SAFE_MALLOC (ptcount, double);
1057 tperm = CC_SAFE_MALLOC (ptcount, int);
1058 if (!tlist || !xtlist || !tperm) {
1059 fprintf (stderr, "out of memory in grow_ghteeth\n");
1060 rval = 1; goto CLEANUP;
1061 }
1062
1063 for (i = 0, ep = pteeth; ep; ep = ep->next, i++) {
1064 tlist[i] = ep->this;
1065 xtlist[i] = -ep->this->x;
1066 tperm[i] = i;
1067 }
1068
1069 CCutil_double_perm_quicksort (tperm, xtlist, ptcount);
1070 if (ptcount % 2 == 0) ptcount--;
1071
1072 z += tlist[tperm[0]]->x;
1073 rval = edgeptr_listadd (&teeth, tlist[tperm[0]], &G->edgeptr_world);
1074 if (rval) goto CLEANUP;
1075 i = 1;
1076
1077 if (z > (double) hcount + (double) ((i - 1)/2) + BLOTOLERANCE) {
1078 rval = work_blossom (G, handle, i, teeth, cuts, cutcount);
1079 if (rval) {
1080 fprintf (stderr, "work_blossom failed\n"); goto CLEANUP;
1081 }
1082 } else {
1083 while (i < ptcount) {
1084 z += tlist[tperm[i]]->x;
1085 rval = edgeptr_listadd (&teeth, tlist[tperm[i]], &G->edgeptr_world);
1086 if (rval) goto CLEANUP;
1087 i++;
1088 z += tlist[tperm[i]]->x;
1089 rval = edgeptr_listadd (&teeth, tlist[tperm[i]], &G->edgeptr_world);
1090 if (rval) goto CLEANUP;
1091 i++;
1092
1093 if (z > (double) hcount + (double) ((i - 1)/2) + BLOTOLERANCE) {
1094 rval = work_blossom (G, handle, i, teeth, cuts, cutcount);
1095 if (rval) {
1096 fprintf (stderr, "work_blossom failed\n"); goto CLEANUP;
1097 }
1098 break;
1099 }
1100 }
1101 }
1102
1103 CLEANUP:
1104
1105 edgeptr_listfree (&G->edgeptr_world, pteeth);
1106 edgeptr_listfree (&G->edgeptr_world, teeth);
1107 CC_IFFREE (tlist, edge *);
1108 CC_IFFREE (xtlist, double);
1109 CC_IFFREE (tperm, int);
1110 return rval;
1111 }
1112
1113