1 /****************************************************************************/
2 /* */
3 /* This file is part of CONCORDE */
4 /* */
5 /* (c) Copyright 1995--1999 by David Applegate, Robert Bixby, */
6 /* Vasek Chvatal, and William Cook */
7 /* */
8 /* Permission is granted for academic research use. For other uses, */
9 /* contact the authors for licensing options. */
10 /* */
11 /* Use at your own risk. We make no guarantees about the */
12 /* correctness or usefulness of this code. */
13 /* */
14 /****************************************************************************/
15
16 /****************************************************************************/
17 /* */
18 /* MIN WEIGHT FRACTIONAL 2-MATCHINGS */
19 /* */
20 /* TSP CODE */
21 /* */
22 /* */
23 /* Written by: Applegate, Bixby, Chvatal, and Cook */
24 /* Date: February 24, 1995 - Modified 10/4/95 (Bico) */
25 /* */
26 /* */
27 /* EXPORTED FUNCTIONS: */
28 /* */
29 /* int CCfmatch_fractional_2match (int ncount, int ecount, int *elist, */
30 /* int *elen, CCdatagroup *dat, double *val, int *thematching, */
31 /* int *thedual, int *thebasis, int wantbasic, int silent, */
32 /* CCrandstate *rstate) */
33 /* int ncount (the number of nodes in the graph) */
34 /* int ecount (the number of edges) */
35 /* int *elist (the edgelist in end1 end2 format) */
36 /* int *elen (the weights on the edges) */
37 /* CCdatagroup *dat (the info to price edges - NULL if no pricing) */
38 /* double *xcoord (the x-coordinates for geometric problems - this */
39 /* field can be NULL) */
40 /* double *ycoord (the y-coordinates) */
41 /* int innorm (the NORM for pricing the complete edgeset) */
42 /* double *val (returns the optimal weight) */
43 /* int *thematching (if non-NULL, then returns the optimal matching */
44 /* in end1 end2 value format, where value is 1 if */
45 /* edge gets assigned 0.5 and value is 2 if edge */
46 /* gets 1.0 - note that the array should be */
47 /* allocated by the calling routine, and should */
48 /* be 6 * ncount + 1 long - it is terminated by a */
49 /* -1) */
50 /* int *thedual (if non-NULL, then returns the optimal dual solution */
51 /* with values twice their actual value (so they will */
52 /* be integers - the array should be alloced by the */
53 /* calling routine, and should be ncount long) */
54 /* int *thebasis (if non-NULL, then returns the edges in the optimal */
55 /* basis in end1 end2 format) */
56 /* int wantbasis (if nonzero, then the optimal basic solution will */
57 /* be found) */
58 /* int silent (if nonzero, will suppress print messages) */
59 /* */
60 /* NOTES: */
61 /* Use to find an optimal basis for the initial tsp LP. By changing */
62 /* MATCHDEGREE to 1, it should find min-wieght fractional matchings. */
63 /* The nodes should be numbered from 0 up to ncount - 1. If dat */
64 /* is specified, then the code will use a price-repair loop to solve */
65 /* the problem over the complete graph. */
66 /* */
67 /****************************************************************************/
68
69 #include "machdefs.h"
70 #include "fmatch.h"
71 #include "kdtree.h"
72 #include "util.h"
73 #include "macrorus.h"
74
75 #define MATCHDEGREE 2 /* Only set up for 1 or 2? */
76 #define MAXWEIGHT 1000000000
77
78 #define ZERO ((unsigned char) 0)
79 #define ONE ((unsigned char) 1)
80 #define TWO ((unsigned char) 2)
81 #ifdef FALSE
82 #undef FALSE
83 #endif
84 #define FALSE ((unsigned char) 0)
85 #ifdef TRUE
86 #undef TRUE
87 #endif
88 #define TRUE ((unsigned char) 1)
89
90 #define OTHEREND(e,n) (e->ends[0] == n ? e->ends[1] : e->ends[0])
91
92 typedef struct edgeptr {
93 struct edge *this;
94 struct node *other;
95 struct edgeptr *next;
96 } edgeptr;
97
98 typedef struct edge {
99 struct edge *next;
100 struct edge *pnext;
101 struct node *ends[2];
102 int weight;
103 int z;
104 unsigned char x;
105 unsigned char basic;
106 } edge;
107
108 typedef struct shortedge {
109 struct node *ends[2];
110 struct shortedge *next;
111 } shortedge;
112
113 typedef struct node {
114 struct node *next;
115 struct node *pnext;
116 edgeptr *adj;
117 edge *parentedge;
118 int name;
119 int y;
120 int label;
121 struct {
122 int order;
123 struct node *next;
124 struct node **prev;
125 } sort;
126 unsigned char flag;
127 unsigned char matchcnt;
128 } node;
129
130 typedef struct graph {
131 int ncount;
132 node *nodelist;
133 edge *edgelist;
134 node **nodenames;
135 int PLUS;
136 int MINUS;
137 CCptrworld node_world;
138 CCptrworld edge_world;
139 CCptrworld edgeptr_world;
140 CCptrworld shortedge_world;
141 } graph;
142
143
144
145 static void
146 init_graph (graph *G),
147 free_graph (graph *G),
148 basic_grab_basic (node *n, int parity, int PLUS, int MINUS),
149 basic_mark_component_as_done (node *n),
150 basic_expand (node *n, int *hit_odd_circuit, int PLUS, int MINUS),
151 basic_minalpha (node *n, node **new, int *alpha, int flip_plus_and_minus,
152 int PLUS, int MINUS),
153 basic_subalpha (node *n, int alpha, int flip_plus_and_minus, int PLUS,
154 int MINUS),
155 initmat (graph *G),
156 initlist (graph *G, node *head, node *tail, node *head2, node *tail2,
157 CCdatagroup *dat),
158 expand (node *n, int *found_one, int PLUS, int MINUS),
159 minalpha (node *n, node **new, int *alpha, int PLUS, int MINUS),
160 subalpha (node *n, int alpha, int PLUS, int MINUS),
161 augmentpath (node *n),
162 augmentpath1 (node *n),
163 augmentplushalf (node *n, edge *e),
164 augmentminushalf (node *n, edge *e),
165 flipcycle (node *n, edge *e, unsigned char v),
166 augmentpair (node *n, node *n1, edge *e),
167 setflag (node *n, unsigned char v),
168 flipseg (node *n, node *n2),
169 stringup (node *n, node *n1, node *n2, edge *e),
170 y_quicksort (node **list, int *y, int l, int u),
171 free_linked_world (graph *G);
172
173 static int
174 buildgraph (graph *G, int ncount, int ecount, int *elist, int *elen),
175 basicrun (graph *G),
176 basicscan (graph *G, node *n),
177 basic_check_scan (graph *G, node *n),
178 basicgrow (graph *G, node *n),
179 basic_grab_ones (node *n, int parity, edge **odd_circuit, int PLUS,
180 int MINUS),
181 basic_checkout_basic (node *n, int parity, edge **odd_circuit, int PLUS,
182 int MINUS),
183 twomatch (graph *G),
184 chkmat (graph *G, double *val),
185 fixmatch (graph *G, int *radded, CCdatagroup *dat, CCrandstate *rstate),
186 kd_fixmatch (graph *G, int *radded, CCdatagroup *dat, CCrandstate *rstate),
187 x_fixmatch (graph *G, int *radded, CCdatagroup *dat),
188 junk_fixmatch (graph *G, int *radded, CCdatagroup *dat),
189 checkoutedge (graph *G, node *n1, node *n2, int *hit, CCdatagroup *dat),
190 precheckoutedge (node *n1, node *n2, shortedge **list, CCdatagroup *dat,
191 CCptrworld *shortedge_world),
192 addbadedge (graph *G, node *n1, node *n2, int w),
193 augment (graph *G, node *n);
194
195 static node
196 *basic_dualchange (node *n, int PLUS, int MINUS),
197 *dualchange (node *n, int PLUS, int MINUS),
198 *findflag (node *n),
199 *findhole (node *n, node *n2);
200
201 static edge
202 *newedge (graph *G, node *n1, node *n2),
203 *findedge (node *n1, node *n2);
204
205
206
207
208 /********** Allocation routines **********/
209
210
CC_PTRWORLD_ROUTINES(node,nodealloc,node_bulk_alloc,nodefree)211 CC_PTRWORLD_ROUTINES (node, nodealloc, node_bulk_alloc, nodefree)
212 CC_PTRWORLD_LEAKS_ROUTINE (node, node_check_leaks, flag, unsigned char)
213
214 CC_PTRWORLD_ROUTINES (edge, edgealloc, edge_bulk_alloc, edgefree)
215 CC_PTRWORLD_LISTFREE_ROUTINE (edge, edge_listfree, edgefree)
216 CC_PTRWORLD_LEAKS_ROUTINE (edge, edge_check_leaks, basic, unsigned char)
217
218 CC_PTRWORLD_ROUTINES (edgeptr, edgeptralloc, edgeptr_bulk_alloc, edgeptrfree)
219 CC_PTRWORLD_LISTFREE_ROUTINE (edgeptr, edgeptr_listfree, edgeptrfree)
220 CC_PTRWORLD_LEAKS_ROUTINE (edgeptr, edgeptr_check_leaks, this, edge *)
221
222 CC_PTRWORLD_ROUTINES (shortedge, shortedgealloc, shortedge_bulk_alloc,
223 shortedgefree)
224 CC_PTRWORLD_LEAKS_ROUTINE (shortedge, shortedge_check_leaks, ends[0], node *)
225
226
227 static void free_linked_world (graph *G)
228 {
229 int ntotal;
230 int nreserve;
231 if (node_check_leaks (&G->node_world, &ntotal, &nreserve)) {
232 fprintf (stderr, "WARNING: Outstanding nodes %d (total %d)\n",
233 ntotal - nreserve, ntotal);
234 }
235 CCptrworld_delete (&G->node_world);
236
237 if (edge_check_leaks (&G->edge_world, &ntotal, &nreserve)) {
238 fprintf (stderr, "WARNING: Outstanding edges %d (total %d)\n",
239 ntotal - nreserve, ntotal);
240 }
241 CCptrworld_delete (&G->edge_world);
242
243 if (edgeptr_check_leaks (&G->edgeptr_world, &ntotal, &nreserve)) {
244 fprintf (stderr, "WARNING: Outstanding edgeptrs %d (total %d)\n",
245 ntotal - nreserve, ntotal);
246 }
247 CCptrworld_delete (&G->edgeptr_world);
248
249 if (shortedge_check_leaks (&G->shortedge_world, &ntotal, &nreserve)) {
250 fprintf (stderr, "WARNING: Outstanding shortedges %d (total %d)\n",
251 ntotal - nreserve, ntotal);
252 }
253 CCptrworld_delete (&G->shortedge_world);
254 }
255
256
257 /********** main routines **********/
258
259
buildgraph(graph * G,int ncount,int ecount,int * elist,int * elen)260 static int buildgraph (graph *G, int ncount, int ecount, int *elist, int *elen)
261 {
262 int i;
263 node *n;
264 edge *e;
265
266 init_graph (G);
267
268 G->ncount = ncount;
269 G->nodenames = CC_SAFE_MALLOC (ncount, node *);
270 if (G->nodenames == (node **) NULL) {
271 fprintf (stderr, "out of memory in buildgraph\n"); return 1;
272 }
273
274 for (i = 0; i < ncount; i++) {
275 n = nodealloc (&G->node_world);
276 if (n == (node *) NULL) {
277 fprintf (stderr, "out of memory in buildgraph\n"); return 1;
278 }
279 n->name = i;
280 n->adj = (edgeptr *) NULL;
281 n->next = G->nodelist;
282 G->nodelist = n;
283 G->nodenames[i] = n;
284 }
285 for (i = 0; i < ecount; i++) {
286 e = newedge (G, G->nodenames[elist[2 * i]],
287 G->nodenames[elist[(2 * i) + 1]]);
288 if (e == (edge *) NULL) {
289 fprintf (stderr, "out of memory in buildgraph\n"); return 1;
290 }
291 e->weight = elen[i] + elen[i];
292 }
293 return 0;
294 }
295
init_graph(graph * G)296 static void init_graph (graph *G)
297 {
298 G->ncount = 0;
299 G->nodenames = (node **) NULL;
300 G->nodelist = (node *) NULL;
301 G->edgelist = (edge *) NULL;
302 G->PLUS = 1;
303 G->MINUS = 2;
304 CCptrworld_init (&G->node_world);
305 CCptrworld_init (&G->edge_world);
306 CCptrworld_init (&G->edgeptr_world);
307 CCptrworld_init (&G->shortedge_world);
308 }
309
free_graph(graph * G)310 static void free_graph (graph *G)
311 {
312 node *n, *nnext;
313
314 for (n = G->nodelist; n; n = nnext) {
315 nnext = n->next;
316 edgeptr_listfree (&G->edgeptr_world, n->adj);
317 nodefree (&G->node_world, n);
318 }
319 G->nodelist = (node *) NULL;
320
321 edge_listfree (&G->edge_world, G->edgelist);
322 G->edgelist = (edge *) NULL;
323
324 CC_IFFREE (G->nodenames, node *);
325 G->ncount = 0;
326 G->PLUS = 0;
327 G->MINUS = 0;
328 }
329
330
newedge(graph * G,node * n1,node * n2)331 static edge *newedge (graph *G, node *n1, node *n2)
332 {
333 edge *e;
334 edgeptr *ep;
335
336 e = edgealloc (&G->edge_world);
337 if (e == (edge *) NULL)
338 return (edge *) NULL;
339 e->ends[0] = n1;
340 e->ends[1] = n2;
341 e->next = G->edgelist;
342 G->edgelist = e;
343
344 ep = edgeptralloc (&G->edgeptr_world);
345 if (ep == (edgeptr *) NULL) {
346 edgefree (&G->edge_world, e);
347 return (edge *) NULL;
348 }
349 ep->this = e;
350 ep->other = n2;
351 ep->next = n1->adj;
352 n1->adj = ep;
353
354 ep = edgeptralloc (&G->edgeptr_world);
355 if (ep == (edgeptr *) NULL) {
356 edgefree (&G->edge_world, e);
357 return (edge *) NULL;
358 }
359 ep->this = e;
360 ep->other = n1;
361 ep->next = n2->adj;
362 n2->adj = ep;
363
364 return e;
365 }
366
CCfmatch_fractional_2match(int ncount,int ecount,int * elist,int * elen,CCdatagroup * dat,double * val,int * thematching,int * thedual,int * thebasis,int wantbasic,int silent,CCrandstate * rstate)367 int CCfmatch_fractional_2match (int ncount, int ecount, int *elist, int *elen,
368 CCdatagroup *dat, double *val, int *thematching, int *thedual,
369 int *thebasis, int wantbasic, int silent, CCrandstate *rstate)
370 {
371 double v, vbasic, szeit, tzeit;
372 int added;
373 int rval = 0;
374 graph G;
375
376 tzeit = CCutil_zeit ();
377
378 rval = buildgraph (&G, ncount, ecount, elist, elen);
379 if (rval) {
380 fprintf (stderr, "buildgraph failed\n"); goto CLEANUP;
381 }
382
383 initmat (&G);
384 rval = twomatch (&G);
385 if (rval) {
386 fprintf (stderr, "twomatch failed\n"); goto CLEANUP;
387 }
388 rval = chkmat (&G, &v);
389 if (rval) {
390 fprintf (stderr, "Chkmat found error in matching\n"); goto CLEANUP;
391 }
392
393 if (!silent) {
394 printf ("Fractional Matching: %.1f\n", v);
395 printf ("Initial Running Time: %.2f (seconds)\n",
396 CCutil_zeit () - tzeit);
397 fflush (stdout);
398 }
399
400 if (wantbasic) {
401 szeit = CCutil_zeit ();
402 if (basicrun (&G)) {
403 fprintf (stderr, "Did not find a basic optimal solution\n");
404 rval = 1;
405 goto CLEANUP;
406 }
407 if (chkmat (&G, &vbasic)) {
408 fprintf (stderr, "Chkmat found error in matching\n");
409 rval = 1;
410 goto CLEANUP;
411 }
412 if (vbasic != v) {
413 fprintf (stderr, "ERROR: Basis routine altered objective\n");
414 rval = 1;
415 goto CLEANUP;
416 }
417 if (!silent) {
418 printf ("Basis Running Time: %.2f (seconds)\n",
419 CCutil_zeit () - szeit);
420 fflush (stdout);
421 }
422 }
423
424 if (dat != (CCdatagroup *) NULL) {
425 if (!silent) {
426 printf ("Price-Repair ...\n"); fflush (stdout);
427 }
428 szeit = CCutil_zeit ();
429 rval = fixmatch (&G, &added, dat, rstate);
430 if (rval) {
431 fprintf (stderr, "fixmatch failed\n"); goto CLEANUP;
432 }
433 if (chkmat (&G, &v)) {
434 fprintf (stderr, "Chkmat found error in matching\n");
435 rval = 1;
436 goto CLEANUP;
437 }
438 if (wantbasic) {
439 do {
440 if (!silent) {
441 printf ("Find basis ...\n"); fflush (stdout);
442 }
443 if (basicrun (&G)) {
444 fprintf (stderr, "Did not find a basic solution\n");
445 rval = 1;
446 goto CLEANUP;
447 }
448 if (chkmat (&G, &vbasic)) {
449 fprintf (stderr, "Chkmat found error in matching\n");
450 rval = 1;
451 goto CLEANUP;
452 }
453 if (vbasic != v) {
454 fprintf (stderr, "ERROR: Basis routine altered obj\n");
455 rval = 1;
456 goto CLEANUP;
457 }
458
459 if (!silent) {
460 printf ("Price-repair basic solution ...\n");
461 fflush (stdout);
462 }
463 rval = fixmatch (&G, &added, dat, rstate);
464 if (rval) {
465 fprintf (stderr, "fixmatch failed\n"); goto CLEANUP;
466 }
467 if (chkmat (&G, &v)) {
468 fprintf (stderr, "Chkmat found error in matching\n");
469 rval = 1;
470 goto CLEANUP;
471 }
472 } while (added);
473 }
474 if (!silent) {
475 printf ("Running Time for Price-Repair: %.2f\n",
476 CCutil_zeit () - szeit);
477 printf ("Fractional Matching on Complete Graph: %.1f\n",v);
478 fflush (stdout);
479 }
480 }
481
482 *val = v;
483
484 if (thematching != (int *) NULL) {
485 int k = 0;
486 edge *e;
487 for (e = G.edgelist; e; e = e->next) {
488 if ((unsigned int) e->x != (unsigned int) ZERO) {
489 thematching[k++] = e->ends[0]->name;
490 thematching[k++] = e->ends[1]->name;
491 thematching[k++] = ((unsigned int) e->x == (unsigned int) ONE ? 1 : 2);
492 }
493 }
494 thematching[k] = -1;
495 }
496 if (thedual != (int *) NULL) {
497 int i;
498 for (i = 0; i < G.ncount; i++)
499 thedual[i] = G.nodenames[i]->y;
500 }
501 if (wantbasic && thebasis != (int *) NULL) {
502 int k = 0;
503 edge *e;
504 for (e = G.edgelist; e; e = e->next) {
505 if ((unsigned int) e->basic) {
506 thebasis[k++] = e->ends[0]->name;
507 thebasis[k++] = e->ends[1]->name;
508 }
509 }
510 }
511
512 CLEANUP:
513
514 free_graph (&G);
515 free_linked_world (&G);
516
517 if (!silent) {
518 printf ("Total fractional matching time: %.2f (seconds)\n",
519 CCutil_zeit () - tzeit);
520 fflush (stdout);
521 }
522
523 return rval;
524 }
525
initmat(graph * G)526 static void initmat (graph *G)
527 {
528 node *n;
529 edge *e;
530
531 for (n = G->nodelist; n; n = n->next) {
532 n->y = 0;
533 n->matchcnt = 2 - MATCHDEGREE;
534 n->parentedge = (edge *) NULL;
535 n->label = 0;
536 }
537 for (e = G->edgelist; e; e = e->next) {
538 e->z = 0;
539 e->x = ZERO;
540 e->pnext = (edge *) NULL;
541 }
542 G->PLUS = 1;
543 G->MINUS = 2;
544 }
545
chkmat(graph * G,double * val)546 static int chkmat (graph *G, double *val) /* val may be 1/2 integer and big */
547 {
548 int k;
549 double v = 0.0, dualval = 0.0;
550 node *n;
551 edge *e;
552 edgeptr *ep;
553
554 for (n = G->nodelist; n; n = n->next) {
555 k = 0;
556 for (ep = n->adj; ep; ep = ep->next)
557 k += (unsigned int) ep->this->x;
558 if (k != 2 * MATCHDEGREE) {
559 fprintf (stderr, "Not a matching, node %d has 2-degree %d\n",
560 n->name, k);
561 return 1;
562 }
563 dualval += (double) n->y;
564 }
565 dualval *= (double) MATCHDEGREE;
566 for (e = G->edgelist; e; e = e->next) {
567 switch ((unsigned int) e->x) {
568 case (unsigned int) TWO:
569 if (e->z < 0 ||
570 e->z != e->ends[0]->y + e->ends[1]->y - e->weight) {
571 fprintf (stderr, "Error in dual solution - 2\n");
572 return 1;
573 }
574 v += (double) e->weight;
575 v += (double) e->weight;
576 dualval -= (double) e->z;
577 break;
578 case (unsigned int) ONE:
579 if (e->z != 0 ||
580 e->ends[0]->y + e->ends[1]->y != e->weight) {
581 fprintf (stderr, "Error in dual solution - 1\n");
582 return 1;
583 }
584 v += (double) e->weight;
585 break;
586 case (unsigned int) ZERO:
587 if (e->z != 0 || e->ends[0]->y + e->ends[1]->y > e->weight) {
588 fprintf (stderr, "Error in dual solution - 0\n");
589 return 1;
590 }
591 break;
592 default:
593 fprintf (stderr, "Error in matching values\n");
594 return 1;
595 }
596 }
597 v /= 4.0;
598 dualval /= 2.0;
599
600 if (v != dualval) {
601 fprintf (stderr, "The primal and dual objective values differ.\n");
602 return 1;
603 }
604 *val = v;
605 return 0;
606 }
607
608
609 /********** Core fractional matching routines **********/
610
611
twomatch(graph * G)612 static int twomatch (graph *G)
613 {
614 node *n;
615 int rval = 0;
616
617 for (n = G->nodelist; n; n = n->next) {
618 while ((unsigned int) n->matchcnt != (unsigned int) TWO) {
619 rval = augment (G, n);
620 if (rval) {
621 fprintf (stderr, "augment failed - probably no fmatching\n");
622 return rval;
623 }
624 }
625 }
626 return 0;
627 }
628
augment(graph * G,node * n)629 static int augment (graph *G, node *n)
630 {
631 node *auglist;
632 int found_augmenting_path = 0;
633
634 G->PLUS += 2;
635 G->MINUS += 2;
636 n->label = G->PLUS;
637 n->parentedge = (edge *) NULL;
638 n->matchcnt = (unsigned int) n->matchcnt + 1;
639 expand (n, &found_augmenting_path, G->PLUS, G->MINUS);
640 if (found_augmenting_path)
641 return 0;
642 while ((auglist = dualchange (n, G->PLUS, G->MINUS)) != (node *) NULL) {
643 while (auglist) {
644 expand (auglist, &found_augmenting_path, G->PLUS, G->MINUS);
645 if (found_augmenting_path)
646 return 0;
647 auglist = auglist->pnext;
648 }
649 }
650 fprintf (stderr, "Error - dual change did not create new edges\n");
651 return 1;
652 }
653
expand(node * n,int * found_one,int PLUS,int MINUS)654 static void expand (node *n, int *found_one, int PLUS, int MINUS)
655 {
656 node *n1;
657 edgeptr *ep;
658 edge *e;
659
660 if (n->label == PLUS) {
661 for (ep = n->adj; ep; ep = ep->next) {
662 if ((unsigned int) ep->this->x == (unsigned int) ONE) {
663 augmentplushalf (n, ep->this);
664 *found_one = 1;
665 return;
666 }
667 }
668 for (ep = n->adj; ep; ep = ep->next) {
669 if ((unsigned int) ep->this->x == (unsigned int) ZERO &&
670 ep->other->y + n->y == ep->this->weight) {
671 e = ep->this;
672 n1 = ep->other;
673 if (n1->label < PLUS) { /* n1 has no label */
674 n1->label = MINUS;
675 n1->parentedge = e;
676 if ((unsigned int) n1->matchcnt != (unsigned int) TWO) {
677 augmentpath (n1);
678 *found_one = 1;
679 return;
680 }
681 expand (n1, found_one, PLUS, MINUS);
682 if (*found_one)
683 return;
684 } else if (n1->label == PLUS) {
685 augmentpair (n, n1, e);
686 *found_one = 1;
687 return;
688 }
689 }
690 }
691 } else { /* MINUS */
692 for (ep = n->adj; ep; ep = ep->next) {
693 if ((unsigned int) ep->this->x == (unsigned int) ONE) {
694 augmentminushalf (n, ep->this);
695 *found_one = 1;
696 return;
697 }
698 }
699 for (ep = n->adj; ep; ep = ep->next) {
700 if ((unsigned int) ep->this->x == (unsigned int) TWO &&
701 ep->this->z == 0) {
702 e = ep->this;
703 n1 = ep->other;
704 if (n1->label < PLUS) { /* n1 has no label */
705 n1->label = PLUS;
706 n1->parentedge = e;
707 expand (n1, found_one, PLUS, MINUS);
708 if (*found_one)
709 return;
710 } else if (n1->label == MINUS) {
711 augmentpair (n, n1, e);
712 *found_one = 1;
713 return;
714 }
715 }
716 }
717 }
718 return;
719 }
720
dualchange(node * n,int PLUS,int MINUS)721 static node *dualchange (node *n, int PLUS, int MINUS)
722 {
723 node *new = (node *) NULL;
724 int alpha = MAXWEIGHT;
725
726 minalpha (n, &new, &alpha, PLUS, MINUS);
727 if (alpha == MAXWEIGHT) {
728 fprintf (stderr, "Dual change required, but no candidate edges\n");
729 return (node *) NULL;
730 }
731 if (alpha & 0x1) {
732 fprintf (stderr, "Whoops, 2 * alpha = %d, not even\n", alpha);
733 return (node *) NULL;
734 }
735 alpha /= 2;
736 subalpha (n, alpha, PLUS, MINUS);
737 return new;
738 }
739
minalpha(node * n,node ** new,int * alpha,int PLUS,int MINUS)740 static void minalpha (node *n, node **new, int *alpha, int PLUS, int MINUS)
741 {
742 int minv = MAXWEIGHT;
743 int thisv;
744 node *n1;
745 edgeptr *ep;
746 edge *e;
747
748 if (n->label == PLUS) {
749 for (ep = n->adj; ep; ep = ep->next) {
750 e = ep->this;
751 if ((unsigned int) e->x == (unsigned int) ZERO) {
752 n1 = ep->other;
753 if (n1->label < PLUS) { /* n1 is unlabeled */
754 thisv = e->weight - n->y - n1->y;
755 thisv += thisv;
756 if (thisv < minv)
757 minv = thisv;
758 } else if (n1->label == PLUS) {
759 thisv = e->weight - n->y - n1->y;
760 if (thisv < minv)
761 minv = thisv;
762 } else { /* n1 has a minus label */
763 if (n1->parentedge == e)
764 minalpha (n1, new, alpha, PLUS, MINUS);
765 }
766 }
767 }
768 } else { /* MINUS case */
769 for (ep = n->adj; ep; ep = ep->next) {
770 e = ep->this;
771 if ((unsigned int) e->x == (unsigned int) TWO) {
772 n1 = ep->other;
773 if (n1->label < PLUS) {
774 thisv = e->z + e->z;
775 if (thisv < minv)
776 minv = thisv;
777 } else if (n1->label == PLUS) {
778 if (n1->parentedge == e)
779 minalpha (n1, new, alpha, PLUS, MINUS);
780 } else { /* n1 has a MINUS label */
781 thisv = e->z;
782 if (thisv < minv)
783 minv = thisv;
784 }
785 }
786 }
787 }
788 if (minv < *alpha) {
789 *new = n;
790 n->pnext = (node *) NULL;
791 *alpha = minv;
792 } else if (minv == *alpha) {
793 n->pnext = *new;
794 *new = n;
795 }
796 }
797
subalpha(node * n,int alpha,int PLUS,int MINUS)798 static void subalpha (node *n, int alpha, int PLUS, int MINUS)
799 {
800 edgeptr *ep;
801 edge *e;
802 node *n1;
803
804 if (n->label == PLUS) {
805 n->y += alpha;
806 for (ep = n->adj; ep; ep = ep->next) {
807 e = ep->this;
808 if ((unsigned int) e->x == (unsigned int) TWO)
809 e->z += alpha;
810 else if ((unsigned int) e->x == (unsigned int) ZERO) {
811 n1 = ep->other;
812 if (n1->parentedge == e && n1->label >= PLUS)
813 subalpha (n1, alpha, PLUS, MINUS);
814 }
815 }
816 } else { /* MINUS */
817 n->y -= alpha;
818 for (ep = n->adj; ep; ep = ep->next) {
819 e = ep->this;
820 if ((unsigned int) e->x == (unsigned int) TWO) {
821 e->z -= alpha;
822 n1 = ep->other;
823 if (n1->parentedge == e && n1->label >= PLUS)
824 subalpha (n1, alpha, PLUS, MINUS);
825 }
826 }
827 }
828 }
829
augmentpath(node * n)830 static void augmentpath (node *n)
831 {
832 n->matchcnt = (unsigned int) n->matchcnt + 1;
833 augmentpath1 (n);
834 }
835
augmentpath1(node * n)836 static void augmentpath1 (node *n)
837 {
838 while (n->parentedge != (edge *) NULL) {
839 n->parentedge->x = (unsigned int) TWO - (unsigned int) n->parentedge->x;
840 n = OTHEREND(n->parentedge, n);
841 }
842 }
843
augmentplushalf(node * n,edge * e)844 static void augmentplushalf (node *n, edge *e)
845 {
846 flipcycle (n, e, TWO);
847 augmentpath1 (n);
848 }
849
augmentminushalf(node * n,edge * e)850 static void augmentminushalf (node *n, edge *e)
851 {
852 flipcycle (n, e, ZERO);
853 augmentpath1 (n);
854 }
855
flipcycle(node * n,edge * e,unsigned char v)856 static void flipcycle (node *n, edge *e, unsigned char v)
857 {
858 edge *e1;
859 edge *e2;
860
861 e1 = e->pnext;
862 if (e1->ends[0] == n || e1->ends[1] == n)
863 e = e1;
864 e1 = e;
865 do {
866 e1->x = v;
867 v = (unsigned int) TWO - (unsigned int) v;
868 e2 = e1->pnext;
869 e1->pnext = (edge *) NULL;
870 e1 = e2;
871 } while (e1 != e);
872 }
873
augmentpair(node * n,node * n1,edge * e)874 static void augmentpair (node *n, node *n1, edge *e)
875 {
876 node *n2, *n3;
877
878 setflag (n, FALSE);
879 setflag (n1, TRUE);
880 n2 = findflag (n);
881 if ((n3 = findhole (n, n2)) != (node *) NULL) {
882 n3->matchcnt = (unsigned int) n3->matchcnt + 1;
883 flipseg (n, n3);
884 e->x = (unsigned int) TWO - (unsigned int) e->x;
885 augmentpath1 (n1);
886 return;
887 }
888 if ((n3 = findhole (n1, n2)) != (node *) NULL) {
889 n3->matchcnt = (unsigned int) n3->matchcnt + 1;
890 flipseg (n1, n3);
891 e->x = (unsigned int) TWO - (unsigned int) e->x;
892 augmentpath1 (n);
893 return;
894 }
895 stringup (n, n1, n2, e);
896 augmentpath1 (n2);
897 }
898
setflag(node * n,unsigned char v)899 static void setflag (node *n, unsigned char v)
900 {
901 n->flag = v;
902 while (n->parentedge != (edge *) NULL) {
903 n = OTHEREND(n->parentedge, n);
904 n->flag = v;
905 }
906 }
907
findflag(node * n)908 static node *findflag (node *n)
909 {
910 while ((unsigned int) n->flag == 0) {
911 n = OTHEREND(n->parentedge, n);
912 }
913 return n;
914 }
915
findhole(node * n,node * n2)916 static node *findhole (node *n, node *n2)
917 {
918 while (n != n2) {
919 if ((unsigned int) n->matchcnt != (unsigned int) TWO)
920 return n;
921 n = OTHEREND(n->parentedge, n);
922 }
923 return (unsigned int) n2->matchcnt != (unsigned int) TWO ? n2
924 : (node *) NULL;
925 }
926
flipseg(node * n,node * n2)927 static void flipseg (node *n, node *n2)
928 {
929 while (n != n2) {
930 n->parentedge->x = (unsigned int) TWO - (unsigned int) n->parentedge->x;
931 n = OTHEREND(n->parentedge, n);
932 }
933 }
934
stringup(node * n,node * n1,node * n2,edge * e)935 static void stringup (node *n, node *n1, node *n2, edge *e)
936 {
937 edge *preve, *savee;
938
939 preve = e;
940 while (n != n2) {
941 n->parentedge->x = ONE;
942 preve->pnext = n->parentedge;
943 preve = n->parentedge;
944 n = OTHEREND(n->parentedge, n);
945 }
946 savee = preve;
947 preve = e;
948 while (n1 != n2) {
949 n1->parentedge->x = ONE;
950 n1->parentedge->pnext = preve;
951 preve = n1->parentedge;
952 n1 = OTHEREND(n1->parentedge, n1);
953 }
954 savee->pnext = preve;
955 e->x = ONE;
956 }
957
findedge(node * n1,node * n2)958 static edge *findedge (node *n1, node *n2)
959 {
960 edgeptr *ep;
961
962 for (ep = n1->adj; ep; ep = ep->next) {
963 if (ep->other == n2)
964 return ep->this;
965 }
966 return (edge *) NULL;
967 }
968
969
970 /********** Basis finding routines **********/
971
972
basicrun(graph * G)973 static int basicrun (graph *G)
974 {
975 node *n;
976 edge *e;
977
978 G->PLUS = 1;
979 G->MINUS = 2;
980
981 for (n = G->nodelist; n; n = n->next) {
982 n->label = 0;
983 n->flag = 0;
984 }
985 for (e = G->edgelist; e; e = e->next)
986 e->basic = 0;
987
988 for (n = G->nodelist; n; n = n->next) {
989 if (n->label == 0) {
990 if (basicscan (G, n))
991 return 1;
992 }
993 }
994 for (n = G->nodelist; n; n = n->next) {
995 if ((unsigned int) n->flag == 0) {
996 if (basicgrow (G, n))
997 return 1;
998 }
999 }
1000
1001 for (n = G->nodelist; n; n = n->next)
1002 n->label = 0;
1003 for (n = G->nodelist; n; n = n->next) {
1004 if (n->label == 0) {
1005 if (basic_check_scan (G, n))
1006 return 1;
1007 }
1008 }
1009
1010 return 0;
1011 }
1012
basicscan(graph * G,node * n)1013 static int basicscan (graph *G, node *n)
1014 {
1015 edge *odd_circuit = (edge *) NULL;
1016
1017 G->PLUS += 2;
1018 G->MINUS += 2;
1019 n->parentedge = (edge *) NULL;
1020 if (basic_grab_ones (n, 0, &odd_circuit, G->PLUS, G->MINUS))
1021 return 1;
1022 if (odd_circuit != (edge *) NULL) {
1023 basic_mark_component_as_done (n);
1024 }
1025 return 0;
1026 }
1027
basic_check_scan(graph * G,node * n)1028 static int basic_check_scan (graph *G, node *n)
1029 {
1030 edge *odd_circuit = (edge *) NULL;
1031
1032 G->PLUS += 2;
1033 G->MINUS += 2;
1034 n->parentedge = (edge *) NULL;
1035 if (basic_checkout_basic (n, 0, &odd_circuit, G->PLUS, G->MINUS))
1036 return 1;
1037 if (odd_circuit == (edge *) NULL) {
1038 printf ("No odd circuit\n");
1039 return 1;
1040 }
1041 return 0;
1042 }
1043
basicgrow(graph * G,node * n)1044 static int basicgrow (graph *G, node *n)
1045 {
1046 int hit_odd_circuit = 0;
1047 node *expandlist;
1048
1049 G->PLUS += 2;
1050 G->MINUS += 2;
1051 basic_grab_basic (n, 0, G->PLUS, G->MINUS);
1052
1053 n->parentedge = (edge *) NULL;
1054 basic_expand (n, &hit_odd_circuit, G->PLUS, G->MINUS);
1055 if (hit_odd_circuit)
1056 return 0;
1057 else {
1058 while ((expandlist = basic_dualchange (n, G->PLUS, G->MINUS))
1059 != (node *) NULL) {
1060 while (expandlist) {
1061 basic_expand (expandlist, &hit_odd_circuit, G->PLUS, G->MINUS);
1062 if (hit_odd_circuit)
1063 return 0;
1064 expandlist = expandlist->pnext;
1065 }
1066 }
1067 fprintf (stderr, "ERROR: No dual change in basis finding code\n");
1068 return 1;
1069 }
1070 }
1071
basic_grab_ones(node * n,int parity,edge ** odd_circuit,int PLUS,int MINUS)1072 static int basic_grab_ones (node *n, int parity, edge **odd_circuit,
1073 int PLUS, int MINUS)
1074 {
1075 edge *e;
1076 edgeptr *ep;
1077 node *n1;
1078
1079 n->label = (parity == 0 ? PLUS : MINUS);
1080 for (ep = n->adj; ep; ep = ep->next) {
1081 e = ep->this;
1082 if ((unsigned int) e->x == (unsigned int) ONE && e != n->parentedge) {
1083 n1 = ep->other;
1084 if (n1->label == 0) {
1085 n1->parentedge = e;
1086 e->basic = 1;
1087 if (basic_grab_ones (n1, 1 - parity, odd_circuit, PLUS, MINUS))
1088 return 1;
1089 } else if (n1->label == n->label) {
1090 if (*odd_circuit == (edge *) NULL) {
1091 *odd_circuit = e;
1092 e->basic = 1;
1093 } else if (*odd_circuit != e) {
1094 fprintf (stderr, "ERROR: Two odd circuits in 1-graph\n");
1095 printf ("Circuit forming edges: %d-%d %d-%d\n",
1096 (*odd_circuit)->ends[0]->name,
1097 (*odd_circuit)->ends[1]->name,
1098 e->ends[0]->name,
1099 e->ends[1]->name);
1100 return 1;
1101 }
1102 } else {
1103 fprintf (stderr, "ERROR: Even circuit in 1-graph\n");
1104 printf ("Circuit forming edge: %d-%d\n",
1105 e->ends[0]->name,
1106 e->ends[1]->name);
1107 return 1;
1108 }
1109 }
1110 }
1111 return 0;
1112 }
1113
basic_checkout_basic(node * n,int parity,edge ** odd_circuit,int PLUS,int MINUS)1114 static int basic_checkout_basic (node *n, int parity, edge **odd_circuit,
1115 int PLUS, int MINUS)
1116 {
1117 edge *e;
1118 edgeptr *ep;
1119 node *n1;
1120
1121 n->label = (parity == 0 ? PLUS : MINUS);
1122 for (ep = n->adj; ep; ep = ep->next) {
1123 e = ep->this;
1124 if ((unsigned int) e->basic && e != n->parentedge) {
1125 n1 = ep->other;
1126 if (n1->label == 0) {
1127 n1->parentedge = e;
1128 if (basic_checkout_basic (n1, 1 - parity, odd_circuit,
1129 PLUS, MINUS))
1130 return 1;
1131 } else if (n1->label == n->label) {
1132 if (*odd_circuit == (edge *) NULL) {
1133 *odd_circuit = e;
1134 } else if (*odd_circuit != e) {
1135 fprintf (stderr, "ERROR: Two odd circuits in basish\n");
1136 printf ("Circuit forming edges: %d-%d %d-%d\n",
1137 (*odd_circuit)->ends[0]->name,
1138 (*odd_circuit)->ends[1]->name,
1139 e->ends[0]->name,
1140 e->ends[1]->name);
1141 return 1;
1142 }
1143 } else {
1144 fprintf (stderr, "ERROR: Even circuit in basis\n");
1145 printf ("Circuit forming edge: %d-%d\n",
1146 e->ends[0]->name,
1147 e->ends[1]->name);
1148 return 1;
1149 }
1150 }
1151 }
1152 return 0;
1153 }
1154
basic_grab_basic(node * n,int parity,int PLUS,int MINUS)1155 static void basic_grab_basic (node *n, int parity, int PLUS, int MINUS)
1156 {
1157 edgeptr *ep;
1158
1159 n->label = (parity == 0 ? PLUS : MINUS);
1160 n->flag = 1;
1161 for (ep = n->adj; ep; ep = ep->next) {
1162 if ((unsigned int) ep->this->basic) {
1163 if ((unsigned int) ep->other->flag == 0)
1164 basic_grab_basic (ep->other, 1 - parity, PLUS, MINUS);
1165 }
1166 }
1167 }
1168
basic_mark_component_as_done(node * n)1169 static void basic_mark_component_as_done (node *n)
1170 {
1171 edgeptr *ep;
1172
1173 n->flag = 1;
1174 for (ep = n->adj; ep; ep = ep->next) {
1175 if ((unsigned int) ep->this->basic) {
1176 if ((unsigned int) ep->other->flag == 0)
1177 basic_mark_component_as_done (ep->other);
1178 }
1179 }
1180 }
1181
basic_expand(node * n,int * hit_odd_circuit,int PLUS,int MINUS)1182 static void basic_expand (node *n, int *hit_odd_circuit, int PLUS, int MINUS)
1183 {
1184 edge *e;
1185 edgeptr *ep;
1186 node *n1;
1187
1188 for (ep = n->adj; ep; ep = ep->next) {
1189 e = ep->this;
1190 if (e != n->parentedge) {
1191 n1 = ep->other;
1192 if ((unsigned int) e->basic) {
1193 n1->parentedge = e;
1194 basic_expand (n1, hit_odd_circuit, PLUS, MINUS);
1195 if (*hit_odd_circuit)
1196 return;
1197 } else if (n->y + n1->y == e->weight) {
1198 if (n1->label < PLUS) {
1199 e->basic = 1;
1200 if ((unsigned int) n1->flag) {
1201 *hit_odd_circuit = 1;
1202 return;
1203 } else {
1204 n1->parentedge = e;
1205 if (n->label == PLUS)
1206 basic_grab_basic (n1, 1, PLUS, MINUS);
1207 else
1208 basic_grab_basic (n1, 0, PLUS, MINUS);
1209 basic_expand (n1, hit_odd_circuit, PLUS, MINUS);
1210 if (*hit_odd_circuit)
1211 return;
1212 }
1213 } else if (n1->label == n->label) {
1214 e->basic = 1;
1215 *hit_odd_circuit = 1;
1216 return;
1217 }
1218 }
1219 }
1220 }
1221 }
1222
basic_dualchange(node * n,int PLUS,int MINUS)1223 static node *basic_dualchange (node *n, int PLUS, int MINUS)
1224 {
1225 node *new = (node *) NULL;
1226 int alpha = MAXWEIGHT;
1227
1228 basic_minalpha (n, &new, &alpha, 0, PLUS, MINUS);
1229 if (alpha != MAXWEIGHT) {
1230 alpha /= 2;
1231 basic_subalpha (n, alpha, 0, PLUS, MINUS);
1232 } else {
1233 /* reverse sense of PLUS and MINUS */
1234 basic_minalpha (n, &new, &alpha, 1, PLUS, MINUS);
1235 if (alpha == MAXWEIGHT) {
1236 printf ("Basic dual change required, but no candidate edges\n");
1237 return (node *) NULL;
1238 }
1239 alpha /= 2;
1240 basic_subalpha (n, alpha, 1, PLUS, MINUS);
1241 }
1242 return new;
1243 }
1244
basic_minalpha(node * n,node ** new,int * alpha,int flip_plus_and_minus,int PLUS,int MINUS)1245 static void basic_minalpha (node *n, node **new, int *alpha,
1246 int flip_plus_and_minus, int PLUS, int MINUS)
1247 {
1248 int minv = MAXWEIGHT;
1249 int thisv;
1250 node *n1;
1251 edge *e;
1252 edgeptr *ep;
1253
1254 if ((n->label == PLUS && !flip_plus_and_minus) ||
1255 (n->label == MINUS && flip_plus_and_minus)) {
1256 for (ep = n->adj; ep; ep = ep->next) {
1257 e = ep->this;
1258 n1 = ep->other;
1259 if ((unsigned int) e->x != (unsigned int) TWO) {
1260 if (n1->label < PLUS) { /* n1 is unlabeled */
1261 thisv = e->weight - n->y - n1->y;
1262 thisv += thisv;
1263 if (thisv < minv)
1264 minv = thisv;
1265 } else if ((n1->label == PLUS && !flip_plus_and_minus) ||
1266 (n1->label == MINUS && flip_plus_and_minus)) {
1267 thisv = e->weight - n->y - n1->y;
1268 if (thisv < minv)
1269 minv = thisv;
1270 } else { /* n1 has a minus label */
1271 if (n1->parentedge == e)
1272 basic_minalpha (n1, new, alpha, flip_plus_and_minus,
1273 PLUS, MINUS);
1274 }
1275 } else if ((unsigned int) e->basic && n1->parentedge == e) {
1276 basic_minalpha (n1, new, alpha, flip_plus_and_minus,
1277 PLUS, MINUS);
1278 }
1279 }
1280 } else { /* MINUS case */
1281 for (ep = n->adj; ep; ep = ep->next) {
1282 e = ep->this;
1283 n1 = ep->other;
1284 if ((unsigned int) e->x == (unsigned int) TWO) {
1285 if (n1->label < PLUS) {
1286 thisv = e->z + e->z;
1287 if (thisv < minv)
1288 minv = thisv;
1289 } else if ((n1->label == PLUS && !flip_plus_and_minus) ||
1290 (n1->label == MINUS && flip_plus_and_minus)) {
1291 if (n1->parentedge == e)
1292 basic_minalpha (n1, new, alpha, flip_plus_and_minus,
1293 PLUS, MINUS);
1294 } else { /* n1 has a MINUS label */
1295 thisv = e->z;
1296 if (thisv < minv)
1297 minv = thisv;
1298 }
1299 } else if ((unsigned int) e->basic && n1->parentedge == e) {
1300 basic_minalpha (n1, new, alpha, flip_plus_and_minus,
1301 PLUS, MINUS);
1302 }
1303 }
1304 }
1305
1306 if (minv < *alpha) {
1307 *new = n;
1308 n->pnext = (node *) NULL;
1309 *alpha = minv;
1310 } else if (minv == *alpha) {
1311 n->pnext = *new;
1312 *new = n;
1313 }
1314 }
1315
basic_subalpha(node * n,int alpha,int flip_plus_and_minus,int PLUS,int MINUS)1316 static void basic_subalpha (node *n, int alpha, int flip_plus_and_minus,
1317 int PLUS, int MINUS)
1318 {
1319 edge *e;
1320 edgeptr *ep;
1321
1322 if ((n->label == PLUS && !flip_plus_and_minus) ||
1323 (n->label == MINUS && flip_plus_and_minus)) {
1324 n->y += alpha;
1325 for (ep = n->adj; ep; ep = ep->next) {
1326 e = ep->this;
1327 if ((unsigned int) e->x == (unsigned int) TWO)
1328 e->z += alpha;
1329 if ((unsigned int) e->basic) {
1330 if (ep->other->parentedge == e)
1331 basic_subalpha (ep->other, alpha, flip_plus_and_minus,
1332 PLUS, MINUS);
1333 }
1334 }
1335 } else { /* MINUS */
1336 n->y -= alpha;
1337 for (ep = n->adj; ep; ep = ep->next) {
1338 e = ep->this;
1339 if ((unsigned int) e->x == (unsigned int) TWO)
1340 e->z -= alpha;
1341 if ((unsigned int) e->basic) {
1342 if (ep->other->parentedge == e)
1343 basic_subalpha (ep->other, alpha, flip_plus_and_minus,
1344 PLUS, MINUS);
1345 }
1346 }
1347 }
1348 }
1349
1350
1351 /********** Price - Repair routines **********/
1352
1353
checkoutedge(graph * G,node * n1,node * n2,int * hit,CCdatagroup * dat)1354 static int checkoutedge (graph *G, node *n1, node *n2, int *hit,
1355 CCdatagroup *dat)
1356 {
1357 int w, wbar;
1358 edge *e;
1359
1360 *hit = 0;
1361 w = CCutil_dat_edgelen (n1->name, n2->name, dat);
1362 w += w;
1363 wbar = w - n1->y - n2->y;
1364 if (wbar < 0) {
1365 if ((e = findedge (n1, n2)) != (edge *) NULL) {
1366 if (e->z != -wbar) {
1367 printf ("Hmmm. edge (%d-%d) has z %d, wbar %d\n",
1368 e->ends[0]->name, e->ends[1]->name, e->z, wbar);
1369 }
1370 } else {
1371 if (addbadedge (G, n1, n2, w)) {
1372 fprintf (stderr, "addbadedge failed\n");
1373 return 1;
1374 }
1375 *hit = 1;
1376 }
1377 }
1378 return 0;
1379 }
1380
precheckoutedge(node * n1,node * n2,shortedge ** list,CCdatagroup * dat,CCptrworld * shortedge_world)1381 static int precheckoutedge (node *n1, node *n2, shortedge **list,
1382 CCdatagroup *dat, CCptrworld *shortedge_world)
1383 {
1384 int w, wbar;
1385 edge *e;
1386 shortedge *s;
1387
1388 w = CCutil_dat_edgelen (n1->name, n2->name, dat);
1389 w += w;
1390 wbar = w - n1->y - n2->y;
1391 if (wbar < 0) {
1392 if ((e = findedge (n1, n2)) != (edge *) NULL) {
1393 if (e->z != -wbar) {
1394 printf ("Hmmm. edge (%d-%d) has z %d, wbar %d\n",
1395 e->ends[0]->name, e->ends[1]->name, e->z, wbar);
1396 }
1397 } else {
1398 s = shortedgealloc (shortedge_world);
1399 s->ends[0] = n1;
1400 s->ends[1] = n2;
1401 s->next = *list;
1402 *list = s;
1403 return 1;
1404 }
1405 }
1406 return 0;
1407 }
1408
fixmatch(graph * G,int * radded,CCdatagroup * dat,CCrandstate * rstate)1409 static int fixmatch (graph *G, int *radded, CCdatagroup *dat,
1410 CCrandstate *rstate)
1411 {
1412 int datnorm;
1413
1414 CCutil_dat_getnorm (dat, &datnorm);
1415 if ((datnorm & CC_NORM_BITS) == CC_KD_NORM_TYPE)
1416 return kd_fixmatch (G, radded, dat, rstate);
1417 else if ((datnorm & CC_NORM_BITS) == CC_X_NORM_TYPE)
1418 return x_fixmatch (G, radded, dat);
1419 else
1420 return junk_fixmatch (G, radded, dat);
1421 }
1422
1423 #define NEAR_TRY_NUM 1 /* The number of nearest (wbar) neighbors */
1424
1425 #define PULL_DIVISOR 100 /* Do not pull more than ncount/PULL_DIVISOR. */
1426 #define PULL_UNIT 100 /* Pull if PULL_UNIT nodes will cut the spread. */
1427 #define PULL_CUT 2 /* A unit must cut at least spread/PULL_CUT. */
1428
kd_fixmatch(graph * G,int * radded,CCdatagroup * dat,CCrandstate * rstate)1429 static int kd_fixmatch (graph *G, int *radded, CCdatagroup *dat,
1430 CCrandstate *rstate)
1431 {
1432 int rval = 0;
1433 int i, j, added, totaladded = 0;
1434 int hit, passcount = 0;
1435 int maxy = -MAXWEIGHT;
1436 int miny = MAXWEIGHT;
1437 double *xcoord = (double *) NULL;
1438 double *ycoord = (double *) NULL;
1439 double *wcoord = (double *) NULL;
1440 CCdatagroup ldat;
1441 node *n;
1442 /*
1443 NEEDED WHEN RADIX SORT IS WORKING
1444 node *ysorted;
1445 */
1446 node **heavy, **light, **order = (node **) NULL;
1447 int top, bottom, nlight, nheavy = 0; /* nheavy should be set to 0 */
1448 int *invnames = (int *) NULL;
1449 double *lbnds = (double *) NULL, *gbnds = (double *) NULL;
1450 int datnorm;
1451
1452 *radded = 0;
1453
1454 CCutil_init_datagroup (&ldat);
1455
1456 CCutil_dat_getnorm (dat, &datnorm);
1457 if (CCutil_dat_setnorm (&ldat, datnorm)) {
1458 rval = 1;
1459 goto CLEANUP;
1460 }
1461
1462 for (n = G->nodelist; n; n = n->next) {
1463 if (n->y < miny)
1464 miny = n->y;
1465 if (n->y > maxy)
1466 maxy = n->y;
1467 }
1468 printf ("Node weight spread: (%d, %d)\n", miny, maxy);
1469 fflush (stdout);
1470
1471 /*
1472 THIS CODE CANNOT BE USED UNDER OS2 WITH CURRENT RADIX
1473 for (n = G->nodelist; n; n = n->next)
1474 n->pnext = n->next;
1475 ysorted = (node *) CCutil_linked_radixsort ((char *) G->nodelist,
1476 (char *) &(G->nodelist->pnext),
1477 (char *) &(G->nodelist->y), sizeof (int));
1478
1479 order = CC_SAFE_MALLOC (G->ncount, node *);
1480 if (!order) {
1481 rval = 1;
1482 goto CLEANUP;
1483 }
1484
1485 THIS IS THE CODE WHEN RADIXSORT WORKS WITH NEGATIVES
1486 for (i = 0, n = ysorted; n; i++, n = n->pnext) {
1487 order[i] = n;
1488 }
1489
1490 INSTEAD, THIS CODE WORKS WITH CURRENT RADIX ON THE RS6000
1491 for (n = ysorted; n; n = n->pnext)
1492 if (n->y < 0)
1493 break;
1494 for (i = 0; n; n = n->pnext, i++)
1495 order[i] = n;
1496 for (n = ysorted; n && n->y >= 0; n = n->pnext, i++)
1497 order[i] = n;
1498 */
1499
1500 {
1501 /* ONLY HERE UNTIL RADIX WORKS */
1502 int *y;
1503
1504 order = CC_SAFE_MALLOC (G->ncount, node *);
1505 if (!order) {
1506 rval = 1;
1507 goto CLEANUP;
1508 }
1509 y = CC_SAFE_MALLOC (G->ncount, int);
1510 if (!y) {
1511 rval = 1;
1512 goto CLEANUP;
1513 }
1514 for (i = 0; i < G->ncount; i++) {
1515 order[i] = G->nodenames[i];
1516 y[i] = G->nodenames[i]->y;
1517 }
1518 y_quicksort (order, y, 0, G->ncount - 1);
1519
1520 CC_FREE (y, int);
1521 }
1522
1523
1524 {
1525 int new, newspread, newtop, newbottom, spread;
1526
1527 newtop = top = -1;
1528 newbottom = bottom = G->ncount;
1529 nheavy = 0;
1530 spread = maxy - miny;
1531
1532 do {
1533 new = 0;
1534 while (new < PULL_UNIT && nheavy + new < G->ncount/ PULL_DIVISOR) {
1535 if (order[newbottom - 1]->y > -2 * order[newtop + 1]->y)
1536 newbottom--;
1537 else
1538 newtop++;
1539 new++;
1540 }
1541 newspread = order[newbottom - 1]->y - order[newtop + 1]->y;
1542 if (spread - newspread > spread / PULL_CUT) {
1543 spread = newspread;
1544 bottom = newbottom;
1545 top = newtop;
1546 nheavy += new;
1547 }
1548 } while (spread == newspread && nheavy < G->ncount/PULL_DIVISOR);
1549 }
1550
1551
1552 printf ("Truncated %d nodes to get spread: (%d, %d)\n",
1553 nheavy, order[top + 1]->y, order[bottom - 1]->y);
1554 fflush (stdout);
1555
1556
1557 if (nheavy) {
1558 heavy = order;
1559 light = order + nheavy;
1560 nlight = G->ncount - nheavy;
1561 /*
1562 THIS IS THE CODE WHEN RADIXSORT WORKS WITH NEGATIVES
1563 for (i = 0, n = ysorted; i <= top; i++, n = n->pnext) {
1564 heavy[i] = n;
1565 }
1566 for (i = 0; i < nlight; i++, n = n->pnext) {
1567 light[i] = n;
1568 }
1569 for (i = top + 1; i < nheavy; i++, n = n->pnext) {
1570 heavy[i] = n;
1571 }
1572 */
1573 {
1574 node **temporder = (node **) NULL;
1575 int k;
1576
1577 temporder = CC_SAFE_MALLOC (G->ncount, node *);
1578 if (!temporder) {
1579 rval = 1;
1580 goto CLEANUP;
1581 }
1582 for (i = 0; i < G->ncount; i++)
1583 temporder[i] = order[i];
1584
1585 for (i = 0, k = 0; i <= top; i++)
1586 heavy[i] = temporder[k++];
1587 for (i = 0; i < nlight; i++)
1588 light[i] = temporder[k++];
1589 for (i = top + 1; i < nheavy; i++)
1590 heavy[i] = temporder[k++];
1591 CC_FREE (temporder, node *);
1592 }
1593
1594 lbnds = CC_SAFE_MALLOC (nheavy, double);
1595 if (!lbnds) {
1596 rval = 1;
1597 goto CLEANUP;
1598 }
1599 gbnds = CC_SAFE_MALLOC (nheavy, double);
1600 if (!gbnds) {
1601 rval = 1;
1602 goto CLEANUP;
1603 }
1604 xcoord = CC_SAFE_MALLOC (nlight, double);
1605 if (!xcoord) {
1606 rval = 1;
1607 goto CLEANUP;
1608 }
1609 ldat.x = xcoord;
1610 ycoord = CC_SAFE_MALLOC (nlight, double);
1611 if (!ycoord) {
1612 rval = 1;
1613 goto CLEANUP;
1614 }
1615 ldat.y = ycoord;
1616 for (i = 0; i < nlight; i++) {
1617 xcoord[i] = dat->x[light[i]->name];
1618 ycoord[i] = dat->y[light[i]->name];
1619 }
1620 } else {
1621 nlight = G->ncount;
1622 light = G->nodenames;
1623 heavy = (node **) NULL;
1624 xcoord = dat->x;
1625 ycoord = dat->y;
1626 ldat.x = xcoord;
1627 ldat.y = ycoord;
1628 CC_FREE (order, node *);
1629 }
1630
1631 wcoord = CC_SAFE_MALLOC (nlight, double);
1632 if (!wcoord) {
1633 rval = 1;
1634 goto CLEANUP;
1635 }
1636 invnames = CC_SAFE_MALLOC (G->ncount, int);
1637 if (!invnames) {
1638 rval = 1;
1639 goto CLEANUP;
1640 }
1641 for (i = 0; i < nlight; i++)
1642 invnames[light[i]->name] = i;
1643 for (i = 0; i < nheavy; i++)
1644 invnames[heavy[i]->name] = -i;
1645
1646
1647 do {
1648 int nodeschecked = 0;
1649 int saver = 0;
1650 int list[NEAR_TRY_NUM];
1651 shortedge *s, *snext, *slist = (shortedge *) NULL;
1652 CCkdtree localkt;
1653 added = 0;
1654 maxy = -MAXWEIGHT;
1655 for (i = 0; i < nlight; i++) {
1656 if (light[i]->y > maxy)
1657 maxy = light[i]->y;
1658 }
1659 for (i = 0; i < nlight; i++)
1660 wcoord[i] = ((double) (maxy - light[i]->y)) * 0.5;
1661 if (CCkdtree_build (&localkt, nlight, &ldat, wcoord, rstate)) {
1662 fprintf (stderr, "Unable to build CCkdtree\n");
1663 rval = 1;
1664 goto CLEANUP;
1665 }
1666 for (i = 0; i < nheavy; i++) { /* 1.0 for floats */
1667 lbnds[i] = dat->x[heavy[i]->name] -
1668 (((double) (heavy[i]->y)) * 0.5) - 1.0;
1669 gbnds[i] = dat->x[heavy[i]->name] +
1670 (((double) (heavy[i]->y)) * 0.5) + 1.0;
1671 }
1672
1673 for (i = 0; i < nlight; i++) {
1674 if (light[i]->label != -1) {
1675 edgeptr *ep;
1676 hit = 0;
1677 for (ep = light[i]->adj; ep; ep = ep->next) {
1678 if ((unsigned int) ep->this->x != (unsigned int) ZERO &&
1679 invnames[ep->other->name] >= 0)
1680 CCkdtree_delete (&localkt, invnames[ep->other->name]);
1681 }
1682 nodeschecked++;
1683 if (CCkdtree_node_k_nearest (&localkt, nlight, i, NEAR_TRY_NUM,
1684 &ldat, wcoord, list, rstate)) {
1685 fprintf (stderr, "node nearest failed\n");
1686 CCkdtree_free (&localkt);
1687 rval = 1;
1688 goto CLEANUP;
1689 }
1690 for (j = NEAR_TRY_NUM - 1; j >= 0; j--) {
1691 if (list[j] != -1)
1692 hit += precheckoutedge (light[i], light[list[j]],
1693 &slist, dat, &G->shortedge_world);
1694 }
1695 for (j = 0; j < nheavy; j++) {
1696 if (heavy[j]->label == -1) {
1697 if (xcoord[i] +
1698 (((double) light[i]->y) * 0.5) > lbnds[j] &&
1699 xcoord[i] -
1700 (((double) light[i]->y) * 0.5) < gbnds[j]) {
1701 hit += precheckoutedge (light[i], heavy[j],
1702 &slist, dat, &G->shortedge_world);
1703 } else {
1704 saver++;
1705 }
1706 }
1707 }
1708 added += hit;
1709 if (hit == 0)
1710 light[i]->label = -1;
1711 for (ep = light[i]->adj; ep; ep = ep->next) {
1712 if ((unsigned int) ep->this->x != (unsigned int) ZERO &&
1713 invnames[ep->other->name] >= 0)
1714 CCkdtree_undelete (&localkt,
1715 invnames[ep->other->name]);
1716 }
1717 }
1718 }
1719 for (j = 0; j < nheavy; j++) {
1720 if (heavy[j]->label != -1) {
1721 hit = 0;
1722 nodeschecked++;
1723 for (i = 0; i < nlight; i++) {
1724 if (xcoord[i]+(((double) light[i]->y) * 0.5) > lbnds[j] &&
1725 xcoord[i]-(((double) light[i]->y) * 0.5) < gbnds[j]) {
1726 hit += precheckoutedge (light[i], heavy[j], &slist,
1727 dat, &G->shortedge_world);
1728 } else {
1729 saver++;
1730 }
1731 }
1732 for (i = 0; i < j; i++) {
1733 if (dat->x[heavy[i]->name] +
1734 (((double) heavy[i]->y) * 0.5) > lbnds[j] &&
1735 dat->x[heavy[i]->name] -
1736 (((double) heavy[i]->y) * 0.5) < gbnds[j]) {
1737 hit += precheckoutedge (heavy[i], heavy[j], &slist,
1738 dat, &G->shortedge_world);
1739 }
1740 }
1741 for (i = j + 1; i < nheavy; i++) {
1742 if (heavy[i]->label == -1) {
1743 if (dat->x[heavy[i]->name] +
1744 (((double) heavy[i]->y) * 0.5) > lbnds[j] &&
1745 dat->x[heavy[i]->name] -
1746 (((double) heavy[i]->y) * 0.5) < gbnds[j]) {
1747 hit += precheckoutedge (heavy[i], heavy[j],
1748 &slist, dat, &G->shortedge_world);
1749 }
1750 }
1751 }
1752 added += hit;
1753 if (hit == 0)
1754 heavy[j]->label = -1;
1755 }
1756 }
1757
1758 printf ("Need to check %d edges (saved %d checks)\n", added, saver);
1759 fflush (stdout);
1760 CCkdtree_free (&localkt);
1761
1762 added = 0;
1763 for (s = slist; s; s = snext) {
1764 snext = s->next;
1765 if (checkoutedge (G, s->ends[0], s->ends[1], &hit, dat)) {
1766 fprintf (stderr, "checkoutedge failed\n");
1767 rval = 1;
1768 goto CLEANUP;
1769 }
1770 added += hit;
1771 shortedgefree (&G->shortedge_world, s);
1772 }
1773 totaladded += added;
1774 printf ("Pass %d: %d edges added (%d total), %d nodes checked\n",
1775 passcount++, added, totaladded, nodeschecked);
1776 fflush (stdout);
1777 } while (added);
1778 *radded = totaladded;
1779
1780 CLEANUP:
1781
1782 CC_IFFREE (invnames, int);
1783 CC_IFFREE (wcoord, double);
1784 CC_IFFREE (order, node *);
1785 if (nheavy) {
1786 CC_IFFREE (xcoord, double);
1787 CC_IFFREE (ycoord, double);
1788 CC_IFFREE (lbnds, double);
1789 CC_IFFREE (gbnds, double);
1790 }
1791 return rval;
1792 }
1793
initlist(graph * G,node * head,node * tail,node * head2,node * tail2,CCdatagroup * dat)1794 static void initlist (graph *G, node *head, node *tail, node *head2,
1795 node *tail2, CCdatagroup *dat)
1796 {
1797 node *n, *p;
1798 int bound;
1799 int datnorm;
1800 double *xcoord = dat->x;
1801 double scale;
1802
1803 CCutil_dat_getnorm (dat, &datnorm);
1804 if (datnorm == CC_GEOGRAPHIC) scale = CC_GEOGRAPHIC_SCALE;
1805 else if (datnorm == CC_GEOM) scale = CC_GEOM;
1806 else if (datnorm == CC_ATT) scale = CC_ATT_SCALE;
1807 else scale = 1.0;
1808
1809 head->sort.order = -MAXWEIGHT;
1810 tail->sort.order = MAXWEIGHT;
1811 head->sort.next = tail;
1812 head->sort.prev = (node **) NULL;
1813 tail->sort.next = (node *) NULL;
1814 tail->sort.prev = &(head->sort.next);
1815 head2->sort.order = MAXWEIGHT;
1816 tail2->sort.order = -MAXWEIGHT;
1817 head2->sort.next = tail2;
1818 head2->sort.prev = (node **) NULL;
1819 tail2->sort.next = (node *) NULL;
1820 tail2->sort.prev = &(head2->sort.next);
1821
1822 for (n = G->nodelist; n; n = n->next) {
1823 bound = (2 * ((int) (scale * xcoord[n->name]))) - n->y;
1824 for (p = head->sort.next; p->sort.order < bound; p = p->sort.next);
1825 n->sort.order = bound;
1826 n->sort.next = p;
1827 n->sort.prev = p->sort.prev;
1828 *(n->sort.prev) = n;
1829 p->sort.prev = &(n->sort.next);
1830 }
1831 }
1832
x_fixmatch(graph * G,int * radded,CCdatagroup * dat)1833 static int x_fixmatch (graph *G, int *radded, CCdatagroup *dat)
1834 {
1835 node *n1, *n2;
1836 int i;
1837 int added, hit;
1838 int nodeschecked;
1839 int edgeschecked;
1840 int totaladded = 0;
1841 node high_fakehead, high_faketail, low_fakehead, low_faketail;
1842 int bound;
1843 double *xcoord = dat->x;
1844 double scale;
1845 int datnorm;
1846
1847 CCutil_dat_getnorm (dat, &datnorm);
1848 if ((datnorm & CC_NORM_BITS) != CC_X_NORM_TYPE &&
1849 (datnorm & CC_NORM_BITS) != CC_KD_NORM_TYPE) {
1850 fprintf (stderr, "Cannot run x_fixmatch with norm %d\n", datnorm);
1851 return 1;
1852 }
1853
1854 if (datnorm == CC_GEOGRAPHIC) scale = CC_GEOGRAPHIC_SCALE;
1855 else if (datnorm == CC_GEOM) scale = CC_GEOM;
1856 else if (datnorm == CC_ATT) scale = CC_ATT_SCALE;
1857 else scale = 1.0;
1858
1859 initlist (G, &high_fakehead, &high_faketail, &low_fakehead, &low_faketail,
1860 dat);
1861
1862 do {
1863 added = 0;
1864 nodeschecked = 0;
1865 edgeschecked = 0;
1866 for (i = 0; i < G->ncount; i++) {
1867 n1 = G->nodenames[i];
1868 *(n1->sort.prev) = n1->sort.next;
1869 n1->sort.next->sort.prev = n1->sort.prev;
1870 if (n1->label != -1) {
1871 nodeschecked++;
1872 n1->label = -1;
1873 bound = (2 * ((int) (scale * xcoord[n1->name]))) + n1->y + 3;
1874 /* Need the +3 to handle floating point data */
1875 for (n2 = high_fakehead.sort.next; n2->sort.order < bound;
1876 n2 = n2->sort.next) {
1877 edgeschecked++;
1878 if (checkoutedge (G, n1, n2, &hit, dat)) {
1879 fprintf (stderr, "checkoutedge failed\n");
1880 return 1;
1881 }
1882 added += hit;
1883 }
1884 bound = (2 * ((int) (scale * xcoord[n1->name]))) - n1->y - 3;
1885 for (n2 = low_fakehead.sort.next; n2->sort.order > bound;
1886 n2 = n2->sort.next) {
1887 edgeschecked++;
1888 if (checkoutedge (G, n1, n2, &hit, dat)) {
1889 fprintf (stderr, "checkoutedge failed\n");
1890 return 1;
1891 }
1892 added += hit;
1893 }
1894 }
1895 bound = (2 * ((int) (scale * xcoord[n1->name]))) + n1->y;
1896 for (n2 = low_fakehead.sort.next; n2->sort.order > bound;
1897 n2 = n2->sort.next);
1898 n1->sort.order = bound;
1899 n1->sort.next = n2;
1900 n1->sort.prev = n2->sort.prev;
1901 *(n1->sort.prev) = n1;
1902 n2->sort.prev = &n1->sort.next;
1903 }
1904 totaladded += added;
1905 printf ("Forward pass completed, %d nodes checked, %d edges checked\n",
1906 nodeschecked, edgeschecked);
1907 printf (" %d edges added, total %d edges added\n",
1908 added, totaladded);
1909 if (added == 0)
1910 break;
1911 added = 0;
1912 nodeschecked = 0;
1913 edgeschecked = 0;
1914 for (i = G->ncount - 1; i >= 0; i--) {
1915 n1 = G->nodenames[i];
1916 *(n1->sort.prev) = n1->sort.next;
1917 n1->sort.next->sort.prev = n1->sort.prev;
1918 if (n1->label != -1) {
1919 nodeschecked++;
1920 n1->label = -1;
1921 bound = (2 * ((int) (scale * xcoord[n1->name]))) + n1->y + 3;
1922 for (n2 = high_fakehead.sort.next; n2->sort.order < bound;
1923 n2 = n2->sort.next) {
1924 edgeschecked++;
1925 if (checkoutedge (G, n1, n2, &hit, dat)) {
1926 fprintf (stderr, "checkoutedge failed\n");
1927 return 1;
1928 }
1929 added += hit;
1930 }
1931 bound = (2 * ((int) (scale * xcoord[n1->name]))) - n1->y - 3;
1932 for (n2 = low_fakehead.sort.next; n2->sort.order > bound;
1933 n2 = n2->sort.next) {
1934 edgeschecked++;
1935 if (checkoutedge (G, n1, n2, &hit, dat)) {
1936 fprintf (stderr, "checkoutedge failed\n");
1937 return 1;
1938 }
1939 added += hit;
1940 }
1941 }
1942 bound = (2 * ((int) (scale * xcoord[n1->name]))) - n1->y;
1943 for (n2 = high_fakehead.sort.next; n2->sort.order < bound;
1944 n2 = n2->sort.next);
1945 n1->sort.order = bound;
1946 n1->sort.next = n2;
1947 n1->sort.prev = n2->sort.prev;
1948 *(n1->sort.prev) = n1;
1949 n2->sort.prev = &n1->sort.next;
1950 }
1951 totaladded += added;
1952 printf ("Backward pass completed, %d nodes checked, %d edges checked\n",
1953 nodeschecked, edgeschecked);
1954 printf (" %d edges added, total %d edges added\n",
1955 added, totaladded);
1956 } while (added);
1957 *radded = totaladded;
1958 return 0;
1959 }
1960
junk_fixmatch(graph * G,int * radded,CCdatagroup * dat)1961 static int junk_fixmatch (graph *G, int *radded, CCdatagroup *dat)
1962 {
1963 int added, hit, totaladded = 0;
1964 node *n1, *n2;
1965
1966 *radded = 0;
1967
1968 do {
1969 added = 0;
1970 for (n1 = G->nodelist; n1; n1 = n1->next) {
1971 for (n2 = n1->next; n2; n2 = n2->next) {
1972 if (checkoutedge (G, n1, n2, &hit, dat)) {
1973 fprintf (stderr, "checkoutedge failed\n");
1974 return 1;
1975 }
1976 added += hit;
1977 }
1978 }
1979 totaladded += added;
1980 printf ("Pass completed: %d edges added, total %d edges added\n",
1981 added, totaladded);
1982 fflush (stdout);
1983 } while (added);
1984
1985 *radded = totaladded;
1986 return 0;
1987 }
1988
addbadedge(graph * G,node * n1,node * n2,int w)1989 static int addbadedge (graph *G, node *n1, node *n2, int w)
1990 {
1991 int wbar = -(w - n1->y - n2->y);
1992 edgeptr *ep;
1993 edge *newe;
1994 node *other1 = 0, *other2 = 0;
1995
1996 for (ep = n1->adj; ep; ep = ep->next) {
1997 switch ((unsigned int) ep->this->x) {
1998 case (unsigned int) ONE:
1999 flipcycle (n1, ep->this, ZERO);
2000 n1->matchcnt = (unsigned int) n1->matchcnt - 1;
2001 break;
2002 case (unsigned int) TWO:
2003 if ((ep->this->z -= wbar) < 0) {
2004 ep->this->z = 0;
2005 ep->this->x = ZERO;
2006 n1->matchcnt = (unsigned int) n1->matchcnt - 1;
2007 ep->other->matchcnt = (unsigned int) ep->other->matchcnt - 1;
2008 other2 = other1;
2009 other1 = ep->other;
2010 }
2011 break;
2012 }
2013 }
2014 n1->y -= wbar;
2015 newe = newedge (G, n1, n2);
2016 if (!newe)
2017 return 1;
2018 newe->weight = w;
2019 newe->z = 0;
2020 newe->x = ZERO;
2021 newe->pnext = (edge *) NULL;
2022 while ((unsigned int) n1->matchcnt != (unsigned int) TWO)
2023 augment (G, n1);
2024 if (other1) {
2025 while ((unsigned int) other1->matchcnt != (unsigned int) TWO)
2026 augment (G, other1);
2027 }
2028 if (other2) {
2029 while ((unsigned int) other2->matchcnt != (unsigned int) TWO)
2030 augment (G, other2);
2031 }
2032 return 0;
2033 }
2034
y_quicksort(node ** list,int * y,int l,int u)2035 static void y_quicksort (node **list, int *y, int l, int u)
2036 {
2037 int i, j, t;
2038 int itemp;
2039 node *ntemp;
2040
2041 if (l >= u)
2042 return;
2043
2044 CC_SWAP (y[l], y[(l+u)/2], itemp);
2045 CC_SWAP (list[l], list[(l+u)/2], ntemp);
2046
2047 i = l;
2048 j = u + 1;
2049 t = y[l];
2050
2051 while (1) {
2052 do i++; while (i <= u && y[i] < t);
2053 do j--; while (y[j] > t);
2054 if (j < i) break;
2055 CC_SWAP (y[i], y[j], itemp);
2056 CC_SWAP (list[i], list[j], ntemp);
2057 }
2058 CC_SWAP (y[l], y[j], itemp);
2059 CC_SWAP (list[l], list[j], ntemp);
2060 y_quicksort (list, y, l, j - 1);
2061 y_quicksort (list, y, i, u);
2062 }
2063