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 /*    EXPORTED FUNCTIONS:                                                   */
19 /*                                                                          */
20 /*  CCtiny_bnb_msp (int nnodes, int nedges, int *elist, int *weight,        */
21 /*      int depot, int *lbound, int *ubound, double *objlimit,              */
22 /*      int objdir, double *objval, int *xsol, int searchlimit)             */
23 /*    MISSING                                                               */
24 /*                                                                          */
25 /****************************************************************************/
26 
27 #include "machdefs.h"
28 #include "util.h"
29 #include "tinytsp.h"
30 
31 /* the node which can have degree > 2 */
32 #define OUTSIDE (0)
33 #define OUTSIDEEDGE(n0,n1) ((n0) == OUTSIDE || (n1) == OUTSIDE)
34 
35 #define BIGDOUBLE (1e30)
36 #undef  USE_DUMBO
37 
38 #ifdef USE_DUMBO
39 #define TINYBRANCH_WEIGHT 10.0
40 #define TINYBRANCH_VAL(x,y) (((x)<(y)) ? TINYBRANCH_VAL2((y),(x)) \
41                                        : TINYBRANCH_VAL2((x),(y)))
42 #define TINYBRANCH_VAL2(x,y) (((x)+(y)*TINYBRANCH_WEIGHT) / \
43                               (TINYBRANCH_WEIGHT+1.0))
44 #endif
45 
46 #define SELECT_INFEASIBLE 1
47 #define SELECT_TOUR 2
48 
49 /* #define DEBUG 1 */
50 #undef TIMINGS
51 
52 #define SEARCHLIMIT (1<<30)
53 
54 typedef struct node {
55     int availdeg;
56     int deg;
57 /* OUTSIDE's pathend == OUTSIDE */
58     int pathend;
59     struct edge *adjvec;
60     int adjlist;
61 } node;
62 
63 typedef struct edge {
64     int next;
65     int prev;
66     int weight;
67     int lo;
68     int hi;
69     int value;
70     int bestvalue;
71 } edge;
72 
73 typedef struct histent {
74     int from;
75     int to;
76     int pfrom;
77     int pto;
78 } histent;
79 
80 typedef struct tspsearch {
81     int nnodes;
82     node *nodes;
83     histent *searchhist;
84     int histtop;
85     double currentval;
86     double bestval;
87     int searchcount;
88     int searchlimit;
89 } tspsearch;
90 
91 static int
92     checkout_node (tspsearch *s, int i),
93     build_paths (tspsearch *s),
94     set_edge(tspsearch *s, int n0, int n1, int v),
95     select_edge (tspsearch *s, int *n0, int *n1);
96 
97 static void
98     build_adj (tspsearch *s),
99     search (tspsearch *s),
100     unset_edge(tspsearch *s, histent *h),
101     savetour (tspsearch *s),
102     unroll_stack (tspsearch *s, int smark);
103 
104 static double
105     lowerbound (tspsearch *s);
106 
CCtiny_bnb_msp(int nnodes,int nedges,int * elist,int * weight,int depot,int * lbound,int * ubound,double * objlimit,int objdir,double * objval,int * xsol,int searchlimit)107 int CCtiny_bnb_msp (int nnodes, int nedges, int *elist, int *weight, int depot,
108         int *lbound, int *ubound, double *objlimit, int objdir,
109         double *objval, int *xsol, int searchlimit)
110 {
111     tspsearch s;
112     int i, j;
113     int rval;
114     int n0, n1;
115     double initlimit;
116 #ifdef TIMINGS
117     double sz = CCutil_zeit();
118 #endif
119 
120     if (depot != OUTSIDE) {
121         fprintf (stderr, "bnbtsp doesn't know how to handle depot %d\n",
122                  depot);
123         return CC_TINYTSP_ERROR;
124     }
125 
126     rval = CC_TINYTSP_ERROR; /* Failures in this part are from out of memory */
127     s.nnodes = nnodes;
128     s.nodes = (node *) NULL;
129     s.searchhist = (histent *) NULL;
130     s.histtop = 0;
131     s.currentval = 0.0;
132     if (objlimit) {
133         initlimit = objdir * (*objlimit);
134     } else {
135         initlimit = BIGDOUBLE;
136     }
137     s.bestval = initlimit;
138 
139     if (searchlimit < 0) {
140         s.searchlimit = SEARCHLIMIT;
141     } else {
142         s.searchlimit = searchlimit;
143     }
144 
145     s.nodes = CC_SAFE_MALLOC (nnodes, node);
146     if (!s.nodes) goto CLEANUP;
147     for (i=0; i<nnodes; i++) {
148         s.nodes[i].adjvec = (edge *) NULL;
149     }
150     s.searchhist = CC_SAFE_MALLOC (nnodes * nnodes, histent);
151     if (!s.searchhist) goto CLEANUP;
152     for (i=0; i<nnodes; i++) {
153         s.nodes[i].adjvec = CC_SAFE_MALLOC (nnodes, edge);
154         if (!s.nodes[i].adjvec) goto CLEANUP;
155     }
156 
157     for (i=0; i<nnodes; i++) {
158         s.nodes[i].pathend = i;
159         for (j=0; j<nnodes; j++) {
160             s.nodes[i].adjvec[j].weight = 0;
161             s.nodes[i].adjvec[j].lo = 0;
162             s.nodes[i].adjvec[j].hi = 0;
163             s.nodes[i].adjvec[j].value = -1;
164             s.nodes[i].adjvec[j].bestvalue = -1;
165         }
166     }
167 
168     for (i=0; i<nedges; i++) {
169         n0 = elist[2*i];
170         n1 = elist[2*i+1];
171         s.nodes[n0].adjvec[n1].weight = objdir * weight[i];
172         s.nodes[n1].adjvec[n0].weight = objdir * weight[i];
173         s.nodes[n0].adjvec[n1].lo = lbound[i];
174         s.nodes[n1].adjvec[n0].lo = lbound[i];
175         s.nodes[n0].adjvec[n1].hi = ubound[i];
176         s.nodes[n1].adjvec[n0].hi = ubound[i];
177     }
178 
179     build_adj (&s);
180 
181 #ifdef DEBUG
182     printf ("Performing initial checks\n");
183 #endif
184     if (build_paths (&s)) {
185         rval = CC_TINYTSP_INFEASIBLE;
186         goto CLEANUP;
187     }
188 
189     for (i=0; i<nnodes; i++) {
190         if (checkout_node (&s, i)) {
191             rval = CC_TINYTSP_INFEASIBLE;
192             goto CLEANUP;
193         }
194         if (s.nodes[i].deg == 1 && i != OUTSIDE &&
195             s.nodes[i].pathend != OUTSIDE &&
196             s.nodes[i].adjvec[s.nodes[i].pathend].value == -1) {
197 #ifdef DEBUG
198 #if DEBUG>1
199             printf ("node %d pathend, forcing %d to 0\n",
200                     i,s.nodes[i].pathend);
201 #endif
202 #endif
203             if (set_edge (&s, i, s.nodes[i].pathend, 0)) {
204                 rval = CC_TINYTSP_INFEASIBLE;
205                 goto CLEANUP;
206             }
207         }
208     }
209 
210 #ifdef DEBUG
211     printf ("Beginning search\n");
212 #endif
213     s.searchcount = 0;
214     search (&s);
215 #ifdef DEBUG
216     printf ("Search finished\n");
217 #endif
218 
219     if (s.searchcount >= s.searchlimit) {
220 #ifdef DEBUG
221         fprintf (stderr, "search node limit %d exceeded\n", s.searchlimit);
222         fprintf (stderr, "SNLE obj:\n");
223         fprintf (stderr, "%d %d\n", nnodes, nedges);
224         for (i=0; i<nedges; i++) {
225             fprintf (stderr, "%d %d %d %d %d\n", elist[2*i], elist[2*i+1],
226                      lbound[i], ubound[i], objdir * weight[i]);
227         }
228 #endif
229 #ifdef TIMINGS
230         printf ("BNB Search Limit Nodes: %d   Time %.3f\n",
231                 s.searchcount, CCutil_zeit() - sz); fflush (stdout);
232 #endif
233         rval = CC_TINYTSP_SEARCHLIMITEXCEEDED;
234         goto CLEANUP;
235     }
236 
237 #ifdef TIMINGS
238     printf ("BNB Search in Nodes: %d   Time %.3f\n",
239             s.searchcount, CCutil_zeit() - sz); fflush (stdout);
240 #endif
241 #ifdef DEBUG
242     printf ("bnbtsp search finished in %d nodes\n", s.searchcount);
243     fflush (stdout);
244 #endif
245 
246     if (objval) *objval = objdir * s.bestval;
247     if (s.bestval >= initlimit) {
248         rval = CC_TINYTSP_INFEASIBLE;
249         goto CLEANUP;
250     }
251 
252     if (xsol) {
253         for (i=0; i<nedges; i++) {
254             xsol[i] = s.nodes[elist[2*i]].adjvec[elist[2*i+1]].bestvalue;
255         }
256     }
257     rval = 0;
258 
259   CLEANUP:
260     if (s.nodes) {
261         for (i=0; i<nnodes; i++) {
262             CC_IFFREE (s.nodes[i].adjvec, edge);
263         }
264         CC_FREE (s.nodes, node);
265     }
266     CC_IFFREE (s.searchhist, histent);
267     return rval;
268 }
269 
checkout_node(tspsearch * s,int i)270 static int checkout_node (tspsearch *s, int i)
271 {
272     node *nodes = s->nodes;
273     int j;
274 
275     if (nodes[i].availdeg < 2 || (i != OUTSIDE && nodes[i].deg > 2)) return 1;
276     if (nodes[i].deg < 2) {
277         if (nodes[i].availdeg == 2) {
278             while ((j = nodes[i].adjlist) != -1) {
279 #ifdef DEBUG
280 #if DEBUG>1
281                 printf ("node %d avail 2, forcing %d to %d\n",
282                         i,j,nodes[i].adjvec[j].hi);
283 #endif
284 #endif
285                 if (set_edge (s, i, j, nodes[i].adjvec[j].hi)) return 1;
286             }
287         }
288     }
289     j = nodes[i].adjlist;
290     if (j != -1) {
291         if (nodes[i].adjvec[j].next == -1) {
292             if (i == OUTSIDE) {
293                 if (nodes[i].deg & 1) {
294 #ifdef DEBUG
295 #if DEBUG>1
296                     printf ("node %d single adj + odd deg, forcing %d to %d\n",
297                              i,j,nodes[i].adjvec[j].lo + 1);
298 #endif
299 #endif
300                     if (set_edge (s, i, j, nodes[i].adjvec[j].lo + 1))
301                         return 1;
302                 } else if (nodes[i].adjvec[j].hi - nodes[i].adjvec[j].lo
303                            == 1) {
304 #ifdef DEBUG
305 #if DEBUG>1
306                     printf ("node %d single adj + even deg, forcing %d to %d\n",
307                              i,j,nodes[i].adjvec[j].lo);
308 #endif
309 #endif
310                     if (set_edge (s, i, j, nodes[i].adjvec[j].lo)) return 1;
311                 }
312             } else {
313 #ifdef DEBUG
314 #if DEBUG>1
315                 printf ("node %d single adj, forcing %d to %d\n",
316                         i,j,2-nodes[i].deg);
317 #endif
318 #endif
319 
320                 if (set_edge (s, i, j, 2 - nodes[i].deg)) return 1;
321             }
322         }
323     }
324     if (i != OUTSIDE && nodes[i].deg == 2) {
325         while ((j = nodes[i].adjlist) != -1) {
326 #ifdef DEBUG
327 #if DEBUG>1
328             printf ("node %d deg 2, forcing %d to %d\n",
329                     i,j,nodes[i].adjvec[j].lo);
330 #endif
331 #endif
332             if (set_edge (s, i, j, nodes[i].adjvec[j].lo)) return 1;
333         }
334     }
335 
336     return 0;
337 }
338 
build_paths(tspsearch * s)339 static int build_paths (tspsearch *s)
340 {
341     int i, j;
342     int pi, pj;
343     int nnodes = s->nnodes;
344     node *nodes = s->nodes;
345 
346     for (i=0; i<nnodes; i++) {
347         nodes[i].pathend = i;
348     }
349     for (i=0; i<nnodes; i++) {
350         for (j=i+1; j<nnodes; j++) {
351             if (s->nodes[i].adjvec[j].lo > 0) {
352                 if (nodes[i].pathend == j) return 1;
353                 pi = nodes[i].pathend;
354                 pj = nodes[j].pathend;
355                 if (i != OUTSIDE) {
356                     nodes[i].pathend = pj;
357                     if (pi != OUTSIDE && pi != i) nodes[pi].pathend = pj;
358                 }
359                 if (j != OUTSIDE) {
360                     nodes[j].pathend = pi;
361                     if (pj != OUTSIDE && pj != j) nodes[pj].pathend = pi;
362                 }
363             }
364         }
365     }
366     return 0;
367 }
368 
build_adj(tspsearch * s)369 static void build_adj (tspsearch *s)
370 {
371     int i, j;
372     int p;
373     int nnodes = s->nnodes;
374     node *nodes = s->nodes;
375 
376     for (i=0; i<nnodes; i++) {
377         p = -1;
378         nodes[i].adjlist = -1;
379         nodes[i].deg = 0;
380         nodes[i].availdeg = 0;
381         for (j=0; j<nnodes; j++) {
382             nodes[i].deg += nodes[i].adjvec[j].lo;
383             nodes[i].availdeg += nodes[i].adjvec[j].hi;
384             if (nodes[i].adjvec[j].lo < nodes[i].adjvec[j].hi) {
385                 nodes[i].adjvec[j].prev = p;
386                 nodes[i].adjvec[j].next = -1;
387                 if (p == -1) {
388                     nodes[i].adjlist = j;
389                 } else {
390                     nodes[i].adjvec[p].next = j;
391                 }
392                 p = j;
393             } else {
394                 nodes[i].adjvec[j].value = nodes[i].adjvec[j].lo;
395                 if (i <= j) { /* only add in each edge once */
396                     s->currentval += nodes[i].adjvec[j].value *
397                         nodes[i].adjvec[j].weight;
398                 }
399             }
400         }
401     }
402 }
403 
404 /* minimization */
search(tspsearch * s)405 static void search (tspsearch *s)
406 {
407     int n0, n1;
408     int smark;
409     int i, lo, hi;
410     int val;
411 
412     s->searchcount++;
413     if (s->searchcount >= s->searchlimit) return;
414 
415     if (s->currentval + lowerbound(s) >= s->bestval) {
416 #ifdef DEBUG
417         printf ("search, currentval %.0f lowerbound %.0f, bestval %.0f, returning\n",
418                 s->currentval, lowerbound(s), s->bestval);
419 #endif
420         return;
421     }
422     val = select_edge (s, &n0, &n1);
423     if (val == SELECT_TOUR) {
424         if (s->currentval + lowerbound(s) < s->bestval) {
425             savetour(s);
426         }
427         return;
428     } else if (val == SELECT_INFEASIBLE) {
429         return;
430     }
431 
432     smark = s->histtop;
433     lo = s->nodes[n0].adjvec[n1].lo;
434     hi = s->nodes[n0].adjvec[n1].hi;
435     if (s->nodes[n0].adjvec[n1].weight > 0) {
436         for (i = lo; i<=hi; i++) {
437 #ifdef DEBUG
438             printf ("branching %d-%d to %d (lo %d hi %d)\n",
439                     n0,n1,i,lo,hi);
440 #endif
441             if (!set_edge (s, n0, n1, i)) {
442                 search (s);
443             }
444             unroll_stack (s, smark);
445         }
446     } else {
447         for (i = hi; i>=lo; i--) {
448 #ifdef DEBUG
449             printf ("branching %d-%d to %d (lo %d hi %d)\n",
450                     n0,n1,i,lo,hi);
451 #endif
452             if (!set_edge (s, n0, n1, i)) {
453                 search (s);
454             }
455             unroll_stack (s, smark);
456         }
457     }
458 #ifdef DEBUG
459     printf ("done branching on %d-%d (lo %d hi %d)\n",
460             n0,n1,lo,hi);
461 #endif
462     return;
463 }
464 
465 /* 1 means contradiction found */
set_edge(tspsearch * s,int n0,int n1,int v)466 static int set_edge(tspsearch *s, int n0, int n1, int v)
467 {
468     int p0, p1;
469     edge *e;
470     node *nodes = s->nodes;
471     int lo, hi;
472 
473 #ifdef DEBUG
474 #if DEBUG>2
475     printf ("set_edge %d-%d = %d\n",n0,n1,v);
476 #endif
477 #endif
478     if (n1 == OUTSIDE) {
479         n1 = n0;
480         n0 = OUTSIDE;
481     }
482     if (nodes[n0].adjvec[n1].value != -1) {
483         if (nodes[n0].adjvec[n1].value == v) return 0;
484         else return 1;
485     }
486     lo = nodes[n0].adjvec[n1].lo;
487     hi = nodes[n0].adjvec[n1].hi;
488     if (n0 != OUTSIDE && nodes[n0].deg - lo + v > 2) return 1;
489     if (nodes[n1].deg - lo + v > 2) return 1;
490     if (nodes[n0].availdeg - hi + v < 2) return 1;
491     if (nodes[n1].availdeg - hi + v < 2) return 1;
492     p0 = nodes[n0].pathend;
493     p1 = nodes[n1].pathend;
494     if (v >= 1 && n0 != OUTSIDE && n1 != OUTSIDE && p0 == n1) return 1;
495 
496     nodes[n0].deg += v - lo;
497     nodes[n1].deg += v - lo;
498     nodes[n0].availdeg += v - hi;
499     nodes[n1].availdeg += v - hi;
500 
501     if (v >= 1) { /* The only situation in which pathends matter */
502         if (p0 != OUTSIDE) nodes[p0].pathend = p1;
503         if (p1 != OUTSIDE) nodes[p1].pathend = p0;
504     }
505 
506     e = &nodes[n0].adjvec[n1];
507     if (e->next != -1) nodes[n0].adjvec[e->next].prev = e->prev;
508     if (e->prev == -1) {
509         nodes[n0].adjlist = e->next;
510     } else {
511         nodes[n0].adjvec[e->prev].next = e->next;
512     }
513     e->value = v;
514 
515     e = &nodes[n1].adjvec[n0];
516     if (e->next != -1) nodes[n1].adjvec[e->next].prev = e->prev;
517     if (e->prev == -1) {
518         nodes[n1].adjlist = e->next;
519     } else {
520         nodes[n1].adjvec[e->prev].next = e->next;
521     }
522     e->value = v;
523 
524     s->currentval += v * e->weight;
525 
526     s->searchhist[s->histtop].from = n0;
527     s->searchhist[s->histtop].to = n1;
528     s->searchhist[s->histtop].pfrom = p0;
529     s->searchhist[s->histtop].pto = p1;
530     s->histtop++;
531 
532     if (v >= 1) {
533         if (p0 != OUTSIDE && p1 != OUTSIDE &&
534             (p0 != n0 || p1 != n1)) {
535 #ifdef DEBUG
536 #if DEBUG>1
537             printf ("Avoiding subtour, forcing %d-%d to 0\n",
538                     p0,p1);
539 #endif
540 #endif
541             if (set_edge (s, p0, p1, 0)) return 1;
542         }
543     }
544     checkout_node (s, n0);
545     checkout_node (s, n1);
546 
547     return 0;
548 }
549 
unset_edge(tspsearch * s,histent * h)550 static void unset_edge(tspsearch *s, histent *h)
551 {
552     int n0 = h->from;
553     int n1 = h->to;
554     edge *e;
555     node *nodes = s->nodes;
556     int v = nodes[n0].adjvec[n1].value;
557     int lo, hi;
558 
559     if (v == -1) return;
560 
561 #ifdef DEBUG
562 #if DEBUG>2
563     printf ("unset_edge %d-%d from %d\n",n0,n1,v);
564 #endif
565 #endif
566     lo = nodes[n0].adjvec[n1].lo;
567     hi = nodes[n0].adjvec[n1].hi;
568     nodes[n0].deg += lo - v;
569     nodes[n1].deg += lo - v;
570     nodes[n0].availdeg += hi - v;
571     nodes[n1].availdeg += hi - v;
572 
573     if (v >= 1) { /* The only situation in which pathends matter */
574         if (n0 != OUTSIDE) nodes[n0].pathend = h->pfrom;
575         if (n1 != OUTSIDE) nodes[n1].pathend = h->pto;
576         if (h->pfrom != OUTSIDE && h->pfrom != n0)
577             nodes[h->pfrom].pathend = n0;
578         if (h->pto != OUTSIDE && h->pto != n1) nodes[h->pto].pathend = n1;
579     }
580 
581     e = &nodes[n0].adjvec[n1];
582     if (e->next != -1) nodes[n0].adjvec[e->next].prev = n1;
583     if (e->prev == -1) {
584         nodes[n0].adjlist = n1;
585     } else {
586         nodes[n0].adjvec[e->prev].next = n1;
587     }
588     e->value = -1;
589 
590     e = &nodes[n1].adjvec[n0];
591     if (e->next != -1) nodes[n1].adjvec[e->next].prev = n0;
592     if (e->prev == -1) {
593         nodes[n1].adjlist = n0;
594     } else {
595         nodes[n1].adjvec[e->prev].next = n0;
596     }
597     e->value = -1;
598 
599     s->currentval -= v * e->weight;
600 }
601 
602 #ifdef USE_DUMBO
603 
select_edge(tspsearch * s,int * n0,int * n1)604 static int select_edge (tspsearch *s, int *n0, int *n1)
605 {
606     int i, j;
607     int bestn0, bestn1;
608     double bestval;
609     node *nodes = s->nodes;
610     int nnodes = s->nnodes;
611     edge *adjvec;
612     int nforced;
613     double val;
614     double val0;
615     double val1;
616     int infeas0;
617     int infeas1;
618     int smark;
619 
620     do {
621         nforced = 0;
622         bestn0 = bestn1 = -1;
623         bestval = 0.0;
624         for (i=1; i<nnodes; i++) {
625             adjvec = nodes[i].adjvec;
626             for (j=nodes[i].adjlist; j != -1; j = adjvec[j].next) {
627                 if (j != OUTSIDE && adjvec[j].value == -1) {
628 /* adjvec[j].value is tested since the loop can change the list */
629 /* since set_edge and unset_edge preserve list order, and don't */
630 /* change the next/prev of the deleted items, this is safe */
631                     smark = s->histtop;
632                     infeas0 = set_edge (s, i, j, 0);
633                     if (infeas0 == 0) {
634                         val0 = s->currentval + lowerbound(s);
635                     }
636                     unroll_stack (s, smark);
637                     if (infeas0 == 0 && val0 < s->bestval) {
638                         infeas1 = set_edge (s, i, j, 1);
639                         if (infeas1 == 0) {
640                             val1 = s->currentval + lowerbound(s);
641                         }
642                         unroll_stack (s, smark);
643                     }
644 
645 #ifdef DEBUG
646 #if DEBUG>1
647                     printf ("eval edge %d-%d =  %.0f %.0f = %.0f\n",
648                             i,j,infeas0?-999.0:val0,infeas1?-999.0:val1,
649                             TINYBRANCH_VAL (val0, val1));
650 #endif
651 #endif
652                     if ((infeas0 || val0 >= s->bestval) &&
653                         (infeas1 || val1 >= s->bestval)) {
654                         return SELECT_INFEASIBLE;
655                     } else if (infeas0 || val0 >= s->bestval) {
656                         if (set_edge (s, i, j, 1)) {
657                             return SELECT_INFEASIBLE;
658                         }
659                         nforced++;
660                     } else if (infeas1 || val1 >= s->bestval) {
661                         if (set_edge (s, i, j, 0)) {
662                             return SELECT_INFEASIBLE;
663                         }
664                         nforced++;
665                     } else {
666                         val = TINYBRANCH_VAL(val0, val1);
667                         if (bestn0 == -1 || val > bestval) {
668                             bestn0 = i;
669                             bestn1 = j;
670                             bestval = val;
671                         }
672                     }
673                 }
674             }
675         }
676 #ifdef DEBUG
677         printf ("%d edges forced (bound %.0f val %.0f)\n",
678                 nforced, s->currentval + lowerbound(s), bestval);
679         fflush (stdout);
680 #endif
681     } while (nforced);
682 
683     if (bestn0 == -1) {
684         return SELECT_TOUR;
685     }
686     *n0 = bestn0;
687     *n1 = bestn1;
688     return 0;
689 }
690 
691 #else /* USE_DUMBO */
692 
select_edge(tspsearch * s,int * n0,int * n1)693 static int select_edge (tspsearch *s, int *n0, int *n1)
694 {
695     int i, j;
696     int bestn0 = -1, bestn1 = -1;
697     int bestval;
698     edge *adjvec;
699     node *nodes = s->nodes;
700     int nnodes = s->nnodes;
701 
702     bestval = -1;
703 
704     for (i=0; i<nnodes; i++) {
705         if (i != OUTSIDE) {
706             adjvec = nodes[i].adjvec;
707             for (j=nodes[i].adjlist; j != -1; j = adjvec[j].next) {
708                 if (j != OUTSIDE) {
709                     if (adjvec[j].weight > bestval) {
710                         bestn0 = i;
711                         bestn1 = j;
712                         bestval = adjvec[j].weight;
713                     } else if (-adjvec[j].weight > bestval) {
714                         bestn0 = i;
715                         bestn1 = j;
716                         bestval = -adjvec[j].weight;
717                     }
718                 }
719             }
720         }
721     }
722     if (bestval >= 0) {
723         *n0 = bestn0;
724         *n1 = bestn1;
725         return 0;
726     }
727 
728     i = OUTSIDE;
729     adjvec = nodes[i].adjvec;
730     for (j=nodes[i].adjlist; j != -1; j = adjvec[j].next) {
731         if (adjvec[j].weight > bestval) {
732             bestn0 = i;
733             bestn1 = j;
734             bestval = adjvec[j].weight;
735         } else if (-adjvec[j].weight > bestval) {
736             bestn0 = i;
737             bestn1 = j;
738             bestval = -adjvec[j].weight;
739         }
740     }
741 
742     if (bestval == -1) {
743         return SELECT_TOUR;
744     } else {
745         *n0 = bestn0;
746         *n1 = bestn1;
747         return 0;
748     }
749 }
750 
751 #endif /* USE_DUMBO */
752 
savetour(tspsearch * s)753 static void savetour (tspsearch *s)
754 {
755     int i, j;
756 
757 #ifdef DEBUG
758     printf ("Solution found, saving tour, value %.0f\n", s->currentval);
759 #endif
760 
761     s->bestval = s->currentval;
762     for (i=0; i<s->nnodes; i++) {
763         for (j=0; j<s->nnodes; j++) {
764             s->nodes[i].adjvec[j].bestvalue = s->nodes[i].adjvec[j].value;
765         }
766     }
767 }
768 
unroll_stack(tspsearch * s,int smark)769 static void unroll_stack (tspsearch *s, int smark)
770 {
771     int i = s->histtop;
772 
773     while (i > smark) {
774         i--;
775         unset_edge (s, &s->searchhist[i]);
776     }
777     s->histtop = smark;
778 }
779 
lowerbound(tspsearch * s)780 static double lowerbound (tspsearch *s)
781 {
782     int i, j;
783     node *nodes = s->nodes;
784     int nnodes = s->nnodes;
785     double lb = 0.0;
786     int min1, min2;
787     int cnt;
788 
789     for (i=0; i<nnodes; i++) {
790         if (i == OUTSIDE) {
791             min1 = INT_MAX;
792             cnt = 0;
793             for (j=nodes[i].adjlist; j != -1; j = nodes[i].adjvec[j].next) {
794                 if (nodes[i].adjvec[j].weight < 0) {
795                     lb += nodes[i].adjvec[j].weight;
796                     cnt++;
797                     if (-nodes[i].adjvec[j].weight < min1) {
798                         min1 = -nodes[i].adjvec[j].weight;
799                     }
800                 } else if (nodes[i].adjvec[j].weight < min1) {
801                     min1 = nodes[i].adjvec[j].weight;
802                 }
803             }
804             if ((nodes[i].deg + cnt) & 1) {
805                 lb += min1;
806             }
807         } else if (nodes[i].deg == 0) {
808             min1 = INT_MAX;
809             min2 = INT_MAX;
810             for (j=nodes[i].adjlist; j != -1; j = nodes[i].adjvec[j].next) {
811                 if (nodes[i].adjvec[j].weight < min2) {
812                     if (nodes[i].adjvec[j].weight < min1) {
813                         if (nodes[i].adjvec[j].hi - nodes[i].adjvec[j].lo
814                             <= 1) {
815                             min2 = min1;
816                             min1 = nodes[i].adjvec[j].weight;
817                         } else {
818                             min1 = nodes[i].adjvec[j].weight;
819                             min2 = min1;
820                         }
821                     } else {
822                         min2 = nodes[i].adjvec[j].weight;
823                     }
824                 }
825             }
826             lb += min1 + min2;
827         } else if (nodes[i].deg == 1) {
828             min1 = INT_MAX;
829             for (j=nodes[i].adjlist; j != -1; j = nodes[i].adjvec[j].next) {
830                 if (nodes[i].adjvec[j].weight < min1) {
831                     min1 = nodes[i].adjvec[j].weight;
832                 }
833             }
834             lb += min1;
835         }
836     }
837     return ceil(lb/2.0);
838 }
839