PROGRAM QAPBB C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IMPLICIT INTEGER (A - Z) parameter (ndim = 33) c auslegung fuer ndim = 33 DIMENSION A(ndim,ndim),B(ndim,ndim),C(ndim,ndim),ZUL(ndim,ndim) dimension UMSPEI(ndim*ndim),LOESG(ndim),U(ndim), V(ndim) * , DD(ndim), PARTPE(ndim), Y(ndim),LAB(ndim),Z1(ndim) * ,H1(ndim),MENG(ndim-2),PHIOFM(ndim-2),MENGE(ndim-2), * VEKSUM(ndim-2), VEKQU(ndim-2), ALTER(ndim-2) dimension ASPEI(11968),BSPEI(11968), CSPEI(12528) LOGICAL BOOL(ndim), BOOL1(ndim), HL1(ndim) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - UNENDL = 1000000000 call ilese(n, a, b, c, ndim) c print *,'0=no feasible sol., else input known value' c read *,olwert olwert = 0 if (olwert .eq. 0) olwert = unendl c WRITE(6,1000) N c DO 12 I=1,N c 12 WRITE(6,1010) (A(I,J),J=1,N) c WRITE(6,1015) c DO 13 I=1,N c 13 WRITE(6,1010) (B(I,J),J=1,N) c write(6,1016) c do 14 i=1,n c14 write(6,1010) (c(i,j),j=1,n) c write( 6,1011) olwert CALL QAP (N,A,B,UNENDL,LOESG,OLWERT,C,UMSPEI,ZUL,U,V,DD, * PARTPE,Y,LAB,Z1,H1,MENG,PHIOFM,MENGE,VEKSUM,VEKQU, * ALTER,ASPEI,BSPEI,CSPEI,BOOL,BOOL1,HL1,NDIM) WRITE(6,2000) DO 200 I=1,N 200 WRITE(6,2005) I, LOESG(I) WRITE(6,2010) OLWERT 102 FORMAT(I2) 103 FORMAT(20I3) 1000 FORMAT(20X,'N =',I3,///, 15X,17HDISTANCE MATRIX A,/) 1010 FORMAT(6X,20I4) 1011 format(6x,'a feasible solution with value=',i10,4x,'was supplied') 1015 FORMAT(//15X,19HCONNECTION MATRIX B,/) 1016 format(//15x,'LINEAR TERM:',/) 2000 FORMAT(///9X,15HEXACT SOLUTION:) 2005 FORMAT(9X,I5,4H -->,I2) 2010 FORMAT(/9X,25HOBJECTIVE FUNCTION VALUE:,I6) END SUBROUTINE QAP (N,A,B,UNENDL,LOESG,OLWERT,C,UMSPEI,ZUL,U,V,DD, * PARTPE,Y,LAB,Z1,H1,MENG,PHIOFM,MENGE,VEKSUM, * VEKQU,ALTER,ASPEI,BSPEI,CSPEI,BOOL,BOOL1,HL1,NDIM) C *** **************************************************************** C * * C * PROCEDURE FOR SOLVING QUADRATIC ASSIGNMENT * C * PROBLEMS * C * * C *** **************************************************************** C * * C * 1. CALL: * C * CALL QAP (N,A,B,UNENDL,LOESG,OLWERT,C,UMSPEI,ZUL,U,V,DD, * C * PARTPE,Y,LAB,Z1,H1,MENG,PHIOFM,MENGE,VEKSUM, * C * VEKQU,ALTER,ASPEI,BSPEI,CSPEI,BOOL,BOOL1,HL1) * C * * C * 2. COMPUTER CODE: * C * FORTRAN IV * C * * C * 3. METHOD: * C * PERTURBATION METHOD FOR THE OPTIMAL SOLUTION OF * C * QUADRATIC ASSIGNMENT PROBLEMS * C * * C * 4. PARAMETERS: * C * INPUT: * C * N DIMENSION OF THE PROBLEM * C * NDIM DIMENSION OF CALLING PROGRAM C * A(I,J) DISTANCE MATRIX (INTEGER) I,J=1,...,N * C * B(I,J) CONNECTION MATRIX (INTEGER) I,J=1,...,N * C * UNENDL LARGE MACHINE NUMBER (INTEGER) * C * OUTPUT: * C * LOESG(I) OPTIMAL ASSIGNMENT (INTEGER) * C * I --> LOESG(I) I=1,...,N * C * OLWERT OPTIMAL OBJECTIVE FUNCTION VALUE (INTEGER) * C * INTEGER ARRAYS OF DIMENSION (N,N) * C * C(I,J), ZUL(I,J) * C * INTEGER ARRAYS OF DIMENSION N * C * U(I), V(I), DD(I), PARTPE(I), Y(I), LAB(I), Z1(I), * C * H1(I) * C * INTEGER ARRAY OF DIMENSION N*N * C * UMSPEI(I) * C * INTEGER ARRAYS OF DIMENSION N-2 * C * MENG(I), PHIOFM(I), MENGE(I), VEKSUM(I), VEKQU(I), * C * ALTER(I) * C * INTEGER ARRAYS OF DIMENSION N*(N+1)*(2*N-2)/6 * C * ASPEI(I), BSPEI(I) * C * INTEGER ARRAY OF DIMENSION (N*(N+1)*(2*N+1)/6)-1 * C * CSPEI(I) * C * LOGICAL ARRAYS OF DIMENSION N * C * BOOL(I), BOOL1(I), HL1(I) * C * * C * 5. EXTERNAL SUBROUTINES: * C * SUBROUTINE ALTKOS * C * SUBROUTINE WEGSPE * C * SUBROUTINE PROGNO * C * SUBROUTINE LSAP * C * SUBROUTINE SSORT * C * * C<< C * 6. COMMENT: * C * THE ELEMENTS OF THE MATRICES A AND B HAVE TO BE * C * NONNEGATIVE INTEGERS. THEY ARE CHANGED DURING THE * C * PROCEDURE. * C * * C * 7. AUTHORS: * C * T. BOENNIGER * C * R.E. BURKARD * C * K.-H. STRATMANN * C * * C *** **************************************************************** C C<< IMPLICIT INTEGER (A - Z) INTEGER A(NDIM,1), B(NDIM,1), C(NDIM,1), UMSPEI(1), U(1), V(1), * DD(1), ZUL(NDIM,1), MENG(1), PHIOFM(1), PARTPE(1), * LOESG(1), Y(1), LAB(1), Z1(1), MENGE(1), ASPEI(1), * BSPEI(1), CSPEI(1), VEKSUM(1), VEKQU(1), ALTER(1), * H1(1) LOGICAL BOOL(1), BOOL1(1), HL1(1), WEITER C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** INITIALIZING OF PARAMETERS C K = 0 c *** OLWERT = UNENDL (here a user defined value can come in) NM2 = N-2 I2 = N + 1 J = 0 J1 = 0 DO 5035 I = 1,NM2 NMK = I2 - I J2 = NMK*NMK J = J + J2 - NMK VEKSUM(I) = J J1 = J1 + J2 5035 VEKQU(I) = J1 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** 1. REDUCTION: C(I,J) = A(I,I) * B(J,J) C DO 5010 I = 1,N BOOL(I) = .FALSE. BOOL1(I) = .FALSE. ZSTERN = A(I,I) DO 5010 J = 1,N ZUL(I,J) = -1 c consider linear term c 5010 C(I,J) = ZSTERN * B(J,J) + c( i, j) DO 5030 J = 1,N A(J,J) = UNENDL 5030 B(J,J) = UNENDL C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** REDUCTION OF MATRICES A AND B C DO 5103 I=1,N RA = A(I,1) RB = B(I,1) DO 5101 J=2,N RAA = A(I,J) RBB = B(I,J) IF(RAA .LT. RA) RA = RAA IF(RBB .LT. RB) RB = RBB 5101 CONTINUE DO 5102 J=1,N A(I,J) = A(I,J)-RA 5102 B(I,J) = B(I,J)-RB U(I) = RA 5103 V(I) = RB C<< C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** DETERMINATION OF THE MATRIX C C ROWWISE REDUCTION OF MATRICES A AND B C DO 5106 I=1,N CA1 = U(I) DO 5106 J=1,N BMJ = V(J) BM = (N-1)*BMJ DO 5104 KK=1,N IF(KK .EQ. J) GOTO 5104 BM = BM+B(J,KK) 5104 CONTINUE CA = CA1*BM AM = 0 DO 5105 KK=1,N IF(KK .EQ. I) GOTO 5105 AM = AM+A(I,KK) 5105 CONTINUE 5106 C(I,J) = CA+BMJ*AM+C(I,J) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COLUMNWISE REDUCTION OF MATRICES A AND B C DO 5109 I=1,N RA = A(1,I) RB = B(1,I) DO 5107 J=2,N RAA = A(J,I) RBB = B(J,I) IF(RAA .LT. RA) RA = RAA IF(RBB .LT. RB) RB = RBB 5107 CONTINUE DO 5108 J=1,N A(J,I) = A(J,I)-RA 5108 B(J,I) = B(J,I)-RB U(I) = RA 5109 V(I) = RB DO 5112 I=1,N A(I,I) = 0 B(I,I) = 0 CA1 = U(I) DO 5112 J=1,N BMJ = V(J) BM = (N-1)*BMJ DO 5110 KK=1,N IF(KK .EQ. J) GOTO 5110 BM = BM+B(KK,J) 5110 CONTINUE CA = CA1*BM AM = 0 DO 5111 KK=1,N IF(KK .EQ. I) GOTO 5111 AM = AM+A(KK,I) 5111 CONTINUE CCC = C(I,J)+CA+BMJ*AM C(I,J) = CCC 5112 CONTINUE ZPART = 0 NMK = N WEITER = .TRUE. C<< C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** IMPROVEMENT OF THE BOUNDS C(I,J) BY ADDING MINIMAL SCALAR C PRODUCTS (GILMORE BOUNDS). THE SUBROUTINES WEGSPE AND PROGNO C COMPUTE THESE MINIMAL SCALAR PRODUCTS. C CALL WEGSPE (N,K,NMK,IZAEHL,JZAEHL,A,B,VEKSUM,U,BOOL,BOOL1, * ASPEI,BSPEI,H1,NDIM) CALL PROGNO (N,K,NMK,UMSPEI,VEKSUM,VEKQU,WEITER,ASPEI,BSPEI, * CSPEI) WEITER = .FALSE. 5210 IZ = 0 I2 = 0 J2 = 0 IKAP = 0 DO 5120 I = 1,N IF (BOOL(I)) GOTO 5120 JJ = IZ*NMK IZ = IZ + 1 JZ = 0 DO 5125 J = 1,N IF (BOOL1(J)) GOTO 5125 JZ = JZ + 1 JJ = JJ+1 IF (ZUL(I,J) .LT. 0) GOTO 5145 UMSPEI(JJ) = UNENDL IKAP = IKAP + 1 IF (IKAP - 2) 5375, 5380, 5125 5375 IHALT = IZ JHALT = JZ GOTO 5125 5380 IF (IZ .EQ. IHALT) I2 = IZ IF (JZ .EQ. JHALT) J2 = JZ GOTO 5125 5145 UMSPEI(JJ) = C(I,J) + UMSPEI(JJ) 5125 CONTINUE 5120 CONTINUE C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF A BOUND BY SOLVING LSAP WITH COSTMATRIX UMSPEI C 5180 CALL LSAP (NMK,UNENDL,UMSPEI,ZSTERN,Y,Z1,LAB,DD,U,V,H1,HL1) JJ = 0 DO 5182 I=1,NMK DO 5181 J=1,NMK JJ = JJ+1 UMSPEI(JJ) = UMSPEI(JJ)-U(I)-V(J) 5181 CONTINUE 5182 CONTINUE C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** ZPART IS THE FIXED FRACTION OF THE OBJECTIVE FUNCTION VALUE C IMPLIED BY THE PRESENT PARTIAL PERMUTATION. C THE PRESENT BOUND IS ZPART+ZSTERN. C IF (ZPART+ZSTERN.GE.OLWERT) GOTO 5250 if (lbd .lt. zpart+zstern) then lbd = zpart + zstern c print 998, lbd endif 998 format(3x,'lb=',i12) IF (WEITER) GOTO 5220 5135 IF (IKAP - 1) 5400, 5410, 5420 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF THE ALTERNATIVE COSTS C 5400 CALL ALTKOS (NMK,UMSPEI,Z1,UNENDL,IZAEHL,JZAEHL,ALTERK) GOTO 5490 C<< C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF THE NEXT SINGLE ASSIGNMENT IN A FIXED ROW C OR COLUMN C 5410 MIN = UNENDL J1 = Z1(IHALT) JJ = (IHALT-1)*NMK DO 5430 J = 1,NMK JJ = JJ+1 T = UMSPEI(JJ) IF (MIN .LE. T .OR. J .EQ. J1) GOTO 5430 MIN = T 5430 CONTINUE MIN1 = MIN MIN = UNENDL JJ = J1 DO 5510 I = 1,NMK T = UMSPEI(JJ) JJ = JJ+NMK IF (MIN .LE. T .OR. I .EQ. IHALT) GOTO 5510 MIN = T 5510 CONTINUE MIN1 = MIN1 + MIN MIN = UNENDL JJ = JHALT DO 5440 I = 1,NMK T = UMSPEI(JJ) JJ = JJ+NMK IF (T .GE. MIN .OR. JHALT .EQ.Z1(I)) GOTO 5440 MIN = T 5440 CONTINUE MIN2 = MIN DO 5530 I = 1,NMK IF (Z1(I) .EQ. JHALT) GOTO 5540 5530 CONTINUE 5540 I1 = I JJ = (I1-1)*NMK MIN = UNENDL DO 5520 J = 1,NMK JJ = JJ+1 T = UMSPEI(JJ) IF (T .GE. MIN .OR. J .EQ. JHALT) GOTO 5520 MIN = T 5520 CONTINUE IF ((MIN+MIN2) .LT. MIN1) GOTO 5450 IZAEHL = I1 JZAEHL = JHALT ALTERK = MIN + MIN2 GOTO 5490 5450 IZAEHL = IHALT JZAEHL = J1 ALTERK = MIN1 GOTO 5490 C<< C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF THE NEXT SINGLE ASSIGNMENT IN THE PREVIOUSLY C FIXED ROW (RESP. PREVIOUSLY FIXED COLUMN) C 5420 IF (I2 .EQ. 0) GOTO 5480 IZAEHL = I2 JZAEHL = Z1(IZAEHL) GOTO 5370 5480 JZAEHL = J2 DO 5470 I = 1,NMK IF (Z1(I) .EQ. JZAEHL) GOTO 5485 5470 CONTINUE 5485 IZAEHL = I 5370 MIN = UNENDL JJ = (IZAEHL-1)*NMK DO 5350 I = 1,NMK JJ = JJ+1 T = UMSPEI(JJ) IF (T .GE. MIN .OR. I .EQ. JZAEHL) GOTO 5350 MIN = T 5350 CONTINUE ALTERK = MIN MIN = UNENDL JJ = JZAEHL DO 5360 J = 1,NMK T = UMSPEI(JJ) JJ = JJ+NMK IF (T .GE. MIN .OR. J .EQ. IZAEHL) GOTO 5360 MIN = T 5360 CONTINUE ALTERK = ALTERK + MIN 5490 IZ = 0 ALTER(K+1) = ALTERK + ZSTERN DO 5155 I = 1,N IF (BOOL(I)) GOTO 5155 IZ = IZ + 1 IF (IZAEHL .EQ. IZ) GOTO 5165 5155 CONTINUE 5165 IZAEHL = I IZ = 0 DO 5175 J = 1,N IF (BOOL1(J)) GOTO 5175 IZ = IZ + 1 IF (JZAEHL .EQ. IZ) GOTO 5185 5175 CONTINUE 5185 JZAEHL = J ZUL(IZAEHL,JZAEHL) =K WEITER = .TRUE. BOOL(IZAEHL) = .TRUE. BOOL1(JZAEHL) = .TRUE. NMK = N - K - 1 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF THE COST MATRIX CORRESPONDING TO THE NEW PARTIAL C PERMUTATION C CALL WEGSPE (N,K,NMK,IZAEHL,JZAEHL,A,B,VEKSUM,U,BOOL,BOOL1, * ASPEI,BSPEI,H1,NDIM) CALL PROGNO (N,K,NMK,UMSPEI,VEKSUM,VEKQU,WEITER,ASPEI,BSPEI, * CSPEI) C<< IZ = 0 DO 5130 I = 1,N IF (BOOL(I)) GOTO 5130 T1 = A(I,IZAEHL) T2 = A(IZAEHL,I) JZ = IZ DO 5140 J = 1,N IF (BOOL1(J)) GOTO 5140 JZ = JZ+1 UMSPEI(JZ) = C(I,J)+T1*B(J,JZAEHL)+T2*B(JZAEHL,J)+UMSPEI(JZ) 5140 CONTINUE IZ = IZ+NMK 5130 CONTINUE ZPART = ZPART + C(IZAEHL,JZAEHL) GOTO 5180 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** THE BOUND FOR THE NEW PARTIAL PERMUTATION IS NOT LESS THAN A C PREVIOUSLY FOUND OBJECTIVE FUNCTION VALUE. C - BACKTRACKING - C 5250 IF (.NOT. WEITER) GOTO 5255 WEITER = .FALSE. K=K+1 GOTO 5230 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** EXIT: THE SOLUTION TREE IS COMPLETELY FATHOMED. C 5255 IF (K. EQ. 0) GOTO 5600 IZAEHL = MENGE(K) JZAEHL = PHIOFM(K) 5220 DO 5150 I = 1,N IF (BOOL(I)) GOTO 5150 T1 = A(IZAEHL,I) T2 = A(I,IZAEHL) DO 5160 J = 1,N IF (BOOL1(J)) GOTO 5160 T = T1 * B(JZAEHL,J) + T2 * B(J,JZAEHL) IF (.NOT. WEITER) T = -T C(I,J) = C(I,J) + T 5160 CONTINUE 5150 CONTINUE IF (.NOT. WEITER) GOTO 5230 PARTPE(IZAEHL) = JZAEHL K = K + 1 MENGE(K) = IZAEHL PHIOFM(K) = JZAEHL IF (K .EQ. (N-2)) GOTO 5270 IKAP = 0 GOTO 5135 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** CANCELLATION OF THE LAST SINGLE ASSIGNMENT C 5230 DO 5260 I = 1,N IF (BOOL(I)) GOTO 5260 DO 5240 J = 1,N IF (.NOT. BOOL1(J) .AND. ZUL(I,J) .EQ. K) ZUL(I,J) = -1 5240 CONTINUE 5260 CONTINUE C<< ZPART = ZPART - C(IZAEHL,JZAEHL) BOOL(IZAEHL) = .FALSE. BOOL1(JZAEHL) = .FALSE. K = K - 1 NMK = N - K IF (ALTER(K+1)+ZPART.GE.OLWERT) GOTO 5255 CALL PROGNO (N,K,NMK,UMSPEI,VEKSUM,VEKQU,WEITER,ASPEI,BSPEI, * CSPEI) GOTO 5210 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF THE OBJECTIVE FUNCTION VALUES FOR THE REMAINING C TWO COMPLETE PERMUTATIONS C 5270 DO 5290 I = 1,N IF (BOOL(I)) GOTO 5290 IZ = I GOTO 5285 5290 CONTINUE 5285 DO 5280 I = 1,N IF (BOOL1(I)) GOTO 5280 J = I GOTO 5295 5280 CONTINUE 5295 BOOL(IZ) = .TRUE. BOOL1(J) = .TRUE. DO 5310 I = 1,N IF (BOOL(I)) GOTO 5310 I1 = I GOTO 5305 5310 CONTINUE 5305 DO 5300 I = 1,N IF(BOOL1(I)) GOTO 5300 J1 = I GOTO 5315 5300 CONTINUE 5315 WEITER = .FALSE. IZE = 0 5330 ZSTERN = C(IZ,J) + C(I1,J1) + A(IZ,I1)*B(J,J1) + A(I1,IZ)*B(J1,J) BOOL(IZ) = .FALSE. BOOL1(J) = .FALSE. IF ((ZSTERN+ZPART) .GE. OLWERT) GOTO 5340 OLWERT=ZSTERN+ZPART c print 999,olwert lbd = 0 999 format(20x,'new solution=',i12) DO 5320 I=1,N IF (BOOL(I)) LOESG(I) = PARTPE(I) 5320 CONTINUE DO 5335 I = 1,NM2 5335 MENG(I) = MENGE(I) LOESG(IZ) = J LOESG(I1) = J1 5340 IF(IZE.NE.0) GOTO 5255 IZE = IZ IZ = I1 I1 = IZE GOTO 5330 5600 CONTINUE RETURN END C<< SUBROUTINE ALTKOS (N,UMSPEI,Z1,UNENDL,IZAEHL,JZAEHL,ALTERK) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER UMSPEI(1), Z1(1) INTEGER A, ALTERK, UNENDL C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** DETERMINATION OF ALTERNATIVE COSTS AND THE RESULTING SINGLE C ASSIGNMENT. C ALTERK = -1 JJ = 0 DO 4110 I = 1,N J = Z1(I) J1 = J-N MIN = UNENDL DO 4120 I1 = 1,N J1 = J1+N IF (I1 .EQ. I) GOTO 4120 A = UMSPEI(J1) IF (A .LT. MIN) MIN = A 4120 CONTINUE MIN1 = MIN MIN = UNENDL DO 4130 I1 = 1,N JJ = JJ+1 IF(I1 .EQ. J) GOTO 4130 A = UMSPEI(JJ) IF (A .LT. MIN) MIN = A 4130 CONTINUE MIN = MIN + MIN1 IF (MIN .LE. ALTERK) GOTO 4110 IZAEHL = I JZAEHL = J ALTERK = MIN 4110 CONTINUE RETURN END C<< SUBROUTINE WEGSPE (N,K,NMK,IZAEHL,JZAEHL,A,B,VEKSUM,VEKT,BOOL, * BOOL1,ASPEI,BSPEI,H1,NDIM) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER A(NDIM,1), B(NDIM,1), VEKT(1), VEKSUM(1), ASPEI(1) INTEGER BSPEI(1), H1(1), T LOGICAL BOOL(1), BOOL1(1), LOGI C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** ORDERING OF THE ELEMENTS IN THE ROWS OF THE MATRIX A C DECREASINGLY AND IN MATRIX B INCREASINGLY. THE RESULTING C MATRICES ARE ROWWISE STORED ON THE VECTORS ASPEI AND BSPEI. C NMKM2 = NMK - 1 IF (NMK .NE. N) GOTO 2900 J1 = 1 DO 2810 I = 1,N IZ = 0 DO 2820 J = 1,N IF (J .EQ. I) GOTO 2820 IZ = IZ + 1 VEKT(IZ) = A(I,J) 2820 CONTINUE CALL SSORT (VEKT,H1,NMKM2) DO 2810 J = 1,NMKM2 NMKJ = NMK-J ASPEI(J1) = VEKT(NMKJ) 2810 J1 = J1 + 1 J1 = 1 DO 2840 I = 1,N IZ = 0 DO 2850 J = 1,N IF (J .EQ. I) GOTO 2850 IZ = IZ + 1 VEKT(IZ) = B(I,J) 2850 CONTINUE CALL SSORT (VEKT,H1,NMKM2) DO 2840 J = 1,NMKM2 BSPEI(J1) = VEKT(J) 2840 J1 = J1 + 1 RETURN 2900 NSUM1 = VEKSUM(K+1) J1 = NSUM1 J2 = NSUM1 - (NMK+1) * NMK DO 2910 I = 1,N IF (BOOL(I)) GOTO 2930 T = A(I,IZAEHL) LOGI = .TRUE. DO 2920 J = 1,NMK J2 = J2 + 1 IF (ASPEI(J2) .EQ. T .AND. LOGI) GOTO 2980 J1 = J1 + 1 ASPEI(J1) = ASPEI(J2) GOTO 2920 2980 LOGI = .FALSE. 2920 CONTINUE GOTO 2910 C<< 2930 IF (I .EQ. IZAEHL) J2 = J2 + NMK 2910 CONTINUE J1 = NSUM1 J2 = NSUM1 - (NMK+1) * NMK DO 2940 I = 1,N IF (BOOL1(I)) GOTO 2960 T = B(I,JZAEHL) LOGI = .TRUE. DO 2950 J = 1,NMK J2 = J2 + 1 IF (BSPEI(J2) .EQ. T .AND. LOGI) GOTO 2970 J1 = J1 + 1 BSPEI(J1) = BSPEI(J2) GOTO 2950 2970 LOGI = .FALSE. 2950 CONTINUE GOTO 2940 2960 IF (I .EQ. JZAEHL) J2 = J2 + NMK 2940 CONTINUE RETURN END C<< SUBROUTINE PROGNO (N,K,NMK,UMSPEI,VEKSUM,VEKQU,WEITER,ASPEI, * BSPEI,CSPEI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER UMSPEI(1), VEKSUM(1), VEKQU(1), ASPEI(1), BSPEI(1) INTEGER CSPEI(1), T LOGICAL WEITER C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF THE FINAL COST MATRIX UMSPEI FOR THE LINEAR C SUBPROBLEM LSAP ON THE K-TH STAGE OF THE DECISION TREE C JU = 0 IF (NMK .NE. N) GOTO 3600 NSUM2 = 0 J1 = 0 IF (WEITER) GOTO 3680 GOTO 3690 3600 IF (WEITER) GOTO 3630 J1 = VEKQU(K) 3690 DO 3640 I = 1,NMK DO 3640 J = 1,NMK J1 = J1 + 1 JU = JU + 1 3640 UMSPEI(JU) = CSPEI(J1) RETURN 3630 NSUM2 = VEKSUM(K+1) J1 = VEKQU(K+1) 3680 J2 = NSUM2 NMKM1 = NMK - 1 DO 3650 I = 1,NMK J3 = NSUM2 DO 3660 J = 1,NMK J1 = J1 + 1 JU = JU + 1 T = 0 DO 3670 IZ = 1,NMKM1 J2IZ = J2+IZ J3IZ = J3+IZ 3670 T = T+ASPEI(J2IZ)*BSPEI(J3IZ) J3 = J3 + NMKM1 CSPEI(J1) = T 3660 UMSPEI(JU) = T 3650 J2 = J2 + NMKM1 RETURN END SUBROUTINE LSAP(N,SUP,C,Z,ZEILE,SPALTE,DMINUS,DPLUS,YS,YT,VOR, 1 LABEL) C *** **************************************************************** C * * C * * C * LINEAR SUM ASSIGNMENT PROBLEM * C * * C *** **************************************************************** C * 1. CALL: * C * CALL LSAP(N,SUP,C,Z,ZEILE,SPALTE,DMINUS,DPLUS,YS,YT,VOR,* C * LABEL) * C * * C * 2. COMPUTER CODE: * C * FORTRAN IV * C * * C * 3. METHOD: * C * SHORTEST AUGMENTING PATH METHOD * C * * C * 4. PARAMETERS: * C * INPUT: * C * N DIMENSION OF THE COST MATRIX C * C * C(I) COST MATRIX (I=1,...N*N) , ROWWISE STORED * C * SUP LARGE MACHINE NUMBER * C * OUTPUT: * C * SPALTE(I) OPTIMAL ASSIGNMENT (I=1,...,N) * C * Z OPTIMAL VALUE * C * YS(I) OPTIMAL DUAL (ROW) VARIABLES (I=1,...,N)* C * YT(I) OPTIMAL DUAL (COLUMN) VARIABLES (I=1,...,N)* C * INTEGER ARRAYS OF LENGHTS N: * C * ZEILE(I),DMINUS(I),DPLUS(I),YS(I),YT(I),VOR(1) * C * LOGICAL ARRAY OF LENGTH N: * C * LABEL(I) * C * * C * 5. EXTERNAL SUBROUTINES: * C * NONE * C * * C * 6. AUTHOR: * C * U.DERIGS * C * * C *** **************************************************************** C C<< IMPLICIT INTEGER (A-Z) DIMENSION YS(1),YT(1),VOR(1) DIMENSION C(1),ZEILE(1),SPALTE(1),DMINUS(1),DPLUS(1) LOGICAL LABEL(1) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** STARTPROCEDURE C CONSTRUCTION OF AN INITIAL (PARTIAL) ASSIGNMENT C DO 50 I=1,N ZEILE(I)=0 SPALTE(I)=0 VOR(I)=0 YS(I)=0 50 YT(I)=0 IK=0 DO 2 I=1,N DO 3 J=1,N IK=IK+1 CC=C(IK) IF(J.EQ.1) GOTO 4 IF((CC-UI).GE.0) GOTO 3 4 UI=CC JO=J 3 CONTINUE YS(I)=UI IF(ZEILE(JO).NE.0) GOTO 2 ZEILE(JO)=I SPALTE(I)=JO 2 CONTINUE DO 5 J=1,N YT(J)=0 IF(ZEILE(J).EQ.0) YT(J)=SUP 5 CONTINUE IK=0 DO 6 I=1,N UI=YS(I) DO 7 J=1,N IK=IK+1 VJ=YT(J) IF(VJ.LE.0) GOTO 7 CC=C(IK)-UI IF(CC.GE.VJ) GOTO 7 YT(J)=CC VOR(J)=I 7 CONTINUE 6 CONTINUE DO 8 J=1,N I=VOR(J) IF(I.EQ.0) GOTO 8 IF(SPALTE(I).NE.0) GOTO 8 SPALTE(I)=J ZEILE(J)=I 8 CONTINUE C<< DO 9 I=1,N IF(SPALTE(I).NE.0) GOTO 9 UI=YS(I) IK=(I-1)*N DO 10 J=1,N IK=IK+1 IF(ZEILE(J).NE.0) GOTO 10 CC=C(IK) IF((CC-UI-YT(J)).GT.0) GOTO 10 SPALTE(I)=J ZEILE(J)=I GOTO 9 10 CONTINUE 9 CONTINUE C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** CONSTRUCTION OF THE OPTIMAL ASSIGNMENT C DO 1000 U=1,N IF(SPALTE(U).GT.0) GOTO 1000 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** SHORTEST PATH COMPUTATION C US=(U-1)*N DO 100 I=1,N VOR(I)=U LABEL(I)=.FALSE. DPLUS(I)=SUP USI=US+I 100 DMINUS(I)=C(USI)-YS(U)-YT(I) DPLUS(U)=0 105 D=SUP DO 110 I=1,N IF(LABEL(I)) GOTO 110 IF(DMINUS(I).GE.D) GOTO 110 D=DMINUS(I) INDEX=I 110 CONTINUE IF(ZEILE(INDEX).LE.0) GOTO 400 LABEL(INDEX)=.TRUE. W=ZEILE(INDEX) WS=(W-1)*N DPLUS(W)=D DO 130 I=1,N IF(LABEL(I)) GOTO 130 WSI=WS+I VGL=D+C(WSI)-YS(W)-YT(I) IF(DMINUS(I).LE.VGL) GOTO 130 DMINUS(I)=VGL VOR(I)=W 130 CONTINUE GOTO 105 C<< C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** AUGMENTATION C 400 W=VOR(INDEX) ZEILE(INDEX)=W IND=SPALTE(W) SPALTE(W)=INDEX IF(W.EQ.U) GOTO 500 INDEX=IND GOTO 400 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** TRANSFORMATION C 500 DO 510 I=1,N IF(DPLUS(I).EQ.SUP) GOTO 505 YS(I)=YS(I)+D-DPLUS(I) 505 IF(DMINUS(I).GE.D) GOTO 510 YT(I)=YT(I)+DMINUS(I)-D 510 CONTINUE 1000 CONTINUE C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C *** COMPUTATION OF THE OPTIMAL VALUE C Z=0 DO 2000 I=1,N IS=(I-1)*N J=SPALTE(I) ISJ=IS+J Z=Z+C(ISJ) 2000 CONTINUE RETURN END SUBROUTINE SSORT (A,B,L) C *** **************************************************************** C * * C * SORT - PROGRAM * C * * C *** **************************************************************** C * * C * 1. CALL: * C * CALL SSORT (A,B,N) * C * * C * 2. COMPUTER CODE: * C * FORTRAN IV * C * * C * 3. METHOD: * C * SORTING OF THE VECTOR A(I) IN INCREASING ORDER BY THE * C * SHELL-ALGORITHM. * C * * C * 4. PARAMETERS: * C * INPUT: * C * N DIMENSION OF THE VECTOR * C * A(I) VECTOR TO SORT (INTEGER) I=1,...,N * C * B(I) = I (INTEGER) I=1,...,N * C * OUTPUT: * C * A(I) THE SORTED VECTOR * C * B(I) PERMUTATION VECTOR OF THE SORTED VECTOR * C * * C * 5. EXTERNAL SUBROUTINES: * C * NONE * C * * C *** **************************************************************** C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IMPLICIT INTEGER (A-Z) DIMENSION A(1), B(1) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - F=1 IF(L.LE.F) RETURN N2 =(L-F+1)/2 S =1023 DO 100 T=1,10 IF(S.GT.N2) GOTO 90 LS=L-S DO 20 I=F,LS IS =I+S AH =A(IS) BH = B(IS) J=I JS =IS 5 IF(AH.GE.A(J)) GOTO 10 A(JS) =A(J) B(JS) = B(J) JS =J J =J-S IF(J.GE.F) GOTO 5 10 A(JS) =AH B(JS) = BH 20 CONTINUE 90 S=S/2 100 CONTINUE RETURN END subroutine ilese(n, a, b, c, ndim) integer a( ndim, ndim), b(ndim, ndim), c(ndim, ndim) read (5,*) n if (n .gt. ndim) then print *,'n zu gross: n, nmax=', n, ndim stop endif do 10 i=1, n 10 read (5,*) (a(i, j), j=1, n) do 20 i=1, n 20 read (5,*) (b(i, j), j=1, n) do 30 i=1, n 30 read (5,*, end = 40) (c(i, j), j=1, n) return 40 continue print *,'no linear term.' do 50 i=1,n do 50 j=1,n 50 c(i,j) = 0 return end