From: dank@blacks.jpl.nasa.gov (Dan Kegel) Newsgroups: alt.sources,sci.math Subject: SUMMARY: Any Traveling Salesman Problem source code out there? Summary: book titles and two source codes Date: 27 May 92 20:30:22 GMT Archive-name: travelling-salesman-kegel In article I asked: >Does anyone have any source code for a TSP solver out there? >I need to rewrite a program that chooses a wirewrapping path, and >would like to avoid reinventing the wheel if possible. >Reply to dank@blacks.jpl.nasa.gov and I'll summarize. I got many helpful answers and one complete Fortran source code. In the end, I implemented the Furthest Insertion Heuristic in C from a Pascal example in one of the books. I enclose first the responses, then my source code. Sorry for taking so long. Thanks, everybody! - Dan Kegel From: trick@kaos.gsia.cmu.edu (Michael Trick) Date: Thu, 12 Mar 92 20:24:24 EST There was one program out a while ago called TRAVEL, by Boyd, Cornuejols, and Pulleyblank, but it was more aligned towards graphical representation of results (it was written in BASIC). ... From: mjs@hubcap.clemson.edu (M. J. Saltzman) Date: Fri, 13 Mar 92 10:04:11 -0500 Organization: Clemson University, Clemson SC Syslo, Deo and Kowalik's book, _Discrete_Optimization_Algorithms_, has some Pascal codes for Euclidean TSP, including a branch-and-bound code, a farthest-insertion heuristic, and 2-opt and 3-opt improvements. _Numerical_Recipes_ has a general-purpose simulated annealing code that could be adapted. By and large, I don't think these are going to be all that good (although they may be OK for problems of modest size). The best improvement heuristic by far is Lin-Kernighan's "variable-depth interchange", carefully tuned. The B&B code does not use any of the state-of-the-art theory for TSP, so it will probably not find optimal solutions in acceptable time for large problems. Unfortunately, I don't know of anyone distributing a state-of-the-art TSP code, although I know of a few being worked on. This is a hot research area, and good codes are also highly commercializable. If you do hear of any, I would like to know, as I am collecting data on optimization codes for eventual publication on the net. -- Matthew Saltzman Clemson University Math Sciences mjs@clemson.edu Date: Fri, 13 Mar 92 17:19:11 EST From: alan@mel.dms.csiro.au (Alan Miller) To: dank@blacks.jpl.nasa.gov Subject: Re: Any Traveling Salesman Problem source code out there? Newsgroups: sci.math,comp.sources.wanted In article you write: >Does anyone have any source code for a TSP solver out there? >I need to rewrite a program that chooses a wirewrapping path, and >would like to avoid reinventing the wheel if possible. >Reply to dank@blacks.jpl.nasa.gov and I'll summarize. >Thanks! SUBROUTINE TSP (N,NN,NN2,DIST,IDIM,BIG,EPS,ISOL, + WK1,WK2,WK3,WK4,WK5,IWK6,IWK7,IWK8,IWK9,IWK10, + IWK11,IWK12,IWK13,IWK14,IWK15,IWK16,IWK17,IWK18) c c From the book: c Lau, H.T. (1986). 'Combinatorial heuristic algorithms with c FORTRAN ', Lecture notes in Economics and Mathematical c systems no. 280, Springer-Verlag. C C HEURISTIC FOR THE TRAVELING SALESMAN PROBLEM C SATISFYING THE TRIANGLE INEQUALITY C REAL DIST(IDIM,1),WK1(NN2),WK2(N),WK3(N),WK4(N),WK5(N) INTEGER ISOL(N),IWK6(NN),IWK7(NN),IWK8(NN),IWK9(NN), + IWK10(NN),IWK11(N),IWK12(N),IWK13(N),IWK14(N), + IWK15(N),IWK16(N),IWK17(N),IWK18(N) C C CONSTRUCT A MINIMUM SPANNING TREE C CALL MINTRE(N,DIST,IDIM,BIG,IWK6,IWK7,WK2,IWK15,IWK16) C C DETERMINE THE SET OF NODES WITH ODD DEGREES IN THE C MINIMUM SPANNING TREE C DO 10 I = 1, N 10 IWK10(I) = 0 NM1 = N - 1 DO 20 I = 1, NM1 II = IWK6(I) JJ = IWK7(I) IWK10(II) = IWK10(II) + 1 IWK10(JJ) = IWK10(JJ) + 1 20 CONTINUE ICT = 0 DO 30 I = 1, N IF (MOD(IWK10(I),2) .NE. 0) THEN ICT = ICT + 1 IWK11(ICT) = I ENDIF 30 CONTINUE C C DETERMINE A MINIMUM WEIGHT PERFECT MATCHING FOR C THE SET OF ODD-DEGREE NODES C K = 0 DO 50 I = 2, ICT I1 = I - 1 DO 40 J = 1, I1 K = K + 1 WK1(K) = DIST(IWK11(J),IWK11(I)) 40 CONTINUE 50 CONTINUE C CALL PMATCH(ICT,NN2,WK1,BIG,EPS,IWK18,WK2,WK3,WK4,WK5,IWK8, + IWK9,IWK10,IWK12,IWK13,IWK14,IWK15,IWK16,IWK17) C C STORE UP THE EDGES IN THE PERFECT MATCHING C DO 60 I = 1, ICT 60 IWK18(I) = IWK11(IWK18(I)) DO 70 I = 1, N 70 IWK10(I) = 0 K = N - 1 DO 80 I = 1, ICT IF (IWK10(IWK11(I)) .EQ. 0) THEN IWK10(IWK11(I)) = 1 IWK10(IWK18(I)) = 1 K = K + 1 IWK6(K) = IWK11(I) IWK7(K) = IWK18(I) ENDIF 80 CONTINUE C C FIND AN EULER CIRCUIT C M = N - 1 + (ICT / 2) CALL EULER(N,M,IWK6,IWK7,IWK10,IWK8, + IWK9,IWK12,IWK13,IWK14,IWK15) C C FORM THE HAMILTONIAN CIRCUIT C ISOL(1) = IWK10(1) J = 2 ISOL(J) = IWK10(J) K = 2 90 J = J + 1 IBASE = IWK10(J) J1 = J - 1 DO 100 I = 1, J1 IF (IBASE .EQ. IWK10(I)) GOTO 90 100 CONTINUE K = K + 1 ISOL(K) = IWK10(J) IF (K .LT. N) GOTO 90 C RETURN END SUBROUTINE EULER (N,M,INODE,JNODE,LOOP, + IWORK1,IWORK2,IWORK3,IWORK4,IWORK5,IWORK6) C INTEGER INODE(M),JNODE(M),LOOP(M),IWORK1(M),IWORK2(M), + IWORK3(N),IWORK4(N),IWORK5(N),IWORK6(N) LOGICAL FOUND,COPYON C DO 10 I = 1, N 10 IWORK3(I) = 0 DO 20 I = 1, M LOOP(I) = 0 IWORK1(I) = 0 20 IWORK2(I) = 0 NUMARC = 1 IWORK2(1) = 1 NNODE = 1 I = INODE(1) IWORK1(NNODE) = I IWORK3(I) = 1 NNODE = NNODE + 1 J = JNODE(1) IWORK1(NNODE) = J IWORK3(J) = 1 IBASE = J NBREAK = 0 C C LOOK FOR THE NEXT ARC C 30 DO 40 I = 2, M IF (IWORK2(I) .EQ. 0) THEN FOUND = .FALSE. IF (IBASE .EQ. INODE(I)) THEN FOUND = .TRUE. IBASE = JNODE(I) ELSE IF (IBASE .EQ. JNODE(I)) THEN FOUND = .TRUE. IBASE = INODE(I) ENDIF ENDIF IF (FOUND) THEN IWORK2(I) = 1 NUMARC = NUMARC + 1 NNODE = NNODE + 1 IF (NNODE .LE. M) IWORK1(NNODE) = IBASE IWORK3(IBASE) = 1 GOTO 30 ENDIF ENDIF 40 CONTINUE C C A CYCLE HAS BEEN FOUND C IF (NBREAK .GT. 0) THEN NNODE = NNODE - 1 IWORK5(NBREAK) = NNODE ENDIF IF (NUMARC .LT. M) THEN C C FIND A NODE IN THE CURRENT EULER CIRCUIT C DO 50 I = 2, M IF (IWORK2(I) .EQ. 0) THEN FOUND = .FALSE. IF (IWORK3(INODE(I)) .NE. 0) THEN FOUND = .TRUE. J = INODE(I) K = JNODE(I) ELSE IF (IWORK3(JNODE(I)) .NE. 0) THEN FOUND = .TRUE. J = JNODE(I) K = INODE(I) ENDIF ENDIF C C IDENTIFY THE PATH WHICH WILL BE ADDED TO THE CIRCUIT C IF (FOUND) THEN NBREAK = NBREAK + 1 IWORK6(NBREAK) = J IBASE = K IWORK3(K) = 1 NNODE = NNODE + 1 IWORK4(NBREAK) = NNODE IWORK1(NNODE) = IBASE IWORK2(I) = 1 NUMARC = NUMARC + 1 GOTO 30 ENDIF ENDIF 50 CONTINUE ENDIF C C FORM THE EULER CIRCUIT C IF (NBREAK .EQ. 0) THEN NNODE = NNODE - 1 DO 60 I = 1, NNODE 60 LOOP(I) = IWORK1(I) RETURN ENDIF INSERT = 1 IPIVOT = IWORK6(INSERT) IFORWD = 0 70 NCOPY = 1 IBASE = IWORK1(1) LOCBAS = 1 LOOP(NCOPY) = IBASE C C A PATH WHICH WAS IDENTIFIED BEFORE IS ADDED TO THE CIRCUIT C 80 IF (IBASE .EQ. IPIVOT) THEN J = IWORK4(INSERT) + IFORWD K = IWORK5(INSERT) + IFORWD DO 90 L = J, K NCOPY = NCOPY + 1 LOOP(NCOPY) = IWORK1(L) IWORK1(L) = 0 90 CONTINUE NCOPY = NCOPY + 1 C C ADD THE INTERSECTING NODE TO THE CIRCUIT C LOOP(NCOPY) = IBASE IFORWD = IFORWD + 1 IF (NCOPY .LT. NNODE) THEN 100 IF (NCOPY .LT. M) THEN LOCBAS = LOCBAS + 1 IF (LOCBAS .LT. M) THEN IBASE = IWORK1(LOCBAS) IF (IBASE .NE. 0) THEN NCOPY = NCOPY + 1 LOOP(NCOPY) = IBASE ENDIF GOTO 100 ENDIF ENDIF ENDIF ELSE NCOPY = NCOPY + 1 IF (NCOPY .LE. NNODE) THEN LOCBAS = LOCBAS + 1 IBASE = IWORK1(LOCBAS) LOOP(NCOPY) = IBASE GOTO 80 ENDIF ENDIF C C CHECK IF MORE PATHS ARE TO BE ADDED TO THE CIRCUIT C COPYON = .FALSE. INSERT = INSERT + 1 IF (INSERT .LE. NBREAK) THEN COPYON = .TRUE. IPIVOT = IWORK6(INSERT) ENDIF IF (COPYON) THEN DO 110 I = 1, M C IF (LOOP(I) .NE. 0) IWORK1(I) = LOOP(I) C LOOP(I) = 0 IWORK1(I) = LOOP(I) 110 CONTINUE GOTO 70 ENDIF C RETURN END SUBROUTINE MINTRE (N,DIST,IDIM,BIG,INODE,JNODE, + WORK,IWORK1,IWORK2) C INTEGER INODE(N),JNODE(N),IWORK1(N),IWORK2(N) REAL DIST(IDIM,1),WORK(N) C DO 10 I = 1, N WORK(I) = BIG IWORK1(I) = 0 IWORK2(I) = 0 10 CONTINUE C C FIND THE FIRST NON-ZERO ARC C DO 20 IJ = 1, N DO 20 KJ = 1, N IF (DIST(IJ,KJ) .LT. BIG) THEN I = IJ GO TO 30 ENDIF 20 CONTINUE 30 WORK(I) = 0 IWORK1(I) = 1 XLEN = 0. KK4 = N - 1 DO 80 JJ = 1, KK4 DO 40 K = 1, N 40 WORK(K) = BIG DO 60 I = 1, N C C FOR EACH FORWARD ARC ORIGINATING AT NODE I CALCULATE C THE LENGTH OF THE PATH TO NODE I. C IF (IWORK1(I) .EQ. 1) THEN DO 50 J = 1, N IF (DIST(I,J).LT.BIG.AND.IWORK1(J).EQ.0) THEN D = XLEN + DIST(I,J) IF (D .LT. WORK(J)) THEN WORK(J) = D IWORK2(J) = I ENDIF ENDIF 50 CONTINUE ENDIF 60 CONTINUE C C FIND THE MINIMUM POTENTIAL. C D = BIG IENT = 0 DO 70 I = 1, N IF (IWORK1(I) .EQ. 0 .AND. WORK(I) .LT. D) THEN D = WORK(I) IENT = I ITR = IWORK2(I) ENDIF 70 CONTINUE C C INCLUDE THE NODE IN THE CURRENT PATH. C IF (D .LT. BIG) THEN IWORK1(IENT) = 1 XLEN = XLEN + DIST(ITR,IENT) INODE(JJ) = ITR JNODE(JJ) = IENT ENDIF 80 CONTINUE C RETURN END SUBROUTINE PMATCH (N,NN2,COST,BIG,EPS,IPAIR, + WORK1,WORK2,WORK3,WORK4,JWK1,JWK2, + JWK3,JWK4,JWK5,JWK6,JWK7,JWK8,JWK9) C INTEGER IPAIR(N),JWK1(N),JWK2(N),JWK3(N),JWK4(N),JWK5(N), + JWK6(N),JWK7(N),JWK8(N),JWK9(N) REAL COST(NN2),WORK1(N),WORK2(N),WORK3(N),WORK4(N) C C INITIALIZATION C JWK1(2) = 0 DO 10 I = 3, N JWK1(I) = JWK1(I-1) + I - 2 10 CONTINUE IHEAD = N + 2 DO 20 I = 1, N JWK2(I) = I JWK3(I) = I JWK4(I) = 0 JWK5(I) = I JWK6(I) = IHEAD JWK7(I) = IHEAD JWK8(I) = IHEAD IPAIR(I) = IHEAD WORK1(I) = BIG WORK2(I) = 0. WORK3(I) = 0. WORK4(I) = BIG 20 CONTINUE C C START PROCEDURE C DO 50 I = 1, N IF (IPAIR(I) .EQ. IHEAD) THEN NN = 0 CWK2 = BIG DO 40 J = 1, N MIN = I MAX = J IF (I .NE. J) THEN IF (I .GT. J) THEN MAX = I MIN = J ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) CSWK = COST(ISUB) - WORK2(J) IF (CSWK .LE. CWK2) THEN IF (CSWK .EQ. CWK2) THEN IF (NN .EQ. 0) GO TO 30 GOTO 40 ENDIF CWK2 = CSWK NN = 0 30 IF (IPAIR(J) .EQ. IHEAD) NN = J ENDIF ENDIF 40 CONTINUE IF (NN .NE. 0) THEN WORK2(I) = CWK2 IPAIR(I) = NN IPAIR(NN) = I ENDIF ENDIF 50 CONTINUE C C INITIAL LABELING C NN = 0 DO 70 I = 1, N IF (IPAIR(I) .EQ. IHEAD) THEN NN = NN + 1 JWK6(I) = 0 WORK4(I) = 0. XWK2 = WORK2(I) DO 60 J = 1, N MIN = I MAX = J IF (I .NE. J) THEN IF (I .GT. J) THEN MAX = I MIN = J ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) CSWK = COST(ISUB) - XWK2 - WORK2(J) IF (CSWK .LT. WORK1(J)) THEN WORK1(J) = CSWK JWK4(J) = I ENDIF ENDIF 60 CONTINUE ENDIF 70 CONTINUE IF (NN .LE. 1) GO TO 340 C C EXAMINE THE LABELING AND PREPARE FOR THE NEXT STEP C 80 CSTLOW = BIG DO 90 I = 1, N IF (JWK2(I) .EQ. I) THEN CST = WORK1(I) IF (JWK6(I) .LT. IHEAD) THEN CST = 0.5 * (CST + WORK4(I)) IF (CST .LE. CSTLOW) THEN INDEX = I CSTLOW = CST ENDIF ELSE IF (JWK7(I) .LT. IHEAD) THEN IF (JWK3(I) .NE. I) THEN CST = CST + WORK2(I) IF (CST .LT. CSTLOW) THEN INDEX = I CSTLOW = CST ENDIF ENDIF ELSE IF (CST .LT. CSTLOW) THEN INDEX = I CSTLOW = CST ENDIF ENDIF ENDIF ENDIF 90 CONTINUE IF (JWK7(INDEX) .LT. IHEAD) GO TO 190 IF (JWK6(INDEX) .LT. IHEAD) THEN LL4 = JWK4(INDEX) LL5 = JWK5(INDEX) KK4 = INDEX KK1 = KK4 KK5 = JWK2(LL4) KK2 = KK5 100 JWK7(KK1) = KK2 MM5 = JWK6(KK1) IF (MM5 .NE. 0) THEN KK2 = JWK2(MM5) KK1 = JWK7(KK2) KK1 = JWK2(KK1) GO TO 100 ENDIF LL2 = KK1 KK1 = KK5 KK2 = KK4 110 IF (JWK7(KK1) .GE. IHEAD) THEN JWK7(KK1) = KK2 MM5 = JWK6(KK1) IF (MM5 .EQ. 0) GO TO 280 KK2 = JWK2(MM5) KK1 = JWK7(KK2) KK1 = JWK2(KK1) GO TO 110 ENDIF 120 IF (KK1 .EQ. LL2) GO TO 130 MM5 = JWK7(LL2) JWK7(LL2) = IHEAD LL1 = IPAIR(MM5) LL2 = JWK2(LL1) GO TO 120 ENDIF C C GROWING AN ALTERNATING TREE, ADD TWO EDGES C JWK7(INDEX) = JWK4(INDEX) JWK8(INDEX) = JWK5(INDEX) LL1 = IPAIR(INDEX) LL3 = JWK2(LL1) WORK4(LL3) = CSTLOW JWK6(LL3) = IPAIR(LL3) CALL SUBB(LL3,N,NN2,BIG,COST,JWK1,JWK2,JWK3,JWK4, + JWK5,JWK7,JWK9,WORK1,WORK2,WORK3,WORK4) GO TO 80 C C SHRINK A BLOSSOM C 130 XWORK = WORK2(LL2) + CSTLOW - WORK4(LL2) WORK2(LL2) = 0. MM1 = LL2 140 WORK3(MM1) = WORK3(MM1) + XWORK MM1 = JWK3(MM1) IF (MM1 .NE. LL2) GO TO 140 MM5 = JWK3(LL2) IF (LL2 .NE. KK5) GO TO 160 150 KK5 = KK4 KK2 = JWK7(LL2) 160 JWK3(MM1) = KK2 LL1 = IPAIR(KK2) JWK6(KK2) = LL1 XWK2 = WORK2(KK2) + WORK1(KK2) - CSTLOW MM1 = KK2 170 MM2 = MM1 WORK3(MM2) = WORK3(MM2) + XWK2 JWK2(MM2) = LL2 MM1 = JWK3(MM2) IF (MM1 .NE. KK2) GO TO 170 JWK5(KK2) = MM2 WORK2(KK2) = XWK2 KK1 = JWK2(LL1) JWK3(MM2) = KK1 XWK2 = WORK2(KK1) + CSTLOW - WORK4(KK1) MM2 = KK1 180 MM1 = MM2 WORK3(MM1) = WORK3(MM1) + XWK2 JWK2(MM1) = LL2 MM2 = JWK3(MM1) IF (MM2 .NE. KK1) GO TO 180 JWK5(KK1) = MM1 WORK2(KK1) = XWK2 IF (KK5 .NE. KK1) THEN KK2 = JWK7(KK1) JWK7(KK1) = JWK8(KK2) JWK8(KK1) = JWK7(KK2) GO TO 160 ENDIF IF (KK5 .NE. INDEX) THEN JWK7(KK5) = LL5 JWK8(KK5) = LL4 IF (LL2 .NE. INDEX) GO TO 150 ELSE JWK7(INDEX) = LL4 JWK8(INDEX) = LL5 ENDIF JWK3(MM1) = MM5 KK4 = JWK3(LL2) JWK4(KK4) = MM5 WORK4(KK4) = XWORK JWK7(LL2) = IHEAD WORK4(LL2) = CSTLOW CALL SUBB(LL2,N,NN2,BIG,COST,JWK1,JWK2,JWK3,JWK4, + JWK5,JWK7,JWK9,WORK1,WORK2,WORK3,WORK4) GO TO 80 C C EXPAND A T-LABELED BLOSSOM C 190 KK4 = JWK3(INDEX) KK3 = KK4 LL4 = JWK4(KK4) MM2 = KK4 200 MM1 = MM2 LL5 = JWK5(MM1) XWK2 = WORK2(MM1) 210 JWK2(MM2) = MM1 WORK3(MM2) = WORK3(MM2) - XWK2 IF (MM2 .NE. LL5) THEN MM2 = JWK3(MM2) GO TO 210 ENDIF MM2 = JWK3(LL5) JWK3(LL5) = MM1 IF (MM2 .NE. LL4) GO TO 200 XWK2 = WORK4(KK4) WORK2(INDEX) = XWK2 JWK3(INDEX) = LL4 MM2 = LL4 220 WORK3(MM2) = WORK3(MM2) - XWK2 IF (MM2 .NE. INDEX) THEN MM2 = JWK3(MM2) GO TO 220 ENDIF MM1 = IPAIR(INDEX) KK1 = JWK2(MM1) MM2 = JWK6(KK1) LL2 = JWK2(MM2) IF (LL2 .NE. INDEX) THEN KK2 = LL2 230 MM5 = JWK7(KK2) KK1 = JWK2(MM5) IF (KK1 .NE. INDEX) THEN KK2 = JWK6(KK1) KK2 = JWK2(KK2) GO TO 230 ENDIF JWK7(LL2) = JWK7(INDEX) JWK7(INDEX) = JWK8(KK2) JWK8(LL2) = JWK8(INDEX) JWK8(INDEX) = MM5 MM3 = JWK6(LL2) KK3 = JWK2(MM3) MM4 = JWK6(KK3) JWK6(LL2) = IHEAD IPAIR(LL2) = MM1 KK1 = KK3 240 MM1 = JWK7(KK1) MM2 = JWK8(KK1) JWK7(KK1) = MM4 JWK8(KK1) = MM3 JWK6(KK1) = MM1 IPAIR(KK1) = MM1 KK2 = JWK2(MM1) IPAIR(KK2) = MM2 MM3 = JWK6(KK2) JWK6(KK2) = MM2 IF (KK2 .NE. INDEX) THEN KK1 = JWK2(MM3) MM4 = JWK6(KK1) JWK7(KK2) = MM3 JWK8(KK2) = MM4 GO TO 240 ENDIF ENDIF MM2 = JWK8(LL2) KK1 = JWK2(MM2) WORK1(KK1) = CSTLOW KK4 = 0 IF (KK1 .NE. LL2) THEN MM1 = JWK7(KK1) KK3 = JWK2(MM1) JWK7(KK1) = JWK7(LL2) JWK8(KK1) = MM2 250 MM5 = JWK6(KK1) JWK6(KK1) = IHEAD KK2 = JWK2(MM5) MM5 = JWK7(KK2) JWK7(KK2) = IHEAD KK5 = JWK8(KK2) JWK8(KK2) = KK4 KK4 = KK2 WORK4(KK2) = CSTLOW KK1 = JWK2(MM5) WORK1(KK1) = CSTLOW IF (KK1 .NE. LL2) GO TO 250 JWK7(LL2) = KK5 JWK8(LL2) = MM5 JWK6(LL2) = IHEAD IF (KK3 .EQ. LL2) GO TO 270 ENDIF KK1 = 0 KK2 = KK3 260 MM5 = JWK6(KK2) JWK6(KK2) = IHEAD JWK7(KK2) = IHEAD JWK8(KK2) = KK1 KK1 = JWK2(MM5) MM5 = JWK7(KK1) JWK6(KK1) = IHEAD JWK7(KK1) = IHEAD JWK8(KK1) = KK2 KK2 = JWK2(MM5) IF (KK2 .NE. LL2) GO TO 260 CALL SUBA(KK1,N,NN2,BIG,COST,JWK1,JWK2,JWK3,JWK4, + JWK5,JWK6,JWK8,WORK1,WORK2,WORK3,WORK4) C 270 IF (KK4 .EQ. 0) GO TO 80 LL2 = KK4 CALL SUBB(LL2,N,NN2,BIG,COST,JWK1,JWK2,JWK3,JWK4, + JWK5,JWK7,JWK9,WORK1,WORK2,WORK3,WORK4) KK4 = JWK8(LL2) JWK8(LL2) = IHEAD GO TO 270 C C AUGMENTATION OF THE MATCHING C EXCHANGE THE MATCHING AND NON-MATCHING EDGES ALONG THE C AUGMENTING PATH C 280 LL2 = KK4 MM5 = LL4 290 KK1 = LL2 300 IPAIR(KK1) = MM5 MM5 = JWK6(KK1) JWK7(KK1) = IHEAD IF (MM5 .NE. 0) THEN KK2 = JWK2(MM5) MM1 = JWK7(KK2) MM5 = JWK8(KK2) KK1 = JWK2(MM1) IPAIR(KK2) = MM1 GO TO 300 ENDIF IF (LL2 .EQ. KK4) THEN LL2 = KK5 MM5 = LL5 GO TO 290 ENDIF C C REMOVE ALL LABELS ON ON-EXPOSED BASE NODES C DO 310 I = 1, N IF (JWK2(I) .EQ. I) THEN IF (JWK6(I) .LT. IHEAD) THEN CST = CSTLOW - WORK4(I) WORK2(I) = WORK2(I) + CST JWK6(I) = IHEAD IF (IPAIR(I) .NE. IHEAD) THEN WORK4(I) = BIG ELSE JWK6(I) = 0 WORK4(I) = 0. ENDIF ELSE IF (JWK7(I) .LT. IHEAD) THEN CST = WORK1(I) - CSTLOW WORK2(I) = WORK2(I) + CST JWK7(I) = IHEAD JWK8(I) = IHEAD ENDIF WORK4(I) = BIG ENDIF WORK1(I) = BIG ENDIF 310 CONTINUE NN = NN - 2 IF (NN .LE. 1) GO TO 340 C C DETERMINE THE NEW WORK1 VALUES C DO 330 I = 1, N KK1 = JWK2(I) IF (JWK6(KK1) .EQ. 0) THEN XWK2 = WORK2(KK1) XWK3 = WORK3(I) DO 320 J = 1, N KK2 = JWK2(J) IF (KK1 .NE. KK2) THEN MIN = I MAX = J IF (I .NE. J) THEN IF (I .GT. J) THEN MAX = I MIN = J ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) CSWK = COST(ISUB) - XWK2 - XWK3 CSWK = CSWK - WORK2(KK2) - WORK3(J) IF (CSWK .LT. WORK1(KK2)) THEN JWK4(KK2) = I JWK5(KK2) = J WORK1(KK2) = CSWK ENDIF ENDIF ENDIF 320 CONTINUE ENDIF 330 CONTINUE GO TO 80 C C GENERATE THE ORIGINAL GRAPH BY EXPANDING ALL SHRUNKEN BLOSSOMS C 340 VALUE = 0. DO 350 I = 1, N IF (JWK2(I) .EQ. I) THEN IF (JWK6(I) .GE. 0) THEN KK5 = IPAIR(I) KK2 = JWK2(KK5) KK4 = IPAIR(KK2) JWK6(I) = -1 JWK6(KK2) = -1 MIN = KK4 MAX = KK5 IF (KK4 .NE. KK5) THEN IF (KK4 .GT. KK5) THEN MAX = KK4 MIN = KK5 ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) VALUE = VALUE + XCST ENDIF ENDIF ENDIF 350 CONTINUE DO 420 I = 1, N 360 LL2 = JWK2(I) IF (LL2 .EQ. I) GO TO 420 MM2 = JWK3(LL2) LL4 = JWK4(MM2) KK3 = MM2 XWORK = WORK4(MM2) 370 MM1 = MM2 LL5 = JWK5(MM1) XWK2 = WORK2(MM1) 380 JWK2(MM2) = MM1 WORK3(MM2) = WORK3(MM2) - XWK2 IF (MM2 .NE. LL5) THEN MM2 = JWK3(MM2) GO TO 380 ENDIF MM2 = JWK3(LL5) JWK3(LL5) = MM1 IF (MM2 .NE. LL4) GO TO 370 WORK2(LL2) = XWORK JWK3(LL2) = LL4 MM2 = LL4 390 WORK3(MM2) = WORK3(MM2) - XWORK IF (MM2 .NE. LL2) THEN MM2 = JWK3(MM2) GO TO 390 ENDIF MM5 = IPAIR(LL2) MM1 = JWK2(MM5) MM1 = IPAIR(MM1) KK1 = JWK2(MM1) IF (LL2 .NE. KK1) THEN IPAIR(KK1) = MM5 KK3 = JWK7(KK1) KK3 = JWK2(KK3) 400 MM3 = JWK6(KK1) KK2 = JWK2(MM3) MM1 = JWK7(KK2) MM2 = JWK8(KK2) KK1 = JWK2(MM1) IPAIR(KK1) = MM2 IPAIR(KK2) = MM1 MIN = MM1 MAX = MM2 IF (MM1 .EQ. MM2) GOTO 360 IF (MM1 .GT. MM2) THEN MAX = MM1 MIN = MM2 ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) VALUE = VALUE + XCST IF (KK1 .NE. LL2) GO TO 400 IF (KK3 .EQ. LL2) GO TO 360 ENDIF 410 KK5 = JWK6(KK3) KK2 = JWK2(KK5) KK6 = JWK6(KK2) MIN = KK5 MAX = KK6 IF (KK5 .EQ. KK6) GOTO 360 IF (KK5 .GT. KK6) THEN MAX = KK5 MIN = KK6 ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) VALUE = VALUE + XCST KK6 = JWK7(KK2) KK3 = JWK2(KK6) IF (KK3 .EQ. LL2) GO TO 360 GO TO 410 420 CONTINUE C RETURN END SUBROUTINE SUBB (KK,N,NN2,BIG,COST,JWK1,JWK2,JWK3,JWK4, + JWK5,JWK7,JWK9,WORK1,WORK2,WORK3,WORK4) C INTEGER JWK1(N),JWK2(N),JWK3(N),JWK4(N),JWK5(N), + JWK7(N),JWK9(N) REAL COST(NN2),WORK1(N),WORK2(N),WORK3(N),WORK4(N) C IHEAD = N + 2 XWK1 = WORK4(KK) - WORK2(KK) WORK1(KK) = BIG XWK2 = XWK1 - WORK3(KK) JWK7(KK) = 0 II = 0 DO 10 I = 1, N JJ3 = JWK2(I) IF (JWK7(JJ3) .GE. IHEAD) THEN II = II + 1 JWK9(II) = I MIN = KK MAX = I IF (KK .NE. I) THEN IF (KK .GT. I) THEN MAX = KK MIN = I ENDIF ISUB = JWK1(MAX) + MIN CSWK = COST(ISUB) + XWK2 CSWK = CSWK - WORK2(JJ3) - WORK3(I) IF (CSWK .LT. WORK1(JJ3)) THEN JWK4(JJ3) = KK JWK5(JJ3) = I WORK1(JJ3) = CSWK ENDIF ENDIF ENDIF 10 CONTINUE JWK7(KK) = IHEAD JJ1 = KK JJ1 = JWK3(JJ1) IF (JJ1 .EQ. KK) RETURN 20 XWK2 = XWK1 - WORK3(JJ1) DO 30 I = 1, II JJ2 = JWK9(I) JJ3 = JWK2(JJ2) MIN = JJ1 MAX = JJ2 IF (JJ1 .NE. JJ2) THEN IF (JJ1 .GT. JJ2) THEN MAX = JJ1 MIN = JJ2 ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) CSWK = COST(ISUB) + XWK2 CSWK = CSWK - WORK2(JJ3) - WORK3(JJ2) IF (CSWK .LT. WORK1(JJ3)) THEN JWK4(JJ3) = JJ1 JWK5(JJ3) = JJ2 WORK1(JJ3) = CSWK ENDIF ENDIF 30 CONTINUE JJ1 = JWK3(JJ1) IF (JJ1 .NE. KK) GO TO 20 C RETURN END SUBROUTINE SUBA (KK,N,NN2,BIG,COST,JWK1,JWK2,JWK3,JWK4, + JWK5,JWK6,JWK8,WORK1,WORK2,WORK3,WORK4) C INTEGER JWK1(N),JWK2(N),JWK3(N),JWK4(N),JWK5(N),JWK6(N), + JWK8(N) REAL COST(NN2),WORK1(N),WORK2(N),WORK3(N),WORK4(N) C IHEAD = N + 2 10 JJ1 = KK KK = JWK8(JJ1) JWK8(JJ1) = IHEAD CSTWK = BIG JJ3 = 0 JJ4 = 0 J = JJ1 XWK2 = WORK2(JJ1) 20 XWK3 = WORK3(J) DO 30 I = 1, N JJ2 = JWK2(I) IF (JWK6(JJ2) .LT. IHEAD) THEN MIN = J MAX = I IF (J .NE. I) THEN IF (J .GT. I) THEN MAX = J MIN = I ENDIF ISUB = JWK1(MAX) + MIN XCST = COST(ISUB) CSWK = COST(ISUB) - XWK2 - XWK3 CSWK = CSWK - WORK2(JJ2) - WORK3(I) CSWK = CSWK + WORK4(JJ2) IF (CSWK .LT. CSTWK) THEN JJ3 = I JJ4 = J CSTWK = CSWK ENDIF ENDIF ENDIF 30 CONTINUE J = JWK3(J) IF (J .NE. JJ1) GO TO 20 JWK4(JJ1) = JJ3 JWK5(JJ1) = JJ4 WORK1(JJ1) = CSTWK IF (KK .NE. 0) GO TO 10 C RETURN END -- Alan Miller, Quality Improvement Project CSIRO Division of Mathematics & Statistics, Melbourne, Australia Phone: +61 3 542-2266 Fax: +61 3 543-6613 E-mail: alan@mel.dms.csiro.au Mail: CSIRO DMS, Private Bag 10, Clayton, Vic. 3168, Australia -- And now, my source. -- :::::::::::::: tsp.h :::::::::::::: /*-------------------------------------------------------------------------- Definitions for users of the travelling salesman problem module. --------------------------------------------------------------------------*/ #ifndef tsp_h #define tsp_h static char SCCSID_tsp_h[]="@(#)tsp.h 1.1 5/20/92"; /*-------------------------------------------------------------------------- Approximate solver using the Farthest Insertion heuristic. Returns total distance of one-way trip generated by deleting longest edge of cycle. If nih is TRUE, uses Nearest Insertion heuristic instead. --------------------------------------------------------------------------*/ int tsp_fi( /* Inputs: */ int n, /* number of nodes */ int s, /* starting node */ int *x, /* X coordinate of each node */ int *y, /* Y coordinate of each node */ int nih, /* if TRUE, use nearest insertion heuristic */ /* Output: */ int *route); /* route[i=0..n-1]= index of ith node in tour */ /*-------------------------------------------------------------------------- Run tsp_fi in several different ways, as directed by try_n, nih, and fih. Returns total distance of one-way trip generated by deleting longest edge of cycle. --------------------------------------------------------------------------*/ int tsp_fin( /* Inputs: */ int n, /* number of nodes */ int *x, /* X coordinate of each node */ int *y, /* Y coordinate of each node */ int try_n, /* if TRUE, call solver n times */ int nih, /* if TRUE, call nearest-insertion solver */ int fih, /* if TRUE, call farthest-insertion solver */ /* Output: */ int *route); /* route[i=0..n-1]= index of ith node in tour */ #endif :::::::::::::: tsp.c :::::::::::::: /*-------------------------------------------------------------------------- Approximate TSP solvers. See p.363, "Discrete Optimization Algorithms," Maciej Syslo, Prentice-Hall --------------------------------------------------------------------------*/ static char SCCSID[] = "@(#)tsp.c 1.1 5/20/92"; #include #include #include #include #include #ifdef ARRAYTEST /* see pg. 367, Syslo */ #define DIST(a,b) DISTa[a][b] #define BIG 0 /* book is wrong */ static int DISTa[6][6] = { {BIG, 3, 93, 13, 33, 9}, {4, BIG, 77, 42, 21, 16}, {45, 17, BIG, 36, 16, 28}, {39, 90, 80, BIG, 56, 7}, {28, 46, 88, 33, BIG, 25}, {3, 88, 18, 46, 92, BIG}}; #else /* warning: this macro cannot be called twice in the same arithmetic * expression */ static int dx, dy; #define DIST(a,b) \ (dx=(x[a]-x[b]),dy=(y[a]-y[b]),(int)(0.5+sqrt((double)(dx*dx+dy*dy)))) #endif /*-------------------------------------------------------------------------- Approximate solver using the Farthest Insertion heuristic. Returns total distance of one-way trip generated by deleting longest edge of cycle. If nih is TRUE, uses Nearest Insertion heuristic instead. --------------------------------------------------------------------------*/ int tsp_fi( /* Inputs: */ int n, /* number of nodes */ int s, /* starting node */ int *x, /* X coordinate of each node */ int *y, /* Y coordinate of each node */ int nih, /* if TRUE, use nearest insertion heuristic */ /* Output: */ int *route) /* route[i=0..n-1]= index of ith node in tour */ { int end1, end2; int i, j, argmaxcost; int index, nextindex; int inscost, newcost, total_cost, maxcost; int *cycle; /* cycle[i]=next node after node i in tour; -1 if not yet in */ int *dist; /* dist[i] =dist from any node in tour to node i not in tour */ assert(n > 0); /* Allocate variable-sized arrays */ cycle = (int *)malloc(n * sizeof(cycle[0])); assert(cycle != NULL); dist = (int *)malloc(n * sizeof(dist[0])); assert(dist != NULL); /* No nodes in tour yet */ for (i=0; i maxcost)) || (nih&&(dist[j]= 0 && argmaxcost < n); /* Farthest node is argmaxcost, distance is maxcost. */ /* Now find cheapest place to insert it into the tour. */ inscost = MAXINT; index = s; for (j=0; j<=i; j++) { nextindex = cycle[index]; assert(nextindex >= 0 && nextindex < n); newcost = DIST(index, argmaxcost); newcost += DIST(argmaxcost, nextindex); newcost -= DIST(index, nextindex); if (newcost < inscost) { inscost = newcost; end1 = index; end2 = nextindex; } index = nextindex; } /* Cheapest place is after end1 and before end2. Insert it there. */ assert(cycle[end1] == end2); cycle[end1] = argmaxcost; cycle[argmaxcost] = end2; total_cost += inscost; #if 0 printf("Inserting node %d between %d and %d\n", argmaxcost, end1, end2); printf("Total_cost %d += (%d+%d-%d) = %d\n", total_cost-inscost, DIST(end1,argmaxcost), DIST(argmaxcost,end2),DIST(end1,end2),total_cost); #endif /* Now that there's a new node in the tour, update the array of * distances from other nodes to the tour. */ for (j=0; j maxcost) { maxcost = newcost; argmaxcost = index; } index = cycle[index]; } #if 0 printf("Longest move in tour is %d,%d, cost %d\n", argmaxcost, cycle[argmaxcost], maxcost); #endif /* Return the tour to the caller as an array of nodes to visit * rather than a map {current node, next node} as it is in dist[]. * Start tour such that longest move in tour is (route[n-1],route[0]). */ #ifdef ARRAYTEST maxcost = 0; index = s; #else index = cycle[argmaxcost]; #endif for (i=0; i