1 /************************************************************************/
2 /* */
3 /* ROUTINE FOR WEIGHTED PERFECT MATCHING PROBLEMS */
4 /* */
5 /* Written by: A. Rohe */
6 /* Date: November 7, 1995 (rohe) */
7 /* November 12, 1995 (rohe) */
8 /* February 7, 1996 (rohe - short look at program) */
9 /* July, 1996 (rohe - change from tree to forest) */
10 /* October, 1996 (bico - change to concorde format) */
11 /* */
12 /* */
13 /* EXPORTED FUNCTION: */
14 /* int perfect_match (int ncount, CCdatagroup *dat, int ecount, */
15 /* int **elist, int **elen, char *blo_filename, */
16 /* char *mat_filename, int just_fractional, int no_fractional, */
17 /* int use_all_trees, int use_tree, int partialprice, */
18 /* double *totalzeit) */
19 /* COMPUTES a minumum weight perfect matching. */
20 /* -ncount is the number of nodes in the graph */
21 /* -dat gives the data to generate edge lengths (if the */
22 /* matching is over the complete graph), it can be NULL */
23 /* -ecount is the number of edges in the initial edge set */
24 /* -elist is a pointer to an array of edges in end end format, */
25 /* this array will be freed by perfect_match to save on */
26 /* memory */
27 /* -elen is a pointer to an array of edge lengths for the edges */
28 /* in elist, it will be freed by the perfect_match */
29 /* -blo_filename is the name of a file (it can be NULL), if it */
30 /* is not NULL then the blossoms and weights will be written */
31 /* to it */
32 /* -mat_filename is the name of a file (it can be NULL), if it */
33 /* is not null, then the matching will be writtien to it */
34 /* -just_fractional should be 0 to compute a matching, and 1 */
35 /* to only compute a fractional matching */
36 /* -if no_fractional is nonzero, then a fractional jumpstart */
37 /* will not be used */
38 /* -if use_all_trees is 1, then a common delta is */
39 /* computed for all trees in the forest; if it is 2, then */
40 /* a single tree is grown */
41 /* -if partialprice is > 0, then the pricing routine will first */
42 /* price the partialprice nearest neighbors before moving to */
43 /* the complete edge set */
44 /* -totalzeit will return the time used by the matching routine */
45 /* (excluding the time needed for a final check that the */
46 /* matching is indeed optimal, if this checking is turned */
47 /* on) */
48 /* */
49 /* NOTES: */
50 /* */
51 /* Returns 0 if it worked and 1 otherwise (for example, when one */
52 /* of the mallocs failed). The nodes in the graph should be named */
53 /* 0 through #nodes - 1. */
54 /* */
55 /************************************************************************/
56 #include "concorde.h"
57 #include "match.h"
58 #define NDEBUG /* remove this to turn assertions on */
59 #include <assert.h>
60
61 #define CARELESS_REPAIRS
62 /* #define BALL_DERIGS_REPAIRS */
63 #define GREEDY_DUAL_CHANGE
64 #define SWITCH_LEVEL 32
65 #define MAX_BAD_PORTION 0.25
66
67 #define PRINT_LEVEL 0
68 #define SHRINKS_SIZE 1000
69 /* infinity */
70 #define INTEGER_MIN -999999999
71 #define INTEGER_MAX 999999999
72 /* for labels */
73 #define NIX 0
74 #define PLUS 1
75 #define MINUS 2
76 /* for status */
77 #define UNMATCHED 0
78 #define HALVES 1
79 #define MATCHED 2
80 /* for x */
81 #define HALF 2
82
83 int FAILED_NODE;
84
85 /*
86 #define OTHEREND(e,n,list) ((list) + (e)->nod1 + (e)->nod2 - ((n) - (list)))
87 */
88 #define OTHEREND(e,n,list) (((e)->nod1 + (list)) == (n) ? \
89 (e)->nod2 + (list) : (e)->nod1 + (list))
90 #define OTHEREND_INT(e,n) ((e)->nod1 + (e)->nod2 - (n))
91
92 #define FIND_SURFACE(n) { \
93 while((n)->blossom_parent != -1) { \
94 (n)=G->nodelist+(n)->blossom_parent; \
95 } \
96 }
97
98 #define SHRINK_SORT_CONST 3
99 #define BADCOUNT_CONST 10000
100 #define PNODE(n) (G->nodelist + (n))
101 #define PEDGE(e) (G->edgelist + (e))
102
103 #define INODE(n) ((n) - G->nodelist)
104 #define IEDGE(e) ((e) - G->edgelist)
105
106 typedef struct edge {
107 int slack;
108 char mark;
109 char x;
110 int ptrs[2];
111 int nod1;
112 int nod2;
113 int orig_nod1;
114 int orig_nod2;
115 } edge;
116
117 typedef struct node {
118 int edg_list;
119 int matched_edge;
120
121 /* These are for the augmenting path tree */
122 int child;
123 int sibling;
124 int parent;
125 int parent_edge;
126
127 /* These are for the blossom_tree */
128 int blossom_next;
129 int blossom_parent;
130 int blossom_root;
131 int blossom_next_edge;
132 int penultimate;
133
134 int pi;
135 int mark; /* unused pseudo nodes are linked by mark */
136 int tree_root;
137
138 char status;
139 char label;
140 char hit;
141
142 int dummy; /* just to pad to multiple of 8 bytes */
143 } node;
144
145 typedef struct nodeptr {
146 int next;
147 int surf;
148 int delta;
149 int dad;
150 int sum;
151 int status;
152 #ifdef GREEDY_DUAL_CHANGE
153 char tree_processed;
154 int gnext; /* used to pad to multiple of 8 bytes and in reordering */
155 #endif
156 } nodeptr;
157
158 typedef struct shrink_ary_T {
159 int node_i;
160 int node_j;
161 int node_k;
162 int edge;
163 int size;
164 int next;
165 } shrink_ary_T;
166
167 typedef struct shrink_T {
168 shrink_ary_T *ary;
169 int length;
170 } shrink_T;
171
172 typedef struct expand_T {
173 int *node;
174 int length;
175 } expand_T;
176
177 typedef struct stats {
178 int expand_count;
179 int shrink_count;
180 int dualchange_count;
181 int dualzero_count;
182 } stats;
183
184 typedef struct graph {
185 edge *edgelist;
186 node *nodelist;
187 int nnodes;
188 int nedges;
189 int max_nedges;
190 int unused;
191 int unmatched;
192 nodeptr *roots;
193 int unmatched_count;
194 }graph;
195
196 typedef struct srkstuff {
197 expand_T expands;
198 shrink_T shrinks;
199 int shrinks_max;
200 int expands_max;
201 int shrinks_[SHRINKS_SIZE];
202 } srkstuff;
203
204 typedef struct stack {
205 int* ary;
206 int count;
207 } stack;
208
209 #ifdef CC_PROTOTYPE_ANSI
210
211 static void
212 adjust_match (graph *G),
213 print_node (graph *G, node *n),
214 print_edges (graph *G),
215 init_tree (graph *G, node *n),
216 clear_tree (graph *G, node *n),
217 make_cycle_halves (graph *G, node *node_i, node *node_j, node *node_k,
218 edge *e),
219 shrink_tree (graph *G, node *newnode),
220 label_penultimate (graph *G, node *n, int label),
221 fix_matching (graph *G, node *oldnode, node *x),
222 fix_tree (graph *G, node *oldnode, node *x, node *y),
223 flip_cycle (graph *G, node *node_i),
224 augment_path (graph *G, node *node_i, int fractional),
225 vereinigung (graph *G, int x, int y),
226 set_vereinigung (graph *G, node *n),
227 unmatch_edge (graph *G, edge *e),
228 dual_match_len (graph *G, int fractional, double *val),
229 fix_match (graph *G, int blossom);
230
231 static int
232 build_graph (graph *G, int ncount, int ecount, int *elist, int *elen),
233 init (graph *G, srkstuff *srk),
234 match_main_frac (graph *G, stats *scount, srkstuff *srk),
235 match_main (graph *G, stats *scount, srkstuff *srk, int use_all),
236 augment_blossom (graph *G, edge *, int fractional, srkstuff *srk),
237 lower_edges (graph *G, node *old),
238 expand_blossom (graph *G, node *oldnode, stats *scount),
239 lift_edges (graph *G, node *newnode),
240 add_to_tree (graph *G, edge *e, node *node_i, node *node_j, int fractional),
241 checkout_node (graph *G, node *n, int fractional, srkstuff *srk),
242 grow_tree_no_shrink (graph *G, node *n, int fractional, srkstuff *srk),
243 grow_tree (graph *G, node *n, int fractional, stats *scount, srkstuff *srk,
244 int *found_aug),
245 apply_dual_change (graph *G, node *n, int delta),
246 find_parity_sum (graph *G, int n),
247 parity_correction (graph *G, stats *scount),
248 find_single_dual_change (graph *G, node *n),
249 find_dual_change_forest (graph *G, node *n),
250 make_dual_change_forest (graph *G, stats *scount),
251 match_main_more_in_forest (graph *G, stats *scount, srkstuff *srk),
252 match_main_forest (graph *G, stats *scount, srkstuff *srk),
253 match_main_tree (graph *G, stats *scount, srkstuff *srk),
254 make_match (graph *G),
255 test_matching (graph *G, CCdatagroup *dat, int *elen, int fractional,
256 int *bad),
257 price_repair (graph *G, int *finished, CCdatagroup *dat, stats *scount,
258 srkstuff *srk, int partialprice),
259 pricing (graph *G, CCdatagroup *dat, stats *scount, srkstuff *srk,
260 int *badcount, int **badlist, int **badlen, CCtsp_edgegenerator *eg,
261 int usehit),
262 bring_to_surface (graph *G, node *p, stats *scount, srkstuff *srk),
263 add_repair_edge (graph *G, int n1, int n2, int len, stats *scount,
264 srkstuff *srk),
265 run_repair (graph *G, int badcount,int* badlist,int* badlen, stats *scount,
266 srkstuff *srk),
267 add_edges (graph *G, CCdatagroup *dat, int badcount,int* badlist,int* badlen),
268 write_match (graph *G, CCdatagroup *dat, int *elen, char* match_file, int *outp),
269 write_blossom_tree (graph *G, char* blossom_tree_file);
270
271 static node
272 *shrink_blossom (graph *G, node *node_i, node *node_j, node *node_k,
273 edge *e, stats *scount),
274 *common_parent (graph *G, node *node_i, node *node_j, int *size),
275 *find_below (graph *G, node *n, int blossom);
276
277 #else
278
279 static void
280 adjust_match (),
281 print_node (),
282 print_edges (),
283 init_tree (),
284 clear_tree (),
285 make_cycle_halves (),
286 shrink_tree (),
287 label_penultimate (),
288 fix_matching (),
289 fix_tree (),
290 flip_cycle (),
291 augment_path (),
292 vereinigung (),
293 set_vereinigung (),
294 unmatch_edge (),
295 dual_match_len (),
296 fix_match ();
297
298 static int
299 build_graph (),
300 init (),
301 match_main_frac (),
302 match_main (),
303 augment_blossom (),
304 lower_edges (),
305 expand_blossom (),
306 lift_edges (),
307 add_to_tree (),
308 checkout_node (),
309 grow_tree_no_shrink (),
310 grow_tree (),
311 apply_dual_change (),
312 find_parity_sum (),
313 parity_correction (),
314 find_single_dual_change (),
315 find_dual_change_forest (),
316 make_dual_change_forest (),
317 make_match (),
318 match_main_more_in_forest (),
319 match_main_forest (),
320 match_main_tree (),
321 test_matching (),
322 price_repair (),
323 pricing (),
324 add_repair_edge (),
325 bring_to_surface (),
326 run_repair (),
327 add_edges (),
328 write_match (),
329 write_blossom_tree ();
330
331 static node
332 *shrink_blossom (),
333 *common_parent (),
334 *find_below ();
335
336 #endif
337
338 #if PRINT_LEVEL
339 static graph *PG = (graph *) NULL;
340 #endif
341
342 #ifdef CC_PROTOTYPE_ANSI
perfect_match(int ncount,CCdatagroup * dat,int ecount,int ** elist,int ** elen,char * blo_filename,char * mat_filename,int just_fractional,int no_fractional,int use_all_trees,int partialprice,double * totalzeit)343 int perfect_match (int ncount, CCdatagroup *dat, int ecount, int **elist,
344 int **elen, char *blo_filename, char *mat_filename,
345 int just_fractional, int no_fractional,
346 int use_all_trees, int partialprice, double *totalzeit)
347 #else
348 int perfect_match (ncount, dat, ecount, elist, elen, blo_filename,
349 mat_filename, just_fractional, no_fractional,
350 use_all_trees, partialprice, totalzeit)
351 int ncount;
352 CCdatagroup *dat;
353 int ecount;
354 int **elist, **elen;
355 char *blo_filename, *mat_filename;
356 int use_all_trees;
357 int partialprice;
358 double *totalzeit;
359 #endif
360 {
361 graph G;
362 srkstuff srk;
363 stats scount;
364 double val, startzeit;
365 int finished = 1;
366 int bad;
367 int sedges = 0;
368 double zero_zeit = CCutil_zeit ();
369
370 #if PRINT_LEVEL
371 PG = &G;
372 #endif
373
374 int status = 0;
375
376 *totalzeit = 0.0;
377
378 G.edgelist = (edge *) NULL;
379 G.nodelist = (node *) NULL;
380 G.unused = -1;
381 G.unmatched = -1;
382 G.roots = (nodeptr *) NULL;
383
384 srk.expands.node = (int *) NULL;
385 srk.shrinks.ary = (shrink_ary_T *) NULL;
386
387 if (build_graph (&G, ncount, ecount, *elist, *elen)) {
388 // fprintf (stderr, "build_graph failed\n");
389 CC_FREE (*elist, int);
390 CC_FREE (*elen, int);
391 status = 1;
392 goto CLEANUP;
393 }
394 CC_FREE (*elist, int);
395 if (dat) {
396 CC_FREE (*elen, int);
397 }
398 sedges = G.nedges;
399
400 do {
401 // printf (" Starting Init...");
402 fflush (stdout);
403 startzeit = CCutil_zeit ();
404 if (init (&G, &srk)) {
405 // fprintf (stderr, "init failed\n");
406 status = 1;
407 goto CLEANUP;
408 }
409 // printf (" ... ready in %.2f sec !! \n", CCutil_zeit () - startzeit);
410 fflush (stdout);
411
412 if (no_fractional == 0) {
413 // printf (" Matching with %i nodes and %i edges ...\n",
414 // G.nnodes, G.nedges);
415 // fflush (stdout);
416 // printf (" Make a fractional Matching...\n");
417 fflush (stdout);
418 startzeit = CCutil_zeit ();
419 scount.expand_count = 0;
420 scount.shrink_count = 0;
421 scount.dualchange_count = 0;
422 scount.dualzero_count = 0;
423 if (match_main_frac (&G, &scount, &srk)) {
424 // fprintf (stderr, "match_main_frac failed\n");
425 status = 1;
426 goto CLEANUP;
427 }
428 // printf (" ... ready in %.2f sec !! \n", CCutil_zeit () - startzeit);
429 // fflush (stdout);
430 }
431
432 if (just_fractional == 0) {
433 // printf (" Make an integral Matching ...\n");
434 // fflush (stdout);
435 startzeit = CCutil_zeit ();
436 if (make_match (&G)) {
437 // fprintf (stderr, "make_match failed\n");
438 // fflush (stdout);
439 }
440 // printf (" ... ready in %.2f sec !! \n", CCutil_zeit () - startzeit);
441 // fflush (stdout);
442
443 // printf (" Starting Integral Matching-Code...\n");
444 // fflush (stdout);
445 startzeit = CCutil_zeit ();
446 scount.expand_count = 0;
447 scount.shrink_count = 0;
448 scount.dualchange_count = 0;
449 scount.dualzero_count = 0;
450 if (match_main (&G, &scount, &srk, use_all_trees)) {
451 // fprintf (stderr, "match_main failed\n");
452 status = 1;
453 goto CLEANUP;
454 }
455 // printf (" ... ready in %.2f sec !! \n", CCutil_zeit () - startzeit);
456 // fflush (stdout);
457 dual_match_len (&G, 0, &val);
458 // printf (" Dual Matching Length: %.1f\n", val);
459 // fflush (stdout);
460 if (dat) {
461 // printf ("Start Price-Repair Phase ...\n");
462 // fflush (stdout);
463 startzeit = CCutil_zeit ();
464 if (price_repair (&G, &finished, dat, &scount, &srk,
465 partialprice)) {
466 // fprintf (stderr, "pricing failed\n");
467 return 1;
468 }
469 if (finished && partialprice > 0) {
470 // printf ("Nearest %d are correct, now try complete graph\n",
471 // partialprice);
472 // fflush (stdout);
473 if (price_repair (&G, &finished, dat, &scount, &srk, 0)) {
474 // fprintf (stderr, "pricing failed\n");
475 return 1;
476 }
477 }
478 // printf (" ... ready in %.2f sec !! \n", CCutil_zeit () - startzeit);
479 // fflush (stdout);
480 }
481 CC_IFFREE (G.roots, nodeptr);
482 }
483 if (blo_filename != (char *) NULL) {
484 if (write_blossom_tree (&G, blo_filename)) {
485 // fprintf (stderr, "write_blossom_tree failed\n");
486 status = 1;
487 goto CLEANUP;
488 }
489 }
490 CC_IFFREE (srk.expands.node, int);
491 CC_IFFREE (srk.shrinks.ary, shrink_ary_T);
492 } while (!finished);
493
494 // printf (" Adjust found Matching...\n");
495 // fflush (stdout);
496 startzeit = CCutil_zeit ();
497 adjust_match (&G);
498 // printf (" ... ready in %.2f sec !! \n", CCutil_zeit () - startzeit);
499 // fflush (stdout);
500
501 *totalzeit = CCutil_zeit () - zero_zeit;
502 // printf ("Matching Time: %.2f Seconds\n", *totalzeit);
503 // fflush (stdout);
504
505 if (dat) {
506 // printf ("Testing the solution...\n");
507 // fflush (stdout);
508 startzeit = CCutil_zeit ();
509 if (test_matching (&G, dat, *elen, just_fractional, &bad)) {
510 // fprintf (stderr, "test matching failed\n");
511 status = 1;
512 goto CLEANUP;
513 }
514 if (bad) {
515 // printf ("ERROR: A problem was detected in the matching\n");
516 status = 1;
517 goto CLEANUP;
518 }
519 // printf ("Testing Time: %.2f sec\n", CCutil_zeit () - startzeit);
520 // fflush (stdout);
521 }
522
523 // printf ("ZZ Edges: %d (starting) %d (total)\n", sedges, G.nedges);
524
525 if (mat_filename != (char *) NULL) {
526 *elist = (int*)malloc((3*ncount+1)*sizeof(int));// over allocated
527 if (write_match (&G, dat, *elen, mat_filename,*elist)) {
528 status = 1;
529 // fprintf (stderr, "write_match failed\n");
530 goto CLEANUP;
531 }
532 }
533
534 CLEANUP:
535
536 CC_IFFREE (G.nodelist, node);
537 CC_IFFREE (G.edgelist, edge);
538 CC_IFFREE (G.roots, nodeptr);
539 CC_IFFREE (srk.expands.node, int);
540 CC_IFFREE (srk.shrinks.ary, shrink_ary_T);
541 if (!dat) {
542 CC_FREE (*elen, int);
543 }
544
545 return status;
546 }
547
548 #ifdef CC_PROTOTYPE_ANSI
print_node(graph * G,node * n)549 static void print_node (graph *G, node *n)
550 #else
551 static void print_node (G, n)
552 graph *G;
553 node *n;
554 #endif
555 {
556 edge *e;
557 int ei;
558 int c;
559 node *nodelist = G->nodelist;
560 edge *edgelist = G->edgelist;
561
562 printf ("Node %d, Pi %d", (int) INODE(n), n->pi);
563 if (n->label == PLUS) printf (" label +");
564 if (n->label == MINUS) printf (" label -");
565 if (n->label == NIX) printf (" label 0");
566 if (n->status == MATCHED)
567 printf (" MATCHED (%d,%d)\n",
568 edgelist[n->matched_edge].nod1,
569 edgelist[n->matched_edge].nod2);
570 if (n->status == UNMATCHED) printf (" UNMATCHED\n");
571 if (n->status == HALVES)
572 printf (" HALVES (%i,%i)\n",
573 edgelist[n->matched_edge].nod1, edgelist[n->matched_edge].nod2);
574 fflush (stdout);
575 if (n->label != NIX) {
576 printf ("TreeRoot %d ", n->tree_root);
577 if (n->parent)
578 printf ("TreeParent = %d ", n->parent);
579 printf ("Children: ");
580 for (c = n->child; c != -1; c = nodelist[c].sibling)
581 printf ("%d ", c);
582 printf ("\n");
583 }
584 if (n->blossom_parent != -1) {
585 printf ("BlossomNode %i, Parent %i, Next %i\n", (int) INODE(n),
586 n->blossom_parent, n->blossom_next);
587 }
588 if (n->blossom_root != -1) {
589 printf ("BlossomParent %i, Root %i\n", (int) INODE(n), n->blossom_root);
590 }
591 for (ei = n->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
592 e = PEDGE(ei/2);
593 printf ("| Edge (%i,%i)", e->nod1, e->nod2);
594 printf ("| %4i is ", e->nod1);
595 if (nodelist[e->nod1].label == PLUS) printf ("+");
596 if (nodelist[e->nod1].label == MINUS) printf ("-");
597 if (nodelist[e->nod1].label == NIX) printf ("0");
598 printf ("| %4i is ", e->nod2);
599 if (nodelist[e->nod2].label == PLUS) printf ("+");
600 if (nodelist[e->nod2].label == MINUS) printf ("-");
601 if (nodelist[e->nod2].label == NIX) printf ("0");
602 printf ("| edg Slack %d, x=%d\n", e->slack, (int) e->x);
603 }
604 }
605
606 #ifdef CC_PROTOTYPE_ANSI
print_edges(graph * G)607 static void print_edges (graph *G)
608 #else
609 static void print_edges (G)
610 graph *G;
611 #endif
612 {
613 int i;
614
615 printf (" %i nodes, %i edges\n", G->nnodes, G->nedges);
616 for (i = 0; i <= 3 * G->nnodes / 2 && i >= 0; i++) {
617 if (G->nodelist[i].edg_list != -1)
618 print_node (G, PNODE(i));
619 }
620 }
621
622 #ifdef CC_PROTOTYPE_ANSI
init(graph * G,srkstuff * srk)623 static int init (graph *G, srkstuff *srk)
624 #else
625 static int init (G, srk)
626 graph *G;
627 srkstuff *srk;
628 #endif
629 {
630 int i, count;
631 node *n;
632 edge *e;
633 int ei;
634 int delta;
635 node *nodelist = G->nodelist;
636 int ncount = G->nnodes;
637 int ecount = G->nedges;
638
639 /* set all the pi s */
640 for (i = 0; i < ncount; i++)
641 nodelist[i].pi = INTEGER_MAX;
642 for (i = ncount; i <= 3 * ncount / 2; i++)
643 nodelist[i].pi = 0;
644
645 for (i = 0, e = G->edgelist; i < ecount; i++, e++) {
646 if (nodelist[e->nod1].pi > e->slack) {
647 nodelist[e->nod1].pi = e->slack;
648 }
649 if (nodelist[e->nod2].pi > e->slack) {
650 nodelist[e->nod2].pi = e->slack;
651 }
652 }
653
654 for (i = 0; i < ncount; i++)
655 nodelist[i].pi /= 2;
656
657 /* calculate the slacks with the pi's and set rest to zero/nothing */
658
659 for (i = 0, e = G->edgelist; i < ecount; i++, e++) {
660 e->slack = e->slack - nodelist[e->nod1].pi - nodelist[e->nod2].pi;
661 e->x = 0;
662 }
663
664 for (i = 0, n = nodelist; i <= 3 * ncount / 2; i++, n++) {
665 n->child = -1;
666 n->sibling = -1;
667 n->parent = -1;
668 n->label = NIX;
669 n->parent_edge = -1;
670 n->matched_edge = -1;
671 n->status = UNMATCHED;
672 n->blossom_root = -1;
673 n->blossom_next = -1;
674 n->blossom_next_edge = -1;
675 n->blossom_parent = -1;
676 n->tree_root = -1;
677 n->hit = 0;
678 }
679 #if 0
680 /* Take all 0-slack edges directly, old version */
681 count = 0;
682 for (i = 0, e = G->edgelist; i < ecount; i++, e++) {
683 if (e->slack == 0) {
684 if ((nodelist[e->nod1].status == UNMATCHED) &&
685 (nodelist[e->nod2].status == UNMATCHED)) {
686 e->x = 1;
687 nodelist[e->nod1].matched_edge = i;
688 nodelist[e->nod2].matched_edge = i;
689 nodelist[e->nod1].status = MATCHED;
690 nodelist[e->nod2].status = MATCHED;
691 count += 2;
692 }
693 }
694 }
695 #else
696 /* Take all 0-slack edges directly & make better pi, new version */
697 count = 0;
698 for (i = 0, n = nodelist; i < ncount; i++, n++) {
699 delta = INTEGER_MAX;
700 if (n->status == UNMATCHED) {
701 for (ei = n->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
702 e = PEDGE(ei/2);
703 if (e->slack <= delta) {
704 if (e->slack == 0) {
705 if ((nodelist[e->nod1].status == UNMATCHED) &&
706 (nodelist[e->nod2].status == UNMATCHED)) {
707 e->x = 1;
708 nodelist[e->nod1].matched_edge = IEDGE(e);
709 nodelist[e->nod2].matched_edge = IEDGE(e);
710 nodelist[e->nod1].status = MATCHED;
711 nodelist[e->nod2].status = MATCHED;
712 count += 2;
713 }
714 }
715 delta = e->slack;
716 }
717 }
718 if (delta != 0) {
719 n->pi += delta;
720 for (ei = n->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
721 e = PEDGE(ei/2);
722 e->slack -= delta;
723 }
724 }
725 }
726 }
727 #endif
728 // printf (" %i nodes ", count); fflush (stdout);
729
730 G->unmatched_count = ncount - count;
731
732 srk->shrinks_max = ncount / 10 + 10; /* was just ncount */
733 srk->shrinks.ary = CC_SAFE_MALLOC (srk->shrinks_max, shrink_ary_T);
734 srk->expands_max = ncount / 10 + 10;
735 srk->expands.node = CC_SAFE_MALLOC (srk->expands_max, int);
736
737 if (!srk->shrinks.ary || !srk->expands.node) {
738 fprintf (stderr, "out of memory in init\n");
739 CC_IFFREE (srk->shrinks.ary, shrink_ary_T);
740 CC_IFFREE (srk->expands.node, int);
741 return 1;
742 }
743 return 0;
744 }
745
746 #ifdef CC_PROTOTYPE_ANSI
init_tree(graph * G,node * n)747 static void init_tree (graph *G, node *n)
748 #else
749 static void init_tree (G, n)
750 graph *G;
751 node *n;
752 #endif
753 {
754 #if PRINT_LEVEL
755 printf (" init_tree %i\n", (int) (n - PG->nodelist)); fflush (stdout);
756 #endif
757 n->child = -1;
758 n->sibling = -1;
759 n->parent = -1;
760 n->parent_edge = -1;
761 n->label = PLUS;
762 }
763
764 #ifdef CC_PROTOTYPE_ANSI
clear_tree(graph * G,node * n)765 static void clear_tree (graph *G, node *n)
766 #else
767 static void clear_tree (G, n)
768 graph *G;
769 node *n;
770 #endif
771 {
772 node *c = n;
773 node *stop = G->nodelist -1;
774
775 while (1) {
776 c->mark = 0;
777 c->label = NIX;
778 c->tree_root = -1; /* ??? */
779
780 /* go to next c */
781 if (c->child!=-1) {
782 c=PNODE(c->child);
783 }
784 else {
785 while (PNODE(c->sibling)==stop) {
786 if (c==n) return ; /* this if is for an n without childs */
787 c=PNODE(c->parent);
788 if (c==n) return ;
789 }
790 c=PNODE(c->sibling);
791 }
792 }
793 }
794
795 /* make_cycle_halves - edge from i to j, k is the "root" */
796 #ifdef CC_PROTOTYPE_ANSI
make_cycle_halves(graph * G,node * node_i,node * node_j,node * node_k,edge * e)797 static void make_cycle_halves (graph *G, node *node_i, node *node_j,
798 node *node_k, edge *e)
799 #else
800 static void make_cycle_halves (G, node_i, node_j, node_k, e)
801 graph *G;
802 node *node_i, *node_j, *node_k;
803 edge *e;
804 #endif
805 {
806 edge *edgelist = G->edgelist;
807 node *parent;
808
809 #if PRINT_LEVEL
810 printf (" Make cycle halves: root %i ends %i and %i:",
811 (int) (node_k - PG->nodelist), (int) (node_i - PG->nodelist),
812 (int) (node_j - PG->nodelist));
813 fflush (stdout);
814 #endif
815
816 /* set matched_edge for node_j */
817 node_j->status = HALVES;
818 node_j->matched_edge = e - edgelist;
819 e->x = HALF;
820
821 /* path from i to k : matched_edges are the parent_edges */
822 for (; node_i != node_k; node_i = PNODE(node_i->parent)) {
823 edgelist[node_i->parent_edge].x = HALF;
824 node_i->matched_edge = node_i->parent_edge;
825 node_i->status = HALVES;
826 }
827
828 /* path from j to k : matched_edges are the "child_edges" */
829 for (; node_j != node_k; node_j = parent) {
830 parent = PNODE(node_j->parent);
831 edgelist[node_j->parent_edge].x = HALF;
832 parent->matched_edge = node_j->parent_edge;
833 parent->status = HALVES;
834 }
835 }
836
837 #ifdef CC_PROTOTYPE_ANSI
lift_edges(graph * G,node * newnode)838 static int lift_edges (graph *G, node *newnode)
839 #else
840 static int lift_edges (G, newnode)
841 graph *G;
842 node *newnode;
843 #endif
844 {
845 node *n, *node_k;
846 node *nodelist = G->nodelist;
847 int nn = INODE(newnode);
848 edge *e;
849 int ei, einext;
850
851 #if PRINT_LEVEL
852 printf ("Lift edges %i\n", nn); fflush (stdout);
853 #endif
854
855 n = node_k = PNODE(newnode->blossom_root);
856 do {
857 ei = n->edg_list;
858 n->edg_list = -1;
859 for (; ei != -1; ei = einext) {
860 e = PEDGE(ei/2);
861 einext = e->ptrs[ei % 2];
862 if ((nodelist[e->nod1].blossom_parent == nn) &&
863 (nodelist[e->nod2].blossom_parent == nn)) {
864 e->ptrs[ei % 2] = n->edg_list;
865 n->edg_list = ei;
866 } else {
867 if (e->nod1 == INODE(n)) { /* nod1 is in new blossom */
868 e->nod1 = nn;
869 } else { /* nod2 is in new blossom */
870 e->nod2 = nn;
871 }
872 e->ptrs[ei % 2] = newnode->edg_list;
873 newnode->edg_list = ei;
874 }
875 }
876 n = PNODE(n->blossom_next);
877 } while (n != node_k);
878 return 0;
879 }
880
881 #ifdef CC_PROTOTYPE_ANSI
shrink_tree(graph * G,node * newnode)882 static void shrink_tree (graph *G, node *newnode)
883 #else
884 static void shrink_tree (G, newnode)
885 graph *G;
886 node *newnode;
887 #endif
888 {
889 node *n, *node_k, *cnode;
890 int c, csibling, p;
891
892 #if PRINT_LEVEL
893 printf (" Shrink-Tree %i\n", (int) (newnode - PG->nodelist));
894 fflush (stdout);
895 #endif
896
897 n = node_k = PNODE(newnode->blossom_root);
898 do {
899 for (c = n->child; c != -1; c = csibling) {
900 cnode = PNODE(c);
901 csibling = cnode->sibling;
902 if (cnode->blossom_parent != INODE(newnode)) {
903 cnode->parent = INODE(newnode);
904 cnode->sibling = newnode->child;
905 newnode->child = c;
906 }
907 }
908 n = PNODE(n->blossom_next);
909 } while (n != node_k);
910
911 /* set other items for newnode */
912 p = node_k->parent;
913
914 newnode->sibling = -1; /* Plus node, has no siblings */
915 newnode->parent = p;
916 newnode->parent_edge = node_k->parent_edge;
917 newnode->matched_edge = node_k->matched_edge;
918 newnode->label = PLUS;
919
920 if (p != -1) /* _p has just one child, so this is ok */
921 G->nodelist[p].child = INODE(newnode);
922 }
923
924 /* shrink blossom - edge from i to j, k is the "root" */
925
926 #ifdef CC_PROTOTYPE_ANSI
shrink_blossom(graph * G,node * node_i,node * node_j,node * node_k,edge * e,stats * scount)927 static node *shrink_blossom (graph *G, node *node_i, node *node_j,
928 node *node_k, edge *e, stats *scount)
929 #else
930 static node *shrink_blossom (G, node_i, node_j, node_k, e, scount)
931 graph *G;
932 node *node_i, *node_j, *node_k;
933 edge *e;
934 stats *scount;
935 #endif
936 {
937 node *newnode;
938 node *parent;
939
940 #if PRINT_LEVEL
941 printf (" Shrink blossom: root %i ends %i and %i:\n",
942 (int) (node_k - PG->nodelist), (int) (node_i - PG->nodelist),
943 (int) (node_j - PG->nodelist));
944 fflush (stdout);
945 #endif
946
947 scount->shrink_count++;
948
949 newnode = PNODE(G->unused);
950 G->unused = newnode->mark;
951 newnode->mark = 0;
952 newnode->edg_list = -1;
953
954 if (node_k->status == UNMATCHED) {
955 G->roots[node_k->tree_root].surf = INODE(newnode);
956 }
957
958 /* set blossom_next,next_edge,parent for node_j */
959 /* blossom_root has to stay the same */
960
961 node_j->blossom_next = INODE(node_i);
962 node_j->blossom_next_edge = IEDGE(e);
963 node_j->blossom_parent = INODE(newnode);
964
965 /* path from i to k : blossom_next_edges are the parent_edges */
966
967 for (; node_i != node_k; node_i = PNODE(node_i->parent)) {
968 node_i->blossom_next_edge = node_i->parent_edge;
969 node_i->blossom_next = node_i->parent;
970 node_i->blossom_parent = INODE(newnode);
971 }
972
973 /* path from j to k : blossom_next_edges are the "child_edges" */
974
975 for (; node_j != node_k; node_j = parent) {
976 parent = PNODE(node_j->parent);
977 parent->blossom_next_edge = node_j->parent_edge;
978 parent->blossom_next = INODE(node_j);
979 parent->blossom_parent = INODE(newnode);
980 }
981
982 newnode->blossom_root = INODE(node_k);
983 newnode->blossom_next = -1;
984 newnode->blossom_next_edge = -1;
985 newnode->blossom_parent = -1;
986
987 newnode->mark = 0;
988 newnode->pi = 0;
989 newnode->status = node_k->status;
990 newnode->tree_root = node_k->tree_root;
991
992 shrink_tree (G, newnode);
993 if (lift_edges (G, newnode)) {
994 fprintf (stderr, "lift_edges failed\n");
995 return (node *) NULL;
996 }
997
998 return newnode;
999 }
1000
1001 #ifdef CC_PROTOTYPE_ANSI
label_penultimate(graph * G,node * n,int label)1002 static void label_penultimate (graph *G, node *n, int label)
1003 #else
1004 static void label_penultimate (G, n, label)
1005 graph *G;
1006 node *n;
1007 int label;
1008 #endif
1009 {
1010 /* all the leafs of the blossom_tree under n get the penultimate label */
1011 node *n2;
1012
1013 if (n->blossom_root == -1) {
1014 n->penultimate = label;
1015 return;
1016 }
1017
1018 n2 = n;
1019 while (1) {
1020 if (n2->blossom_root != -1) {
1021 n2 = PNODE(n2->blossom_root);
1022 } else {
1023 n2->penultimate = label;
1024 while (n2->blossom_next ==
1025 PNODE(n2->blossom_parent)->blossom_root) {
1026 n2 = PNODE(n2->blossom_parent);
1027 if (n2 == n) return;
1028 }
1029 n2 = PNODE(n2->blossom_next);
1030 }
1031 }
1032 }
1033
1034 #ifdef CC_PROTOTYPE_ANSI
lower_edges(graph * G,node * old)1035 static int lower_edges (graph *G, node *old)
1036 #else
1037 static int lower_edges (G, old)
1038 graph *G;
1039 node *old;
1040 #endif
1041 {
1042 edge *e;
1043 int ei, einext;
1044 node *x;
1045
1046 #if PRINT_LEVEL
1047 printf ("Lower edges %i !\n", (int) (old - PG->nodelist));
1048 fflush (stdout);
1049 #endif
1050
1051 for (ei = old->edg_list; ei != -1; ei = einext) {
1052 e = PEDGE(ei/2);
1053 einext = e->ptrs[ei % 2];
1054 if (e->nod1 == INODE(old)) {
1055 x = PNODE(e->orig_nod1);
1056 e->nod1 = x->penultimate;
1057 } else {
1058 x = PNODE(e->orig_nod2);
1059 e->nod2 = x->penultimate;
1060 }
1061 e->ptrs[ei % 2] = PNODE(x->penultimate)->edg_list;
1062 PNODE(x->penultimate)->edg_list = ei;
1063 }
1064 return 0;
1065 }
1066
1067 #ifdef CC_PROTOTYPE_ANSI
fix_matching(graph * G,node * oldnode,node * x)1068 static void fix_matching (graph *G, node *oldnode, node *x)
1069 #else
1070 static void fix_matching (G, oldnode, x)
1071 graph *G;
1072 node *oldnode, *x;
1073 #endif
1074 {
1075 node *n, *memo;
1076 node *nodelist = G->nodelist;
1077 edge *e;
1078
1079 #if PRINT_LEVEL
1080 printf (" Fix Matching %i; x = %i\n",
1081 (int) (oldnode - PG->nodelist), (int) (x - PG->nodelist));
1082 fflush (stdout);
1083 #endif
1084
1085 if (x->blossom_next_edge == x->matched_edge) {
1086 n = x;
1087 memo = PNODE(oldnode->blossom_root);
1088 } else {
1089 n = PNODE(oldnode->blossom_root);
1090 memo = x;
1091 }
1092
1093 for (; n != memo; n = PNODE(n->blossom_next)) {
1094 e = PEDGE(n->blossom_next_edge);
1095 e->x = 1 - e->x;
1096 if (e->x == 1) {
1097 nodelist[e->nod1].status = MATCHED;
1098 nodelist[e->nod2].status = MATCHED;
1099 nodelist[e->nod1].matched_edge = IEDGE(e);
1100 nodelist[e->nod2].matched_edge = IEDGE(e);
1101 }
1102 }
1103
1104 oldnode->blossom_root = INODE(x);
1105 x->matched_edge = oldnode->matched_edge;
1106 x->status = MATCHED;
1107 }
1108
1109 #ifdef CC_PROTOTYPE_ANSI
fix_tree(graph * G,node * oldnode,node * x,node * y)1110 static void fix_tree (graph *G, node *oldnode, node *x, node *y)
1111 #else
1112 static void fix_tree (G, oldnode, x, y)
1113 graph *G;
1114 node *oldnode, *x, *y;
1115 #endif
1116 {
1117 node *p, *n, *cnode;
1118 int c, label_help;
1119
1120 p = PNODE(oldnode->parent);
1121
1122 #if PRINT_LEVEL
1123 printf (" Fix Tree %i; x = %i; y = %i\n",
1124 (int) (oldnode - PG->nodelist), (int) (x - PG->nodelist),
1125 (int) (y - PG->nodelist));
1126 fflush (stdout);
1127 #endif
1128
1129 /* Restore child-sibling list for p */
1130
1131 if (p->child == INODE(oldnode)) {
1132 y->sibling = oldnode->sibling;
1133 p->child = INODE(y);
1134 } else {
1135 for (c = p->child; c != -1; c = cnode->sibling) {
1136 cnode = PNODE(c);
1137 if (cnode->sibling == INODE(oldnode)) {
1138 cnode->sibling = INODE(y);
1139 y->sibling = oldnode->sibling;
1140 break;
1141 }
1142 }
1143 assert (c != -1);
1144 }
1145
1146 y->parent = INODE(p);
1147 y->parent_edge = oldnode->parent_edge;
1148 G->nodelist[oldnode->child].parent = INODE(x);
1149 x->child = oldnode->child;
1150
1151 label_help = MINUS;
1152 if (y->blossom_next_edge == y->matched_edge) {
1153 /* From y to x */
1154 for (n = y; n != x; n = PNODE(n->blossom_next)) {
1155 n->child = n->blossom_next;
1156 G->nodelist[n->blossom_next].parent = INODE(n);
1157 G->nodelist[n->blossom_next].parent_edge = n->blossom_next_edge;
1158 n->label = label_help;
1159 if (label_help == MINUS) {
1160 label_help = PLUS;
1161 } else {
1162 label_help = MINUS;
1163 }
1164 }
1165 x->label = MINUS; /* This is not set in for-loop */
1166 } else {
1167 for (n = x; n != y; n = PNODE(n->blossom_next)) {
1168 /* From x to y */
1169 n->parent = n->blossom_next;
1170 n->parent_edge = n->blossom_next_edge;
1171 G->nodelist[n->blossom_next].child = INODE(n);
1172 n->label = label_help;
1173 if (label_help == MINUS) {
1174 label_help = PLUS;
1175 } else {
1176 label_help = MINUS;
1177 }
1178 }
1179 y->label = MINUS; /* This is not set in for-loop */
1180 }
1181 }
1182
1183 #ifdef CC_PROTOTYPE_ANSI
expand_blossom(graph * G,node * oldnode,stats * scount)1184 static int expand_blossom (graph *G, node *oldnode, stats *scount)
1185 #else
1186 static int expand_blossom (G, oldnode, scount)
1187 graph *G;
1188 node *oldnode;
1189 stats *scount;
1190 #endif
1191 {
1192 node *memo, *x, *y, *n;
1193 int child, parent;
1194
1195 scount->expand_count++;
1196
1197 #if PRINT_LEVEL
1198 printf (" Expand blossom %i \n", (int) (oldnode - PG->nodelist));
1199 fflush (stdout);
1200 #endif
1201
1202 child = oldnode->child;
1203 parent = oldnode->parent;
1204
1205 n = memo = PNODE(oldnode->blossom_root);
1206 do {
1207 label_penultimate (G, n, INODE(n));
1208 n = PNODE(n->blossom_next);
1209 } while (n != memo);
1210
1211 if (G->edgelist[oldnode->matched_edge].nod1 == INODE(oldnode))
1212 x = PNODE(G->edgelist[oldnode->matched_edge].orig_nod1);
1213 else
1214 x = PNODE(G->edgelist[oldnode->matched_edge].orig_nod2);
1215
1216 if (G->edgelist[oldnode->parent_edge].nod1 == INODE(oldnode))
1217 y = PNODE(G->edgelist[oldnode->parent_edge].orig_nod1);
1218 else
1219 y = PNODE(G->edgelist[oldnode->parent_edge].orig_nod2);
1220
1221 if (lower_edges (G, oldnode)) {
1222 fprintf (stderr, "lower_edges failed\n");
1223 return 1;
1224 }
1225 fix_matching (G, oldnode, PNODE(x->penultimate));
1226
1227 n = memo;
1228 do {
1229 /* Clear tree structure of nodes in oldnodes blossom */
1230 n->blossom_parent = -1;
1231 n->child = -1;
1232 n->parent = -1;
1233 n->sibling = -1;
1234 n->parent_edge = -1;
1235 n->label = NIX;
1236
1237 n = PNODE(n->blossom_next);
1238 n->hit = n->hit | oldnode->hit;
1239 } while (n != memo);
1240
1241 fix_tree (G, oldnode, PNODE(x->penultimate), PNODE(y->penultimate));
1242
1243 /* update tree roots */
1244
1245 for (child = G->nodelist[child].parent; child != parent;
1246 child = G->nodelist[child].parent) {
1247 G->nodelist[child].tree_root = G->nodelist[parent].tree_root;
1248 }
1249
1250 /* clear oldnode */
1251
1252 oldnode->edg_list = -1;
1253 oldnode->pi = 0;
1254 oldnode->status = UNMATCHED;
1255 oldnode->matched_edge = -1;
1256 oldnode->child = -1;
1257 oldnode->sibling = -1;
1258 oldnode->parent = -1;
1259 oldnode->parent_edge = -1;
1260 oldnode->label = NIX;
1261 oldnode->blossom_root = -1;
1262 oldnode->blossom_next = -1;
1263 oldnode->blossom_next_edge = -1;
1264 oldnode->blossom_parent = -1;
1265 oldnode->penultimate = -1;
1266 oldnode->tree_root = -1;
1267 oldnode->hit = 0;
1268
1269 /* put oldnode on unused list */
1270
1271 oldnode->mark = G->unused;
1272 G->unused = oldnode - G->nodelist;
1273
1274 return 0;
1275 }
1276
1277 #ifdef CC_PROTOTYPE_ANSI
common_parent(graph * G,node * p,node * q,int * size)1278 static node *common_parent (graph *G, node *p, node *q, int *size)
1279 #else
1280 static node *common_parent (G, p, q, size)
1281 graph *G;
1282 node *p, *q;
1283 int *size;
1284 #endif
1285 {
1286 int count;
1287 node *n;
1288 node *stop = G->nodelist - 1;
1289
1290 #if PRINT_LEVEL
1291 printf (" Common-parent from %i and %i is ",
1292 (int) (p - PG->nodelist), (int) (q - PG->nodelist));
1293 fflush (stdout);
1294 #endif
1295
1296 for (count = 0, n = p; n != stop; n = PNODE(n->parent)) {
1297 count++;
1298 n->mark = count;
1299 }
1300
1301 for (count = 0; q != stop; q = PNODE(q->parent)) {
1302 if (q->mark) {
1303 *size = (count + q->mark);
1304 for (n = p; n != stop; n = PNODE(n->parent))
1305 n->mark = 0;
1306 return q;
1307 }
1308 count++;
1309 }
1310 for (n = p; n != stop; n = PNODE(n->parent))
1311 n->mark = 0;
1312
1313 *size = 0;
1314 return (node *) NULL;
1315 }
1316
1317 #ifdef CC_PROTOTYPE_ANSI
flip_cycle(graph * G,node * node_i)1318 static void flip_cycle (graph *G, node *node_i)
1319 #else
1320 static void flip_cycle (G, node_i)
1321 graph *G;
1322 node *node_i;
1323 #endif
1324 {
1325 int count = 0, ok = 1;
1326 node *node_j, *node_k;
1327 edge *e, *memo;
1328
1329 #if PRINT_LEVEL
1330 printf (" Flip cycle with root %i: ", (int) (node_i - PG->nodelist));
1331 fflush (stdout);
1332 #endif
1333
1334 /* init start (node_j, node_k, e) */
1335
1336 e = PEDGE(node_i->matched_edge);
1337 node_j = node_i;
1338 node_k = OTHEREND (e, node_i, G->nodelist);
1339
1340 while (ok) {
1341 /* test if last edge */
1342 if (node_k == node_i)
1343 ok = 0;
1344
1345 e->x = (count % 2);
1346 memo = PEDGE(node_k->matched_edge);
1347
1348 if (count % 2) {
1349 node_j->matched_edge = node_k->matched_edge = IEDGE(e);
1350 node_j->status = MATCHED;
1351 node_k->status = MATCHED;
1352 }
1353
1354 e = memo;
1355 node_j = node_k;
1356 node_k = OTHEREND (e, node_j, G->nodelist);
1357 count++;
1358 }
1359 }
1360
1361 #ifdef CC_PROTOTYPE_ANSI
augment_path(graph * G,node * n,int fractional)1362 static void augment_path (graph *G, node *n, int fractional)
1363 #else
1364 static void augment_path (G, n, fractional)
1365 graph *G;
1366 node *n;
1367 int fractional;
1368 #endif
1369 {
1370 edge *edgelist = G->edgelist;
1371 node *nodelist = G->nodelist;
1372
1373 #if PRINT_LEVEL
1374 printf (" Augment path %i", (int) (n - PG->nodelist));
1375 fflush (stdout);
1376 #endif
1377
1378 if (n->parent != -1) {
1379 for (; n->parent_edge != -1; n = PNODE(n->parent)) {
1380 edgelist[n->parent_edge].x = 1 - edgelist[n->parent_edge].x;
1381 if (edgelist[n->parent_edge].x == 1) {
1382 nodelist[edgelist[n->parent_edge].nod1].matched_edge =
1383 n->parent_edge;
1384 nodelist[edgelist[n->parent_edge].nod2].matched_edge =
1385 n->parent_edge;
1386 }
1387 }
1388 }
1389 n->status = MATCHED; /* This will be corrected in HALVES case */
1390 if (!fractional)
1391 G->roots[n->tree_root].status = MATCHED;
1392 clear_tree (G, n); /* n is the root */
1393 }
1394
1395 #ifdef CC_PROTOTYPE_ANSI
augment_blossom(graph * G,edge * e,int fractional,srkstuff * srk)1396 static int augment_blossom (graph *G, edge *e, int fractional, srkstuff *srk)
1397 #else
1398 static int augment_blossom (G, e, fractional, srk)
1399 graph *G;
1400 edge *e;
1401 int fractional;
1402 srkstuff *srk;
1403 #endif
1404 {
1405 node *node_i = PNODE(e->nod1);
1406 node *node_j = PNODE(e->nod2);
1407 node *node_k;
1408 int size;
1409 shrink_ary_T *ary = srk->shrinks.ary;
1410
1411 #if PRINT_LEVEL
1412 printf (" Augment blossom with ends %i & %i \n", e->nod1, e->nod2);
1413 fflush (stdout);
1414 #endif
1415
1416 node_k = common_parent (G, node_i, node_j, &size);
1417 if (node_k == (node *) NULL) { /* more then one tree */
1418 e->x = 1;
1419 node_i->matched_edge = e -G->edgelist;
1420 node_j->matched_edge = e -G->edgelist;
1421 G->unmatched_count -= 2;
1422 augment_path (G, node_i, fractional);
1423 augment_path (G, node_j, fractional);
1424 return 1;
1425 } else {
1426 if (fractional) {
1427 make_cycle_halves (G, node_i, node_j, node_k, e);
1428 augment_path (G, node_k, fractional);
1429 node_k->status = HALVES;
1430 G->unmatched_count--;
1431 return 1;
1432 } else {
1433 if (srk->shrinks.length < srk->shrinks_max) {
1434 ary[srk->shrinks.length].node_i = node_i - G->nodelist;
1435 ary[srk->shrinks.length].node_j = node_j - G->nodelist;
1436 ary[srk->shrinks.length].node_k = node_k - G->nodelist;
1437 ary[srk->shrinks.length].edge = e - G->edgelist;
1438 ary[srk->shrinks.length].size = size;
1439 srk->shrinks.length++;
1440 if (srk->shrinks.length == srk->shrinks_max) {
1441 printf (" WARNING: shrinks.length==shrinks_max=%i\n",
1442 srk->shrinks_max);
1443 fflush (stdout);
1444 }
1445 }
1446 return 0;
1447 }
1448 }
1449 }
1450
1451 #ifdef CC_PROTOTYPE_ANSI
add_to_tree(graph * G,edge * e,node * node_i,node * node_j,int fractional)1452 static int add_to_tree (graph *G, edge *e, node *node_i, node *node_j,
1453 int fractional)
1454 #else
1455 static int add_to_tree (G, e, node_i, node_j, fractional)
1456 graph *G;
1457 edge *e;
1458 node *node_i, *node_j;
1459 int fractional;
1460 #endif
1461 {
1462 node *node_k;
1463 edge *f;
1464
1465 #if PRINT_LEVEL
1466 printf (" Add (%i,%i)=(%i,%i) to tree \n", e->nod1,
1467 e->nod2, (int) (node_i - PG->nodelist), (int) (node_j - PG->nodelist));
1468 fflush (stdout);
1469 #endif
1470
1471 if (node_j->status == UNMATCHED) {
1472 e->x = 1;
1473 node_j->status = MATCHED;
1474 node_i->matched_edge = e - G->edgelist;
1475 node_j->matched_edge = e - G->edgelist;
1476 G->unmatched_count -= 2;
1477 if (!fractional)
1478 G->roots[node_j->tree_root].status = MATCHED;
1479 augment_path (G, node_i, fractional);
1480 return 1;
1481 } else if (node_j->status == HALVES) {
1482 flip_cycle (G, node_j);
1483 e->x = 1;
1484 node_i->status = MATCHED;
1485 node_j->status = MATCHED;
1486 node_i->matched_edge = e - G->edgelist;
1487 node_j->matched_edge = e - G->edgelist;
1488 augment_path (G, node_i, fractional);
1489 G->unmatched_count--;
1490 return 1;
1491 } else {
1492 /* set node_j in tree */
1493 node_j->sibling = node_i->child;
1494 node_j->parent = INODE(node_i);
1495 node_j->tree_root = node_i->tree_root;
1496 node_j->parent_edge = IEDGE(e);
1497 node_j->label = MINUS;
1498
1499 /* set child for node_i */
1500 node_i->child = INODE(node_j);
1501
1502 /* set child for node_j */
1503 f = PEDGE(node_j->matched_edge);
1504 node_k = OTHEREND (f, node_j, G->nodelist);
1505 node_j->child = INODE(node_k);
1506
1507 /* set node_k in tree */
1508 node_k->child = -1;
1509 node_k->sibling = -1;
1510 node_k->parent = INODE(node_j);
1511 node_k->parent_edge = IEDGE(f);
1512 node_k->label = PLUS;
1513 node_k->tree_root = node_j->tree_root;
1514
1515 return 0;
1516 }
1517 }
1518
1519 #ifdef CC_PROTOTYPE_ANSI
checkout_node(graph * G,node * n,int fractional,srkstuff * srk)1520 static int checkout_node (graph *G, node *n, int fractional, srkstuff *srk)
1521 #else
1522 static int checkout_node (G, n, fractional, srk)
1523 graph *G;
1524 node *n;
1525 int fractional;
1526 srkstuff *srk;
1527 #endif
1528 {
1529 node *m;
1530 int augmented;
1531 edge *e;
1532 int ei;
1533
1534 #if PRINT_LEVEL
1535 printf (" Checkout node %i ...\n", (int) (n - PG->nodelist));
1536 fflush (stdout);
1537 #endif
1538
1539 for (ei = n->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
1540 e = PEDGE(ei/2);
1541 if (e->slack == 0) {
1542 m = OTHEREND (e, n, G->nodelist);
1543 if (m->label == PLUS) {
1544 return augment_blossom (G, e, fractional, srk);
1545 }
1546 if (m->label == NIX) {
1547 if ((augmented = add_to_tree (G, e, n, m, fractional)) != 0) {
1548 return augmented;
1549 }
1550 }
1551 }
1552 }
1553
1554 return 0;
1555 }
1556
1557 #ifdef CC_PROTOTYPE_ANSI
grow_tree_no_shrink(graph * G,node * n,int fractional,srkstuff * srk)1558 static int grow_tree_no_shrink (graph *G, node *n, int fractional,
1559 srkstuff *srk)
1560 #else
1561 static int grow_tree_no_shrink (G, n, fractional, srk)
1562 graph *G;
1563 node *n;
1564 int fractional;
1565 srkstuff *srk;
1566 #endif
1567 {
1568 int augmented;
1569 expand_T *expands = &(srk->expands);
1570 int expands_max = srk->expands_max;
1571 node *c=n;
1572 node *stop = G->nodelist - 1;
1573
1574 #if PRINT_LEVEL
1575 printf (" growtree no shrink %i\n", (int) (c - PG->nodelist));
1576 fflush (stdout);
1577 #endif
1578
1579 while (1) {
1580 if (c->label == PLUS) {
1581 if ((augmented = checkout_node (G, c, fractional, srk)) > 0) {
1582 return augmented;
1583 }
1584 } else { /* MINUS node */
1585 if (c->blossom_root != -1 && c->pi == 0) {
1586 if (expands->length < expands_max) {
1587 expands->node[expands->length] = c - G->nodelist;
1588 expands->length++;
1589 if (expands->length == expands_max) {
1590 printf (" WARNING: expands = expands_max = %i\n",
1591 expands_max);
1592 fflush (stdout);
1593 }
1594 }
1595 }
1596 }
1597 /* go to next c */
1598 if (c->child!=-1) {
1599 c=PNODE(c->child);
1600 } else {
1601 while (PNODE(c->sibling) == stop) {
1602 if (c == n) return 0; /* this if is for childless n */
1603 c = PNODE(c->parent);
1604 if (c == n) return 0;
1605 }
1606 c = PNODE(c->sibling);
1607 }
1608 }
1609 }
1610
1611 #ifdef CC_PROTOTYPE_ANSI
grow_tree(graph * G,node * n,int fractional,stats * scount,srkstuff * srk,int * found_aug)1612 static int grow_tree (graph *G, node *n, int fractional, stats *scount,
1613 srkstuff *srk, int *found_aug)
1614 #else
1615 static int grow_tree (G, n, fractional, scount, srk, found_aug)
1616 graph *G;
1617 node *n;
1618 int fractional;
1619 stats *scount;
1620 srkstuff *srk;
1621 int *found_aug;
1622 #endif
1623 {
1624 int i, size;
1625 node *node_i, *node_j, *node_k, *newnode;
1626 edge *e;
1627 shrink_ary_T *ary = srk->shrinks.ary;
1628 expand_T *expands = &(srk->expands);
1629 int *shrinks_ = srk->shrinks_;
1630 shrink_T *shrinks = &(srk->shrinks);
1631 int old_shrinks_len;
1632
1633 *found_aug = 0;
1634
1635 do {
1636 FIND_SURFACE (n);
1637 shrinks->length = 0;
1638 expands->length = 0;
1639
1640 if (grow_tree_no_shrink (G, n, fractional, srk)) {
1641 *found_aug = 1;
1642 return 0;
1643 }
1644
1645 /* found no new match edge, so do shrinks */
1646
1647 while (shrinks->length) {
1648 old_shrinks_len = shrinks->length;
1649 if (shrinks->length > 4) {
1650 /* build buckets for all different sizes */
1651 for (i = 0; i < SHRINKS_SIZE; i++)
1652 shrinks_[i] = -1;
1653 for (i = 0; i < shrinks->length; i++) {
1654 size = ary[i].size;
1655 if (size >= SHRINKS_SIZE)
1656 size = SHRINKS_SIZE - 1;
1657 if (shrinks_[size] == -1) {
1658 shrinks_[size] = i;
1659 ary[i].next = -1;
1660 } else {
1661 ary[i].next = shrinks_[size];
1662 shrinks_[size] = i;
1663 }
1664 }
1665 for (size = SHRINKS_SIZE - 1; size >= 3; size -= 2) {
1666 for (i = shrinks_[size]; i != -1; i = ary[i].next) {
1667 node_i = PNODE(ary[i].node_i);
1668 node_j = PNODE(ary[i].node_j);
1669 node_k = PNODE(ary[i].node_k);
1670 e = PEDGE(ary[i].edge);
1671 FIND_SURFACE (node_i);
1672 FIND_SURFACE (node_j);
1673 FIND_SURFACE (node_k);
1674 if (node_i != node_j) {
1675 newnode = shrink_blossom (G, node_i, node_j,
1676 node_k, e, scount);
1677 if ((checkout_node(G,newnode,fractional,srk)) >0) {
1678 *found_aug = 1;
1679 return 0;
1680 }
1681 }
1682 }
1683 }
1684 } else { /* just do all of the stupid shrinks */
1685 for (i = 0; i < shrinks->length; i++) {
1686 node_i = PNODE(ary[i].node_i);
1687 node_j = PNODE(ary[i].node_j);
1688 node_k = PNODE(ary[i].node_k);
1689 e = PEDGE(ary[i].edge);
1690 FIND_SURFACE (node_i);
1691 FIND_SURFACE (node_j);
1692 FIND_SURFACE (node_k);
1693 if (node_i != node_j) {
1694 newnode = shrink_blossom (G, node_i, node_j,
1695 node_k, e, scount);
1696 if ((checkout_node (G, newnode, fractional, srk)) > 0) {
1697 *found_aug = 1;
1698 return 0;
1699 }
1700 }
1701 }
1702 }
1703 if (shrinks->length > old_shrinks_len) {
1704 size = shrinks->length - old_shrinks_len;
1705 for (i = 0; i < size; i++) {
1706 ary[i].node_i = ary[i + old_shrinks_len].node_i;
1707 ary[i].node_j = ary[i + old_shrinks_len].node_j;
1708 ary[i].node_k = ary[i + old_shrinks_len].node_k;
1709 ary[i].edge = ary[i + old_shrinks_len].edge;
1710 ary[i].size = ary[i + old_shrinks_len].size;
1711 }
1712 shrinks->length = size;
1713 } else {
1714 shrinks->length = 0;
1715 }
1716 }
1717 for (i = 0; i < expands->length; i++) {
1718 /* blossom_root != -1 is necessary because expand
1719 * of this node might have happened here before */
1720 /* blossom_parent == -1 is necessary because
1721 * this node might be shrunk already */
1722 if (G->nodelist[expands->node[i]].blossom_root != -1 &&
1723 G->nodelist[expands->node[i]].blossom_parent == -1) {
1724 if (expand_blossom (G, G->nodelist+expands->node[i], scount)) {
1725 fprintf (stderr, "expand_blossom failed\n");
1726 return 1;
1727 }
1728 }
1729 }
1730 } while (expands->length > 0);
1731 return 0;
1732 }
1733
1734 #ifdef CC_PROTOTYPE_ANSI
apply_dual_change(graph * G,node * n,int delta)1735 static int apply_dual_change (graph *G, node *n, int delta)
1736 #else
1737 static int apply_dual_change (G, n, delta)
1738 graph *G;
1739 node *n;
1740 int delta;
1741 #endif
1742 {
1743 int ei;
1744 edge *e;
1745 node *c=n;
1746 node *stop = G->nodelist - 1;
1747
1748 if (delta == INTEGER_MAX) {
1749 // fprintf (stderr, "\ndelta=Infinity, node=%i\n", (int) (n-G->nodelist));
1750 // fprintf (stderr, "There seems to be no perfect matching\n");
1751 // FAILED_NODE = (int)INODE(n);
1752 // print_node (G, n);
1753 return 1;
1754 }
1755
1756 while (1) {
1757 if (c->label == PLUS) {
1758 c->hit = 1;
1759 c->pi += delta;
1760 for (ei = c->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
1761 e = PEDGE(ei/2);
1762 e->slack -= delta;
1763 }
1764 } else if (c->label == MINUS) {
1765 c->pi -= delta;
1766 for (ei = c->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
1767 e = PEDGE(ei/2);
1768 e->slack += delta;
1769 }
1770 }
1771
1772 /* go to next c */
1773
1774 if (c->child != -1) {
1775 c = PNODE(c->child);
1776 } else {
1777 while (PNODE(c->sibling) == stop) {
1778 if (c == n) return 0; /* this if is for childless n */
1779 c = PNODE(c->parent);
1780 if (c == n) return 0;
1781 }
1782 c = PNODE(c->sibling);
1783 }
1784 }
1785 }
1786
1787 #ifdef CC_PROTOTYPE_ANSI
find_single_dual_change(graph * G,node * n)1788 static int find_single_dual_change (graph *G, node *n)
1789 #else
1790 static int find_single_dual_change (G, n)
1791 graph *G;
1792 node *n;
1793 #endif
1794 {
1795 edge *e;
1796 int ei;
1797 int lab;
1798 int delta;
1799 node *c = n;
1800 node *stop = G->nodelist - 1;
1801
1802 #if PRINT_LEVEL
1803 printf (" %i", (int) (n - PG->nodelist));
1804 fflush (stdout);
1805 #endif
1806
1807 delta = INTEGER_MAX;
1808
1809 while (1) {
1810 if (c->label == PLUS) {
1811 for (ei = c->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
1812 e = PEDGE(ei/2);
1813 if (e->slack < 2 * delta) {
1814 lab = OTHEREND (e, c, G->nodelist)->label;
1815 if (lab == PLUS) {
1816 delta = e->slack / 2;
1817 } else if (lab == NIX) {
1818 if (e->slack < delta) {
1819 delta = e->slack;
1820 }
1821 }
1822 }
1823 }
1824 } else if (c->label == MINUS) {
1825 if (c->blossom_root != -1 && c->pi < delta) {
1826 delta = c->pi;
1827 }
1828 }
1829
1830 /* go to next c */
1831
1832 if (c->child != -1) {
1833 c = PNODE(c->child);
1834 } else {
1835 while (PNODE(c->sibling) == stop) {
1836 if (c == n) return delta ; /* this if is for childless n */
1837 c = PNODE(c->parent);
1838 if (c == n) return delta;
1839 }
1840 c = PNODE(c->sibling);
1841 }
1842 }
1843 }
1844
1845 #ifdef CC_PROTOTYPE_ANSI
find_dual_change_forest(graph * G,node * n)1846 static int find_dual_change_forest (graph *G, node *n)
1847 #else
1848 static int find_dual_change_forest (G, n)
1849 graph *G;
1850 node *n;
1851 #endif
1852 {
1853 edge *e;
1854 int ei;
1855 int delta;
1856 nodeptr *roots = G->roots;
1857 node *c = n, *o;
1858 node *stop = G->nodelist - 1;
1859
1860 #if PRINT_LEVEL
1861 printf (" %i", (int) (c - PG->nodelist));
1862 fflush (stdout);
1863 #endif
1864
1865 delta = INTEGER_MAX;
1866
1867 while (1) {
1868 if (c->label == PLUS) {
1869 for (ei = c->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
1870 e = PEDGE(ei/2);
1871 o = OTHEREND (e, c, G->nodelist);
1872 switch (o->label) {
1873 #ifdef GREEDY_DUAL_CHANGE
1874 case PLUS:
1875 if (roots[c->tree_root].dad != roots[o->tree_root].dad ) {
1876 /* two trees, not connected with vereinigung */
1877 if (roots[o->tree_root].tree_processed) {
1878 if (delta > e->slack -
1879 roots[roots[o->tree_root].dad].delta) {
1880 delta = e->slack -
1881 roots[roots[o->tree_root].dad].delta;
1882 }
1883 } else {
1884 if (e->slack < delta) {
1885 delta = e->slack;
1886 }
1887 }
1888 } else {
1889 if (e->slack < 2 * delta) {
1890 delta = e->slack / 2;
1891 }
1892 }
1893 break;
1894 case NIX:
1895 if (e->slack < delta) {
1896 delta = e->slack;
1897 }
1898 break;
1899 case MINUS:
1900 if (roots[c->tree_root].dad != roots[o->tree_root].dad) {
1901 /* two trees, not connected with vereinigung */
1902 if (roots[o->tree_root].tree_processed) {
1903 if (delta > e->slack +
1904 roots[roots[o->tree_root].dad].delta) {
1905 delta = e->slack +
1906 roots[roots[o->tree_root].dad].delta;
1907 }
1908 } else {
1909 if (e->slack < delta ) {
1910 delta = e->slack;
1911 }
1912 }
1913 }
1914 }
1915 #else /* GREEDY_DUAL_CHANGE */
1916 case PLUS:
1917 if (e->slack < 2 * delta) {
1918 delta = e->slack / 2;
1919 }
1920 break;
1921 case NIX:
1922 if (e->slack < delta) {
1923 delta = e->slack;
1924 }
1925 break;
1926 case MINUS:
1927 if (roots[c->tree_root].dad != roots[o->tree_root].dad) {
1928 if (e->slack < delta ) {
1929 delta = e->slack;
1930 }
1931 }
1932 }
1933 #endif /* GREEDY_DUAL_CHANGE */
1934 }
1935 } else if (c->label == MINUS) {
1936 if (c->blossom_root != -1 && c->pi <= delta) {
1937 delta = c->pi;
1938 }
1939 }
1940
1941 /* go to next c */
1942
1943 if (c->child != -1) {
1944 c = PNODE(c->child);
1945 } else {
1946 while (PNODE(c->sibling) == stop) {
1947 if (c == n) return delta ; /* this is for childless n */
1948 c = PNODE(c->parent);
1949 if (c == n) return delta;
1950 }
1951 c = PNODE(c->sibling);
1952 }
1953 }
1954 }
1955
1956 #ifdef CC_PROTOTYPE_ANSI
vereinigung(graph * G,int x,int y)1957 static void vereinigung (graph *G, int x, int y)
1958 #else
1959 static void vereinigung (G, x, y)
1960 graph *G;
1961 int x, y;
1962 #endif
1963 {
1964 int xroot, yroot, dad;
1965 nodeptr *roots = G->roots;
1966
1967 #if PRINT_LEVEL
1968 printf (" vereinigung (%i,%i)", x, y); fflush (stdout);
1969 #endif
1970
1971 xroot = x;
1972 while (roots[xroot].dad >= 0)
1973 xroot = roots[xroot].dad;
1974 while (roots[x].dad >= 0) {
1975 dad = roots[x].dad;
1976 roots[x].dad = xroot;
1977 x = dad;
1978 }
1979
1980 yroot = y;
1981 while (roots[yroot].dad >= 0)
1982 yroot = roots[yroot].dad;
1983 while (roots[y].dad >= 0) {
1984 dad = roots[y].dad;
1985 roots[y].dad = yroot;
1986 y = dad;
1987 }
1988
1989 if (xroot != yroot) {
1990 if (roots[yroot].dad < roots[xroot].dad) {
1991 roots[yroot].dad += roots[xroot].dad;
1992 roots[xroot].dad = yroot;
1993 } else {
1994 roots[xroot].dad += roots[yroot].dad;
1995 roots[yroot].dad = xroot;
1996 }
1997 }
1998 }
1999
2000 #ifdef CC_PROTOTYPE_ANSI
set_vereinigung(graph * G,node * n)2001 static void set_vereinigung (graph *G, node *n)
2002 #else
2003 static void set_vereinigung (G, n)
2004 graph *G;
2005 node *n;
2006 #endif
2007 {
2008 edge *e;
2009 int ei;
2010 node *nod = G->nodelist;
2011 node *c = n;
2012 node *stop = G->nodelist - 1;
2013
2014 #if PRINT_LEVEL
2015 printf (" set_vereinigung (%i)", (int) (n - PG->nodelist));
2016 fflush (stdout);
2017 #endif
2018
2019 /* this is called only for PLUS nodes */
2020
2021 while (1) {
2022 for (ei = c->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
2023 e = PEDGE(ei/2);
2024 if (nod[e->nod1].tree_root != nod[e->nod2].tree_root) {
2025 if (e->slack == 0) {
2026 if (OTHEREND (e, c, nod)->label == MINUS) {
2027 vereinigung (G, nod[e->nod1].tree_root,
2028 nod[e->nod2].tree_root);
2029 }
2030 }
2031 }
2032 }
2033
2034 /* now check the chilren of c's children (these are PLUS nodes) */
2035
2036 if (c->child != -1) {
2037 c = PNODE(PNODE(c->child)->child);
2038 } else {
2039 if (c == n) return; /* this if is for childless n */
2040 c = PNODE(c->parent);
2041 while (PNODE(c->sibling) == stop) {
2042 c = PNODE(c->parent);
2043 if (c == n) return;
2044 c = PNODE(c->parent);
2045 }
2046 c = PNODE(PNODE(c->sibling)->child);
2047 }
2048 }
2049 }
2050
2051 #ifdef CC_PROTOTYPE_ANSI
find_parity_sum(graph * G,int n)2052 static int find_parity_sum (graph *G, int n)
2053 #else
2054 static int find_parity_sum (G, n)
2055 graph *G;
2056 int n;
2057 #endif
2058 {
2059 node *v;
2060 int sum, ei;
2061 edge *e;
2062
2063 ei = PNODE(n)->edg_list;
2064 e = PEDGE(ei/2);
2065
2066 if (e->nod1 == n) {
2067 v = PNODE(e->orig_nod1);
2068 } else {
2069 v = PNODE(e->orig_nod2);
2070 }
2071
2072 sum = v->pi;
2073 while (v->blossom_parent != -1) {
2074 v = PNODE(v->blossom_parent);
2075 sum += v->pi;
2076 }
2077
2078 return sum;
2079 }
2080
2081 #ifdef CC_PROTOTYPE_ANSI
parity_correction(graph * G,stats * scount)2082 static int parity_correction (graph *G, stats *scount)
2083 #else
2084 static int parity_correction (G, scount)
2085 graph *G;
2086 stats *scount;
2087 #endif
2088 {
2089 nodeptr *np, *npprev = (nodeptr *) NULL;
2090 int hitme = 0;
2091 nodeptr *roots = G->roots;
2092 nodeptr *rstop = G->roots - 1;
2093
2094 #if PRINT_LEVEL
2095 printf (" parity_correction ...\n");
2096 fflush (stdout);
2097 #endif
2098
2099 np = roots + G->unmatched;
2100 while (np->status != UNMATCHED)
2101 np = roots + np->next;
2102 G->unmatched = np - roots;
2103
2104 while (np != rstop) {
2105 if (np->status != UNMATCHED) {
2106 npprev->next = np->next;
2107 } else {
2108 if (np->sum % 2) {
2109 hitme = 1;
2110 scount->dualchange_count++;
2111 if (apply_dual_change (G, PNODE(np->surf), 1)) {
2112 // fprintf (stderr, "apply_dual_change failed\n");
2113 return -1;
2114 }
2115 np->sum += 1;
2116 }
2117 npprev = np;
2118 }
2119 np = roots + np->next;
2120 }
2121 return hitme;
2122 }
2123
2124 #ifdef CC_PROTOTYPE_ANSI
make_dual_change_forest(graph * G,stats * scount)2125 static int make_dual_change_forest (graph *G, stats *scount)
2126 #else
2127 static int make_dual_change_forest (G, scount)
2128 graph *G;
2129 stats *scount;
2130 #endif
2131 {
2132 nodeptr *np;
2133 int t;
2134 int delta;
2135 nodeptr *roots = G->roots;
2136 nodeptr *rstop = G->roots - 1;
2137
2138 #if PRINT_LEVEL
2139 printf ("\n make_dual_change_forest ...\n");
2140 fflush (stdout);
2141 #endif
2142
2143 if (parity_correction (G, scount) == -1) {
2144 fprintf (stderr, "parity_correction failed\n");
2145 return -1;
2146 }
2147
2148 for (np = roots + G->unmatched; np != rstop; np = roots + np->next)
2149 np->dad = -1; /* Set dad for all trees to -1 */
2150
2151 for (np = roots + G->unmatched; np != rstop; np = roots + np->next)
2152 set_vereinigung (G, PNODE(np->surf));
2153
2154 /* Set roots to point to theirselves */
2155 for (np = roots + G->unmatched; np != rstop; np = roots + np->next) {
2156 if (np->dad < 0) {
2157 np->dad = np - roots;
2158 np->delta = INTEGER_MAX;
2159 }
2160 }
2161
2162 /* Set others to point to root */
2163 for (np = roots + G->unmatched; np != rstop; np = roots + np->next) {
2164 t = np - roots;
2165 while (roots[t].dad != t)
2166 t = roots[t].dad;
2167 np->dad = t;
2168 }
2169
2170 #ifdef GREEDY_DUAL_CHANGE
2171 {
2172 int npnext;
2173 int glist = -1;
2174
2175 /* order the unmatched nodes so that each union appears consecutively */
2176
2177 for (np = roots + G->unmatched; np != rstop; np = roots + np->next) {
2178 np->tree_processed = 0;
2179 if (np->dad == np - roots) {
2180 np->gnext = glist;
2181 glist = np - roots;
2182 }
2183 }
2184
2185 for (np = roots + G->unmatched; np != rstop; np = roots + np->next) {
2186 if (np->dad != np - roots) {
2187 np->gnext = roots[np->dad].gnext;
2188 roots[np->dad].gnext = np - roots;
2189 }
2190 }
2191
2192 for (np = roots + G->unmatched; np != rstop; np = roots + npnext) {
2193 npnext = np->next;
2194 np->next = np->gnext;
2195 }
2196 G->unmatched = glist;
2197 }
2198 #endif
2199
2200 /* Set delta for all vereinigung trees */
2201 for (np = roots + G->unmatched; np != rstop; np = roots + np->next) {
2202 delta = find_dual_change_forest (G, PNODE(np->surf));
2203 if (delta < roots[np->dad].delta)
2204 roots[np->dad].delta = delta;
2205 #ifdef GREEDY_DUAL_CHANGE
2206 np->tree_processed = 1;
2207 #endif
2208 }
2209
2210 /* Apply_dual_change for all vereinigung trees */
2211 for (np = roots + G->unmatched; np != rstop; np = roots + np->next) {
2212 delta = roots[np->dad].delta;
2213 scount->dualchange_count++;
2214 if (delta == 0) {
2215 scount->dualzero_count++;
2216 } else {
2217 if (apply_dual_change (G, PNODE(np->surf), delta)) {
2218 // fprintf (stderr, "apply_dual_change failed\n");
2219 return -1;
2220 }
2221 }
2222 np->sum += delta;
2223 }
2224
2225 return 0;
2226 }
2227
2228 #ifdef CC_PROTOTYPE_ANSI
match_main_frac(graph * G,stats * scount,srkstuff * srk)2229 static int match_main_frac (graph *G, stats *scount, srkstuff *srk)
2230 #else
2231 static int match_main_frac (G, scount, srk)
2232 graph *G;
2233 stats *scount;
2234 srkstuff *srk;
2235 #endif
2236 {
2237 node *n;
2238 int i, delta;
2239 int found_aug = 0;
2240
2241 /* Just work with one tree */
2242
2243 for (i = 0, n = G->nodelist; i < G->nnodes; i++, n++) {
2244 if (n->status == UNMATCHED) {
2245 init_tree (G, n);
2246 if (grow_tree (G, n, 1, scount, srk, &found_aug)) {
2247 fprintf (stderr, "grow_tree failed\n");
2248 return 1;
2249 }
2250 while (!found_aug) {
2251 scount->dualchange_count++;
2252 delta = find_single_dual_change (G, n);
2253 if (delta != 0) {
2254 if (apply_dual_change (G, n, delta)) {
2255 // fprintf (stderr, "apply_dual_change failed\n");
2256 return 1;
2257 }
2258 } else {
2259 scount->dualzero_count++;
2260 }
2261 if (grow_tree (G, n, 1, scount, srk, &found_aug)) {
2262 fprintf (stderr, "grow_tree failed\n");
2263 return 1;
2264 }
2265 }
2266 #if PRINT_LEVEL
2267 printf ("."); fflush (stdout);
2268 #endif
2269 }
2270 }
2271 // printf (" %i Dual Changes, %i with delta=0 ",
2272 // scount->dualchange_count, scount->dualzero_count);
2273 // fflush (stdout);
2274 return 0;
2275 }
2276
2277 #ifdef CC_PROTOTYPE_ANSI
match_main_more_in_forest(graph * G,stats * scount,srkstuff * srk)2278 static int match_main_more_in_forest (graph *G, stats *scount, srkstuff *srk)
2279 #else
2280 static int match_main_more_in_forest (G, scount, srk)
2281 graph *G;
2282 stats *scount;
2283 srkstuff *srk;
2284 #endif
2285 {
2286 nodeptr *np;
2287 node *n;
2288 int i, ni, ninext;
2289 int grow_status_forest;
2290 int found_aug = 0;
2291 nodeptr *roots, *rstop;
2292
2293 /* The unmatched surface nodes should be linked using the mark */
2294 /* fields, with G->unmatched pointing to the start of the list. */
2295 /* G->umatched is assumed to be accurate. */
2296
2297 if (G->unmatched_count == 0)
2298 return 0;
2299
2300 /* ********************************* */
2301 /* First init all forests */
2302 /* ********************************* */
2303
2304 G->roots = CC_SAFE_MALLOC (G->unmatched_count, nodeptr);
2305 if (!G->roots) {
2306 // fprintf (stderr, "out of memory in match_main_more_in_forest\n");
2307 return 1;
2308 }
2309 roots = G->roots;
2310 rstop = G->roots - 1;
2311
2312 for (i = 0, ni = G->unmatched; ni != -1; ni = ninext, i++) {
2313 n = PNODE(ni);
2314 ninext = n->mark;
2315 n->mark = 0;
2316 roots[i].status = UNMATCHED;
2317 roots[i].surf = ni;
2318 roots[i].sum = find_parity_sum (G, ni);
2319 roots[i].next = i+1;
2320 init_tree (G, n);
2321 n->tree_root = i;
2322 }
2323 roots[i - 1].next = -1;
2324 G->unmatched = 0;
2325
2326 // printf ("\nTry to grow the trees in the forest directly\n");
2327 // fflush (stdout);
2328
2329 /* ************************************ */
2330 /* Try to grow all forests directly */
2331 /* ************************************ */
2332
2333 for (np = roots + G->unmatched; np != rstop; np = roots + np->next) {
2334 n = PNODE(np->surf);
2335 if (np->status == UNMATCHED) {
2336 if (grow_tree (G, n, 0, scount, srk, &found_aug)) {
2337 // fprintf (stderr, "grow_tree failed\n");
2338 return 1;
2339 }
2340 if (found_aug) {
2341 // printf ("."); fflush (stdout);
2342 if (G->unmatched_count == 0) {
2343 return 0;
2344 }
2345 }
2346 }
2347 }
2348
2349 /* *************************** */
2350 /* Work with tree-vereinigung */
2351 /* *************************** */
2352
2353 // printf ("\nNow work on tree-vereinigung in a forest (%i points)\n",
2354 // G->unmatched_count / 2);
2355 // fflush (stdout);
2356
2357 /* first make sure that all trees have the correct parity */
2358
2359 if (parity_correction (G, scount) == -1) {
2360 // fprintf (stderr, "parity_correction failed\n");
2361 return 1;
2362 }
2363
2364 while (G->unmatched_count > 0) {
2365 if (make_dual_change_forest (G, scount) == -1) {
2366 // fprintf (stderr, "make_dual_change_forest failed\n");
2367 return 1;
2368 }
2369 do {
2370 grow_status_forest = 0;
2371 for (np = roots + G->unmatched; np != rstop; np = roots+np->next) {
2372 n = PNODE(np->surf);
2373 if (np->status == UNMATCHED) {
2374 if (grow_tree (G, n, 0, scount, srk, &found_aug)) {
2375 // fprintf (stderr, "grow_tree failed\n");
2376 return 1;
2377 }
2378 if (found_aug) {
2379 // printf ("."); fflush (stdout);
2380 if (G->unmatched_count == 0) {
2381 goto DONE;
2382 }
2383 grow_status_forest = 1;
2384 }
2385 }
2386 }
2387 } while (grow_status_forest);
2388 }
2389
2390 DONE:
2391
2392 G->unmatched = -1;
2393 CC_IFFREE (G->roots, nodeptr);
2394
2395 // printf ("\n %i Dual Changes, %i with delta=0 ",
2396 // scount->dualchange_count, scount->dualzero_count);
2397 // printf ("| %i Expands ", scount->expand_count);
2398 // printf ("| %i Shrinks ", scount->shrink_count);
2399 // fflush (stdout);
2400
2401 return 0;
2402 }
2403
2404 #ifdef CC_PROTOTYPE_ANSI
match_main_forest(graph * G,stats * scount,srkstuff * srk)2405 static int match_main_forest (graph *G, stats *scount, srkstuff *srk)
2406 #else
2407 static int match_main_forest (G, scount, srk)
2408 graph *G;
2409 stats *scount;
2410 srkstuff *srk;
2411 #endif
2412 {
2413 nodeptr *np;
2414 node *n;
2415 int i, ni, ninext;
2416 int grow_status_forest;
2417 int found_aug = 0;
2418 int delta, delta2;
2419 int do_parity = 1;
2420 nodeptr *roots, *rstop;
2421
2422 // printf ("match_main_forest (%d points)\n", G->unmatched_count / 2);
2423 // fflush (stdout);
2424
2425 if (G->unmatched_count == 0) {
2426 return 0;
2427 }
2428
2429 /* ********************************* */
2430 /* First init all forests */
2431 /* ********************************* */
2432
2433 G->roots = CC_SAFE_MALLOC (G->unmatched_count, nodeptr);
2434 if (!G->roots) {
2435 fprintf (stderr, "out of memory in match_main_more_in_forest\n");
2436 return 1;
2437 }
2438 roots = G->roots;
2439 rstop = G->roots - 1;
2440
2441 for (i = 0, ni = G->unmatched; ni != -1; ni = ninext, i++) {
2442 n = PNODE(ni);
2443 ninext = n->mark;
2444 n->mark = 0;
2445 roots[i].status = UNMATCHED;
2446 roots[i].surf = ni;
2447 roots[i].sum = find_parity_sum (G, ni);
2448 roots[i].next = i+1;
2449 init_tree (G, n);
2450 n->tree_root = i;
2451 }
2452 roots[i - 1].next = -1;
2453 G->unmatched = 0;
2454
2455 /* ********************************* */
2456 /* Work with a forest */
2457 /* ********************************* */
2458
2459 while (G->unmatched_count > 0) {
2460 do {
2461 grow_status_forest = 0;
2462 for (np = roots + G->unmatched; np != rstop; np = roots+np->next) {
2463 n = PNODE(np->surf);
2464 if (np->status == UNMATCHED) {
2465 if (grow_tree (G, n, 0, scount, srk, &found_aug)) {
2466 // fprintf (stderr, "grow_tree failed\n");
2467 return 1;
2468 }
2469 if (found_aug) {
2470 printf ("."); fflush (stdout);
2471 if (G->unmatched_count == 0) {
2472 goto DONE;
2473 }
2474 grow_status_forest = 1;
2475 }
2476 }
2477 }
2478 } while (grow_status_forest);
2479
2480 if (G->unmatched_count > 0 && do_parity) {
2481 if (parity_correction (G, scount) == -1) {
2482 // fprintf (stderr, "parity_correction failed\n");
2483 return 1;
2484 }
2485 do_parity = 0;
2486 }
2487
2488 if (G->unmatched_count > 0 && !do_parity) {
2489 delta = INTEGER_MAX;
2490 scount->dualchange_count++;
2491 for (np = roots + G->unmatched; np != rstop; np = roots+np->next) {
2492 n = PNODE(np->surf);
2493 if (np->status == UNMATCHED) {
2494 delta2 = find_single_dual_change (G, n);
2495 if (delta2 < delta) {
2496 delta = delta2;
2497 }
2498 }
2499 }
2500
2501 if (delta == 0) {
2502 scount->dualzero_count++;
2503 } else {
2504 for (np = roots+G->unmatched; np!=rstop; np = roots+np->next) {
2505 n = PNODE(np->surf);
2506 if (np->status == UNMATCHED) {
2507 if (apply_dual_change (G, n, delta)) {
2508 // fprintf (stderr, "apply_dual_change failed\n");
2509 return 1;
2510 }
2511 np->sum += delta;
2512 }
2513 }
2514 }
2515 }
2516 }
2517
2518 DONE:
2519
2520 G->unmatched = -1;
2521 return 0;
2522 }
2523
2524 #ifdef CC_PROTOTYPE_ANSI
match_main_tree(graph * G,stats * scount,srkstuff * srk)2525 static int match_main_tree (graph *G, stats *scount, srkstuff *srk)
2526 #else
2527 static int match_main_tree (G, scount, srk)
2528 graph *G;
2529 stats *scount;
2530 srkstuff *srk;
2531 #endif
2532 {
2533 nodeptr *np;
2534 node *n;
2535 int i, ni, ninext;
2536 int found_aug = 0;
2537 int delta;
2538 nodeptr *roots, *rstop;
2539
2540 // printf ("match_main_tree (%d points)\n", G->unmatched_count / 2);
2541 // fflush (stdout);
2542
2543 if (G->unmatched_count == 0) {
2544 return 0;
2545 }
2546
2547 /* ********************************* */
2548 /* First init all forests */
2549 /* ********************************* */
2550
2551 G->roots = CC_SAFE_MALLOC (G->unmatched_count, nodeptr);
2552 if (!G->roots) {
2553 // fprintf (stderr, "out of memory in match_main_more_in_forest\n");
2554 return 1;
2555 }
2556 roots = G->roots;
2557 rstop = G->roots - 1;
2558
2559 for (i = 0, ni = G->unmatched; ni != -1; ni = ninext, i++) {
2560 n = PNODE(ni);
2561 ninext = n->mark;
2562 n->mark = 0;
2563 roots[i].status = UNMATCHED;
2564 roots[i].surf = ni;
2565 roots[i].sum = find_parity_sum (G, ni);
2566 roots[i].next = i+1;
2567 /* init_tree (G, n); */
2568 n->tree_root = i;
2569 }
2570 roots[i - 1].next = -1;
2571 G->unmatched = 0;
2572
2573 /* ********************************* */
2574 /* Work with a tree */
2575 /* ********************************* */
2576
2577 while (G->unmatched_count > 0) {
2578 np = roots + G->unmatched;
2579 while (np->status != UNMATCHED) {
2580 np = roots + np->next;
2581 }
2582 G->unmatched = np - roots;
2583 np = roots + G->unmatched;
2584 n = PNODE(np->surf);
2585 init_tree (G, n);
2586 do {
2587 if (grow_tree (G, n, 0, scount, srk, &found_aug)) {
2588 // fprintf (stderr, "grow_tree failed\n");
2589 return 1;
2590 }
2591 if (found_aug) {
2592 printf ("."); fflush (stdout);
2593 if (G->unmatched_count == 0) {
2594 goto DONE;
2595 }
2596 } else {
2597 scount->dualchange_count++;
2598 n = PNODE(np->surf);
2599 delta = find_single_dual_change (G, n);
2600
2601 if (delta == 0) {
2602 scount->dualzero_count++;
2603 } else {
2604 if (apply_dual_change (G, n, delta)) {
2605 // fprintf (stderr, "apply_dual_change failed\n");
2606 return 1;
2607 }
2608 }
2609 }
2610 } while (!found_aug);
2611 }
2612
2613 DONE:
2614
2615 G->unmatched = -1;
2616 return 0;
2617 }
2618
2619 #ifdef CC_PROTOTYPE_ANSI
match_main(graph * G,stats * scount,srkstuff * srk,int use_all)2620 static int match_main (graph *G, stats *scount, srkstuff *srk, int use_all)
2621 #else
2622 static int match_main (G, scount, srk, use_all)
2623 graph *G;
2624 stats *scount;
2625 srkstuff *srk;
2626 int use_all;
2627 #endif
2628 {
2629 double szeit;
2630
2631 szeit = CCutil_zeit ();
2632 // printf (" Need to find %i edges ...\n", G->unmatched_count / 2);
2633 // fflush (stdout);
2634
2635 if (!G->unmatched_count)
2636 return 0;
2637
2638 if (use_all == 2) {
2639 if (match_main_tree (G, scount, srk)) {
2640 // fprintf (stderr, "match_main_tree failed\n");
2641 return 1;
2642 }
2643 } else if (use_all == 1) {
2644 if (match_main_forest (G, scount, srk)) {
2645 // fprintf (stderr, "match_main_forest failed\n");
2646 return 1;
2647 }
2648 } else {
2649 if (match_main_more_in_forest (G, scount, srk)) {
2650 // fprintf (stderr, "match_main_more_in_forest failed\n");
2651 return 1;
2652 }
2653 }
2654 // printf ("\n %i Dual Changes, %i with delta=0 \n",
2655 // scount->dualchange_count, scount->dualzero_count);
2656 // printf (" %i Expands \n", scount->expand_count);
2657 // printf (" %i Shrinks \n", scount->shrink_count);
2658 // printf (" Running Time in match_main: %.2f\n", CCutil_zeit () - szeit);
2659 // fflush (stdout);
2660
2661 return 0;
2662 }
2663
2664 #ifdef CC_PROTOTYPE_ANSI
find_below(graph * G,node * n,int blossom)2665 static node *find_below (graph *G, node *n, int blossom)
2666 #else
2667 static node *find_below (G, n, blossom)
2668 graph *G;
2669 node *n;
2670 int blossom;
2671 #endif
2672 {
2673 if (!n)
2674 return (node *) NULL;
2675
2676 for (; n->blossom_parent != blossom; n = PNODE(n->blossom_parent)) {
2677 if (n->blossom_parent == -1)
2678 return (node *) NULL;
2679 }
2680 return n;
2681 }
2682
2683 #ifdef CC_PROTOTYPE_ANSI
fix_match(graph * G,int blossom)2684 static void fix_match (graph *G, int blossom)
2685 #else
2686 static void fix_match (G, blossom)
2687 graph *G;
2688 int blossom;
2689 #endif
2690 {
2691 node *x, *n;
2692 node *b = blossom + G->nodelist;
2693
2694 #if PRINT_LEVEL
2695 printf (" Fixing blossom %i ", blossom);
2696 fflush (stdout);
2697 #endif
2698
2699 x = find_below (G, PNODE(G->edgelist[b->matched_edge].orig_nod1),
2700 blossom);
2701 if (x == (node *) NULL) {
2702 x = find_below (G, PNODE(G->edgelist[b->matched_edge].orig_nod2),
2703 blossom);
2704 }
2705 /* x is now something like the match_edge->nod(1 or 2)->penultimate */
2706 fix_matching (G, b, x);
2707 b->mark = 1;
2708
2709 n = PNODE(b->blossom_root);
2710 do {
2711 if (n->blossom_root != -1) /* if n is a real blossom */
2712 fix_match (G, INODE(n)); /* call fix match again */
2713 n = PNODE(n->blossom_next);
2714 } while (n != PNODE(b->blossom_root));
2715 }
2716
2717 #ifdef CC_PROTOTYPE_ANSI
adjust_match(graph * G)2718 static void adjust_match (graph *G)
2719 #else
2720 static void adjust_match (G)
2721 graph *G;
2722 #endif
2723 {
2724 node *n, *b;
2725 int i;
2726 int ncount = G->nnodes;
2727 node *nodelist = G->nodelist;
2728
2729 for (i = 0; i < ncount; i++) {
2730 n = PNODE(i);
2731 if (n->blossom_parent != -1 && nodelist[n->blossom_parent].mark == 0) {
2732 for (b = n; b->blossom_parent != -1 &&
2733 nodelist[b->blossom_parent].mark == 0;
2734 b = PNODE(b->blossom_parent));
2735 fix_match (G, INODE(b));
2736 }
2737 }
2738 }
2739
2740 #ifdef CARELESS_REPAIRS
2741 #ifdef CC_PROTOTYPE_ANSI
bring_to_surface(graph * G,node * p,stats * scount,srkstuff * srk)2742 static int bring_to_surface (graph *G, node *p, stats *scount, srkstuff *srk)
2743 #else
2744 static int bring_to_surface (G, p, scount, srk)
2745 graph *G;
2746 node *p;
2747 stats *scount;
2748 srkstuff *srk;
2749 #endif
2750 {
2751 node *psurf, *n, *memo, *x;
2752 edge *e;
2753 int ei;
2754
2755 if (p->blossom_parent == -1)
2756 return 0;
2757
2758 psurf = p;
2759 while (psurf->blossom_parent != -1) {
2760 psurf = PNODE(psurf->blossom_parent);
2761 }
2762
2763 while (p != psurf) {
2764 scount->expand_count++;
2765
2766 /* Mark edge to be unmatched later */
2767 G->edgelist[psurf->matched_edge].mark = 1;
2768
2769 for (ei = psurf->edg_list; ei != -1; ei = e->ptrs[ei % 2]) {
2770 e = PEDGE(ei/2);
2771 e->slack += psurf->pi;
2772 }
2773
2774 n = memo = PNODE(psurf->blossom_root);
2775 do {
2776 label_penultimate (G, n, INODE(n));
2777 n = PNODE(n->blossom_next);
2778 } while (n != memo);
2779
2780 if (G->edgelist[psurf->matched_edge].nod1 == (psurf - G->nodelist))
2781 x = PNODE(G->edgelist[psurf->matched_edge].orig_nod1);
2782 else
2783 x = PNODE(G->edgelist[psurf->matched_edge].orig_nod2);
2784
2785 if (lower_edges (G, psurf)) {
2786 fprintf (stderr, "lower_edges failed\n");
2787 return 1;
2788 }
2789 fix_matching (G, psurf, PNODE(x->penultimate));
2790
2791 n = memo;
2792 do {
2793 /* Clear tree structure of nodes in psurf's blossom */
2794 n->blossom_parent = -1;
2795 n->child = -1;
2796 n->parent = -1;
2797 n->sibling = -1;
2798 n->parent_edge = -1;
2799 n->label = NIX;
2800 n = PNODE(n->blossom_next);
2801 n->hit = n->hit | psurf->hit;
2802 } while (n != memo);
2803
2804 /* clear oldnode */
2805
2806 psurf->edg_list = -1;
2807 psurf->pi = 0;
2808 psurf->status = UNMATCHED;
2809 psurf->matched_edge = -1;
2810 psurf->child = -1;
2811 psurf->sibling = -1;
2812 psurf->parent = -1;
2813 psurf->parent_edge = -1;
2814 psurf->label = NIX;
2815 psurf->blossom_root = -1;
2816 psurf->blossom_next = -1;
2817 psurf->blossom_next_edge = -1;
2818 psurf->blossom_parent = -1;
2819 psurf->penultimate = -1;
2820 psurf->tree_root = -1;
2821 psurf->hit = 0;
2822
2823 /* add psurf to unused list */
2824
2825 psurf->mark = G->unused;
2826 G->unused = psurf - G->nodelist;
2827
2828 psurf = PNODE(p->penultimate);
2829 }
2830
2831 return 0;
2832 }
2833 #endif /* CARELESS_REPAIRS */
2834
2835 #ifdef BALL_DERIGS_REPAIRS
2836 #ifdef CC_PROTOTYPE_ANSI
bring_to_surface(graph * G,node * p,stats * scount,srkstuff * srk)2837 static int bring_to_surface (graph *G, node *p, stats *scount, srkstuff *srk)
2838 #else
2839 static int bring_to_surface (G, p, scount, srk)
2840 graph *G;
2841 node *p;
2842 stats *scount;
2843 srkstuff *srk;
2844 #endif
2845 {
2846 node *n, *newnode;
2847 node *stop = G->nodelist - 1;
2848 edge *e, *f;
2849 int sum, delta;
2850 int found_aug = 0;
2851
2852 printf ("bring to surface (%d) ...\n", p - G->nodelist);
2853 fflush (stdout);
2854
2855 if (p->blossom_parent == -1) {
2856 return 0;
2857 }
2858
2859 /* Take newnode from unused list */
2860
2861 newnode = PNODE(G->unused);
2862 G->unused = newnode->mark;
2863 newnode->mark = 0;
2864
2865 /* Initialize newnode */
2866
2867 sum = 0;
2868 for (n = p; n != stop; n = PNODE(n->blossom_parent))
2869 sum += n->pi;
2870
2871 newnode->matched_edge = -1;
2872 newnode->status = UNMATCHED;
2873 newnode->pi = sum % 2;
2874
2875 /* Add the edge (newnode, p) */
2876
2877 if (G->nedges >= G->max_nedges) {
2878 fprintf (stderr, "number of edges exceeds max_nedges\n");
2879 return 1;
2880 }
2881 e = PEDGE(G->nedges);
2882 e->slack = 0;
2883 e->x = 0;
2884 e->orig_nod1 = e->nod1 = newnode - G->nodelist;
2885 e->orig_nod2 = p - G->nodelist;
2886 n = p;
2887 while (n->blossom_parent != -1)
2888 n = PNODE(n->blossom_parent);
2889 e->nod2 = INODE(n);
2890
2891 /* Attach the new edge to the edg_lists of newnode and n */
2892
2893 e->ptrs[0] = n->edg_list;
2894 n->edg_list = 2*G->nedges;
2895 e->ptrs[1] = -1;
2896 newnode->edg_list = 2*G->nedges + 1;
2897
2898 G->nedges++;
2899
2900 init_tree (G, newnode);
2901 newnode->tree_root = INODE(newnode);
2902
2903 if (grow_tree (G, newnode, 0, scount, srk, &found_aug)) {
2904 fprintf (stderr, "grow_tree failed\n");
2905 return 1;
2906 }
2907 if (found_aug) {
2908 fprintf (stderr, "found an augmentation in bring_to_surface\n");
2909 return 1;
2910 }
2911
2912 while (p->blossom_parent != -1) {
2913 delta = find_single_dual_change (G, newnode);
2914 fflush (stdout);
2915 if (delta == 0) {
2916 printf ("O"); fflush (stdout);
2917 }
2918 if (apply_dual_change (G, newnode, delta)) {
2919 // fprintf (stderr, "apply_dual_change failed\n");
2920 return 1;
2921 }
2922 if (grow_tree (G, newnode, 0, scount, srk, &found_aug)) {
2923 fprintf (stderr, "grow_tree failed\n");
2924 return 1;
2925 }
2926 if (found_aug) {
2927 fprintf (stderr, "found an augmentation in bring_to_surface\n");
2928 return 1;
2929 }
2930 }
2931
2932 clear_tree (G, newnode);
2933
2934 /* reset newnode */
2935
2936 newnode->edg_list = -1;
2937 newnode->pi = 0;
2938 newnode->status = UNMATCHED;
2939 newnode->matched_edge = -1;
2940 newnode->child = -1;
2941 newnode->sibling = -1;
2942 newnode->parent = -1;
2943 newnode->parent_edge = -1;
2944 newnode->label = NIX;
2945 newnode->blossom_root = -1;
2946 newnode->blossom_next = -1;
2947 newnode->blossom_next_edge = -1;
2948 newnode->blossom_parent = -1;
2949 newnode->penultimate = -1;
2950 newnode->tree_root = -1;
2951 newnode->hit = 0;
2952
2953 /* Delete the added edge */
2954
2955 f = PEDGE (p->edg_list/2);
2956 if (f == e) {
2957 p->edg_list = e->ptrs[p->edg_list % 2];
2958 } else {
2959 int ei, eiprev = p->edg_list;
2960
2961 for (ei = f->ptrs[p->edg_list % 2]; ei != -1; ei = f->ptrs[ei % 2]) {
2962 f = PEDGE(ei/2);
2963 if (f == e) {
2964 G->edgelist[eiprev/2].ptrs[eiprev % 2] = f->ptrs[ei % 2];
2965 break;
2966 }
2967 eiprev = ei;
2968 }
2969 if (ei == -1) {
2970 fprintf (stderr, "Screw up in bring_to_surface\n");
2971 return 1;
2972 }
2973 }
2974 G->nedges--;
2975
2976 /* Return newnode from unused list */
2977
2978 newnode->mark = G->unused;
2979 G->unused = newnode - G->nodelist;
2980
2981 return 0;
2982 }
2983 #endif /* BALL_DERIGS_REPAIRS */
2984
2985 #ifdef CC_PROTOTYPE_ANSI
add_repair_edge(graph * G,int n1,int n2,int len,stats * scount,srkstuff * srk)2986 static int add_repair_edge (graph *G, int n1, int n2, int len,
2987 stats *scount, srkstuff *srk)
2988 #else
2989 static int add_repair_edge (G, n1, n2, len, scount, srk)
2990 graph *G;
2991 int n1, n2, len;
2992 stats *scount;
2993 srkstuff *srk;
2994 #endif
2995 {
2996 edge *e, *f = (edge *) NULL;
2997 node *n1surf, *n2surf;
2998 int ei, eiprev = -1, sum2;
2999
3000 n1surf = PNODE(n1);
3001 n2surf = PNODE(n2);
3002 sum2 = n2surf->pi;
3003 while (n2surf->blossom_parent != -1) {
3004 n2surf = PNODE(n2surf->blossom_parent);
3005 sum2 += n2surf->pi;
3006 }
3007
3008 if (G->nedges >= G->max_nedges) {
3009 fprintf (stderr, "number of edges exceeds max_nedges\n");
3010 return 1;
3011 }
3012 e = PEDGE(G->nedges);
3013 e->slack = len - n1surf->pi - sum2;
3014 e->mark = 0;
3015 e->x = 0;
3016 e->orig_nod1 = n1;
3017 e->orig_nod2 = n2;
3018 e->nod1 = n1surf - G->nodelist;
3019 e->nod2 = n2surf - G->nodelist;
3020
3021 for (ei = n1surf->edg_list; ei != -1; ei = f->ptrs[ei % 2]) {
3022 f = PEDGE(ei/2);
3023 eiprev = ei;
3024 }
3025 f->ptrs[eiprev % 2] = 2*G->nedges;
3026 e->ptrs[0] = -1;
3027
3028 for (ei = n2surf->edg_list; ei != -1; ei = f->ptrs[ei % 2]) {
3029 f = PEDGE(ei/2);
3030 eiprev = ei;
3031 }
3032 f->ptrs[eiprev % 2] = 2*G->nedges + 1;
3033 e->ptrs[1] = -1;
3034
3035 G->nedges++;
3036
3037 if (e->slack < 0) {
3038 n1surf->pi += e->slack;
3039 for (ei = n1surf->edg_list; ei != -1; ei = f->ptrs[ei % 2]) {
3040 f = PEDGE(ei/2);
3041 f->slack -= e->slack;
3042 }
3043 }
3044
3045 #ifdef BALL_DERIGS_REPAIRS
3046 if (n1surf->matched_edge != -1) {
3047 unmatch_edge (G, PEDGE(n1surf->matched_edge));
3048 if (match_main_more_in_forest (G, scount, srk)) {
3049 /* if (match_main_forest (G, scount, srk)) { */
3050 fprintf (stderr, "match_main_more_in_forest failed\n");
3051 return 1;
3052 }
3053 }
3054 #endif
3055 #ifdef CARELESS_REPAIRS
3056 if (n1surf->matched_edge != -1) {
3057 /* mark edges to be unmatched later */
3058 G->edgelist[n1surf->matched_edge].mark = 1;
3059 }
3060 #endif
3061
3062 return 0;
3063 }
3064
3065 #ifdef CC_PROTOTYPE_ANSI
unmatch_edge(graph * G,edge * e)3066 static void unmatch_edge (graph *G, edge *e)
3067 #else
3068 static void unmatch_edge (G, e)
3069 graph *G;
3070 edge *e;
3071 #endif
3072 {
3073 G->nodelist[e->nod1].status = UNMATCHED;
3074 G->nodelist[e->nod2].status = UNMATCHED;
3075 e->x = 0;
3076 G->nodelist[e->nod1].matched_edge = -1;
3077 G->nodelist[e->nod2].matched_edge = -1;
3078 G->unmatched_count += 2;
3079
3080 G->nodelist[e->nod1].mark = G->unmatched;
3081 G->nodelist[e->nod2].mark = e->nod1;
3082 G->unmatched = e->nod2;
3083 }
3084
3085 #ifdef CC_PROTOTYPE_ANSI
run_repair(graph * G,int badcount,int * badlist,int * badlen,stats * scount,srkstuff * srk)3086 static int run_repair (graph *G, int badcount, int *badlist, int *badlen,
3087 stats *scount, srkstuff *srk)
3088 #else
3089 static int run_repair (G, badcount, badlist, badlen, scount, srk)
3090 graph *G;
3091 int badcount;
3092 int *badlist;
3093 int *badlen;
3094 stats *scount;
3095 srkstuff *srk;
3096 #endif
3097 {
3098 int i, sum1, sum2;
3099 int n1, n2;
3100 node *n;
3101 edge *e, *lastedge;
3102
3103 for (i = 0; i < badcount; i++) {
3104 #ifdef PRINT_REPAIRS
3105 printf ("Bad Edge [%d, %d]\n", badlist[2*i], badlist[2*i + 1]);
3106 fflush (stdout);
3107 #endif
3108 sum1 = 0;
3109 n = PNODE(badlist[2*i]);
3110 while (n->blossom_parent != -1) {
3111 n = PNODE(n->blossom_parent);
3112 sum1 += n->pi;
3113 }
3114 sum2 = 0;
3115 n = PNODE(badlist[2*i + 1]);
3116 while (n->blossom_parent != -1) {
3117 n = PNODE(n->blossom_parent);
3118 sum2 += n->pi;
3119 }
3120 if (sum1 < sum2) {
3121 n1 = badlist[2*i];
3122 n2 = badlist[2*i + 1];
3123 } else {
3124 n1 = badlist[2*i + 1];
3125 n2 = badlist[2*i];
3126 }
3127 if (bring_to_surface (G, PNODE(n1), scount, srk)){
3128 fprintf (stderr, "bring_to_surface failed\n");
3129 return 1;
3130 }
3131 if (add_repair_edge (G, n1, n2, 2 * badlen[i], scount, srk)) {
3132 fprintf (stderr, "add_repair_edge failed\n");
3133 return 1;
3134 }
3135 printf ("+"); fflush (stdout);
3136 }
3137 printf ("\n");
3138
3139 lastedge = PEDGE(G->nedges);
3140 for (e = G->edgelist; e != lastedge; e++) {
3141 if (e->mark) {
3142 e->mark = 0;
3143 if (e->x == 1) {
3144 unmatch_edge (G, e);
3145 }
3146 }
3147 }
3148
3149 if (G->unmatched_count > 0) {
3150 /*
3151 if (match_main_forest (G, scount, srk)) {
3152 fprintf (stderr, "match_main_forest failed\n");
3153 return 1;
3154 }
3155 */
3156 if (match_main_more_in_forest (G, scount, srk)) {
3157 fprintf (stderr, "match_main_more_in_forest failed\n");
3158 return 1;
3159 }
3160 }
3161
3162 return 0;
3163 }
3164
3165 #ifdef CC_PROTOTYPE_ANSI
add_edges(graph * G,CCdatagroup * dat,int badcount,int * badlist,int * badlen)3166 static int add_edges (graph *G, CCdatagroup *dat, int badcount, int *badlist,
3167 int *badlen)
3168 #else
3169 static int add_edges (G, dat, badcount, badlist, badlen)
3170 graph *G;
3171 CCdatagroup *dat;
3172 int badcount;
3173 int *badlist;
3174 int *badlen;
3175 #endif
3176 {
3177 int i,j;
3178 int ncount = G->nnodes;
3179 int ecount = G->nedges;
3180 edge *e;
3181 edge *edgelist = G->edgelist;
3182 node *nodelist = G->nodelist;
3183 int maxbad = G->nnodes * MAX_BAD_PORTION;
3184
3185 /* init unused */
3186 j = 3*ncount/2;
3187 for (i = ncount; i < j; i++) {
3188 nodelist[i].mark = i + 1;
3189 }
3190 nodelist[j].mark = -1;
3191 G->unused = ncount;
3192
3193 /* Init edg_list */
3194 for (i = 0; i < 3*ncount/2; i++) {
3195 nodelist[i].edg_list = -1;
3196 }
3197
3198 /* correct old edges */
3199 for (i = 0, e = edgelist; i < ecount; i++, e++) {
3200 e->nod1 = e->orig_nod1;
3201 e->nod2 = e->orig_nod2;
3202 e->slack = 2 * CCutil_dat_edgelen (e->nod1, e->nod2, dat);
3203 }
3204
3205 if (badcount <= maxbad) {
3206 for (i = 0; i < badcount; i++) {
3207 if (ecount + 1 < G->max_nedges) {
3208 edgelist[ecount].orig_nod1 =
3209 edgelist[ecount].nod1 = badlist[2*i];
3210 edgelist[ecount].orig_nod2 =
3211 edgelist[ecount].nod2 = badlist[2*i+1];
3212 edgelist[ecount].slack = 2 * badlen[i];
3213 ecount++;
3214 } else {
3215 fprintf(stderr, "Exceeded limit on the number of edges\n");
3216 G->nedges = ecount;
3217 return 1;
3218 }
3219 }
3220 } else {
3221 printf ("Using only %.2f of the bad edges\n",
3222 ((double) maxbad) / ((double) badcount));
3223 fflush (stdout);
3224 for (i = 0; i < badcount; i++) {
3225 if (CCutil_lprand () % badcount < maxbad) {
3226 if (ecount + 1 < G->max_nedges) {
3227 edgelist[ecount].orig_nod1 =
3228 edgelist[ecount].nod1 = badlist[2*i];
3229 edgelist[ecount].orig_nod2 =
3230 edgelist[ecount].nod2 = badlist[2*i+1];
3231 edgelist[ecount].slack = 2 * badlen[i];
3232 ecount++;
3233 } else {
3234 fprintf(stderr, "Exceeded limit on the number of edges\n");
3235 G->nedges = ecount;
3236 return 1;
3237 }
3238 }
3239 }
3240 }
3241
3242 printf(" Added %i edges ", ecount - G->nedges);
3243 fflush (stdout);
3244
3245 G->nedges = ecount;
3246
3247 /* build the edg_lists */
3248
3249 for (i = 0, e = edgelist; i < ecount; i++, e++) {
3250 e->ptrs[0] = nodelist[e->nod1].edg_list;
3251 nodelist[e->nod1].edg_list = 2*i;
3252 e->ptrs[1] = nodelist[e->nod2].edg_list;
3253 nodelist[e->nod2].edg_list = 2*i + 1;
3254 }
3255 for (i = 0; i < ncount; i++) {
3256 G->nodelist[i].mark = 0;
3257 }
3258
3259 return 0;
3260 }
3261
3262 #ifdef CC_PROTOTYPE_ANSI
price_repair(graph * G,int * finished,CCdatagroup * dat,stats * scount,srkstuff * srk,int partialprice)3263 static int price_repair (graph *G, int *finished, CCdatagroup *dat, stats *scount,
3264 srkstuff *srk, int partialprice)
3265 #else
3266 static int price_repair (G, finished, dat, scount, srk, partialprice)
3267 graph *G;
3268 int *finished;
3269 CCdatagroup *dat;
3270 stats *scount;
3271 srkstuff *srk;
3272 int partialprice;
3273 #endif
3274 {
3275 int badcount = 0;
3276 int *badlist = (int *) NULL;
3277 int *badlen = (int *) NULL;
3278 int i;
3279 CCtsp_edgegenerator eg;
3280 CCtsp_edgegenerator *myeg = (CCtsp_edgegenerator *) NULL;
3281
3282 *finished = 0;
3283
3284 if (partialprice > 0) {
3285 printf ("Price using the nearest %d neighbor graph\n", partialprice);
3286 fflush (stdout);
3287 if (CCtsp_init_edgegenerator (&eg, G->nnodes, dat, (CCtsp_genadj *) NULL,
3288 partialprice)) {
3289 fprintf (stderr, "init_edgegenerator failed\n");
3290 return 1;
3291 }
3292 myeg = ⪚
3293 }
3294
3295 if (pricing (G, dat, scount, srk, &badcount, &badlist, &badlen, myeg, 0)) {
3296 fprintf (stderr, "pricing failed\n");
3297 if (myeg)
3298 CCtsp_free_edgegenerator (myeg);
3299 return 1;
3300 }
3301
3302 while (badcount > 0) {
3303 if (badcount > G->nnodes / SWITCH_LEVEL /* 10000 */) {
3304 if (myeg)
3305 CCtsp_free_edgegenerator (myeg);
3306 if (add_edges (G, dat, badcount, badlist, badlen)) {
3307 fprintf (stderr, "add_edges failed\n");
3308 CC_IFFREE (badlist, int);
3309 CC_IFFREE (badlen, int);
3310 return 1;
3311 } else {
3312 CC_IFFREE (badlist, int);
3313 CC_IFFREE (badlen, int);
3314 return 0;
3315 }
3316 } else {
3317 for (i = 0; i < 3 * G->nnodes / 2; i++)
3318 G->nodelist[i].hit = 0;
3319 if (run_repair (G, badcount, badlist, badlen, scount, srk)) {
3320 fprintf (stderr, "run_repair failed\n");
3321 CC_IFFREE (badlist, int);
3322 CC_IFFREE (badlen, int);
3323 if (myeg)
3324 CCtsp_free_edgegenerator (myeg);
3325 return 1;
3326 }
3327 }
3328 CC_IFFREE (badlist, int);
3329 CC_IFFREE (badlen, int);
3330
3331 if (pricing(G,dat,scount,srk, &badcount, &badlist, &badlen, myeg, 1)) {
3332 fprintf (stderr, "pricing failed\n");
3333 if (myeg)
3334 CCtsp_free_edgegenerator (myeg);
3335 return 1;
3336 }
3337 }
3338
3339 *finished = 1;
3340 if (myeg)
3341 CCtsp_free_edgegenerator (myeg);
3342
3343 return 0;
3344 }
3345
3346 #ifdef CC_PROTOTYPE_ANSI
pricing(graph * G,CCdatagroup * dat,stats * scount,srkstuff * srk,int * badcount,int ** badlist,int ** badlen,CCtsp_edgegenerator * eg,int usehit)3347 static int pricing (graph *G, CCdatagroup *dat, stats *scount, srkstuff *srk,
3348 int *badcount, int **badlist, int **badlen,
3349 CCtsp_edgegenerator *eg, int usehit)
3350 #else
3351 static int pricing (G, dat, scount, srk, badcount, badlist, badlen, eg, usehit)
3352 graph *G;
3353 CCdatagroup *dat;
3354 stats *scount;
3355 srkstuff *srk;
3356 int *badcount;
3357 int **badlist, **badlen;
3358 CCtsp_edgegenerator *eg;
3359 int usehit;
3360 #endif
3361 {
3362 int i, k;
3363 double szeit, penalty;
3364 double *orig_pi = (double *) NULL;
3365 int *orig_parent = (int *) NULL;
3366 int rval = 0;
3367 int ncount = G->nnodes;
3368 node *nodelist = G->nodelist;
3369 char *hit = (char *) NULL;
3370
3371 szeit = CCutil_zeit();
3372 printf ("Pricing ...\n");
3373 fflush (stdout);
3374
3375 orig_pi = CC_SAFE_MALLOC (ncount*3/2+1, double);
3376 orig_parent = CC_SAFE_MALLOC (ncount*3/2+1, int);
3377 if (!orig_pi || !orig_parent) {
3378 fprintf (stderr, "out of memory in pricing\n");
3379 rval = 1; goto CLEANUP;
3380 }
3381
3382 if (usehit) {
3383 hit = CC_SAFE_MALLOC (ncount*3/2 + 1, char);
3384 if (!hit) {
3385 fprintf (stderr, "out of memory in pricing\n");
3386 rval = 1; goto CLEANUP;
3387 }
3388 }
3389
3390 for (i = 0; i < 3*ncount/2; i++) {
3391 k = nodelist[i].blossom_parent;
3392 orig_parent[i] = k;
3393 orig_pi[i] = ((double) nodelist[i].pi) / 2.0;
3394 if (hit) {
3395 hit[i] = nodelist[i].hit;
3396 }
3397 }
3398
3399 if (matching_price (ncount, dat, orig_pi, orig_parent, badcount,
3400 badlist, badlen, &penalty, eg, hit) ) {
3401 fprintf (stderr, "matching_price failed\n");
3402 rval = 1; goto CLEANUP;
3403 }
3404
3405 printf("\n%i Edges are wrong in the extra graph, penalty=%.4f\n",
3406 *badcount, penalty);
3407
3408 printf (" Total time in pricing: %.2f seconds\n", CCutil_zeit () - szeit);
3409 fflush(stdout);
3410
3411 CLEANUP:
3412
3413 CC_IFFREE (orig_pi, double);
3414 CC_IFFREE (orig_parent, int);
3415 CC_IFFREE (hit, char);
3416
3417 return rval;
3418 }
3419
3420 #ifdef CC_PROTOTYPE_ANSI
test_matching(graph * G,CCdatagroup * dat,int * elen,int fractional,int * bad)3421 static int test_matching (graph *G, CCdatagroup *dat, int *elen, int fractional,
3422 int *bad)
3423 #else
3424 static int test_matching (G, dat, elen, fractional, bad)
3425 graph *G;
3426 CCdatagroup *dat;
3427 int *elen;
3428 int fractional;
3429 int *bad;
3430 #endif
3431 {
3432 int i, c, len;
3433 double a = 0.0, b = 0.0;
3434 edge *e;
3435 int ncount = G->nnodes;
3436 node *nodelist = G->nodelist;
3437
3438 *bad = 0;
3439
3440 if (dat) {
3441 if (fractional) {
3442 for (i = 0; i < G->nedges; i++) {
3443 e = G->edgelist + i;
3444 if (e->x == 1) {
3445 a += (double) CCutil_dat_edgelen (e->orig_nod1, e->orig_nod2, dat);
3446 a += (double) CCutil_dat_edgelen (e->orig_nod1, e->orig_nod2, dat);
3447 } else if (e->x == HALF) {
3448 a += (double) CCutil_dat_edgelen (e->orig_nod1, e->orig_nod2, dat);
3449 }
3450 }
3451 for (i = 0; i < ncount; i++)
3452 b += (double) nodelist[i].pi;
3453 } else {
3454 for (i = 0; i < ncount; i++) {
3455 e = nodelist[i].matched_edge + G->edgelist;
3456 a += (double) CCutil_dat_edgelen (e->orig_nod1, e->orig_nod2, dat);
3457 }
3458 for (i = 0; i <= 3*ncount/2; i++)
3459 b += (double) nodelist[i].pi;
3460 }
3461 } else {
3462 if (fractional) {
3463 for (i = 0; i < ncount; i++) {
3464 a += (double) elen[nodelist[i].matched_edge];
3465 b += (double) nodelist[i].pi;
3466 }
3467 } else {
3468 for (i = 0; i < ncount; i++) {
3469 a += (double) elen[nodelist[i].matched_edge];
3470 }
3471 for (i = 0; i <= 3*ncount/2; i++)
3472 b += (double) nodelist[i].pi;
3473 }
3474 }
3475 b /= 2.0;
3476 a /= 2.0;
3477 printf(" Matching Length: %.1f, Dual %.1f\n", a, b);
3478 fflush (stdout);
3479
3480 if (a != b) {
3481 printf ("ERROR: the primal and dual values do not agree\n");
3482 fflush (stdout);
3483 *bad = 1;
3484 return 0;
3485 }
3486
3487 if (fractional) {
3488 for (i = 0, e = G->edgelist; i < G->nedges; i++, e++) {
3489 if (dat) {
3490 len = 2 * CCutil_dat_edgelen (e->orig_nod1, e->orig_nod2, dat);
3491 } else {
3492 len = 2 * elen[i];
3493 }
3494 if ((c = len - nodelist[e->nod1].pi-nodelist[e->nod2].pi) < 0) {
3495 printf("ERROR: Edge %i %i has slack < 0\n", e->nod1, e->nod2);
3496 fflush (stdout);
3497 *bad = 1;
3498 return 0;
3499 }
3500 if (e->x != 0 && c > 0) {
3501 printf("ERROR: Edge %i %i has x != 0 and slack > 0\n",
3502 e->nod1, e->nod2);
3503 fflush (stdout);
3504 *bad = 1;
3505 return 0;
3506 }
3507 }
3508 } else {
3509 int *mlist = (int *) NULL;
3510 int *mlen = (int *) NULL;
3511 double *pi = (double *) NULL;
3512 int *parent = (int *) NULL;
3513 int csbad, k;
3514 int *badlist = (int *) NULL;
3515 int *badlen = (int *) NULL;
3516 int badcount = 0;
3517 double penalty;
3518
3519 mlist = CC_SAFE_MALLOC (ncount, int);
3520 mlen = CC_SAFE_MALLOC (ncount/2, int);
3521 pi = CC_SAFE_MALLOC (ncount*3/2+1, double);
3522 parent = CC_SAFE_MALLOC (ncount*3/2+1, int);
3523 if (!mlist || !mlen || !pi || !parent) {
3524 fprintf (stderr, "out of memory in test_matching\n");
3525 CC_IFFREE (mlist, int);
3526 CC_IFFREE (mlen, int);
3527 CC_IFFREE (pi, double);
3528 CC_IFFREE (parent, int);
3529 return 1;
3530 }
3531
3532 for (i = 0, e = G->edgelist, k = 0; i < G->nedges; i++, e++) {
3533 if (e->x == 1) {
3534 mlist[2*k] = e->orig_nod1;
3535 mlist[2*k + 1] = e->orig_nod2;
3536 if (dat) {
3537 mlen[k] = CCutil_dat_edgelen (e->orig_nod1, e->orig_nod2, dat);
3538 } else {
3539 mlen[k] = elen[i];
3540 }
3541 k++;
3542 }
3543 }
3544 for (i = 0; i < 3*ncount/2; i++) {
3545 parent[i] = nodelist[i].blossom_parent;
3546 pi[i] = ((double) nodelist[i].pi) / 2.0;
3547 }
3548
3549 if (matching_check (ncount, pi, parent, mlist, mlen, &csbad)) {
3550 fprintf (stderr, "matching_check failed\n");
3551 CC_FREE (mlist, int);
3552 CC_FREE (mlen, int);
3553 CC_FREE (pi, double);
3554 CC_FREE (parent, int);
3555 return 1;
3556 }
3557
3558 if (csbad) {
3559 printf ("matching_check found an error\n");
3560 fflush (stdout);
3561 *bad = 1;
3562 CC_FREE (mlist, int);
3563 CC_FREE (mlen, int);
3564 CC_FREE (pi, double);
3565 CC_FREE (parent, int);
3566 return 0;
3567 }
3568
3569 if (matching_price (ncount, dat, pi, parent, &badcount,
3570 &badlist, &badlen, &penalty, (CCtsp_edgegenerator *) NULL,
3571 (char *) NULL)) {
3572 fprintf (stderr, "matching_price failed\n");
3573 CC_FREE (mlist, int);
3574 CC_FREE (mlen, int);
3575 CC_FREE (pi, double);
3576 CC_FREE (parent, int);
3577 return 1;
3578 }
3579
3580 if (badcount) {
3581 printf ("pricing on complete graph found %d bad edges\n", badcount);
3582 fflush (stdout);
3583 *bad = 1;
3584 CC_IFFREE (badlist, int);
3585 CC_IFFREE (badlen, int);
3586 CC_FREE (mlist, int);
3587 CC_FREE (mlen, int);
3588 CC_FREE (pi, double);
3589 CC_FREE (parent, int);
3590 return 0;
3591 }
3592
3593 CC_FREE (mlist, int);
3594 CC_FREE (mlen, int);
3595 CC_FREE (pi, double);
3596 CC_FREE (parent, int);
3597 }
3598
3599 return 0;
3600 }
3601
3602 #ifdef CC_PROTOTYPE_ANSI
dual_match_len(graph * G,int fractional,double * val)3603 static void dual_match_len (graph *G, int fractional, double *val)
3604 #else
3605 static void dual_match_len (G, fractional, val)
3606 graph *G;
3607 int fractional;
3608 double *val;
3609 #endif
3610 {
3611 int i;
3612 double b = 0.0;
3613 int ncount = G->nnodes;
3614 node *nodelist = G->nodelist;
3615
3616 if (fractional) {
3617 for (i = 0; i < ncount; i++) {
3618 b += (double) nodelist[i].pi;
3619 }
3620 } else {
3621 for (i = 0; i <= 3*ncount/2; i++)
3622 b += (double) nodelist[i].pi;
3623 }
3624 *val = b / 2.0;
3625 }
3626
3627 #ifdef CC_PROTOTYPE_ANSI
make_match(graph * G)3628 static int make_match (graph *G)
3629 #else
3630 static int make_match (G)
3631 graph *G;
3632 #endif
3633 {
3634 int i, i2, j, k, start, counter=0;
3635 int ncount = G->nnodes;
3636 edge *e;
3637 node *nodelist = G->nodelist;
3638 edge *edgelist = G->edgelist;
3639
3640 for (i = 0; i < ncount; i++) {
3641 if (nodelist[i].status == MATCHED)
3642 counter++;
3643 if (nodelist[i].status == HALVES) {
3644 k = 1;
3645 start = i;
3646 i2 = i;
3647 j = OTHEREND_INT (PEDGE(nodelist[i2].matched_edge), i2);
3648 do {
3649 e = PEDGE(nodelist[j].matched_edge);
3650 if (k % 2) {
3651 nodelist[i2].status = MATCHED;
3652 nodelist[j].status = MATCHED;
3653 edgelist[nodelist[j].matched_edge].x = 0;
3654 nodelist[j].matched_edge = nodelist[i2].matched_edge;
3655 edgelist[nodelist[j].matched_edge].x = 1;
3656 }
3657 i2 = j;
3658 j = OTHEREND_INT (e, i2);
3659 k++;
3660 } while (j!=start);
3661 nodelist[i2].status = UNMATCHED;
3662 edgelist[nodelist[i2].matched_edge].x = 0;
3663 nodelist[i2].matched_edge = -1;
3664 counter++;
3665 }
3666 nodelist[i].label = NIX;
3667 nodelist[i].child = -1;
3668 nodelist[i].sibling = -1;
3669 nodelist[i].parent = -1;
3670 nodelist[i].parent_edge = -1;
3671 nodelist[i].mark = 0;
3672 nodelist[i].tree_root = -1;
3673 }
3674
3675 counter=0;
3676 for (i = 0; i < ncount; i++)
3677 if (nodelist[i].status == MATCHED)
3678 counter++;
3679 // printf(" %i nodes are matched now !\n",counter);
3680
3681 /* build list of unmatched nodes */
3682
3683 G->unmatched_count = 0;
3684 for (i = 0; i < ncount; i++) {
3685 if (nodelist[i].status == UNMATCHED) {
3686 G->unmatched_count++;
3687 }
3688 }
3689
3690 j = -1;
3691 for (i = 0; i < ncount; i++) {
3692 if (nodelist[i].status == UNMATCHED) {
3693 if (j != -1)
3694 nodelist[j].mark = i;
3695 j = i;
3696 }
3697 }
3698 if (j != -1)
3699 nodelist[j].mark = -1;
3700
3701
3702 if (G->unmatched_count == 0) {
3703 G->unmatched = -1;
3704 } else {
3705 for (i = 0; nodelist[i].status != UNMATCHED; i++);
3706 G->unmatched = i;
3707 }
3708 // printf(" unmatched_count = %i\n", G->unmatched_count);
3709 // fflush (stdout);
3710
3711 /* hit is used to identify the blossoms that receive a label + during */
3712 /* the matching run; this will help the pricing code */
3713
3714 for (i = 0; i < 3 * ncount / 2; i++)
3715 nodelist[i].hit = 0;
3716
3717 return 0;
3718 }
3719
3720 #ifdef CC_PROTOTYPE_ANSI
build_graph(graph * G,int ncount,int ecount,int * elist,int * elen)3721 static int build_graph (graph *G, int ncount, int ecount, int *elist, int *elen)
3722 #else
3723 static int build_graph (G, ncount, ecount, elist, elen)
3724 graph *G;
3725 int ncount, ecount;
3726 int *elist, *elen;
3727 #endif
3728 {
3729 int i, j;
3730 edge *e;
3731
3732 if (ncount % 2 != 0) {
3733 fprintf (stderr, "problem has an odd number of nodes\n");
3734 return 1;
3735 }
3736
3737 G->nnodes = ncount;
3738 G->nodelist = CC_SAFE_MALLOC ((3*ncount/2+1), node);
3739 if (!G->nodelist) {
3740 fprintf (stderr, "out of memory in build_graph\n");
3741 return 1;
3742 }
3743 for (i = 0; i < 3*ncount/2 + 1; i++) {
3744 G->nodelist[i].blossom_parent = -1;
3745 G->nodelist[i].matched_edge = -1;
3746 G->nodelist[i].edg_list = -1;
3747 }
3748
3749 /* init unused */
3750
3751 j = 3*ncount/2;
3752 for (i = ncount; i < j; i++) {
3753 G->nodelist[i].mark = i + 1;
3754 }
3755 G->nodelist[j].mark = -1;
3756 G->unused = ncount;
3757
3758 G->nedges = ecount;
3759 G->max_nedges = ecount + (1.5 * ncount);
3760
3761 G->edgelist = CC_SAFE_MALLOC (G->max_nedges, edge);
3762 if (!G->edgelist) {
3763 fprintf (stderr, "out of memory in build_graph\n");
3764 CC_FREE (G->nodelist, node);
3765 return 1;
3766 }
3767
3768 for (i = 0; i < ncount; i++) {
3769 G->nodelist[i].mark = 0;
3770 }
3771
3772 for (i = 0, e = G->edgelist; i < ecount; i++, e++) {
3773 e->nod1 = elist[2*i];
3774 e->nod2 = elist[2*i+1];
3775 e->orig_nod1 = e->nod1;
3776 e->orig_nod2 = e->nod2;
3777 e->slack = 2 * elen[i]; /* to work with ints */
3778 e->mark = 0;
3779 e->ptrs[0] = G->nodelist[e->nod1].edg_list;
3780 G->nodelist[e->nod1].edg_list = 2*i;
3781 e->ptrs[1] = G->nodelist[e->nod2].edg_list;
3782 G->nodelist[e->nod2].edg_list = 2*i + 1;
3783 }
3784 for (; i < G->max_nedges; i++)
3785 G->edgelist[i].mark = 0;
3786
3787 return 0;
3788 }
3789
3790 #ifdef CC_PROTOTYPE_ANSI
write_match(graph * G,CCdatagroup * dat,int * elen,char * fname,int * outp)3791 static int write_match (graph *G, CCdatagroup *dat, int *elen, char* fname, int *outp)
3792 #else
3793 static int write_match (G, dat, elen, fname, outp)
3794 graph *G;
3795 CCdatagroup *dat;
3796 int *elen;
3797 char *fname;
3798 int *outp;
3799 #endif
3800 {
3801 int i, c, len = 0;
3802 edge *e;
3803 FILE* fp;
3804 int KOUNT=0;
3805 double szeit;
3806
3807 szeit = CCutil_zeit();
3808 // printf(" Writing %s ...", fname);
3809 // fflush(stdout);
3810
3811 // if ((fp = fopen (fname, "w")) == (FILE *) NULL) {
3812 // fprintf(stderr," Error: Can't open file %s\n",fname);
3813 // return 1;
3814 // }
3815
3816 c = 0;
3817 for (i = 0; i < G->nnodes; i++) {
3818 if (i == G->edgelist[G->nodelist[i].matched_edge].orig_nod2 ||
3819 G->edgelist[G->nodelist[i].matched_edge].x == HALF) {
3820 c++;
3821 }
3822 }
3823 // printf(" %i nodes, %i edges ", G->nnodes, c);
3824 // fprintf (fp, "%i %i\n", G->nnodes, c);
3825 outp[0] = c;
3826
3827 for (i = 0; i < G->nnodes; i++) {
3828 e = PEDGE(G->nodelist[i].matched_edge);
3829 if (i == e->orig_nod2 || e->x == HALF) {
3830 if (dat)
3831 len = CCutil_dat_edgelen (e->orig_nod1, e->orig_nod2, dat);
3832 else
3833 len = elen[e - G->edgelist];
3834 outp[1+3*KOUNT] = G->edgelist[G->nodelist[i].matched_edge].orig_nod1;
3835 outp[1+3*KOUNT+1] = G->edgelist[G->nodelist[i].matched_edge].orig_nod2;
3836 outp[1+3*KOUNT+2] = len;
3837 KOUNT++;
3838 // fprintf (fp,"%i %i %i\n",
3839 // G->edgelist[G->nodelist[i].matched_edge].orig_nod1,
3840 // G->edgelist[G->nodelist[i].matched_edge].orig_nod2, len);
3841 }
3842 }
3843 // fclose(fp);
3844
3845 // printf(" ... in %.2f sec !! \n", CCutil_zeit () - szeit);
3846 // fflush(stdout);
3847
3848 return 0;
3849 }
3850
3851 #ifdef CC_PROTOTYPE_ANSI
write_blossom_tree(graph * G,char * blossom_tree_file)3852 static int write_blossom_tree (graph *G, char* blossom_tree_file)
3853 #else
3854 static int write_blossom_tree (G, blossom_tree_file)
3855 graph *G;
3856 char *blossom_tree_file;
3857 #endif
3858 {
3859 int i,j;
3860 FILE *fp;
3861 double szeit;
3862 double t;
3863
3864 szeit = CCutil_zeit();
3865 printf(" Writing %s ...",blossom_tree_file);
3866 fflush(stdout);
3867
3868 if ((fp = fopen (blossom_tree_file, "w")) == (FILE *) NULL) {
3869 fprintf(stderr," cannot open file %s\n",blossom_tree_file);
3870 return 1;
3871 }
3872
3873 for (i = 0; i < 3*G->nnodes/2; i++) {
3874 j = G->nodelist[i].blossom_parent;
3875 t = ((double) G->nodelist[i].pi) / 2.0;
3876 fwrite (&j, sizeof(int), 1, fp);
3877 fwrite (&t, sizeof(double), 1, fp);
3878 }
3879 fclose(fp);
3880
3881 printf(" ... in %.2f sec !! \n", CCutil_zeit () - szeit);
3882 fflush(stdout);
3883
3884 return 0;
3885 }
3886