C ****************************************************************** C C FLOATING POINT AUCTION CODE FOR N BY N ASSIGNMENT PROBLEMS C C WRITTEN BY DIMITRI P. BERTSEKAS C (MODIFICATIONS BY DAVID CASTANON AND PAUL TSENG) C C VERSION 1.2, OCT. 1991 C C C MODIFIED BY DAVID CASTANON (OCT. 1991) TO WORK WITH ARBITRARILY C LARGE COST RANGE (UP TO THE INTEGER RANGE OF THE MACHINE); C TO ACCOMPLISH THIS, THE OBJECT PRICES AND EPSILON ARE C DOUBLE PRECISION VARIABLES, SO THE PRICES CANNOT OVERFLOW THE C MACHINE'S INTEGER RANGE. C C THIS CODE IMPLEMENTS THE AUCTION ALGORITHM WITH E-SCALING. C IT SOLVES A SEQUENCE OF SUBPROBLEMS AND DECREASES C EPSILON BY A CONSTANT FACTOR BETWEEN SUBPROBLEMS. C THIS VERSION CORRESPONDS TO A GAUSS-SEIDEL MODE C AND SOLVES EPSILON SUBPROBLEMS INEXACTLY. C C THE CODE IS AN IMPROVED VERSION OF AN EARLIER (SEPT. 1985) C AUCTION CODE WITH E-SCALING WRITTEN BY DIMITRI P. BERTSEKAS C C THE CODE TREATS THE PROBLEM AS A MAXIMIZATION PROBLEM. C TO SOLVE A MINIMIZATION PROBLEM, REVERSE THE SIGN OF THE C ARC COSTS PRIOR TO CALLING AUCTION, AND REVERSE AGAIN C THE SIGN OF THE OPTIMAL COST UPON RETURN FROM AUCTION. C THIS CODE ALLOWS MULTIPLE ARCS BETWEEN A ROW AND A COLUMN. C C THIS VERSION OF THE AUCTION ALGORITHM IS ADAPTIVE. HERE THE MINIMAL C BIDDING INCREMENT (STORED IN THE VARIABLE INCR) MAY BE SMALLER THAN C EPSILON. FOR EVERY SUBPROBLEM, INCR STARTS AT THE PARAMETER VALUE C STARTINCR (WHICH IS PASSED TO THE AUCTION ROUTINE) AND IS INCREASED C BY A FACTOR OF 2 AT THE END OF EACH CYCLE, UP TO A MAXIMUM VALUE OF C EPSILON. THIS ADAPTIVE FEATURE IS PARTICULARLY EFFECTIVE FOR C DENSE PROBLEMS. IT CAN BE DEFEATED BY SELECTING STARTINCR=BEGEPS C C ******************************************************************** C C THE USER MUST SUPPLY THE FOLLOWING PROBLEM DATA IN FORWARD STAR FORMAT C (THAT IS, ALL ARCS OF THE SAME ROW ARE NUMBERED CONSECUTIVELY): C N=NUMBER OF ROWS (EQUALS NUMBER OF COLUMNS) C A=NUMBER OF ARCS C FOUT(ROW)=FIRST ARC COMING OUT OF ROW C COST(ARC)=COST OF ARC C END(ARC)=COLUMN CORRESPONDING TO ARC C C AND THE FOLLOWING PARAMETERS FOR THE AUCTION ALGORITHM: C BEGEPS=STARTING VALUE OF EPSILON (MUST BE NO LESS THAN 1/(N+1)) C ENDEPS=FINAL VALUE OF EPSILON BEFORE IT IS SET TO 1/(N+1) C FACTOR=FACTOR BY WHICH EPSILON IS DECREASED BETWEEN SUBPROBLEMS C STARTINCR=THE STARTING VALUE OF THE BIDDING INCREMENT C ENDEPS SHOULD NOT EXCEED BEGEPS. C FACTOR MUST BE GREATER THAN 1. C C FOUT(.) IS AN ARRAY OF LENGTH N. C COST(.),END(.) ARE ARRAYS OF LENGTH A. C C THE SOLUTION IS CONTAINED IN THE ARRAY ASSIGN(.) WHERE C ASSIGN(COL) GIVES THE ROW ASSIGNED TO COL. C C THIS ALGORITHM DOES NOT CHECK FOR INFEASIBILITY OF THE PROBLEM. C TO MAKE SURE THE PROBLEM IS FEASIBLE THE USER MAY ADD C ADDITIONAL VERY SMALL COST ARCS. C C ****************************************************************** C ALL PROBLEM DATA ARE INTEGER C ****************************************************************** SUBROUTINE AUCTION(BEGEPS,FACTOR,ENDEPS,STARTINCR) PARAMETER(MAXNODES=10000, MAXARCS=100000) IMPLICIT NONE INTEGER NONEWLIST,THRESH INTEGER A,K,N,I,J,CURARC,CURCOL INTEGER ROW,FSTARC,FSTCOL,SNDARC,SNDCOL,BSTCOL INTEGER TRDARC INTEGER M,ISMALL,ILARGE,LARGEINCR,CYCLES INTEGER NUMPHASES,LSTARC,NOLIST,OLDROW INTEGER FOUT(MAXNODES),LIST(MAXNODES) INTEGER ASSIGN(MAXNODES) INTEGER COST(MAXARCS),END(MAXARCS) REAL*8 BEGEPS,ENDEPS,FACTOR,STARTINCR,INCRFACTOR REAL*8 PCOL(MAXNODES) REAL*8 MAX1,MAX2,TMAX,EPSILON,INCR REAL*8 AVERAGE COMMON/ARRAYC/COST/ARRAYS/END/ARRAYF/FOUT/BK1/N,A,ISMALL $ /BK2/CYCLES,AVERAGE,NUMPHASES/ARRAYA/ASSIGN/PCOLA/PCOL C *************************************************************** C ******* CHECK VALIDITY OF PARAMETERS PASSED ******* IF (BEGEPS.LE.1.0/(N+2)) THEN PRINT*,'STARTING VALUE OF EPSILON IS LESS THAN 1/(N+1)' PRINT*,'EXECUTION ABORTED' STOP END IF IF (ENDEPS.GT.BEGEPS) THEN PRINT*,'PARAMETER ENDEPS IS GREATER THAN PARAMETER BEGEPS' PRINT*,'ENDEPS IS SET AT THE DEFAULT VALUE OF 1/(N+1)' ENDEPS=1.0/(N+1) END IF IF ((FACTOR.LE.1).AND.(BEGEPS.GE.1.0/N)) THEN PRINT*,'EPSILON REDUCTION FACTOR IS NOT GREATER THAN 1' PRINT*,'EXECUTION ABORTED' STOP END IF IF (STARTINCR.LE.1.0/(N+2)) THEN PRINT*,'MIN BIDDING INCREMENT IS LESS THAN 1/(N+1)' PRINT*,'STARTINCR IS SET AT THE DEFAULT VALUE OF 1/(N+1)' STARTINCR=1.0/(N+1) END IF C ******* INITIALIZATION ******* EPSILON=BEGEPS ILARGE=-ISMALL LARGEINCR=INT(ILARGE/10) THRESH=INT(0.2*N) INCRFACTOR=2.0 IF (THRESH.GT.100) THRESH=100 X CYCLES=1 X AVERAGE=N NUMPHASES=1 DO 10 J=1,N ASSIGN(J)=0 PCOL(J)=ISMALL 10 CONTINUE FOUT(N+1)=A+1 NOLIST=N DO 20 I=1,N LIST(I)=I 20 CONTINUE C **************************************************************** C C THIS IMPLEMENTATION OF THE AUCTION ALGORITHM OPERATES IN CYCLES. C EACH CYCLE CONSISTS OF ONE BID BY EACH OF THE ROWS THAT ARE C UNASSIGNED AT THE START OF THE CYCLE (THESE ROWS ARE STORED IN C THE ARRAY LIST(.)). AS THE CYCLE PROGRESSES NEW C ROWS BECOME UNASSIGNED; THESE ARE STORED IN LIST(.) C AND WILL SUBMIT A BID AT THE NEXT CYCLE. C NOTE THAT ONLY ONE ROW SUBMITS A BID AT A TIME; THIS IS KNOWN AS C THE GAUSS-SEIDEL VERSION OF THE ALGORITHM. THE VERSION WHERE ALL C UNASSIGNED ROWS SUBMIT A BID AT THE SAME TIME IS KNOWN AS THE JACOBI C VERSION OF THE ALGORITHM, AND HAS NOT BEEN IMPLEMENTED. IT GENERALLY C TENDS TO RUN SOMEWHAT SLOWER THAN THE GAUSS-SEIDEL VERSION, BUT C IT ADMITS A HIGHER DEGREE OF PARALLELIZATION. A JACOBI VERSION C MAY BE PREFERABLE ON A PARALLEL MACHINE, BUT IS USUALLY INFERIOR C TO THE GAUSS-SEIDEL VERSION ON A SERIAL MACHINE. C C **************************************************************** C C START SUBPROBLEM (SCALING PHASE) W/ NEW EPSILON C C **************************************************************** 12 CONTINUE IF (EPSILON.LE.1.0/(N+1)) THRESH=0 INCR=STARTINCR IF (INCR.GT.EPSILON) INCR=EPSILON C **************************************************************** C C START AUCTION CYCLE WITH NEW LIST C C **************************************************************** 15 CONTINUE C INITIALIZE COUNT OF NEXT LIST OF UNASSIGNED ROWS NONEWLIST=0 C CYCLE THROUGH THE CURRENT LIST OF UNASSIGNED ROWS DO 100 I=1,NOLIST ROW=LIST(I) FSTARC=FOUT(ROW) LSTARC=FOUT(ROW+1)-1 FSTCOL=END(FSTARC) C FIRST TAKE CARE OF THE EXCEPTIONAL CASE WHERE ROW HAS ONLY ONE ARC IF (FSTARC.EQ.LSTARC) THEN PCOL(FSTCOL)=PCOL(FSTCOL)+LARGEINCR OLDROW=ASSIGN(FSTCOL) ASSIGN(FSTCOL)=ROW IF (OLDROW.GT.0) THEN NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF GO TO 100 END IF C NEXT TAKE CARE OF THE REGULAR CASE WHERE ROW HAS MULTIPLE ARCS SNDARC=FSTARC+1 SNDCOL=END(SNDARC) MAX1=COST(FSTARC)-PCOL(FSTCOL) MAX2=COST(SNDARC)-PCOL(SNDCOL) IF (MAX1.GE.MAX2) THEN BSTCOL=FSTCOL ELSE TMAX=MAX1 MAX1=MAX2 MAX2=TMAX BSTCOL=SNDCOL END IF IF (SNDARC.LT.LSTARC) THEN TRDARC=SNDARC+1 DO 40 CURARC=TRDARC,LSTARC CURCOL=END(CURARC) TMAX=COST(CURARC)-PCOL(CURCOL) IF (TMAX.GT.MAX2) THEN IF (TMAX.GT.MAX1) THEN MAX2=MAX1 MAX1=TMAX BSTCOL=CURCOL ELSE MAX2=TMAX END IF END IF 40 CONTINUE END IF C ROW BIDS FOR BSTCOL INCREASING ITS PRICE, AND GETS ASSIGNED C TO BSTCOL, WHILE ANY ROW ASSIGNED TO BSTCOL BECOMES UNASSIGNED PCOL(BSTCOL)=PCOL(BSTCOL)+MAX1-MAX2+INCR OLDROW=ASSIGN(BSTCOL) ASSIGN(BSTCOL)=ROW IF (OLDROW.GT.0) THEN NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF 100 CONTINUE C ************* END OF AN AUCTION CYCLE *************** C OPTIONALLY COLLECT STATISTICS X AVERAGE=(CYCLES*AVERAGE+NOLIST)/(CYCLES+1) X CYCLES=CYCLES+1 C CHECK IF THERE ARE STILL 'MANY' UNASSIGNED ROWS, THAT IS, IF THE C NUMBER OF UNASSIGNED ROWS IS GREATER THAN C THE PARAMETER THRESH. IF NOT, REPLACE CURRENT LIST WITH THE NEW LIST, C AND GO FOR ANOTHER CYCLE. OTHERWISE, IF EPSILON > 1, REDUCE EPSILON, C RESET THE ASSIGNMENT TO EMPTY AND RESTART AUCTION; C IF EPSILON = 1 TERMINATE. C ALSO INCREASE THE MINIMAL BIDDING INCREMENT UP TO A MAXIMUM C VALUE OF EPSILON (THIS IS THE ADAPTIVE FEATURE) INCR=INCR*INCRFACTOR IF (INCR.GT.EPSILON) INCR=EPSILON IF (NONEWLIST.GT.THRESH) THEN NOLIST=NONEWLIST GO TO 15 END IF C *************************************************************** C C END OF SUBPROBLEM (SCALING PHASE) C C *************************************************************** C ****** IF EPSILON IS 1 TERMINATE ****** IF (EPSILON.LE.1.0/N) THEN RETURN ELSE C ELSE REDUCE EPSILON AND RESET THE ASSIGNMENT TO EMPTY NUMPHASES=NUMPHASES+1 EPSILON=EPSILON/FACTOR IF (EPSILON.GT.INCR) EPSILON=EPSILON/FACTOR IF ((EPSILON.LT.1.0/N).OR.(EPSILON.LT.ENDEPS)) THEN EPSILON=1.0/(N+1) END IF THRESH=INT(THRESH/FACTOR) X PRINT*,'*** END OF A SCALING PHASE; NEW EPSILON=',EPSILON DO 200 J=1,N IF (ASSIGN(J).GT.0) THEN NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=ASSIGN(J) ASSIGN(J)=0 END IF 200 CONTINUE C FINAL PARAMETER UPDATES BEFORE RETURNING FOR ANOTHER SCALING PHASE NOLIST=NONEWLIST IF (STARTINCR.LT.EPSILON) STARTINCR=FACTOR*STARTINCR GO TO 12 END IF END SUBROUTINE SETRAN(ISEED) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C PORTABLE CONGRUENTIAL (UNIFORM) RANDOM NUMBER GENERATOR: C NEXT_VALUE = [(7**5) * PREVIOUS_VALUE] MODULO[(2**31)-1] C C THIS GENERATOR CONSISTS OF TWO ROUTINES: C (1) SETRAN - INITIALIZES CONSTANTS AND SEED C (2) RRAN - GENERATES A REAL RANDOM NUMBER C C THE GENERATOR REQUIRES A MACHINE WITH AT LEAST 32 BITS OF PRECISION. C THE SEED (ISEED) MUST BE IN THE RANGE (1,(2**31)-1). C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IF(ISEED.LT.1) STOP 77 MULT=16807 MODUL=2147483647 I15=2**15 I16=2**16 JRAN=ISEED RETURN END REAL FUNCTION RAN() IMPLICIT REAL*4 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C RAN GENERATES A REAL RANDOM NUMBER BETWEEN 0 AND 1 C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IXHI=JRAN/I16 IXLO=JRAN-IXHI*I16 IXALO=IXLO*MULT LEFTLO=IXALO/I16 IXAHI=IXHI*MULT IFULHI=IXAHI+LEFTLO IRTLO=IXALO-LEFTLO*I16 IOVER=IFULHI/I15 IRTHI=IFULHI-IOVER*I15 JRAN=((IRTLO-MODUL)+IRTHI*I16)+IOVER IF(JRAN.LT.0) JRAN=JRAN+MODUL RAN = FLOAT(JRAN)/FLOAT(MODUL) RETURN END AUCTION-AS This code implements the forward/reverse auction algorithm with $\e$-scaling for the asymmetric assignment problem; cf.\ Section 4.2.1. It can also solve as a special case the symmetric assignment problem. In this case as well as in the case where $\e$-scaling is bypassed ($\e=1$ throughout the algorithm), only the forward part of the algorithm is used. Note also that this code uses the ``third best'' object implementation described in Exercise 1.7 of Section 4.1. C ****************************************************************** C C SAMPLE CALLING PROGRAM FOR AUCTION ALGORITHM C FOR THE ASYMMETRIC ASSIGNMENT PROBLEM C C THIS DRIVER CREATES AN ASYMMETRIC ASSIGNMENT PROBLEM C WITH LESS OR EQUAL NUMBER OF ROWS THAN COLUMNS, C AND CALLS THE AUCTION_AS SUBROUTINE TO FIND AN C ASSIGNMENT OF MAXIMAL VALUE. C C THIS VERSION USES A THIRD BEST VALUE C C ****************************************************************** PARAMETER(MAXNODES=10000, MAXARCS=100000) IMPLICIT NONE INTEGER M,N,NA,A,IA,K,ILARGE,BEGEPS,ENDEPS,CYCLES INTEGER NUMPHASES,STARTINCR,LAMBDA INTEGER I,J,ARC,NOASS,ICOST,ABSCOST,CURARC INTEGER CURCOL,FSTARC,LSTARC,MINCOST,MAXCOST,MCOST INTEGER EXTRA,REMAINDER,INDEX,COUNT INTEGER COL_ASSIGNED_TO(MAXNODES),ROW_ASSIGNED_TO(MAXNODES) INTEGER FOUT(MAXNODES),FIN(MAXNODES),NXTIN(MAXARCS) INTEGER COST(MAXARCS),START(MAXARCS),END(MAXARCS) INTEGER PCOL(MAXNODES),PROW(MAXNODES) INTEGER PRDARC(MAXNODES) REAL*8 FACTOR,TT1,TT2,TCOST,AVERAGE COMMON/ARRAYC/COST/ARRAYS/START/ARRAYE/END COMMON/ARRAYFO/FOUT/ARRAYFI/FIN/ARRAYNI/NXTIN COMMON/ARRAYRA/ROW_ASSIGNED_TO/ARRAYCA/COL_ASSIGNED_TO COMMON/ARRAYPC/PCOL/ARRAYPR/PROW COMMON/BK1/M,N,A,ILARGE COMMON/BK2/CYCLES,AVERAGE,NUMPHASES,LAMBDA C **************************************************************** C C PROBLEM GENERATION CODE STARTS HERE C THE USER MAY REPLACE THIS CODE WITH A CODE THAT READS C HIS/HER PROBLEM FROM A FILE C C **************************************************************** C THIS CODE INCLUDES A UNIFORM RANDOM NUMBER GENERATOR C WHICH RETURNS A VALUE IN (0,1) C INITIALIZE RANDOM GENERATOR INTEGER MULT,MODUL,I15,I16,JRAN,ISEED REAL RAN COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN ISEED=13502460 CALL SETRAN(ISEED) PRINT*,'GENERATING AN ASYMMETRIC ASSIGNMENT PROBLEM' PRINT*,'********************************************' C **** READ THE NUMBER OF ROWS M, COLUMNS N, AND ARCS A **** PRINT*,'ENTER THE NUMBER OF ROWS' READ*,M 2 PRINT*,'ENTER THE NUMBER OF COLUMNS (>= # OF ROWS)' READ*,N IF (N.LT.M) GO TO 2 5 PRINT*,'ENTER THE NUMBER OF ARCS PER ROW (>1)' READ*,NA IF (NA.LT.2) GOTO 5 PRINT*,'ENTER THE MINIMUM AND THE MAXIMUM COST' READ*,MINCOST,MAXCOST C THE NUMBER OF ARCS IS M*NA+N-M A=M*NA+N-M PRINT*,'THE NUMBER OF ARCS IS ',A C THE ARCS INCIDENT TO ROW I ARE FOUT(I) TO FOUT(I+1)-1 C ALSO, FOR FEASIBILITY EACH ROW IS DIRECTLY CONNECTED C WITH THE CORRESPONDING COLUMN; C ALSO EACH COLUMN HAS AT LEAST ONE ARC EXTRA=INT((N-M)/M) DO 20 I=1,M FOUT(I)=1+(I-1)*(NA+EXTRA) 20 CONTINUE FOUT(M+1)=A+1 DO 22 I=1,M DO 23 ARC=FOUT(I),FOUT(I+1)-1 START(ARC)=I 23 CONTINUE 22 CONTINUE C GENERATE THE END(ARC) AND COST(ARC) WHICH ARE THE COLUMN C AND THE COST COEFFICIENT ASSOCIATED WITH ARC DO 25 ARC=1,A END(ARC)=1+RAN()*N IF ((END(ARC).GT.N).OR.(END(ARC).LT.1)) THEN PRINT*,'ERROR IN PROBLEM GENERATION' PAUSE STOP END IF COST(ARC)=MINCOST+RAN()*(MAXCOST-MINCOST) 25 CONTINUE C MODIFY THE END OF THE LAST ARC OUT OF EACH ROW FOR FEASIBILITY C AND SET ITS COST TO -MAXCOST DO 30 I=1,M END(FOUT(I+1)-1)=I COST(FOUT(I+1)-1)=-MAXCOST 30 CONTINUE C MODIFY THE END OF SOME ARCS SO THAT EACH COLUMN HAS AT LEAST ONE ARC IF (EXTRA.GE.1) THEN DO 32 I=1,M DO 33 K=1,EXTRA END(FOUT(I)+K-1)=K*M+I 33 CONTINUE 32 CONTINUE END IF REMAINDER=N-M*(1+EXTRA) IF (REMAINDER.GE.1) THEN INDEX=FOUT(M)+EXTRA-1 DO 35 J=1,REMAINDER END(INDEX+J)=(1+EXTRA)*M+J 35 CONTINUE END IF C CONSTRUCT THE FIN() AND NXTIN() ARRAYS DO 42 J=1,N FIN(J)=0 PRDARC(J)=0 42 CONTINUE DO 43 ARC=1,A NXTIN(ARC)=0 J=END(ARC) IF (FIN(J).NE.0) THEN NXTIN(PRDARC(J))=ARC ELSE FIN(J)=ARC END IF PRDARC(J)=ARC 43 CONTINUE C **************************************************************** C C PROBLEM GENERATION CODE ENDS HERE C C **************************************************************** C SCALE THE COST TO WORK WITH INTEGER EPSILON MAXCOST=0 DO 45 IA=1,A ABSCOST=IABS(COST(IA)) IF (ABSCOST.GT.MAXCOST) MAXCOST=ABSCOST COST(IA)=COST(IA)*(M+1) 45 CONTINUE C *** ILARGE IS A VERY LARGE INTEGER FOR YOUR MACHINE *** ILARGE=2000000000 IF (MAXCOST.GT.INT(ILARGE/(M+1))) THEN PRINT*,'THE COST RANGE IS TOO LARGE FOR INTEGER ARITHMETIC' PAUSE STOP END IF MAXCOST=MAXCOST*(M+1) LAMBDA=-ILARGE C THE FOLLOWING PARAMETERS BEGEPS, FACTOR, ENDEPS, AND STARTINCR C ARE PASSED TO THE AUCTION ALGORITHM. VALUES BETWEEN C (A) MAXCOST/5 AND MAXCOST/2 FOR BEGEPS C (B) 4 AND 6 FOR FACTOR C (C) M/10 AND 1 FOR ENDEPS C (D) 1 AND BEGEPS FOR STARTINCR C HAVE WORKED WELL FOR LARGE SPARSE PROBLEMS. C FOR DENSE PROBLEMS AND FOR VERY ASYMMETRIC PROBLEMS C IT IS RECOMMENDED THAT C BEGEPS BE SET TO A SMALLER VALUE (POSSIBLY 1), C ENDEPS BE SET TO 1, C STARTINCR BE SET TO 1. C FOR VERY ASYMMETRIC PROBLEMS, BEGEPS=1, CORRESPONDING TO C NO E-SCALING, MAY WORK BEST AND SHOULD BE AT LEAST TRIED. PRINT*,'MAXIMUM COST IS ',MAXCOST PRINT*,'ENTER THE STARTING EPSILON' READ*,BEGEPS IF (BEGEPS.LT.1) BEGEPS=1 PRINT*,'ENTER THE EPSILON REDUCTION FACTOR' READ*,FACTOR ENDEPS=M/10 IF (ENDEPS.LT.1) ENDEPS=1 IF (ENDEPS.GT.BEGEPS) ENDEPS=BEGEPS STARTINCR=BEGEPS/10 IF (STARTINCR.LT.1) STARTINCR=1 PRINT*,'********************************************' PRINT*,'STARTING EPSILON = ',BEGEPS PRINT*,'EPSILON REDUCTION FACTOR = ',FACTOR PRINT*,'THRESHOLD EPSILON BEFORE IT IS SET TO 1 = ',ENDEPS PRINT*,'STARTING MIN BIDDING INCREMENT = ',STARTINCR PRINT*,'CALLING ASYMMETRIC ASSIGNMENT AUCTION' PRINT*,'********************************************' C GET STARTING TIME FOR THE MAC II TT1 = LONG(362)/60.0 CALL AUCTION_AS(BEGEPS,FACTOR,ENDEPS,STARTINCR) C GET ENDING TIME FOR THE MAC II TT2 = LONG(362)/60.0 - TT1 PRINT *,'FINISHED --- TOTAL CPU TIME', TT2,' SECS' PRINT*,'********************************************' C *** DISPLAY RESULTS *** X WRITE(9,2010) CYCLES X2010 FORMAT(' NO OF AUCTION CYCLES',I7) X WRITE(9,2020) AVERAGE X2020 FORMAT(' AVERAGE NUMBER OF BIDS PER CYCLE',F9.3) WRITE(9,2030) NUMPHASES 2030 FORMAT('NO OF EPSILON SUBPROBLEMS SOLVED =',I7) C CHECK OPTIMALITY & CALCULATE COST DO 48 J=1,N I=ROW_ASSIGNED_TO(J) IF (N.GT.M) THEN ROW_ASSIGNED_TO(J)=0 ELSE IF (I.GT.0) THEN COL_ASSIGNED_TO(I)=J END IF END IF 48 CONTINUE TCOST=0 DO 50 I=1,M J=COL_ASSIGNED_TO(I) ROW_ASSIGNED_TO(J)=I IF (J.EQ.0) THEN PRINT*,'ROW ',I,' IS UNASSIGNED' END IF IF (PCOL(J).LT.LAMBDA) THEN PRINT*,'CS VIOLATION AT ASSIGNED COLUMN ',J END IF FSTARC=FOUT(I) LSTARC=FOUT(I+1)-1 MCOST=-ILARGE DO 55 ARC=FSTARC,LSTARC CURCOL=END(ARC) IF (PROW(I)+PCOL(CURCOL).LT.COST(ARC)-1) THEN PRINT*,'1-CS VIOLATED AT ARC ',ARC END IF IF (CURCOL.EQ.J) THEN IF (MCOST.LT.COST(ARC)) THEN MCOST=COST(ARC) END IF END IF 55 CONTINUE TCOST=TCOST+MCOST/(M+1) IF (PROW(I)+PCOL(J).NE.MCOST) THEN PRINT*,'1-CS VIOLATED AT ROW ',I END IF 50 CONTINUE COUNT=0 DO 60 J=1,N IF (ROW_ASSIGNED_TO(J).EQ.0) THEN IF (PCOL(J).GT.LAMBDA) THEN PRINT*,'CS VIOLATION AT UNASSIGNED COLUMN ',J END IF ELSE COUNT=COUNT+1 END IF 60 CONTINUE IF (COUNT.LT.M) THEN PRINT*,'THE NUMBER OF ASSIGNED COLUMNS IS WRONG' END IF WRITE(9,2100) TCOST 2100 FORMAT(' ASSIGNMENT COST=',F18.2) PRINT *, ' PROGRAM ENDED; TO EXIT ' PAUSE END C ****************************************************************** C C AUCTION CODE FOR ASYMMETRIC M BY N ASSIGNMENT PROBLEMS C C WRITTEN BY DIMITRI P. BERTSEKAS C C DEC. 1990 C C THIS CODE IMPLEMENTS A FORWARD/REVERSE AUCTION ALGORITHM C WITH E-SCALING FOR THE ASYMMETRIC ASSIGNMENT PROBLEM. C IT SOLVES A SEQUENCE OF SUBPROBLEMS AND DECREASES C EPSILON BY A CONSTANT FACTOR BETWEEN SUBPROBLEMS. C THIS VERSION CORRESPONDS TO A GAUSS-SEIDEL MODE. C C THE CODE TREATS THE PROBLEM AS A MAXIMIZATION PROBLEM. C TO SOLVE A MINIMIZATION PROBLEM, REVERSE THE SIGN OF THE C ARC COSTS PRIOR TO CALLING AUCTION, AND REVERSE AGAIN C THE SIGN OF THE OPTIMAL COST UPON RETURN FROM AUCTION. C THIS CODE ALLOWS MULTIPLE ARCS BETWEEN A ROW AND A COLUMN. C C THIS VERSION OF THE AUCTION ALGORITHM IS ADAPTIVE. HERE THE MINIMAL C BIDDING INCREMENT (STORED IN THE VARIABLE INCR) MAY BE SMALLER THAN C EPSILON. FOR EVERY SUBPROBLEM, INCR STARTS AT THE PARAMETER VALUE C STARTINCR (WHICH IS PASSED TO THE AUCTION ROUTINE) AND IS INCREASED C BY A FACTOR OF 2 AT THE END OF EACH CYCLE, UP TO A MAXIMUM VALUE OF C EPSILON. THIS ADAPTIVE FEATURE IS PARTICULARLY EFFECTIVE FOR C DENSE PROBLEMS. IT CAN BE DEFEATED BY SELECTING STARTINCR=BEGEPS C C ******************************************************************** C C THE USER MUST SUPPLY THE FOLLOWING PROBLEM DATA IN FORWARD STAR FORMAT C (THAT IS, ALL ARCS OF THE SAME ROW ARE NUMBERED CONSECUTIVELY): C M=NUMBER OF ROWS C N=NUMBER OF COLUMNS C A=NUMBER OF ARCS C FOUT(ROW)=FIRST ARC COMING OUT OF ROW C FIN(COL)=FIRST ARC COMING INTO COL C NXTIN(ARC)=NEXT ARC INCIDENT TO THE SAME COLUMN AS ARC C COST(ARC)=COST OF ARC C START(ARC)=ROW CORRESPONDING TO ARC C END(ARC)=COLUMN CORRESPONDING TO ARC C C AND THE FOLLOWING PARAMETERS FOR THE AUCTION ALGORITHM: C BEGEPS=STARTING VALUE OF EPSILON (MUST BE NO LESS THAN 1) C ENDEPS=FINAL VALUE OF EPSILON BEFORE IT IS SET TO 1 C FACTOR=FACTOR BY WHICH EPSILON IS DECREASED BETWEEN SUBPROBLEMS C STARTINCR=THE STARTING VALUE OF THE BIDDING INCREMENT C ENDEPS SHOULD NOT EXCEED BEGEPS. C FACTOR MUST BE GREATER THAN 1. C C FOUT(.) IS AN ARRAY OF LENGTH N. C COST(.),END(.) ARE ARRAYS OF LENGTH A. C C THE SOLUTION IS CONTAINED IN THE ARRAY ROW_ASSIGNED_TO(.) WHERE C ROW_ASSIGNED_TO(COL) GIVES THE ROW ASSIGNED TO COL. C ALSO COL_ASSIGNED_TO(ROW) GIVES THE COLUMN ASSIGNED TO ROW. C C THIS ALGORITHM DOES NOT CHECK FOR INFEASIBILITY OF THE PROBLEM. C TO MAKE SURE THE PROBLEM IS FEASIBLE THE USER MAY ADD C ADDITIONAL VERY SMALL COST ARCS. C C THIS CODE MAY FAIL DUE TO INTEGER OVERFLOW IF THE NUMBER OF NODES C OR THE COST RANGE (OR BOTH) ARE LARGE. TO CORRECT THIS SITUATION, C THE PRICES AND OTHER RELATED VARIABLES (MAX1,MAX2,TMAX ETC) SHOULD C BE DECLARED AS DOUBLE PRECISION REALS. C ****************************************************************** C ALL PROBLEM DATA ARE INTEGER C ****************************************************************** SUBROUTINE AUCTION_AS(BEGEPS,FACTOR,ENDEPS,STARTINCR) PARAMETER(MAXNODES=10000, MAXARCS=100000) IMPLICIT NONE INTEGER A,K,N,I,J,M,CURARC,CURCOL,CURROW,LAMBDA,DELTA INTEGER THRESH,INCR,STARTINCR,INCRFACTOR,EPS_THRESH INTEGER NOLIST,RNOLIST,NONEWLIST,RNONEWLIST INTEGER ROW,COLUMN, BSTROW,BSTCOL INTEGER FSTARC,FSTCOL,NXTARC,NXTROW,LSTARC,SNDARC,SNDCOL INTEGER MAX1,MAX2,TMAX,TMIN,TRDARC,EPSILON,BEGEPS,ENDEPS INTEGER ISMALL,ILARGE,LARGEINCR,CYCLES INTEGER NUMPHASES,OLDROW,OLDCOL INTEGER MAX3,BSTARC,SBSTARC INTEGER CUR_BARC,TRDVAL(MAXNODES) LOGICAL INIT INTEGER BEST_ARC(MAXNODES),SECD_ARC(MAXNODES) INTEGER COL_ASSIGNED_TO(MAXNODES),ROW_ASSIGNED_TO(MAXNODES) INTEGER FOUT(MAXNODES),FIN(MAXNODES),NXTIN(MAXARCS) INTEGER COST(MAXARCS),START(MAXARCS),END(MAXARCS) INTEGER LIST(MAXNODES),REV_LIST(MAXNODES) INTEGER PCOL(MAXNODES),PROW(MAXNODES) REAL*8 AVERAGE,FACTOR COMMON/ARRAYC/COST/ARRAYS/START/ARRAYE/END COMMON/ARRAYFO/FOUT/ARRAYFI/FIN/ARRAYNI/NXTIN COMMON/ARRAYRA/ROW_ASSIGNED_TO/ARRAYCA/COL_ASSIGNED_TO COMMON/ARRAYPC/PCOL/ARRAYPR/PROW COMMON/BK1/M,N,A,ILARGE COMMON/BK2/CYCLES,AVERAGE,NUMPHASES,LAMBDA C *************************************************************** C ******* CHECK VALIDITY OF PARAMETERS PASSED ******* IF (BEGEPS.LT.1) THEN PRINT*,'STARTING VALUE OF EPSILON IS LESS THAN 1' PRINT*,'EXECUTION ABORTED' STOP END IF IF (ENDEPS.GT.BEGEPS) THEN PRINT*,'PARAMETER ENDEPS IS GREATER THAN PARAMETER BEGEPS' PRINT*,'ENDEPS IS SET AT THE DEFAULT VALUE OF 1' ENDEPS=1 END IF IF ((FACTOR.LE.1).AND.(BEGEPS.GT.1)) THEN PRINT*,'EPSILON REDUCTION FACTOR IS NOT GREATER THAN 1' PRINT*,'EXECUTION ABORTED' STOP END IF IF (STARTINCR.LT.1) THEN PRINT*,'MIN BIDDING INCREMENT IS LESS THAN 1' PRINT*,'STARTINCR IS SET AT THE DEFAULT VALUE OF 1' STARTINCR=1 END IF C ******* INITIALIZATION ******* EPSILON=BEGEPS ISMALL=-ILARGE LARGEINCR=INT(ILARGE/10) THRESH=INT(0.2*M) INCRFACTOR=2 EPS_THRESH=BEGEPS/(FACTOR**3) IF (EPS_THRESH.LT.2) EPS_THRESH=2 INIT=.FALSE. IF (THRESH.GT.100) THRESH=100 X CYCLES=1 X AVERAGE=N NUMPHASES=1 DO 10 I=1,N PCOL(I)=ISMALL 10 CONTINUE FOUT(M+1)=A+1 C **************************************************************** C C THIS IMPLEMENTATION OF THE AUCTION ALGORITHM OPERATES IN CYCLES. C EACH CYCLE CONSISTS OF ONE BID BY EACH OF THE ROWS THAT ARE C UNASSIGNED AT THE START OF THE CYCLE (THESE ROWS ARE STORED IN C THE ARRAY LIST(.)). AS THE CYCLE PROGRESSES NEW C ROWS BECOME UNASSIGNED; THESE ARE STORED IN LIST(.) C AND WILL SUBMIT A BID AT THE NEXT CYCLE. C C THE ALGORITHM FIRST FINDS A FEASIBLE ASSIGNMENT, WHERE ALL PERSONS C ARE ASSIGNED, USING A COMBINED FORWARD/REVERSE AUCTION. C THEN THE ALGORITHM USES A MODIFIED REVERSE AUCTION TO TO SATISFY C THE REMAINING E-CS CONDITIONS. C C **************************************************************** C C START SUBPROBLEM (SCALING PHASE) W/ NEW EPSILON C C **************************************************************** 12 CONTINUE C INITIALIZE ROW ASSIGNMENT LISTS NOLIST=M DO 20 I=1,M LIST(I)=I 20 CONTINUE DO 22 J=1,N ROW_ASSIGNED_TO(J)=0 22 CONTINUE INCR=STARTINCR IF (INCR.GT.EPSILON) INCR=EPSILON IF (EPSILON.EQ.1) THRESH=0 C **************************************************************** C C START FORWARD AUCTION CYCLE C C **************************************************************** IF (EPSILON.LT.EPS_THRESH) THEN 17 CONTINUE C INITIALIZE COUNT OF NEXT LIST OF UNASSIGNED ROWS NONEWLIST=0 C CYCLE THROUGH THE CURRENT LIST OF UNASSIGNED ROWS DO 103 I=1,NOLIST ROW=LIST(I) FSTARC=FOUT(ROW) LSTARC=FOUT(ROW+1)-1 FSTCOL=END(FSTARC) C FIRST TAKE CARE OF THE EXCEPTIONAL CASE WHERE ROW HAS ONLY ONE ARC IF (FSTARC.EQ.LSTARC) THEN PCOL(FSTCOL)=PCOL(FSTCOL)+LARGEINCR PROW(ROW)=COST(FSTARC)-PCOL(FSTCOL) OLDROW=ROW_ASSIGNED_TO(FSTCOL) ROW_ASSIGNED_TO(FSTCOL)=ROW IF (OLDROW.GT.0) THEN NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF GO TO 103 END IF C NEXT TAKE CARE OF THE REGULAR CASE WHERE ROW HAS MULTIPLE ARCS IF (((FSTARC+2).LE.LSTARC).AND.(INIT)) THEN BSTARC=BEST_ARC(ROW) SBSTARC=SECD_ARC(ROW) BSTCOL=END(BSTARC) SNDCOL=END(SECD_ARC(ROW)) MAX1=COST(BSTARC)-PCOL(BSTCOL) MAX2=COST(SBSTARC)-PCOL(SNDCOL) IF((MAX1.GE.TRDVAL(ROW)).AND.(MAX2.GE.TRDVAL(ROW))) THEN IF (MAX1.LT.MAX2) THEN TMAX=MAX1 MAX1=MAX2 MAX2=TMAX BSTCOL=SNDCOL END IF GO TO 70 END IF END IF SNDARC=FSTARC+1 SNDCOL=END(SNDARC) MAX1=COST(FSTARC)-PCOL(FSTCOL) MAX2=COST(SNDARC)-PCOL(SNDCOL) IF (MAX1.GE.MAX2) THEN BSTARC=FSTARC SBSTARC=SNDARC ELSE TMAX=MAX1 MAX1=MAX2 MAX2=TMAX BSTARC=SNDARC SBSTARC=FSTARC END IF IF (SNDARC.LT.LSTARC) THEN TRDARC=SNDARC+1 MAX3=ISMALL DO 43 CURARC=TRDARC,LSTARC CURCOL=END(CURARC) TMAX=COST(CURARC)-PCOL(CURCOL) IF (TMAX.GT.MAX2) THEN IF (TMAX.GT.MAX1) THEN SBSTARC=BSTARC BSTARC=CURARC MAX3=MAX2 MAX2=MAX1 MAX1=TMAX ELSE SBSTARC=CURARC MAX3=MAX2 MAX2=TMAX END IF ELSE IF (MAX3.LT.TMAX) MAX3=TMAX END IF 43 CONTINUE END IF BEST_ARC(ROW)=BSTARC BSTCOL=END(BSTARC) SECD_ARC(ROW)=SBSTARC TRDVAL(ROW)=MAX3 70 CONTINUE C ROW BIDS FOR BSTCOL INCREASING ITS PRICE, AND GETS ASSIGNED C TO BSTCOL, WHILE ANY ROW ASSIGNED TO BSTCOL BECOMES UNASSIGNED PCOL(BSTCOL)=PCOL(BSTCOL)+MAX1-MAX2+INCR PROW(ROW)=MAX2-INCR OLDROW=ROW_ASSIGNED_TO(BSTCOL) ROW_ASSIGNED_TO(BSTCOL)=ROW IF (OLDROW.GT.0) THEN NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF 103 CONTINUE C ************* END OF A FORWARD AUCTION CYCLE *************** C OPTIONALLY COLLECT STATISTICS X AVERAGE=(CYCLES*AVERAGE+NOLIST)/(CYCLES+1) X CYCLES=CYCLES+1 C CHECK IF A SWITCH TO MODIFIED REVERSE AUCTION C SHOULD BE MADE (IF NONEWLIST EQUALS ZERO). C ALSO INCREASE THE MINIMAL BIDDING INCREMENT UP TO A MAXIMUM C VALUE OF EPSILON (THIS IS THE ADAPTIVE FEATURE) INCR=INCR*INCRFACTOR IF (INCR.GT.EPSILON) INCR=EPSILON IF (NONEWLIST.GT.THRESH) THEN NOLIST=NONEWLIST INIT=.TRUE. GO TO 17 END IF ELSE C **************************************************************** C C START FORWARD AUCTION CYCLE C C **************************************************************** 15 CONTINUE C INITIALIZE COUNT OF NEXT LIST OF UNASSIGNED ROWS NONEWLIST=0 C CYCLE THROUGH THE CURRENT LIST OF UNASSIGNED ROWS DO 100 I=1,NOLIST ROW=LIST(I) FSTARC=FOUT(ROW) LSTARC=FOUT(ROW+1)-1 FSTCOL=END(FSTARC) C FIRST TAKE CARE OF THE EXCEPTIONAL CASE WHERE ROW HAS ONLY ONE ARC IF (FSTARC.EQ.LSTARC) THEN PCOL(FSTCOL)=PCOL(FSTCOL)+LARGEINCR PROW(ROW)=COST(FSTARC)-PCOL(FSTCOL) OLDROW=ROW_ASSIGNED_TO(FSTCOL) ROW_ASSIGNED_TO(FSTCOL)=ROW IF (OLDROW.GT.0) THEN NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF GO TO 100 END IF C NEXT TAKE CARE OF THE REGULAR CASE WHERE ROW HAS MULTIPLE ARCS SNDARC=FSTARC+1 SNDCOL=END(SNDARC) MAX1=COST(FSTARC)-PCOL(FSTCOL) MAX2=COST(SNDARC)-PCOL(SNDCOL) IF (MAX1.GE.MAX2) THEN BSTCOL=FSTCOL ELSE TMAX=MAX1 MAX1=MAX2 MAX2=TMAX BSTCOL=SNDCOL END IF IF (SNDARC.LT.LSTARC) THEN TRDARC=SNDARC+1 DO 40 CURARC=TRDARC,LSTARC CURCOL=END(CURARC) TMAX=COST(CURARC)-PCOL(CURCOL) IF (TMAX.GT.MAX2) THEN IF (TMAX.GT.MAX1) THEN MAX2=MAX1 MAX1=TMAX BSTCOL=CURCOL ELSE MAX2=TMAX END IF END IF 40 CONTINUE END IF C ROW BIDS FOR BSTCOL INCREASING ITS PRICE, AND GETS ASSIGNED C TO BSTCOL, WHILE ANY ROW ASSIGNED TO BSTCOL BECOMES UNASSIGNED PCOL(BSTCOL)=PCOL(BSTCOL)+MAX1-MAX2+INCR PROW(ROW)=MAX2-INCR OLDROW=ROW_ASSIGNED_TO(BSTCOL) ROW_ASSIGNED_TO(BSTCOL)=ROW IF (OLDROW.GT.0) THEN NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF 100 CONTINUE C ************* END OF A FORWARD AUCTION CYCLE *************** C OPTIONALLY COLLECT STATISTICS X AVERAGE=(CYCLES*AVERAGE+NOLIST)/(CYCLES+1) X CYCLES=CYCLES+1 C CHECK IF A SWITCH TO MODIFIED REVERSE AUCTION C SHOULD BE MADE (IF NONEWLIST IS NO MORE THAN THE THRESHOLD). C ALSO INCREASE THE MINIMAL BIDDING INCREMENT UP TO A MAXIMUM C VALUE OF EPSILON (THIS IS THE ADAPTIVE FEATURE) INCR=INCR*INCRFACTOR IF (INCR.GT.EPSILON) INCR=EPSILON IF (NONEWLIST.GT.THRESH) THEN NOLIST=NONEWLIST GO TO 15 END IF END IF C *************************************************************** C C START OF MODIFIED REVERSE AUCTION C C *************************************************************** IF (EPSILON.GT.1) GOTO 400 IF (N.EQ.M) GOTO 400 INCR=EPSILON C COMPUTE LAMBDA WHICH IS THE MINIMUM COLUMN PRICE OVER ALL C ASSIGNED COLUMNS, AND COMPILE THE LIST OF UNASSIGNED COLUMNS LAMBDA=ILARGE RNOLIST=0 DO 310 J=1,N ROW=ROW_ASSIGNED_TO(J) IF (ROW.GT.0) THEN COL_ASSIGNED_TO(ROW)=J IF (LAMBDA.GT.PCOL(J)) LAMBDA=PCOL(J) ELSE RNOLIST=RNOLIST+1 REV_LIST(RNOLIST)=J END IF 310 CONTINUE IF (RNOLIST.NE.N-M) THEN PRINT*,'NUMBER OF UNASSIGNED ROWS IS WRONG' PAUSE STOP END IF C START OF A NEW MODIFIED REVERSE AUCTION CYCLE 315 CONTINUE C INITIALIZE COUNT OF NEXT LIST OF UNASSIGNED COLUMNS RNONEWLIST=0 C CYCLE THROUGH THE CURRENT LIST OF UNASSIGNED COLUMNS DO 340 J=1,RNOLIST COLUMN=REV_LIST(J) IF (PCOL(COLUMN).LE.LAMBDA) GO TO 340 CURARC=FIN(COLUMN) NXTARC=NXTIN(CURARC) CURROW=START(CURARC) C FIRST TAKE CARE OF THE EXCEPTIONAL CASE WHERE ROW HAS ONLY ONE ARC IF (NXTARC.EQ.0) THEN MAX1=COST(CURARC)-PROW(CURROW) IF (LAMBDA.GE.MAX1-INCR) THEN PCOL(COLUMN)=LAMBDA GO TO 340 END IF PCOL(COLUMN)=LAMBDA PROW(CURROW)=COST(CURARC)-LAMBDA OLDCOL=COL_ASSIGNED_TO(CURROW) COL_ASSIGNED_TO(CURROW)=COLUMN IF (PCOL(OLDCOL).GT.LAMBDA) THEN RNONEWLIST=RNONEWLIST+1 REV_LIST(RNONEWLIST)=OLDCOL END IF GO TO 340 END IF C NEXT TAKE CARE OF THE REGULAR CASE WHERE ROW HAS MULTIPLE ARCS NXTROW=START(NXTARC) MAX1=COST(CURARC)-PROW(CURROW) MAX2=COST(NXTARC)-PROW(NXTROW) IF (MAX1.GE.MAX2) THEN BSTROW=CURROW ELSE TMAX=MAX1 MAX1=MAX2 MAX2=TMAX BSTROW=NXTROW END IF CURARC=NXTIN(NXTARC) 330 IF (CURARC.GT.0) THEN CURROW=START(CURARC) TMAX=COST(CURARC)-PROW(CURROW) IF (TMAX.GT.MAX2) THEN IF (TMAX.GT.MAX1) THEN MAX2=MAX1 MAX1=TMAX BSTROW=CURROW ELSE MAX2=TMAX END IF END IF CURARC=NXTIN(CURARC) GO TO 330 END IF C COLUMN BIDS FOR BSTROW INCREASING ITS PRICE, AND GETS ASSIGNED C TO BSTROW, WHILE ANY COLUMN ASSIGNED TO BSTROW BECOMES UNASSIGNED IF (LAMBDA.GE.MAX1-INCR) THEN PCOL(COLUMN)=LAMBDA GO TO 340 END IF DELTA=MAX1-MAX2+INCR IF (DELTA.GT.MAX1-LAMBDA) DELTA=MAX1-LAMBDA PROW(BSTROW)=PROW(BSTROW)+DELTA PCOL(COLUMN)=MAX1-DELTA OLDCOL=COL_ASSIGNED_TO(BSTROW) COL_ASSIGNED_TO(BSTROW)=COLUMN IF (PCOL(OLDCOL).GT.LAMBDA) THEN RNONEWLIST=RNONEWLIST+1 REV_LIST(RNONEWLIST)=OLDCOL END IF 340 CONTINUE C ************* END OF A MODIFIED REVERSE AUCTION CYCLE *************** C OPTIONALLY COLLECT STATISTICS X AVERAGE=(CYCLES*AVERAGE+RNOLIST)/(CYCLES+1) X CYCLES=CYCLES+1 C CHECK IF THERE ARE STILL 'MANY' UNASSIGNED ELIGIBLE COLUMNS, THAT IS, C IF THE NUMBER OF ELIGIBLE UNASSIGNED COLUMNS IS GREATER THAN C THE PARAMETER THRESH. IF NOT, REPLACE CURRENT LIST WITH THE NEW LIST, C AND GO FOR ANOTHER CYCLE. OTHERWISE, IF EPSILON > 1, REDUCE EPSILON, C RESET THE ASSIGNMENT TO EMPTY AND RESTART AUCTION; C IF EPSILON = 1 TERMINATE. IF (RNONEWLIST.GT.THRESH) THEN RNOLIST=RNONEWLIST GO TO 315 END IF C *************************************************************** C C END OF SUBPROBLEM (SCALING PHASE) C C *************************************************************** 400 CONTINUE C ****** IF EPSILON IS 1 TERMINATE ****** IF (EPSILON.EQ.1) THEN RETURN ELSE C ELSE REDUCE EPSILON AND UPDATE PARAMETERS FOR THE NEXT SCALING PHASE NUMPHASES=NUMPHASES+1 EPSILON=INT(EPSILON/FACTOR) IF (EPSILON.GT.INCR) EPSILON=INT(EPSILON/FACTOR) IF ((EPSILON.LT.1).OR.(EPSILON.LT.ENDEPS)) EPSILON=1 IF (STARTINCR.LT.EPSILON) STARTINCR=FACTOR*STARTINCR THRESH=INT(THRESH/FACTOR) X PRINT*,'*** END OF A SCALING PHASE; NEW EPSILON=',EPSILON GO TO 12 END IF END SUBROUTINE SETRAN(ISEED) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C PORTABLE CONGRUENTIAL (UNIFORM) RANDOM NUMBER GENERATOR: C NEXT_VALUE = [(7**5) * PREVIOUS_VALUE] MODULO[(2**31)-1] C C THIS GENERATOR CONSISTS OF TWO ROUTINES: C (1) SETRAN - INITIALIZES CONSTANTS AND SEED C (2) RRAN - GENERATES A REAL RANDOM NUMBER C C THE GENERATOR REQUIRES A MACHINE WITH AT LEAST 32 BITS OF PRECISION. C THE SEED (ISEED) MUST BE IN THE RANGE (1,(2**31)-1). C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IF(ISEED.LT.1) STOP 77 MULT=16807 MODUL=2147483647 I15=2**15 I16=2**16 JRAN=ISEED RETURN END C REAL FUNCTION RAN() IMPLICIT REAL*4 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C RAN GENERATES A REAL RANDOM NUMBER BETWEEN 0 AND 1 C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IXHI=JRAN/I16 IXLO=JRAN-IXHI*I16 IXALO=IXLO*MULT LEFTLO=IXALO/I16 IXAHI=IXHI*MULT IFULHI=IXAHI+LEFTLO IRTLO=IXALO-LEFTLO*I16 IOVER=IFULHI/I15 IRTHI=IFULHI-IOVER*I15 JRAN=((IRTLO-MODUL)+IRTHI*I16)+IOVER IF(JRAN.LT.0) JRAN=JRAN+MODUL RAN = FLOAT(JRAN)/FLOAT(MODUL) RETURN END AUCTION-FR This code implements the combined forward/reverse auction algorithm with $\e$-scaling for the symmetric assignment problem; cf.\ Section 4.2. For good performance, it is frequently not essential to use $\e$-scaling in this code. It is therefore recommended that the code be tried without $\e$-scaling by setting the starting $\e$ to 1. C ****************************************************************** C C SAMPLE CALLING PROGRAM FOR COMBINED FORWARD/REVERSE C AUCTION ALGORITHM C C THIS DRIVER CREATES A SYMMETRIC ASSIGNMENT PROBLEM C WITH EQUAL NUMBER OF ROWS AND COLUMNS, C AND CALLS THE AUCTION_FR SUBROUTINE TO FIND AN C ASSIGNMENT OF MAXIMAL VALUE. C C ****************************************************************** PARAMETER(MAXNODES=10000, MAXARCS=100000) IMPLICIT NONE INTEGER N,NA,A,IA,ILARGE,BEGEPS,ENDEPS,CYCLES INTEGER NUMPHASES,STARTINCR INTEGER I,J,ARC,NOASS,ICOST,ABSCOST,CURARC INTEGER CURCOL,FSTARC,LSTARC,MINCOST,MAXCOST,MCOST INTEGER COL_ASSIGNED_TO(MAXNODES),ROW_ASSIGNED_TO(MAXNODES) INTEGER FOUT(MAXNODES),FIN(MAXNODES),NXTIN(MAXARCS) INTEGER COST(MAXARCS),START(MAXARCS),END(MAXARCS) INTEGER PCOL(MAXNODES),PROW(MAXNODES) INTEGER PRDARC(MAXNODES) REAL*8 FACTOR,TT1,TT2,TCOST,AVERAGE COMMON/ARRAYC/COST/ARRAYS/START/ARRAYE/END COMMON/ARRAYFO/FOUT/ARRAYFI/FIN/ARRAYNI/NXTIN COMMON/ARRAYRA/ROW_ASSIGNED_TO/ARRAYCA/COL_ASSIGNED_TO COMMON/ARRAYPC/PCOL/ARRAYPR/PROW COMMON/BK1/N,A,ILARGE/BK2/CYCLES,AVERAGE,NUMPHASES C **************************************************************** C C PROBLEM GENERATION CODE STARTS HERE C THE USER MAY REPLACE THIS CODE WITH A CODE THAT READS C HIS/HER PROBLEM FROM A FILE C C **************************************************************** C THIS CODE INCLUDES A UNIFORM RANDOM NUMBER GENERATOR C WHICH RETURNS A VALUE IN (0,1) C INITIALIZE RANDOM GENERATOR INTEGER MULT,MODUL,I15,I16,JRAN,ISEED REAL RAN COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN ISEED=13502460 CALL SETRAN(ISEED) PRINT*,'GENERATING A SYMMETRIC ASSIGNMENT PROBLEM' PRINT*,'********************************************' C **** READ THE NUMBER OF ROWS N & THE NUMBER OF ARCS A **** PRINT*,'ENTER THE NUMBER OF ROWS (AND COLUMNS)' READ*,N 5 PRINT*,'ENTER THE NUMBER OF ARCS PER ROW (>1)' READ*,NA IF (NA.LT.2) GOTO 5 PRINT*,'ENTER THE MINIMUM AND THE MAXIMUM COST' READ*,MINCOST,MAXCOST C THE NUMBER OF ARCS IS N*NA A=N*NA C THE ARCS INCIDENT TO ROW I ARE FOUT(I) TO FOUT(I+1)-1 C ALSO, FOR FEASIBILITY EACH ROW IS DIRECTLY CONNECTED C WITH THE CORRESPONDING COLUMN DO 20 I=1,N FOUT(I)=1+(I-1)*NA 20 CONTINUE FOUT(N+1)=A+1 DO 22 I=1,N DO 23 ARC=FOUT(I),FOUT(I+1)-1 START(ARC)=I 23 CONTINUE 22 CONTINUE C GENERATE THE END(ARC) AND COST(ARC) WHICH ARE THE COLUMN C AND THE COST COEFFICIENT ASSOCIATED WITH ARC DO 25 ARC=1,A END(ARC)=1+RAN()*N IF ((END(ARC).GT.N).OR.(END(ARC).LT.1)) THEN PRINT*,'ERROR IN PROBLEM GENERATION' PAUSE STOP END IF COST(ARC)=MINCOST+RAN()*(MAXCOST-MINCOST) 25 CONTINUE C MODIFY THE END OF THE LAST ARC OUT OF EACH ROW FOR FEASIBILITY C AND SET ITS COST TO -MAXCOST DO 30 I=1,N END(FOUT(I+1)-1)=I COST(FOUT(I+1)-1)=-MAXCOST 30 CONTINUE C CONSTRUCT THE FIN() AND NXTIN() ARRAYS DO 32 J=1,N FIN(J)=0 PRDARC(J)=0 32 CONTINUE DO 33 ARC=1,A NXTIN(ARC)=0 J=END(ARC) IF (FIN(J).NE.0) THEN NXTIN(PRDARC(J))=ARC ELSE FIN(J)=ARC END IF PRDARC(J)=ARC 33 CONTINUE C **************************************************************** C C PROBLEM GENERATION CODE ENDS HERE C C **************************************************************** C SCALE THE COST TO WORK WITH INTEGER EPSILON MAXCOST=0 DO 35 IA=1,A ABSCOST=IABS(COST(IA)) IF (ABSCOST.GT.MAXCOST) MAXCOST=ABSCOST COST(IA)=COST(IA)*(N+1) 35 CONTINUE C *** ILARGE IS A VERY LARGE INTEGER FOR YOUR MACHINE *** C THE SCALED COSTS SHOULD BE SIGNIFICANTLY SMALLER THAN C ILARGE AND SIGNIFICANTLY LARGER THAN -ILARGE ILARGE=2000000000 IF (MAXCOST.GT.INT(ILARGE/(N+1))) THEN PRINT*,'THE COST RANGE IS TOO LARGE FOR INTEGER ARITHMETIC' PAUSE STOP END IF MAXCOST=MAXCOST*(N+1) C THE FOLLOWING PARAMETERS BEGEPS, FACTOR, ENDEPS, AND STARTINCR C ARE PASSED TO THE AUCTION ALGORITHM. VALUES BETWEEN C (A) MAXCOST/2 AND 1 FOR BEGEPS C (B) 4 AND 6 FOR FACTOR C (C) N AND 1 FOR ENDEPS C (D) 1 AND BEGEPS FOR STARTINCR C HAVE WORKED WELL FOR LARGE SPARSE PROBLEMS. C FOR DENSE PROBLEMS IT IS RECOMMENDED THAT C BEGEPS BE SET TO A SMALLER VALUE, C ENDEPS BE SET TO 1, C STARTINCR BE SET TO 1. C GENERALLY, THE COMBINED FORWARD/REVERSE ALGORITHM C WORKS BEST WITH A SMALLER VALUE OF BEGEPS THAN THE C FORWARD ALGORITHM. C IT IS WORTH TRYING BEGEPS=1, IN WHICH CASE THERE IS C ONLY ONE SCALING PHASE (THAT IS, NO E-SCALING IS USED). PRINT*,'********************************************' PRINT*,'MAXIMUM COST IS ',MAXCOST PRINT*,'ENTER THE STARTING EPSILON' READ*,BEGEPS IF (BEGEPS.LT.1) BEGEPS=1 PRINT*,'ENTER THE EPSILON REDUCTION FACTOR' READ*,FACTOR ENDEPS=N/10 IF (ENDEPS.LT.1) ENDEPS=1 IF (ENDEPS.GT.BEGEPS) ENDEPS=BEGEPS STARTINCR=1 IF (STARTINCR.LT.1) STARTINCR=1 PRINT*,'********************************************' PRINT*,'STARTING EPSILON = ',BEGEPS PRINT*,'EPSILON REDUCTION FACTOR = ',FACTOR PRINT*,'THRESHOLD EPSILON BEFORE IT IS SET TO 1 = ',ENDEPS PRINT*,'STARTING MIN BIDDING INCREMENT = ',STARTINCR PRINT*,'CALLING FORWARD/REVERSE AUCTION' PRINT*,'********************************************' C GET STARTING TIME FOR THE MAC II TT1 = LONG(362)/60.0 CALL AUCTION_FR(BEGEPS,FACTOR,ENDEPS,STARTINCR) C GET ENDING TIME FOR THE MAC II TT2 = LONG(362)/60.0 - TT1 PRINT *,'FINISHED --- TOTAL CPU TIME', TT2,' SECS' PRINT*,'********************************************' C *** DISPLAY RESULTS *** X WRITE(9,2010) CYCLES X2010 FORMAT(' NO OF AUCTION CYCLES',I7) X WRITE(9,2020) AVERAGE X2020 FORMAT(' AVERAGE NUMBER OF BIDS PER CYCLE',F9.3) WRITE(9,2030) NUMPHASES 2030 FORMAT('NO OF EPSILON SUBPROBLEMS SOLVED =',I7) C CHECK OPTIMALITY & CALCULATE COST TCOST=0 DO 40 I=1,N J=COL_ASSIGNED_TO(I) IF (J.EQ.0) THEN PRINT*,'ROW ',I,' IS UNASSIGNED' END IF IF (ROW_ASSIGNED_TO(J).NE.I) THEN PRINT*,'ASSIGNMENT MIXUP: ROW ',I,' COLUMN ',J END IF FSTARC=FOUT(I) LSTARC=FOUT(I+1)-1 MCOST=-ILARGE DO 45 ARC=FSTARC,LSTARC CURCOL=END(ARC) IF (PROW(I)+PCOL(CURCOL).LT.COST(ARC)-1) THEN PRINT*,'1-CS VIOLATED AT ARC ',ARC END IF IF (CURCOL.EQ.J) THEN IF (MCOST.LT.COST(ARC)) THEN MCOST=COST(ARC) END IF END IF 45 CONTINUE TCOST=TCOST+MCOST/(N+1) IF (PROW(I)+PCOL(J).NE.MCOST) THEN PRINT*,'1-CS VIOLATED AT ROW ',I END IF 40 CONTINUE DO 50 J=1,N IF (ROW_ASSIGNED_TO(J).EQ.0) THEN PRINT*,'COLUMN ',J,' IS UNASSIGNED' END IF 50 CONTINUE WRITE(9,2100) TCOST 2100 FORMAT(' ASSIGNMENT COST=',F18.2) PRINT *, ' PROGRAM ENDED; TO EXIT ' PAUSE END C ****************************************************************** C C FORWARD/REVERSE AUCTION CODE FOR N BY N ASSIGNMENT PROBLEMS C C WRITTEN BY DIMITRI P. BERTSEKAS C C DEC. 1990 C C THIS CODE IMPLEMENTS THE FORWARD/REVERSE AUCTION ALGORITHM C WITH E-SCALING FOR SYMMETRIC N BY N ASSIGNMENT PROBLEMS. C IT SOLVES A SEQUENCE OF SUBPROBLEMS AND DECREASES C EPSILON BY A CONSTANT FACTOR BETWEEN SUBPROBLEMS. C THIS VERSION CORRESPONDS TO A GAUSS-SEIDEL MODE C AND SOLVES EPSILON SUBPROBLEMS INEXACTLY. C C THE CODE IS AN IMPROVED VERSION OF AN EARLIER (SEPT. 1985) C AUCTION CODE WITH E-SCALING WRITTEN BY DIMITRI P. BERTSEKAS C C THE CODE TREATS THE PROBLEM AS A MAXIMIZATION PROBLEM. C TO SOLVE A MINIMIZATION PROBLEM, REVERSE THE SIGN OF THE C ARC COSTS PRIOR TO CALLING AUCTION, AND REVERSE AGAIN C THE SIGN OF THE OPTIMAL COST UPON RETURN FROM AUCTION. C THIS CODE ALLOWS MULTIPLE ARCS BETWEEN A ROW AND A COLUMN. C C THIS VERSION OF THE AUCTION ALGORITHM IS ADAPTIVE. HERE THE MINIMAL C BIDDING INCREMENT (STORED IN THE VARIABLE INCR) MAY BE SMALLER THAN C EPSILON. FOR EVERY SUBPROBLEM, INCR STARTS AT THE PARAMETER VALUE C STARTINCR (WHICH IS PASSED TO THE AUCTION ROUTINE) AND IS INCREASED C BY A FACTOR OF 2 AT THE END OF EACH CYCLE, UP TO A MAXIMUM VALUE OF C EPSILON. THIS ADAPTIVE FEATURE IS PARTICULARLY EFFECTIVE FOR C DENSE PROBLEMS. IT CAN BE DEFEATED BY SELECTING STARTINCR=BEGEPS C C THIS CODE MAY FAIL DUE TO INTEGER OVERFLOW IF THE NUMBER OF NODES C OR THE COST RANGE (OR BOTH) ARE LARGE. TO CORRECT THIS SITUATION, C THE PRICES AND OTHER RELATED VARIABLES (MAX1,MAX2,TMAX ETC) SHOULD C BE DECLARED AS DOUBLE PRECISION REALS. C ******************************************************************** C C THE USER MUST SUPPLY THE FOLLOWING PROBLEM DATA IN FORWARD STAR FORMAT C (THAT IS, ALL ARCS OF THE SAME ROW ARE NUMBERED CONSECUTIVELY): C N=NUMBER OF ROWS (EQUALS NUMBER OF COLUMNS) C A=NUMBER OF ARCS C FOUT(ROW)=FIRST ARC COMING OUT OF ROW C FIN(COL)=FIRST ARC COMING INTO COL C NXTIN(ARC)=NEXT ARC INCIDENT TO THE SAME COLUMN AS ARC C COST(ARC)=COST OF ARC C START(ARC)=ROW CORRESPONDING TO ARC C END(ARC)=COLUMN CORRESPONDING TO ARC C C AND THE FOLLOWING PARAMETERS FOR THE AUCTION ALGORITHM: C BEGEPS=STARTING VALUE OF EPSILON (MUST BE NO LESS THAN 1) C ENDEPS=FINAL VALUE OF EPSILON BEFORE IT IS SET TO 1 C FACTOR=FACTOR BY WHICH EPSILON IS DECREASED BETWEEN SUBPROBLEMS C STARTINCR=THE STARTING VALUE OF THE BIDDING INCREMENT C ENDEPS SHOULD NOT EXCEED BEGEPS. C FACTOR MUST BE GREATER THAN 1 (UNLESS BEGEPS=1). C C FOUT(.) IS AN ARRAY OF LENGTH N. C COST(.),END(.) ARE ARRAYS OF LENGTH A. C C THE SOLUTION IS CONTAINED IN THE ARRAY ROW_ASSIGNED_TO(.) WHERE C ROW_ASSIGNED_TO(COL) GIVES THE ROW ASSIGNED TO COL. C ALSO COL_ASSIGNED_TO(ROW) GIVES THE COLUMN ASSIGNED TO ROW. C C THIS ALGORITHM DOES NOT CHECK FOR INFEASIBILITY OF THE PROBLEM. C TO MAKE SURE THE PROBLEM IS FEASIBLE THE USER MAY ADD C ADDITIONAL VERY SMALL COST ARCS. C C ****************************************************************** C ALL PROBLEM DATA ARE INTEGER C ****************************************************************** SUBROUTINE AUCTION_FR(BEGEPS,FACTOR,ENDEPS,STARTINCR) PARAMETER(MAXNODES=10000, MAXARCS=100000) IMPLICIT NONE INTEGER A,K,N,I,J,M,CURARC,CURCOL,CURROW INTEGER THRESH,INCR,STARTINCR,INCRFACTOR INTEGER NOLIST,RNOLIST,NONEWLIST,RNONEWLIST INTEGER ROW,COLUMN, BSTROW,BSTCOL INTEGER FSTARC,FSTCOL,NXTARC,NXTROW,LSTARC,SNDARC,SNDCOL INTEGER MAX1,MAX2,TMAX,TMIN,TRDARC,EPSILON,BEGEPS,ENDEPS INTEGER ISMALL,ILARGE,LARGEINCR,CYCLES INTEGER NUMPHASES,OLDROW,OLDCOL INTEGER COL_ASSIGNED_TO(MAXNODES),ROW_ASSIGNED_TO(MAXNODES) INTEGER FOUT(MAXNODES),FIN(MAXNODES),NXTIN(MAXARCS) INTEGER COST(MAXARCS),START(MAXARCS),END(MAXARCS) INTEGER LIST(MAXNODES),REV_LIST(MAXNODES) INTEGER PCOL(MAXNODES),PROW(MAXNODES) REAL*8 AVERAGE,FACTOR LOGICAL READY_SWITCH,SWITCH COMMON/ARRAYC/COST/ARRAYS/START/ARRAYE/END COMMON/ARRAYFO/FOUT/ARRAYFI/FIN/ARRAYNI/NXTIN COMMON/ARRAYRA/ROW_ASSIGNED_TO/ARRAYCA/COL_ASSIGNED_TO COMMON/ARRAYPC/PCOL/ARRAYPR/PROW COMMON/BK1/N,A,ILARGE/BK2/CYCLES,AVERAGE,NUMPHASES C *************************************************************** C ******* CHECK VALIDITY OF PARAMETERS PASSED ******* IF (BEGEPS.LT.1) THEN PRINT*,'STARTING VALUE OF EPSILON IS LESS THAN 1' PRINT*,'EXECUTION ABORTED' STOP END IF IF (ENDEPS.GT.BEGEPS) THEN PRINT*,'PARAMETER ENDEPS IS GREATER THAN PARAMETER BEGEPS' PRINT*,'ENDEPS IS SET AT THE DEFAULT VALUE OF 1' ENDEPS=1 END IF IF ((FACTOR.LE.1).AND.(BEGEPS.GT.1)) THEN PRINT*,'EPSILON REDUCTION FACTOR IS NOT GREATER THAN 1' PRINT*,'EXECUTION ABORTED' STOP END IF IF (STARTINCR.LT.1) THEN PRINT*,'MIN BIDDING INCREMENT IS LESS THAN 1' PRINT*,'STARTINCR IS SET AT THE DEFAULT VALUE OF 1' STARTINCR=1 END IF C ******* INITIALIZATION ******* EPSILON=BEGEPS ISMALL=-ILARGE LARGEINCR=INT(ILARGE/10) THRESH=INT(0.2*N) INCRFACTOR=2 IF (THRESH.GT.100) THRESH=100 X CYCLES=1 X AVERAGE=N NUMPHASES=1 C INITIALIZE FORWARD/REVERSE SWITCH PARAMETERS READY_SWITCH=.FALSE. SWITCH=.TRUE. DO 10 J=1,N PCOL(J)=0 10 CONTINUE FOUT(N+1)=A+1 C **************************************************************** C C THIS IMPLEMENTATION OF THE AUCTION ALGORITHM OPERATES IN CYCLES. C EACH FORWARD AUCTION CYCLE CONSISTS OF ONE BID BY EACH OF THE ROWS C THAT ARE UNASSIGNED AT THE START OF THE CYCLE (THESE ROWS ARE C STORED IN THE ARRAY LIST(.)). AS THE CYCLE PROGRESSES NEW C ROWS BECOME UNASSIGNED; THESE ARE STORED IN LIST(.) C AND WILL SUBMIT A BID AT THE NEXT CYCLE. C REVERSE AUCTION CYCLES ARE STRUCTURED SIMILARLY. C C **************************************************************** C C START SUBPROBLEM (SCALING PHASE) W/ NEW EPSILON C C **************************************************************** C 12 CONTINUE C INITIALIZE ROW AND COLUMN ASSIGNMENT LISTS NOLIST=N DO 20 I=1,N COL_ASSIGNED_TO(I)=0 LIST(I)=I 20 CONTINUE RNOLIST=N DO 22 J=1,N ROW_ASSIGNED_TO(J)=0 REV_LIST(J)=J 22 CONTINUE INCR=STARTINCR IF (INCR.GT.EPSILON) INCR=EPSILON IF (EPSILON.EQ.1) THRESH=0 C **************************************************************** C C START FORWARD AUCTION CYCLE C C **************************************************************** 15 CONTINUE C INITIALIZE COUNT OF NEXT LIST OF UNASSIGNED ROWS NONEWLIST=0 C CYCLE THROUGH THE CURRENT LIST OF UNASSIGNED ROWS DO 100 I=1,NOLIST ROW=LIST(I) IF (COL_ASSIGNED_TO(ROW).GT.0) GO TO 100 FSTARC=FOUT(ROW) LSTARC=FOUT(ROW+1)-1 FSTCOL=END(FSTARC) C FIRST TAKE CARE OF THE EXCEPTIONAL CASE WHERE ROW HAS ONLY ONE ARC IF (FSTARC.EQ.LSTARC) THEN PCOL(FSTCOL)=PCOL(FSTCOL)+LARGEINCR PROW(ROW)=COST(FSTARC)-PCOL(FSTCOL) OLDROW=ROW_ASSIGNED_TO(FSTCOL) ROW_ASSIGNED_TO(FSTCOL)=ROW COL_ASSIGNED_TO(ROW)=FSTCOL IF (OLDROW.GT.0) THEN COL_ASSIGNED_TO(OLDROW)=0 NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF GO TO 100 END IF C NEXT TAKE CARE OF THE REGULAR CASE WHERE ROW HAS MULTIPLE ARCS SNDARC=FSTARC+1 SNDCOL=END(SNDARC) MAX1=COST(FSTARC)-PCOL(FSTCOL) MAX2=COST(SNDARC)-PCOL(SNDCOL) IF (MAX1.GE.MAX2) THEN BSTCOL=FSTCOL ELSE TMAX=MAX1 MAX1=MAX2 MAX2=TMAX BSTCOL=SNDCOL END IF IF (SNDARC.LT.LSTARC) THEN TRDARC=SNDARC+1 DO 40 CURARC=TRDARC,LSTARC CURCOL=END(CURARC) TMAX=COST(CURARC)-PCOL(CURCOL) IF (TMAX.GT.MAX2) THEN IF (TMAX.GT.MAX1) THEN MAX2=MAX1 MAX1=TMAX BSTCOL=CURCOL ELSE MAX2=TMAX END IF END IF 40 CONTINUE END IF C ROW BIDS FOR BSTCOL INCREASING ITS PRICE, AND GETS ASSIGNED C TO BSTCOL, WHILE ANY ROW ASSIGNED TO BSTCOL BECOMES UNASSIGNED PCOL(BSTCOL)=PCOL(BSTCOL)+MAX1-MAX2+INCR PROW(ROW)=MAX2-INCR COL_ASSIGNED_TO(ROW)=BSTCOL OLDROW=ROW_ASSIGNED_TO(BSTCOL) ROW_ASSIGNED_TO(BSTCOL)=ROW IF (OLDROW.GT.0) THEN COL_ASSIGNED_TO(OLDROW)=0 NONEWLIST=NONEWLIST+1 LIST(NONEWLIST)=OLDROW END IF 100 CONTINUE C ************* END OF A FORWARD AUCTION CYCLE *************** C OPTIONALLY COLLECT STATISTICS X AVERAGE=(CYCLES*AVERAGE+NOLIST)/(CYCLES+1) X CYCLES=CYCLES+1 C CHECK IF THERE ARE STILL 'MANY' UNASSIGNED ROWS, THAT IS, IF THE C NUMBER OF UNASSIGNED ROWS IS GREATER THAN C THE PARAMETER THRESH. IF NOT, REPLACE CURRENT LIST WITH THE NEW LIST, C AND GO FOR ANOTHER CYCLE. OTHERWISE, IF EPSILON > 1, REDUCE EPSILON, C RESET THE ASSIGNMENT TO EMPTY AND RESTART AUCTION; C IF EPSILON = 1 TERMINATE. C ALSO INCREASE THE MINIMAL BIDDING INCREMENT UP TO A MAXIMUM C VALUE OF EPSILON (THIS IS THE ADAPTIVE FEATURE) INCR=INCR*INCRFACTOR IF (INCR.GT.EPSILON) INCR=EPSILON IF (NONEWLIST.GT.THRESH) THEN IF (SWITCH) THEN NOLIST=NONEWLIST GO TO 115 END IF IF (NONEWLIST.LT.NOLIST) READY_SWITCH=.TRUE. IF ((NONEWLIST.EQ.NOLIST).AND.(READY_SWITCH)) THEN READY_SWITCH=.FALSE. GO TO 115 ELSE NOLIST=NONEWLIST GO TO 15 END IF ELSE GO TO 300 END IF C **************************************************************** C C START REVERSE AUCTION CYCLE C C **************************************************************** 115 CONTINUE C INITIALIZE COUNT OF NEXT LIST OF UNASSIGNED COLUMNS RNONEWLIST=0 C CYCLE THROUGH THE CURRENT LIST OF UNASSIGNED COLUMNS DO 200 J=1,RNOLIST COLUMN=REV_LIST(J) IF (ROW_ASSIGNED_TO(COLUMN).GT.0) GO TO 200 CURARC=FIN(COLUMN) NXTARC=NXTIN(CURARC) CURROW=START(CURARC) C FIRST TAKE CARE OF THE EXCEPTIONAL CASE WHERE ROW HAS ONLY ONE ARC IF (NXTARC.EQ.0) THEN PROW(CURROW)=PROW(CURROW)+LARGEINCR PCOL(COLUMN)=COST(CURARC)-PROW(CURROW) OLDCOL=COL_ASSIGNED_TO(CURROW) COL_ASSIGNED_TO(CURROW)=COLUMN ROW_ASSIGNED_TO(COLUMN)=CURROW IF (OLDCOL.GT.0) THEN ROW_ASSIGNED_TO(OLDCOL)=0 RNONEWLIST=RNONEWLIST+1 REV_LIST(RNONEWLIST)=OLDCOL END IF GO TO 200 END IF C NEXT TAKE CARE OF THE REGULAR CASE WHERE ROW HAS MULTIPLE ARCS NXTROW=START(NXTARC) MAX1=COST(CURARC)-PROW(CURROW) MAX2=COST(NXTARC)-PROW(NXTROW) IF (MAX1.GE.MAX2) THEN BSTROW=CURROW ELSE TMAX=MAX1 MAX1=MAX2 MAX2=TMAX BSTROW=NXTROW END IF CURARC=NXTIN(NXTARC) 130 IF (CURARC.GT.0) THEN CURROW=START(CURARC) TMAX=COST(CURARC)-PROW(CURROW) IF (TMAX.GT.MAX2) THEN IF (TMAX.GT.MAX1) THEN MAX2=MAX1 MAX1=TMAX BSTROW=CURROW ELSE MAX2=TMAX END IF END IF CURARC=NXTIN(CURARC) GO TO 130 END IF C COLUMN BIDS FOR BSTROW INCREASING ITS PRICE, AND GETS ASSIGNED C TO BSTROW, WHILE ANY COLUMN ASSIGNED TO BSTROW BECOMES UNASSIGNED PROW(BSTROW)=PROW(BSTROW)+MAX1-MAX2+INCR PCOL(COLUMN)=MAX2-INCR ROW_ASSIGNED_TO(COLUMN)=BSTROW OLDCOL=COL_ASSIGNED_TO(BSTROW) COL_ASSIGNED_TO(BSTROW)=COLUMN IF (OLDCOL.GT.0) THEN ROW_ASSIGNED_TO(OLDCOL)=0 RNONEWLIST=RNONEWLIST+1 REV_LIST(RNONEWLIST)=OLDCOL END IF 200 CONTINUE C ************* END OF A REVERSE AUCTION CYCLE *************** C OPTIONALLY COLLECT STATISTICS X AVERAGE=(CYCLES*AVERAGE+RNOLIST)/(CYCLES+1) X CYCLES=CYCLES+1 C CHECK IF THERE ARE STILL 'MANY' UNASSIGNED COLUMNS, THAT IS, IF THE C NUMBER OF UNASSIGNED COLUMNS IS GREATER THAN C THE PARAMETER THRESH. IF NOT, REPLACE CURRENT LIST WITH THE NEW LIST, C AND GO FOR ANOTHER CYCLE. OTHERWISE, IF EPSILON > 1, REDUCE EPSILON, C RESET THE ASSIGNMENT TO EMPTY AND RESTART AUCTION; C IF EPSILON = 1 TERMINATE. C ALSO INCREASE THE MINIMAL BIDDING INCREMENT UP TO A MAXIMUM C VALUE OF EPSILON (THIS IS THE ADAPTIVE FEATURE) INCR=INCR*INCRFACTOR IF (INCR.GT.EPSILON) INCR=EPSILON IF (RNONEWLIST.GT.THRESH) THEN IF (SWITCH) THEN SWITCH=.FALSE. RNOLIST=RNONEWLIST GO TO 15 END IF IF (RNONEWLIST.LT.RNOLIST) READY_SWITCH=.TRUE. IF ((RNONEWLIST.EQ.RNOLIST).AND.(READY_SWITCH)) THEN READY_SWITCH=.FALSE. GO TO 15 ELSE RNOLIST=RNONEWLIST GO TO 115 END IF ELSE GO TO 300 END IF C *************************************************************** C C END OF SUBPROBLEM (SCALING PHASE) C C *************************************************************** 300 CONTINUE C ****** IF EPSILON IS 1 TERMINATE ****** IF (EPSILON.EQ.1) THEN RETURN ELSE C ELSE REDUCE EPSILON AND UPDATE PARAMETERS FOR THE NEXT SCALING PHASE NUMPHASES=NUMPHASES+1 EPSILON=INT(EPSILON/FACTOR) IF (EPSILON.GT.INCR) EPSILON=INT(EPSILON/FACTOR) IF ((EPSILON.LT.1).OR.(EPSILON.LT.ENDEPS)) EPSILON=1 IF (STARTINCR.LT.EPSILON) STARTINCR=FACTOR*STARTINCR THRESH=INT(THRESH/FACTOR) X PRINT*,'*** END OF A SCALING PHASE; NEW EPSILON=',EPSILON GO TO 12 END IF END SUBROUTINE SETRAN(ISEED) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C PORTABLE CONGRUENTIAL (UNIFORM) RANDOM NUMBER GENERATOR: C NEXT_VALUE = [(7**5) * PREVIOUS_VALUE] MODULO[(2**31)-1] C C THIS GENERATOR CONSISTS OF TWO ROUTINES: C (1) SETRAN - INITIALIZES CONSTANTS AND SEED C (2) RRAN - GENERATES A REAL RANDOM NUMBER C C THE GENERATOR REQUIRES A MACHINE WITH AT LEAST 32 BITS OF PRECISION. C THE SEED (ISEED) MUST BE IN THE RANGE (1,(2**31)-1). C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IF(ISEED.LT.1) STOP 77 MULT=16807 MODUL=2147483647 I15=2**15 I16=2**16 JRAN=ISEED RETURN END C REAL FUNCTION RAN() IMPLICIT REAL*4 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C RAN GENERATES A REAL RANDOM NUMBER BETWEEN 0 AND 1 C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IXHI=JRAN/I16 IXLO=JRAN-IXHI*I16 IXALO=IXLO*MULT LEFTLO=IXALO/I16 IXAHI=IXHI*MULT IFULHI=IXAHI+LEFTLO IRTLO=IXALO-LEFTLO*I16 IOVER=IFULHI/I15 IRTHI=IFULHI-IOVER*I15 JRAN=((IRTLO-MODUL)+IRTHI*I16)+IOVER IF(JRAN.LT.0) JRAN=JRAN+MODUL RAN = FLOAT(JRAN)/FLOAT(MODUL) RETURN END *************************************************************************** NAUCTION_SP: Combined Naive Auction and Sequential Shortest Path Code *************************************************************************** This code implements the sequential shortest path method for the assignment problem, preceded by an extensive initialization using the naive auction algorithm (cf.\ Section 1.2.4). The code is quite similar in structure and performance to a code of the author [Ber81] and to the code of [JoV86] and [JoV87]. These codes also combined a naive auction initialization with the sequential shortest path method. C ******************************************************************** C C COMBINED NAIVE AUCTION AND SEQUENTIAL SHORTEST PATH METHOD C FOR SYMMETRIC N X N ASSIGNMENT PROBLEMS C FINDS AN OPTIMAL ASSIGNMENT C C WRITTEN BY DIMITRI P. BERTSEKAS C C ******************************************************************** C C USER MUST SUPPLY THE FOLLOWING PROBLEM DATA C N=NUMBER OF ROWS AND NUMBER OF COLUMNS C A=NUMBER OF ARCS C FOUT(ROW)=1ST ARC COMING OUT OF ROW C COST(ARC)=COST OF ARC C END(ARC)=COLUMN CORRESPONDING TO ARC C C FOUT(.) IS AN N - LENGTH ARRAY C COST(.),END(.),NXTEND(.) ARE ARRAYS OF LENGTH A C THE OPTIMAL ASSIGNMENT IS CONTAINED IN THE ARRAY ASSIGN(.) WHERE C ASSIGN(COL) IS THE ROW ASSIGNED TO COL C THIS ALGORITHM DOES NOT CHECK FOR INFEASIBILITY OF THE PROBLEM C TO MAKE SURE THE PROBLEM IS FEASIBLE THE USER MAY ADD C ADDITIONAL VERY HIGH COST ARCS C C THIS CODE ALLOWS MULTIPLE ARCS BETWEEN A ROW AND A COLUMN. C C ****************************************************************** C ALL PROBLEM DATA ARE INTEGER C ****************************************************************** C IMPLICIT INTEGER (A-Z) INTEGER AUCTNUM,A,STARC,ENDARC,CURROW,ARC,COUNT INTEGER HCOUNT,ROW,FSTARC,FSTCOL,SNDARC,SNDCOL,TMAX,BSTCOL INTEGER TMARG,TA,ROWPR,OLMARG,PRICE,TPRICE,CURCOL,DUMMY INTEGER FOUT(6000),PCOL(6000),PROW(6000),MARG(6000) INTEGER ASSIGN(6000),COLLAB(6000),ROWLAB(6000) INTEGER LIST(6000),SCAN(6000) INTEGER COST(70000),END(70000) LOGICAL MAXSET REAL THRESH,RAND,TT,TIMER,TCOST C C **************************************************************** C C PROBLEM GENERATION CODE STARTS HERE C THE USER MAY REPLACE THIS CODE WITH A CODE THAT READS C HIS/HER PROBLEM FROM A FILE C C **************************************************************** C C RANDOM GENERATOR STUFF C INTEGER MULT,MODUL,I15,I16,JRAN,ISEED REAL RAN COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN C UNIFORM RANDOM NUMBER GENERATOR WHICH RETURNS A VALUE IN (0,1) C C INITIALIZE RANDOM GENERATOR C ISEED=13502460 CALL SETRAN(ISEED) PRINT*,'GENERATING A SYMMETRIC ASSIGNMENT PROBLEM' PRINT*,'********************************************' C **** READ THE NUMBER OF ROWS N & THE NUMBER OF ARCS A **** PRINT*,'ENTER THE NUMBER OF ROWS (AND COLUMNS)' READ*,N 5 PRINT*,'ENTER THE NUMBER OF ARCS PER ROW (>1)' READ*,NA IF (NA.LT.2) GOTO 5 PRINT*,'ENTER THE MINIMUM AND THE MAXIMUM COST' READ*,MINCOST,MAXCOST C THE NUMBER OF ARCS IS N*NA A=N*NA C THE ARCS INCIDENT TO ROW I ARE FOUT(I) TO FOUT(I+1)-1 C FOR FEASIBILITY EACH ROW IS DIRECTLY CONNECTED C WITH THE CORRESPONDING COLUMN DO 20 I=1,N FOUT(I)=1+(I-1)*NA 20 CONTINUE C **** GENERATE THE END(IA) AND COST(IA) WHICH ARE THE COLUMN C AND THE COST COEFFICIENT ASSOCIATED WITH ARC IA **** DO 25 IA=1,A END(IA)=1+RAN()*N IF ((END(IA).GT.N).OR.(END(IA).LT.1)) THEN PRINT*,'ERROR IN PROBLEM GENERATION' PAUSE STOP END IF COST(IA)=MINCOST+RAN()*(MAXCOST-MINCOST) 25 CONTINUE C MODIFY THE END OF THE LAST ARC OUT OF EACH ROW FOR FEASIBILITY C AND SET ITS COST TO -MAXCOST DO 30 I=1,N END(FOUT(I+1)-1)=I COST(FOUT(I+1)-1)=-MAXCOST 30 CONTINUE C C **************************************************************** C C PROBLEM GENERATION CODE ENDS HERE C C **************************************************************** C C MAIN ALGORITHM BEGINS HERE C C THE ALGORITHM INTERNALLY TREATS THE PROBLEM AS A C MINIMIZATION PROBLEM. C TO SOLVE AN ASSIGNMENT MAXIMIZATION PROBLEM, SET C THE FOLLOWING PARAMETER MAXSET TO TRUE, ELSE SET IT TO FALSE MAXSET=.TRUE. IF (MAXSET) THEN DO 35 IA=1,A COST(IA)=-COST(IA) 35 CONTINUE END IF PRINT*,'NO OF ROWS (AND COLUMNS):',N PRINT*,'NO OF ARCS:',A PRINT*,'COMBINED NAIVE AUCTION AND SEQ. SH. PATH METHOD' PRINT *,'**********************************************' TIMER = LONG(362) ISIMPL=0 IPRC=0 ISMALL=-20000000 ILARGE=-ISMALL COUNT=0 HCOUNT=0 FOUT(N+1)=A+1 C THE FOLLOWING INITIALIZATION UP TO THE START OF THE C NAIVE AUCTION CYCLES WAS SUGGESTED BY JONKER AND VOLGENANT C INITIALIZE COLUMN PRICES DO 105 J=1,N PCOL(J)=ILARGE 105 CONTINUE C FIND BEST ROW FOR EACH COLUMN DO 120 I=1,N PROW(I)=0 ROWLAB(I)=0 FSTARC=FOUT(I) LSTARC=FOUT(I+1)-1 DO 110 ARC=FSTARC,LSTARC J=END(ARC) IF (COST(ARC).LT.PCOL(J)) THEN PCOL(J)=COST(ARC) ASSIGN(J)=I END IF 110 CONTINUE 120 CONTINUE C CONSTRUCT ROW ASSIGNMENT AND DEASSIGN MULTIPLY ASSIGNED ROWS DO 130 J=1,N J0=N-J+1 I=ASSIGN(J0) IF (ROWLAB(I).NE.0) THEN ROWLAB(I)=-ABS(ROWLAB(I)) ASSIGN(J0)=0 ELSE ROWLAB(I)=J0 END IF 130 CONTINUE C CONSTRUCT CANDIDATE LIST AND CORRECT COLUMN PRICES NEWNOL=0 DO 150 I=1,N CURCOL=ROWLAB(I) IF (CURCOL.EQ.0) THEN NEWNOL=NEWNOL+1 LIST(NEWNOL)=I GOTO 150 END IF IF (CURCOL.LT.0) THEN ROWLAB(I)=-CURCOL ELSE MAX2=ISMALL FSTARC=FOUT(I) LSTARC=FOUT(I+1)-1 DO 140 ARC=FSTARC,LSTARC J=END(ARC) IF (J.NE.CURCOL) THEN IF (PCOL(J)-COST(ARC).GT.MAX2) THEN MAX2=PCOL(J)-COST(ARC) END IF END IF 140 CONTINUE PROW(I)=MAX2 PCOL(CURCOL)=PCOL(CURCOL)+MAX2 END IF 150 CONTINUE IF (NEWNOL.EQ.0) GO TO 2000 C ************************************************************ C END OF INITIALIZATION; DO A NUMBER OF AUCTION CYCLES C WHICH IS SET BELOW IN THE PARAMETER AUCTNUM BASED ON THE C SPARSITY OF THE PROBLEM IF (N*N.LT.10*A) AUCTNUM=2 IF ((N*N.GE.10*A).AND.(N*N.LT.25*A)) AUCTNUM=3 IF ((N*N.GE.25*A).AND.(N*N.LT.50*A)) AUCTNUM=4 IF ((N*N.GE.50*A).AND.(N*N.LT.100*A)) AUCTNUM=5 IF (N*N.GE.100*A) AUCTNUM=6 C TO RUN A PURE SEQUENTIAL SHORTEST PATH METHOD SET C AUCTNUM=0 IF (AUCTNUM.EQ.0) GOTO 1000 C START OF A NEW CYCLE 80 NOLIST=NEWNOL I=1 NEWNOL=0 C TAKE UP A ROW FOR ITERATION 100 ROW=LIST(I) I=I+1 FSTARC=FOUT(ROW) LSTARC=FOUT(ROW+1)-1 BSTCOL=END(FSTARC) IF (FSTARC.EQ.LSTARC) THEN OLDROW=ASSIGN(BSTCOL) ASSIGN(BSTCOL)=ROW ROWLAB(ROW)=BSTCOL PROW(ROW)=ISMALL PCOL(BSTCOL)=ISMALL+COST(FSTARC) IF (OLDROW.GT.0) THEN I=I-1 LIST(I)=OLDROW ROWLAB(OLDROW)=0 END IF ISIMPL=ISIMPL+1 GO TO 100 END IF SNDARC=FSTARC+1 SNDCOL=END(SNDARC) MAX1=PCOL(BSTCOL)-COST(FSTARC) MAX2=PCOL(SNDCOL)-COST(SNDARC) IF (MAX1.LT.MAX2) THEN TMAX=MAX1 MAX1=MAX2 MAX2=TMAX FSTCOL=BSTCOL BSTCOL=SNDCOL SNDCOL=FSTCOL END IF IF (SNDARC.LT.LSTARC) THEN TRDARC=SNDARC+1 DO 108 ARC=TRDARC,LSTARC CURCOL=END(ARC) TMAX=PCOL(CURCOL)-COST(ARC) IF (TMAX.GT.MAX2) THEN IF (TMAX.GT.MAX1) THEN MAX2=MAX1 MAX1=TMAX SNDCOL=BSTCOL BSTCOL=CURCOL ELSE MAX2=TMAX SNDCOL=CURCOL END IF END IF 108 CONTINUE END IF PROW(ROW)=MAX2 OLDROW=ASSIGN(BSTCOL) INCR=MAX1-MAX2 IF (INCR.GT.0) THEN PCOL(BSTCOL)=PCOL(BSTCOL)-INCR ASSIGN(BSTCOL)=ROW ROWLAB(ROW)=BSTCOL IF (OLDROW.GT.0) THEN I=I-1 LIST(I)=OLDROW ROWLAB(OLDROW)=0 ISIMPL=ISIMPL+1 GOTO 100 END IF ELSE IF (OLDROW.GT.0) THEN BSTCOL=SNDCOL OLDROW=ASSIGN(BSTCOL) END IF IF (OLDROW.EQ.0) THEN ASSIGN(BSTCOL)=ROW ROWLAB(ROW)=BSTCOL ELSE NEWNOL=NEWNOL+1 LIST(NEWNOL)=ROW END IF END IF ISIMPL=ISIMPL+1 IF (I.LE.NOLIST) GOTO 100 C END OF A NAIVE AUCTION CYCLE COUNT=COUNT+1 IF (NEWNOL.EQ.0) THEN NASSIH=NEWNOL GOTO 2000 END IF IF (COUNT.LT.AUCTNUM) GOTO 80 C **************************************************************** C C ******* END OF NAIVE AUCTION PART ******* C C C ******* START OF SEQ. SH. PATH METHOD ******* C C **************************************************************** 1000 NASSIH=NEWNOL X IHPRC=0 1020 DO 180 I=1,NEWNOL SCAN(I)=LIST(I) 180 CONTINUE DO 190 J=1,N MARG(J)=1 190 CONTINUE NOSCAN=NEWNOL IL1=1 IL2=NOSCAN 1030 DO 1060 I=IL1,IL2 CURROW=SCAN(I) ROWPR=PROW(CURROW) FSTARC=FOUT(CURROW) LSTARC=FOUT(CURROW+1)-1 DO 1050 ARC=FSTARC,LSTARC CURCOL=END(ARC) OLMARG=MARG(CURCOL) IF (OLMARG.EQ.0) GO TO 1050 TMARG=PCOL(CURCOL)-ROWPR-COST(ARC) IF (TMARG.EQ.0) THEN IF (ASSIGN(CURCOL).EQ.0) THEN GO TO 1500 C PERFORM AUGMENTATION ELSE MARG(CURCOL)=0 NOSCAN=NOSCAN+1 SCAN(NOSCAN)=ASSIGN(CURCOL) COLLAB(CURCOL)=CURROW END IF ELSE IF ((OLMARG.GT.0).OR.(TMARG.GT.OLMARG)) THEN MARG(CURCOL)=TMARG COLLAB(CURCOL)=CURROW END IF END IF 1050 CONTINUE 1060 CONTINUE C CURRENT SCAN PHASE COMPLETE; CHECK FOR MORE ROWS TO SCAN IF (IL2.LT.NOSCAN) THEN IL1=IL2+1 IL2=NOSCAN GO TO 1030 END IF C PRICE CHANGE IPRC=IPRC+1 X IHPRC=IHPRC+1 INCR=ISMALL DO 200 J=1,N TMARG=MARG(J) IF ((TMARG.LT.0).AND.(TMARG.GT.INCR)) INCR=TMARG 200 CONTINUE DO 210 I=1,NOSCAN PROW(SCAN(I))=PROW(SCAN(I))+INCR 210 CONTINUE DO 220 J=1,N IF (MARG(J).EQ.0) PCOL(J)=PCOL(J)+INCR 220 CONTINUE DO 1400 J=1,N IF (MARG(J).GE.0) GO TO 1400 MARG(J)=MARG(J)-INCR IF (MARG(J).EQ.0) THEN IF (ASSIGN(J).EQ.0) THEN CURCOL=J CURROW=COLLAB(CURCOL) C PERFORM AUGMENTATION GO TO 1500 ELSE NOSCAN=NOSCAN+1 SCAN(NOSCAN)=ASSIGN(J) END IF END IF 1400 CONTINUE IL1=IL2+1 IL2=NOSCAN GO TO 1030 C END PRICE CHANGE C AUGMENTATION STARTS HERE 1500 IF (ROWLAB(CURROW).EQ.0) THEN ASSIGN(CURCOL)=CURROW ROWLAB(CURROW)=CURCOL GO TO 1700 END IF TA=ROWLAB(CURROW) ASSIGN(CURCOL)=CURROW ROWLAB(CURROW)=CURCOL CURCOL=TA CURROW=COLLAB(CURCOL) GO TO 1500 1700 IF (NEWNOL.EQ.1) GO TO 2000 DO 240 I=1,NEWNOL-1 IF (LIST(I).EQ.CURROW) THEN DO 250 K=I+1,NEWNOL LIST(K-I)=LIST(K) 250 CONTINUE DO 260 K=1,I-1 LIST(NEWNOL-I+K)=SCAN(K) 260 CONTINUE NEWNOL=NEWNOL-1 GO TO 1020 END IF 240 CONTINUE NEWNOL=NEWNOL-1 GO TO 1020 C END AUGMENTATION C ******************************************************* C C EXIT ROUTINE. CHECKS FOR OPTIMALITY OF THE SOLUTION C CALCULATES THE OPTIMAL COST, AND COMPILES SOLUTION C STATISTICS. THE USER MAY REPLACE THIS BY CODE THAT C WRITES AT THE APPROPRIATE DEVICE THE OPTIMAL SOLUTION C CONTAINED IN THE ARRAY ASSIGN(.). C C ******************************************************* 2000 CONTINUE TT =(LONG(362) - TIMER)/60 PRINT*, 'TOTAL TIME = ',TT, ' secs.' PRINT*,'NO OF NAIVE AUCTION ITERATIONS:',ISIMPL PRINT*,'NO OF SH. PATH ITERATIONS:',NASSIH X PRINT*,'NO OF PRICE CHANGES:',IHPRC C CHECK FEASIBILITY OF SOLUTION & CALCULATE COST DO 265 I=1,N ROWLAB(I)=0 265 CONTINUE NOASS=0 DO 280 J=1,N IF (ASSIGN(J).GT.0) THEN ROWLAB(ASSIGN(J))=J NOASS=NOASS+1 END IF 280 CONTINUE IF (NOASS.LT.N) THEN PRINT*,'THE NUMBER OF ASSIGNED COLUMNS IS',NOASS,'(TOO SMALL)' END IF TCOST=0 DO 300 I=1,N J=ROWLAB(I) IF (J.EQ.0) THEN PRINT*,'ROW ',I,' IS UNASSIGNED' END IF FSTARC=FOUT(I) LSTARC=FOUT(I+1)-1 MINCOST=ILARGE DO 310 ARC=FSTARC,LSTARC CURCOL=END(ARC) TMARG=PROW(I)-PCOL(CURCOL)+COST(ARC) IF (TMARG.LT.0) THEN PRINT*,'COMPL. SLACKNESS VIOLATION: ROW',I,' COLUMN',CURCOL END IF IF (CURCOL.EQ.J) THEN IF (MINCOST.GT.COST(ARC)) THEN MINCOST=COST(ARC) ROWMARG=TMARG END IF END IF 310 CONTINUE IF (ROWMARG.NE.0) THEN PRINT*,'COMPL. SLACK. VIOLATION AT ASSIGNED ARC OF ROW',I END IF IF (MAXSET) THEN TCOST=TCOST-MINCOST ELSE TCOST=TCOST+MINCOST END IF 300 CONTINUE WRITE(9,2100) TCOST 2100 FORMAT(' ASSIGNMENT COST=',F14.2) PRINT *,'**********************************************' PRINT *, 'PROGRAM ENDED; PRESS ' PAUSE STOP END SUBROUTINE SETRAN(ISEED) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C PORTABLE CONGRUENTIAL (UNIFORM) RANDOM NUMBER GENERATOR: C NEXT_VALUE = [(7**5) * PREVIOUS_VALUE] MODULO[(2**31)-1] C C THIS GENERATOR CONSISTS OF TWO ROUTINES: C (1) SETRAN - INITIALIZES CONSTANTS AND SEED C (2) RRAN - GENERATES A REAL RANDOM NUMBER C C THE GENERATOR REQUIRES A MACHINE WITH AT LEAST 32 BITS OF PRECISION. C THE SEED (ISEED) MUST BE IN THE RANGE (1,(2**31)-1). C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IF(ISEED.LT.1) STOP 77 MULT=16807 MODUL=2147483647 I15=2**15 I16=2**16 JRAN=ISEED RETURN END C REAL FUNCTION RAN() IMPLICIT REAL*4 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C RAN GENERATES A REAL RANDOM NUMBER BETWEEN 0 AND 1 C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IXHI=JRAN/I16 IXLO=JRAN-IXHI*I16 IXALO=IXLO*MULT LEFTLO=IXALO/I16 IXAHI=IXHI*MULT IFULHI=IXAHI+LEFTLO IRTLO=IXALO-LEFTLO*I16 IOVER=IFULHI/I15 IRTHI=IFULHI-IOVER*I15 JRAN=((IRTLO-MODUL)+IRTHI*I16)+IOVER IF(JRAN.LT.0) JRAN=JRAN+MODUL RAN = FLOAT(JRAN)/FLOAT(MODUL) RETURN END *************************************************************************** Max-Flow Codes *************************************************************************** FORD-FULKERSON This code implements the Ford-Fulkerson method for solving the max-flow problem (cf.\ Section 1.2.2). C ******* SAMPLE CALLING PROGRAM FOR FORD-FULKERSON ******* PARAMETER (MAXNODES=10000, MAXARCS=40000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,SOURCE,SINK COMMON /STATS/ NITER COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK3/ PRDCSR COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /LABELS/ LABEL COMMON /MARKS/ MARK INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER U(MAXARCS) INTEGER F(MAXARCS) INTEGER PRDCSR(MAXNODES) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER LABEL(MAXNODES) LOGICAL MARK(MAXNODES) REAL*8 TT,TIMER PRINT*,'FORD-FULKERSON METHOD FOR MAX-FLOW' PRINT*,'**********************************' PRINT *,'READING PROBLEM DATA' OPEN(13,FILE='FOR013.DAT',STATUS='OLD') REWIND(13) LARGE=500000000 C READ NUMBER OF NODES AND ARCS READ(13,1010) N,NA C READ START, END, COST, AND CAPACITY OF EACH ARC DO 20 I=1,NA READ(13,1020) STARTN(I),ENDN(I),DUMMY,U(I) 20 CONTINUE ENDFILE(13) REWIND(13) 1000 FORMAT(1I8) 1010 FORMAT(2I8) 1020 FORMAT(4I8) PRINT*,'THE NUMBER OF NODES IS = ',N 41 PRINT*,'ENTER THE SUPERSOURCE NODE' READ*,SOURCE IF ((SOURCE.LE.0).OR.(SOURCE.GT.N)) GOTO 41 42 PRINT*,'ENTER THE SUPERSINK NODE' READ*,SINK IF ((SINK.LE.0).OR.(SINK.GT.N)) GOTO 42 IF (SINK.EQ.R) GOTO 42 C THE SOURCE AND SINK NODES WILL NOT BE ITERATED ON PRINT *, 'RESTRUCTURING THE DATA ' CALL INIDAT C SET FLOWS TO ZERO DO 50 ARC=1,NA F(ARC)=0 50 CONTINUE PRINT *,'**********************************' TIMER = LONG(362) CALL FORD_FULK TT = FLOAT(LONG(362) - TIMER)/60 PRINT*, 'TOTAL TIME = ',TT, ' secs.' C ***** COMPUTE MAX-FLOW ***** MAX_FLOW=0 DO 60 ARC=1,NA IF (STARTN(ARC).EQ.SOURCE) MAX_FLOW=MAX_FLOW+F(ARC) 60 CONTINUE MAX_FLOW2=0 DO 70 ARC=1,NA IF (ENDN(ARC).EQ.SINK) MAX_FLOW2=MAX_FLOW2+F(ARC) 70 CONTINUE PRINT *,'# OF ITERATIONS = ',NITER PRINT *, 'MAX-FLOW = ',MAX_FLOW2 PRINT *,'**********************************' IF (MAXFLOW.NE.MAXFLOW2) THEN PRINT*,'SOURCE FLOW NOT EQUAL TO SINK FLOW' ENDIF PRINT *, 'PROGRAM ENDED; PRESS TO EXIT' PAUSE END C **************************************************** SUBROUTINE INIDAT C THIS SUBROUTINE USES THE DATA ARRAYS STARTN AND ENDN C TO CONSTRUCT AUXILIARY DATA ARRAYS FOUT, NXTOU, FIN, AND C NXTIN. IN THIS SUBROUTINE WE C ARBITRARILY ORDER THE ARCS LEAVING EACH NODE AND STORE C THIS INFORMATION IN FOUT AND NXTOU. SIMILARLY, WE ARBITRA- C RILLY ORDER THE ARCS ENTERING EACH NODE AND STORE THIS C INFORMATION IN FIN AND NXTIN. AT THE COMPLETION OF THE C CONSTRUCTION, WE HAVE C C FOUT(I) = FIRST ARC LEAVING NODE I. C NXTOU(J) = NEXT ARC LEAVING THE HEAD NODE OF ARC J. C FIN(I) = FIRST ARC ENTERING NODE I. C NXTIN(J) = NEXT ARC ENTERING THE TAIL NODE OF ARC J. PARAMETER (MAXNODES=10000, MAXARCS=40000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER FINALIN(MAXNODES) INTEGER FINALOU(MAXNODES) DO 20 NODE=1,N FIN(NODE)=0 FOUT(NODE)=0 FINALIN(NODE)=0 FINALOU(NODE)=0 20 CONTINUE DO 30 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) IF (FOUT(START).NE.0) THEN NXTOU(FINALOU(START))=ARC ELSE FOUT(START)=ARC END IF IF (FIN(END).NE.0) THEN NXTIN(FINALIN(END))=ARC ELSE FIN(END)=ARC END IF FINALOU(START)=ARC FINALIN(END)=ARC NXTIN(ARC)=0 NXTOU(ARC)=0 30 CONTINUE C RETURN END C **************************************************************** C C BASIC FORD-FULKERSON METHOD FOR MAX-FLOW. C C **************************************************************** SUBROUTINE FORD_FULK IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,SOURCE,SINK COMMON /STATS/ NITER COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK3/ PRDCSR COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /MARKS/ MARK COMMON /LABELS/ LABEL INTEGER STARTN(1),ENDN(1),U(1),F(1),FIN(1),FOUT(1) INTEGER NXTIN(1),NXTOU(1),PRDCSR(1),LABEL(1) LOGICAL MARK(1) NITER=0 DO 10 I=1,N MARK(I)=.FALSE. 10 CONTINUE C START OF NEW ITERATION 15 NLABEL=1 NSCAN=1 MARK(SOURCE)=.TRUE. LABEL(1)=SOURCE 20 CONTINUE C SCAN A NEW NODE NODE=LABEL(NSCAN) C SCAN OUTGOING ARCS OF NODE ARC=FOUT(NODE) 30 IF (ARC.GT.0) THEN NODE2=ENDN(ARC) IF ((.NOT.MARK(NODE2)).AND.(F(ARC).LT.U(ARC))) THEN PRDCSR(NODE2)=ARC IF (NODE2.EQ.SINK) THEN CALL AUGMENT NITER=NITER+1 DO 40 I=1,NLABEL MARK(LABEL(I))=.FALSE. 40 CONTINUE GOTO 15 ELSE MARK(NODE2)=.TRUE. NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 END IF END IF ARC=NXTOU(ARC) GOTO 30 END IF C SCAN INCOMING ARCS OF NODE ARC=FIN(NODE) 50 IF (ARC.GT.0) THEN NODE2=STARTN(ARC) IF ((.NOT.MARK(NODE2)).AND.(F(ARC).GT.0)) THEN PRDCSR(NODE2)=-ARC IF (NODE2.EQ.SINK) THEN CALL AUGMENT NITER=NITER+1 DO 60 I=1,NLABEL MARK(LABEL(I))=.FALSE. 60 CONTINUE GOTO 15 ELSE MARK(NODE2)=.TRUE. NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 END IF END IF ARC=NXTIN(ARC) GOTO 50 END IF C CHECK FOR TERMINATION; SCAN A NEW NODE IF (NSCAN.EQ.NLABEL) THEN RETURN END IF NSCAN=NSCAN+1 GOTO 20 END SUBROUTINE AUGMENT IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,SOURCE,SINK COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK3/ PRDCSR INTEGER STARTN(1),ENDN(1),U(1),F(1),PRDCSR(1) DX=LARGE CURNODE=SINK 10 IF (CURNODE.NE.SOURCE) THEN ARC=PRDCSR(CURNODE) IF (ARC.GT.0) THEN INCR=U(ARC)-F(ARC) IF (DX.GT.INCR) DX=INCR CURNODE=STARTN(ARC) ELSE ARC=-ARC INCR=F(ARC) IF (DX.GT.INCR) DX=INCR CURNODE=ENDN(ARC) END IF GOTO 10 END IF CURNODE=SINK 20 IF (CURNODE.NE.SOURCE) THEN ARC=PRDCSR(CURNODE) IF (ARC.GT.0) THEN F(ARC)=F(ARC)+DX CURNODE=STARTN(ARC) ELSE ARC=-ARC F(ARC)=F(ARC)-DX CURNODE=ENDN(ARC) END IF GOTO 20 END IF RETURN END *************************************************************************** $\e$-RELAX-MF *************************************************************************** This code implements the $\e$-relaxation method for the max-flow problem (cf.\ Section 4.5). Here, $\e=1$ throughout the algorithm. C ******* SAMPLE CALLING PROGRAM FOR E-RELAX/MAX-FLOW ******* C C THIS VERSION USES A CYCLIC QUEUE, A SHORTEST PATH C (BREADTH FIRST) INITIALIZATION C AND A ROUTINE THAT DOES GLOBAL RELABELING AND ALSO C DETECTS A SATURATED CUT AND TERMINATES EARLY C IT THEN FINDS A MAXIMUM FLOW USING A "REVERSE'ALGORITHM C C THIS CODE IS A SOMEWHAT DIFFERENT VERSION OF THE ONE C IN THE BOOK "LINEAR NETWORK OPTIMIZATION: ALGORITHMS AND C CODES, M.I.T. PRESS, 1991, BY DIMITRI P. BERTSEKAS C C *************************************************** PARAMETER (MAXNODES=20000, MAXARCS=120000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,SOURCE,SINK COMMON /STATS/ NITER,CUTEXIT,CUTSEARCH,CUTCOUNT,NAUG COMMON /CHECK/ CHECKTIME,TOT_TIME COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /PRICES/ P COMMON /BLK3/ SURPLUS COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /BLK11/ FPUSHF COMMON /BLK12/ NXTPUSHF COMMON /BLK13/ FPUSHB COMMON /BLK14/ NXTPUSHB COMMON /LABELS/ LABEL COMMON /MARKS/ MARK COMMON /QUEUE/ NXTQUEUE COMMON /PRED/ PRED C RANDOM GENERATOR STUFF COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN REAL RAN INTEGER NXTQUEUE(MAXNODES) INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER U(MAXARCS) INTEGER F(MAXARCS) INTEGER SURPLUS(MAXNODES) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER P(MAXNODES) INTEGER FPUSHF(MAXNODES) INTEGER NXTPUSHF(MAXARCS) INTEGER FPUSHB(MAXNODES) INTEGER NXTPUSHB(MAXARCS) INTEGER LABEL(MAXNODES) INTEGER PRED(MAXNODES) LOGICAL MARK(MAXNODES) REAL*8 TT,TOT_TIME,TIMER,TOTAL,TOTAL_C,CHECKTIME PRINT*,'E-RELAXATION METHOD FOR MAX-FLOW' PRINT*,'USES FIFO QUEUE AND GLOBAL RELABELING' PRINT*,'********************************************' LARGE=1000000000 C ************************************************************ C C READ THE PROBLEM DATA C C THE PROBLEM IS SPECIFIED BY THE FOLLOWING VARIABLES: C C N: THE NUMBER OF NODES C NA: THE NUMBER OF ARCS C SOURCE: THE SOURCE NODE OF THE MAX-FLOW PROBLEM C SINK: THE SINK NODE OF THE MAX-FLOW PROBLEM C STARTN(ARC): THE START NODE OF ARC C ENDN(ARC): THE END NODE OF ARC C U(ARC): THE CAPACITY OF ARC (IT IS LATER CHANGED TO RESIDUAL CAPACITY) C C PRINT*,'READING PROBLEM DATA' OPEN(13,FILE='FOR013.DAT',STATUS='OLD') REWIND(13) C READ NUMBER OF NODES AND ARCS READ(13,*) N,NA DO 5 I=1,NA READ(13,*) STARTN(I),ENDN(I),DUMMY,U(I) 5 CONTINUE ENDFILE(13) REWIND(13) SOURCE=1 SINK=N PRINT*,'********************************************' PRINT*,'NUMBER OF NODES =',N,' NUMBER OF ARCS =',NA C C C END OF READING THE PROBLEM DATA C C INITIALIZE PRICES & SURPLUSES DO 20 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 P(NODE)=N SURPLUS(NODE)=0 20 CONTINUE C SET FLOWS TO UPPER BOUND OR LOWER BOUND TO SATISFY E-CS DO 50 ARC=1,NA IF (STARTN(ARC).EQ.SOURCE) THEN F(ARC)=U(ARC) SURPLUS(ENDN(ARC))=SURPLUS(ENDN(ARC))+F(ARC) SURPLUS(SOURCE)=SURPLUS(SOURCE)-F(ARC) ELSE F(ARC)=0 END IF 50 CONTINUE C CREATE NECESSARY DATA STRUCTURES CALL INIDAT TIMER = LONG(362) CALL E_RELAX_MF TT = FLOAT(LONG(362) - TIMER)/60 PRINT*,'TOTAL TIME = ',TT, ' secs' PRINT*,'TIME TO FIND A MIN CUT = ',TOT_TIME PRINT*,'TIME TO CHECK FOR A MIN CUT = ',CHECKTIME C ***** COMPUTE MAX-FLOW ***** MAX_FLOW=0 ARC=FIN(SINK) 70 IF (ARC.GT.0) THEN MAX_FLOW=MAX_FLOW+F(ARC) ARC=NXTIN(ARC) GOTO 70 END IF IF (MAX_FLOW.NE.SURPLUS(SINK)) THEN PRINT*,'ERROR IN THE CALCULATION OF THE MAX-FLOW' END IF PRINT *, 'MAX-FLOW = ',MAX_FLOW PRINT *,'# OF FLOW CHANGES PER ARC = ',FLOAT(NAUG)/FLOAT(NA) PRINT *,'# OF PRICE RISES PER NODE = ',FLOAT(NITER)/FLOAT(N) PRINT *,'# OF SEARCHES FOR A SATURATED CUT = ',CUTCOUNT PRINT *,'************************************************' C ***** CHECK THAT ALL SURPLUSES ARE ZERO **** DO 80 NODE=1,N IF ((NODE.NE.SOURCE).AND.(NODE.NE.SINK)) THEN IF (SURPLUS(NODE).NE.0) THEN PRINT*,'SURPLUS OF NODE ',NODE,' IS NONZERO' END IF END IF 80 CONTINUE PRINT *, 'PROGRAM ENDED; PRESS TO EXIT' PAUSE END C **************************************************************** C C BASIC E-RELAXATION FOR MAX-FLOW. C THIS VERSION USES A QUEUE TO SELECT NODES TO ITERATE ON, C AND MAINTAINS PUSH LISTS. NODES JOIN THE QUEUE AT THE BOTTOM. C THIS VERSION USES A SHORTEST PATH (BREADTH FIRST) INITIALIZATION C AND A ROUTINE THAT DETECTS A SATURATED CUT AND C PERFORMS GLOBAL RELABELING C C THIS CODE IS A SOMEWHAT DIFFERENT VERSION OF THE ONE C IN THE BOOK "LINEAR NETWORK OPTIMIZATION: ALGORITHMS AND C CODES, M.I.T. PRESS, 1991, BY DIMITRI P. BERTSEKAS C C **************************************************************** SUBROUTINE E_RELAX_MF PARAMETER (MAXNODES=20000, MAXARCS=120000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,SOURCE,SINK COMMON /STATS/ NITER,CUTEXIT,CUTSEARCH,CUTCOUNT,NAUG COMMON /CHECK/ CHECKTIME,TOT_TIME COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /PRICES/ P COMMON /BLK3/ SURPLUS COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /BLK11/ FPUSHF COMMON /BLK12/ NXTPUSHF COMMON /BLK13/ FPUSHB COMMON /BLK14/ NXTPUSHB COMMON /MARKS/ MARK COMMON /LABELS/ LABEL COMMON /QUEUE/ NXTQUEUE COMMON /PRED/ PRED INTEGER NXTQUEUE(MAXNODES) INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER U(MAXARCS) INTEGER F(MAXARCS) INTEGER SURPLUS(MAXNODES) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER P(MAXNODES) INTEGER FPUSHF(MAXNODES) INTEGER NXTPUSHF(MAXARCS) INTEGER FPUSHB(MAXNODES) INTEGER NXTPUSHB(MAXARCS) INTEGER LABEL(MAXNODES) INTEGER PRED(MAXNODES) LOGICAL MARK(MAXNODES) REAL*8 TIMER,TOT_TIME,TIMER1,CHECKTIME LOGICAL OUTFLAG,RELABEL TIMER1=LONG(362) C INITIALIZE COUNT OF NUMBER OF UP ITERATIONS PERFORMED NITER = 0 NAUG=0 COUNT=0 CUTCOUNT=0 CUTEXIT=0 THRESH=N CHECKTIME=0 TOT_TIME=0 RELABEL=.TRUE. C INITIALIZE PRICES BY A BREADTH FIRST SEARCH METHOD DO 51 NODE=1,N MARK(NODE)=.FALSE. 51 CONTINUE MARK(SOURCE)=.TRUE. MARK(SINK)=.TRUE. P(SOURCE)=N P(SINK)=0 NLABEL=1 LABEL(1)=SINK NSCAN=0 55 CONTINUE NSCAN=NSCAN+1 NODE=LABEL(NSCAN) PRICE=P(NODE)+1 ARC=FIN(NODE) 58 IF (ARC.GT.0) THEN START=STARTN(ARC) IF (.NOT.MARK(START)) THEN NLABEL=NLABEL+1 LABEL(NLABEL)=START P(START)=PRICE MARK(START)=.TRUE. END IF ARC=NXTIN(ARC) GOTO 58 END IF IF (NLABEL.GT.NSCAN) GOTO 55 C ***** QUEUE INITIALIZATION ***** 70 CONTINUE FIRST=0 DO 82 I=1,NLABEL NODE=LABEL(I) IF ((NODE.NE.SOURCE).AND.(NODE.NE.SINK)) THEN IF (SURPLUS(NODE).GT.0) THEN IF (FIRST.EQ.0) THEN FIRST=NODE ELSE NXTQUEUE(LAST)=NODE END IF LAST=NODE ELSE NXTQUEUE(NODE)=0 END IF ELSE NXTQUEUE(NODE)=-1 END IF 82 CONTINUE NXTQUEUE(LAST)=FIRST IF (FIRST.EQ.0) RETURN I=FIRST PREVNODE=LAST C *********************************************************** C C START OF THE MAIN ALGORITHM. C DO SINGLE NODE ITERATIONS UNTIL TERMINATION. C C *********************************************************** 100 CONTINUE C TAKE UP NEXT NODE (NODE I) FOR ITERATION SURPI=SURPLUS(I) PRICE = P(I) IF ((SURPI .GT. 0).AND.(PRICE.LT.N)) THEN C PRINT*,'NODE ',I,' START PRICE = ',PRICE ARCF = FPUSHF(I) ARCB = FPUSHB(I) C ********************** D - PUSHES ************************ C START BY TRYING TO PUSH AWAY FLOW ON ARCS THAT WERE C ADMISSIBLE AT THE END OF THE LAST ITERATION (IF ANY) C AT THIS NODE. WE MUST CHECK THAT THEY ARE STILL ADMISSIBLE. 120 IF ((SURPI .GT. 0) .AND. (ARCF .GT. 0)) THEN RESID = U(ARCF) - F(ARCF) J = ENDN(ARCF) IF ((PRICE-P(J).EQ.1).AND.(RESID.GT.0)) THEN NAUG=NAUG+1 IF (SURPI .GE. RESID) THEN F(ARCF) = U(ARCF) SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + RESID IF (NXTQUEUE(J).EQ.0) THEN IF (SURPLUS(J).GT.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ARCF = NXTPUSHF(ARCF) ELSE F(ARCF) = F(ARCF) + SURPI SURPLUS(J) = SURPLUS(J) + SURPI IF (NXTQUEUE(J).EQ.0) THEN IF (SURPLUS(J).GT.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF SURPI = 0 ENDIF ELSE ARCF = NXTPUSHF(ARCF) ENDIF GOTO 120 ENDIF 121 IF ((SURPI .GT. 0) .AND. (ARCB .GT. 0)) THEN J = STARTN(ARCB) IF ((PRICE-P(J).EQ.1).AND.(F(ARCB).GT.0)) THEN NAUG=NAUG+1 IF (SURPI .GE. F(ARCB)) THEN SURPI = SURPI - F(ARCB) SURPLUS(J) = SURPLUS(J) + F(ARCB) IF (NXTQUEUE(J).EQ.0) THEN IF (SURPLUS(J).GT.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF F(ARCB) = 0 ARCB = NXTPUSHB(ARCB) ELSE F(ARCB) = F(ARCB) - SURPI SURPLUS(J) = SURPLUS(J) + SURPI IF (NXTQUEUE(J).EQ.0) THEN IF (SURPLUS(J).GT.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF SURPI = 0 ENDIF ELSE ARCB = NXTPUSHB(ARCB) ENDIF GOTO 121 ENDIF C NOW WE ENTER A REPEAT...UNTIL SORT OF LOOP UNTIL WE C HAVE THE SURPLUS REDUCED TO ZERO. C ********************* PRICE RISE ********************* 129 CONTINUE C FIRST, TRY A (POSSIBLY DEGENERATE) PRICE RISE. 130 IF ((ARCF .EQ. 0) .AND. (ARCB .EQ. 0)) THEN NITER = NITER + 1 PRICE = LARGE ARC = FOUT(I) 131 IF (ARC .GT. 0) THEN IF (F(ARC) .LT. U(ARC)) THEN XP = P(ENDN(ARC)) + 1 IF (XP .LT. PRICE) THEN PRICE = XP ARCF = ARC NXTPUSHF(ARC)=0 ELSE IF (XP .EQ. PRICE) THEN NXTPUSHF(ARC) = ARCF ARCF = ARC END IF ENDIF ENDIF ARC = NXTOU(ARC) GOTO 131 ENDIF ARC = FIN(I) 132 IF (ARC .GT. 0) THEN IF (F(ARC) .GT. 0) THEN XP = P(STARTN(ARC))+1 IF (XP .LT. PRICE) THEN PRICE = XP ARCB = ARC ARCF = 0 NXTPUSHB(ARC) = 0 ELSE IF (XP .EQ. PRICE) THEN NXTPUSHB(ARC) = ARCB ARCB = ARC END IF ENDIF ENDIF ARC = NXTIN(ARC) GOTO 132 ENDIF ENDIF C ********************** D - PUSHES ******************** C IF THE STEP WAS NOT DEGENERATE, TRY AND DO SOME C PUSHING. WE DO NOT REUSE THE CODE ABOVE BECAUSE IT IS C NOT NECESSARY TO CHECK IF THE ARCS USED ARE ADMISSIBLE. IF (SURPI .GT. 0) THEN 141 IF (ARCF .GT. 0) THEN RESID = U(ARCF) - F(ARCF) NAUG=NAUG+1 J = ENDN(ARCF) IF (SURPI .GE. RESID) THEN F(ARCF) = U(ARCF) SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + RESID IF (NXTQUEUE(J).EQ.0) THEN IF (SURPLUS(J).GT.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ARCF = NXTPUSHF(ARCF) IF (SURPI .GT. 0) GOTO 141 ELSE F(ARCF) = F(ARCF) + SURPI SURPLUS(J) = SURPLUS(J) + SURPI IF (NXTQUEUE(J).EQ.0) THEN IF (SURPLUS(J).GT.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF SURPI = 0 ENDIF END IF 142 IF ((SURPI .GT. 0).AND.(ARCB .GT. 0)) THEN J = STARTN(ARCB) NAUG=NAUG+1 IF (SURPI .GE. F(ARCB)) THEN SURPI = SURPI - F(ARCB) SURPLUS(J) = SURPLUS(J) + F(ARCB) IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF F(ARCB) = 0 ARCB = NXTPUSHB(ARCB) ELSE F(ARCB) = F(ARCB) - SURPI SURPLUS(J) = SURPLUS(J) + SURPI IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF SURPI = 0 ENDIF GOTO 142 ENDIF C WE'VE DONE ALL THE PUSHING WE CAN AT THIS PRICE C LEVEL. TRY TO DO A (POSSIBLY DEGENERATE) PRICE C INCREASE. GOTO 129 ENDIF C THIS IS THE END OF THE UP ITERATION. IF WE GET HERE, THE C SURPLUS OF I IS 0, AND WE HAVE DONE OUR LAST PRICE RISE. SURPLUS(I) = 0 P(I) = PRICE FPUSHF(I) = ARCF FPUSHB(I) = ARCB ENDIF C ADVANCE THE QUEUE. C IF THE QUEUE IS EMPTY RETURN ALL SURPLUS TO THE SOURCE AND EXIT NXTNODE=NXTQUEUE(I) IF (I.EQ.NXTNODE) THEN IF (CUTEXIT.EQ.1) THEN RETURN ELSE TOT_TIME= TOT_TIME+FLOAT(LONG(362) - TIMER1)/60 GOTO 600 END IF ELSE NXTQUEUE(PREVNODE)=NXTNODE NXTQUEUE(I)=0 I=NXTNODE END IF IF (COUNT.LT.THRESH) THEN COUNT=COUNT+1 GO TO 100 END IF C **************** PERIODICALLY LOOK FOR A SATURATING CUT *************** C AND RECALCULATE THE PRICES TO BE EQUAL TO THE SHORTEST DISTANCES 500 CONTINUE X PRINT*,'** CHECKING FOR EARLY TERMINATION **' TIMER = LONG(362) C UPDATE THE THRESHOLD THRESH=INT(1.1*THRESH) C REINITIALIZE THE COUNT OF ITERATIONS TO THE NEXT CHECK FOR A MIN CUT COUNT=0 CUTCOUNT=CUTCOUNT+1 DO 510 L=1,NLABEL MARK(LABEL(L))=.FALSE. 510 CONTINUE OUTFLAG=.FALSE. NLABEL=1 LABEL(1)=SINK NSCAN=0 550 CONTINUE NSCAN=NSCAN+1 NODE=LABEL(NSCAN) MARK(NODE)=.TRUE. ARC=FIN(NODE) 580 IF (ARC.GT.0) THEN START=STARTN(ARC) C IF THIS IS AN ARC ALONG WHICH FLOW CAN BE PUSHED TO NODE C CHECK IF THE OPP. NODE HAS 0 SURPLUS AND IF SO LABEL THE NODE IF (F(ARC).LT.U(ARC)) THEN IF (.NOT.MARK(START)) THEN IF (RELABEL) THEN IF (SURPLUS(START).GT.0) THEN OUTFLAG=.TRUE. END IF P(START)=P(NODE)+1 ELSE IF (SURPLUS(START).GT.0) THEN CHECKTIME= CHECKTIME+FLOAT(LONG(362) - TIMER)/60 GOTO 100 END IF END IF NLABEL=NLABEL+1 LABEL(NLABEL)=START MARK(START)=.TRUE. END IF END IF ARC=NXTIN(ARC) GOTO 580 END IF ARC=FOUT(NODE) 585 IF (ARC.GT.0) THEN END=ENDN(ARC) C IF THIS IS AN ARC ALONG WHICH FLOW CAN BE PUSHED TO NODE C CHECK IF THE OPP. NODE HAS 0 SURPLUS AND IF SO LABEL THE NODE IF (F(ARC).GT.0) THEN IF (.NOT.MARK(END)) THEN IF (RELABEL) THEN IF (SURPLUS(END).GT.0) THEN OUTFLAG=.TRUE. END IF P(END)=P(NODE)+1 ELSE IF (SURPLUS(END).GT.0) THEN CHECKTIME= CHECKTIME+FLOAT(LONG(362) - TIMER)/60 GOTO 100 END IF END IF NLABEL=NLABEL+1 LABEL(NLABEL)=END MARK(END)=.TRUE. END IF END IF ARC=NXTOU(ARC) GOTO 585 END IF IF (NLABEL.GT.NSCAN) GOTO 550 C IF A NODE W/ POSITIVE SURPLUS IS REACHABLE FROM THE SINK C RELABEL THE UNREACHABLE NODES TO N AND GO BACK FOR MORE ITERATION IF (OUTFLAG) THEN IF (RELABEL) THEN DO 595 NODE=1,N IF (.NOT.MARK(NODE)) THEN IF (P(NODE).LT.N) P(NODE)=N END IF 595 CONTINUE END IF X PRINT*,'OUT OF CHECK FOR EARLY TERMINATION' CHECKTIME= CHECKTIME+FLOAT(LONG(362) - TIMER)/60 GOTO 100 ELSE X PRINT*,'SATURATED CUT FOUND' X PRINT*,'# OF NODES ON SINK SIDE OF THE CUT =',NLABEL CHECKTIME= CHECKTIME+FLOAT(LONG(362) - TIMER)/60 TOT_TIME= TOT_TIME+FLOAT(LONG(362) - TIMER1)/60 GOTO 600 END IF C C ********************** FIND A MAX-FLOW ************************** C C SATURATED CUT HAS BEEN FOUND BUT SOME NODES MAY HAVE NONZERO SURPLUS. C TO FIND A MAX-FLOW, WE RETURN THE EXCESS FLOW TO THE SOURCE BY C USING THE SAME ALGORITHM C 600 CONTINUE CUTEXIT=1 RELABEL=.FALSE. COUNT=0 THRESH=LARGE DO 620 NODE=1,N MARK(NODE)=.FALSE. 620 CONTINUE MARK(SOURCE)=.TRUE. MARK(SINK)=.TRUE. P(SOURCE)=0 NLABEL=1 LABEL(1)=SOURCE NSCAN=0 650 CONTINUE NSCAN=NSCAN+1 NODE=LABEL(NSCAN) PRICE=P(NODE)+1 ARC=FIN(NODE) 660 IF (ARC.GT.0) THEN IF (U(ARC).GT.0) THEN START=STARTN(ARC) IF (.NOT.MARK(START)) THEN NLABEL=NLABEL+1 LABEL(NLABEL)=START P(START)=PRICE MARK(START)=.TRUE. END IF END IF ARC=NXTIN(ARC) GOTO 660 END IF ARC=FOUT(NODE) 670 IF (ARC.GT.0) THEN IF (F(ARC).GT.0) THEN END=ENDN(ARC) IF (.NOT.MARK(END)) THEN NLABEL=NLABEL+1 LABEL(NLABEL)=END P(END)=PRICE MARK(END)=.TRUE. END IF END IF ARC=NXTOU(ARC) GOTO 670 END IF IF (NLABEL.GT.NSCAN) GOTO 650 GOTO 70 END C **************************************************** SUBROUTINE INIDAT C THIS SUBROUTINE USES THE DATA ARRAYS STARTN AND ENDN C TO CONSTRUCT AUXILIARY DATA ARRAYS FOUT, NXTOU, FIN, AND C NXTIN THAT ARE REQUIRED BY E-RELAX-MF. IN THIS SUBROUTINE WE C ARBITRARILY ORDER THE ARCS LEAVING EACH NODE AND STORE C THIS INFORMATION IN FOUT AND NXTOU. SIMILARLY, WE ARBITRA- C RILLY ORDER THE ARCS ENTERING EACH NODE AND STORE THIS C INFORMATION IN FIN AND NXTIN. AT THE COMPLETION OF THE C CONSTRUCTION, WE HAVE C C FOUT(I) = FIRST ARC LEAVING NODE I. C NXTOU(J) = NEXT ARC LEAVING THE HEAD NODE OF ARC J. C FIN(I) = FIRST ARC ENTERING NODE I. C NXTIN(J) = NEXT ARC ENTERING THE TAIL NODE OF ARC J. PARAMETER (MAXNODES=20000, MAXARCS=120000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,SOURCE,SINK COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER FINALIN(MAXNODES),FINALOU(MAXNODES) C DO 20 NODE=1,N FIN(NODE)=0 FOUT(NODE)=0 FINALIN(NODE)=0 FINALOU(NODE)=0 20 CONTINUE C DO 30 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) IF (FOUT(START).NE.0) THEN NXTOU(FINALOU(START))=ARC ELSE FOUT(START)=ARC END IF IF (FIN(END).NE.0) THEN NXTIN(FINALIN(END))=ARC ELSE FIN(END)=ARC END IF FINALOU(START)=ARC FINALIN(END)=ARC NXTIN(ARC)=0 NXTOU(ARC)=0 30 CONTINUE C RETURN END C **************************************************** SUBROUTINE SETRAN(ISEED) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C PORTABLE CONGRUENTIAL (UNIFORM) RANDOM NUMBER GENERATOR: C NEXT_VALUE = [(7**5) * PREVIOUS_VALUE] MODULO[(2**31)-1] C C THIS GENERATOR CONSISTS OF TWO ROUTINES: C (1) SETRAN - INITIALIZES CONSTANTS AND SEED C (2) RRAN - GENERATES A REAL RANDOM NUMBER C C THE GENERATOR REQUIRES A MACHINE WITH AT LEAST 32 BITS OF PRECISION. C THE SEED (ISEED) MUST BE IN THE RANGE (1,(2**31)-1). C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IF(ISEED.LT.1) STOP 77 MULT=16807 MODUL=2147483647 I15=2**15 I16=2**16 JRAN=ISEED RETURN END C REAL FUNCTION RAN() IMPLICIT REAL*4 (A-H,O-Z) , INTEGER*4 (I-N) C*********************************************************************** C RAN GENERATES A REAL RANDOM NUMBER BETWEEN 0 AND 1 C*********************************************************************** COMMON /RANDM/ MULT,MODUL,I15,I16,JRAN IXHI=JRAN/I16 IXLO=JRAN-IXHI*I16 IXALO=IXLO*MULT LEFTLO=IXALO/I16 IXAHI=IXHI*MULT IFULHI=IXAHI+LEFTLO IRTLO=IXALO-LEFTLO*I16 IOVER=IFULHI/I15 IRTHI=IFULHI-IOVER*I15 JRAN=((IRTLO-MODUL)+IRTHI*I16)+IOVER IF(JRAN.LT.0) JRAN=JRAN+MODUL RAN = FLOAT(JRAN)/FLOAT(MODUL) RETURN END *************************************************************************** $\e$-Relaxation Codes *************************************************************************** $\e$-RELAX This code implements the $\e$-relaxation method with $\e$-scaling for the minimum cost flow problem (cf.\ Section 4.5). C ** SAMPLE CALLING PROGRAM FOR E-RELAX ** C C THIS PROGRAM WILL READ A PROBLEM FILE IN STANDARD FORMAT C AND SOLVE IT USING E-RELAX. C C THIS IS AN E-SCALED VERSION THAT USES A CYCLIC QUEUE C C *************************************************** C PARAMETER (MAXNODES=8000, MAXARCS=25000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,FACTOR,EPS, &THRESHSURPLUS,SFACTOR COMMON /STATS/ NCYC,NITER,TOTALITER COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /PRICES/ P COMMON /BLK3/ SURPLUS COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /BLK11/ FPUSHF COMMON /BLK12/ NXTPUSHF COMMON /BLK13/ FPUSHB COMMON /BLK14/ NXTPUSHB COMMON /COST/ COST COMMON /QUEUE/ NXTQUEUE INTEGER NXTQUEUE(MAXNODES) INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER U(MAXARCS) INTEGER F(MAXARCS) INTEGER SURPLUS(MAXNODES) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER P(MAXNODES) INTEGER FPUSHF(MAXNODES) INTEGER NXTPUSHF(MAXARCS) INTEGER FPUSHB(MAXNODES) INTEGER NXTPUSHB(MAXARCS) INTEGER COST(MAXARCS) REAL*8 TT,TIMER DOUBLE PRECISION TCOST,PRODUCT PRINT*,'E-RELAXATION METHOD FOR MIN COST FLOW' PRINT*,'*************************************' PRINT *,'READING PROBLEM DATA' OPEN(13,FILE='FOR013.DAT',STATUS='OLD') REWIND(13) C READ NUMBER OF NODES AND ARCS READ(13,1010) N,NA C READ START, END, COST, AND CAPACITY OF EACH ARC C AND GENERATE MAX COST VALUE MAXCOST=0 DO 5 I=1,NA READ(13,1020) STARTN(I),ENDN(I),COST(I),U(I) COST(I)=COST(I)*(N+1) IF (ABS(COST(I)).GT.MAXCOST) MAXCOST=ABS(COST(I)) 5 CONTINUE C READ SUPPLY OF EACH NODE DO 8 I=1,N READ(13,1000) SURPLUS(I) 8 CONTINUE ENDFILE(13) REWIND(13) PRINT*,'END OF READING' 1000 FORMAT(1I8) 1010 FORMAT(2I8) 1020 FORMAT(4I8) C LARGE=1500000000 C INITIALIZE PRICES & PUSH LISTS DO 10 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 P(NODE)=-INT(LARGE/10) 10 CONTINUE C IN THE FOLLOWING STATEMENTS, THE E-SCALING PARAMETERS ARE SET. C THE FOLLOWING RANGES ARE RECOMMENDED: C STARTING EPSILON: BETWEEN MAXCOST/20 TO MAXCOST/2 C SCALING FACTOR: BETWEEN FOUR AND TEN 25 CONTINUE PRINT*,'ENTER STARTING EPSILON' PRINT*,'SHOULD BE BETWEEN 1 & ',MAXCOST READ*,EPS IF (EPS.LT.1) GO TO 25 PRINT*,'ENTER SCALING FACTOR' 30 CONTINUE READ*,FACTOR IF ((EPS.GT.1).AND.(FACTOR.LE.1)) THEN PRINT*,'ENTER SCALING FACTOR; SHOULD BE GREATER THAN 1' GO TO 30 END IF C MAXSURPLUS=-LARGE DO 920 NODE=1,N IF (SURPLUS(NODE).GT.MAXSURPLUS) THEN MAXSURPLUS=SURPLUS(NODE) END IF 920 CONTINUE C SET THE PARAMETER SSCALE TO 1 TO ALLOW SURPLUS SCALING C A HEURISTIC SCHEME IS USED TO SET THE SFACTOR AND C THRESHSURPLUS PARAMETERS THAT CONTROL SURPLUS SCALING SSCALE=1 IF (SSCALE.EQ.1) THEN SFACTOR=2+INT(FACTOR*MAXSURPLUS/MAXCOST) IF (SFACTOR.LT.4) SFACTOR=4 THRESHSURPLUS=INT(MAXSURPLUS/SFACTOR) IF (EPS.EQ.1) THRESHSURPLUS=0 ELSE THRESHSURPLUS=0 END IF PRINT*,'INITIALIZING DATA STRUCTURES' CALL INIDAT PRINT *,'**********************************' PRINT *,'CALLING E-RELAX FOR MCF (+ SURPLUS NODES ITERATED ONLY)' TIMER = LONG(362) CALL EPS_RELAX TT =(LONG(362) - TIMER)/60 PRINT*, 'TOTAL TIME = ',TT, ' secs.' TCOST=0 DO 330 I=1,NA COST(I)=COST(I)/(N+1) PRODUCT=FLOAT(F(I)*COST(I)) TCOST=TCOST+PRODUCT 330 CONTINUE WRITE(9,1100) TCOST 1100 FORMAT(' ','OPTIMAL COST =',F14.2) PRINT *,'**********************************' PRINT *,'# OF ITERATIONS = ',TOTALITER PRINT *,'# OF S. N. PRICE RISES = ',NITER PRINT *,'**********************************' C CHECK CORRECTNESS OF THE ANSWER IF (EPS.NE.1) THEN PRINT*,'* CAUTION * THE FINAL EPSILON IS EQUAL TO ',EPS END IF DO 80 NODE=1,N IF (SURPLUS(NODE).NE.0) THEN PRINT*,'NONZERO SURPLUS AT NODE ',NODE ENDIF 80 CONTINUE DO 90 ARC=1,NA COST(ARC)=COST(ARC)*(N+1) IF (F(ARC).GT.0) THEN IF (P(STARTN(ARC))-P(ENDN(ARC)).LT.-EPS+COST(ARC)) THEN PRINT*,'E-CS VIOLATED AT ARC ',ARC ENDIF ENDIF IF (F(ARC).LT.U(ARC)) THEN IF (P(STARTN(ARC))-P(ENDN(ARC)).GT.EPS+COST(ARC)) THEN PRINT*,'E-CS VIOLATED AT ARC ',ARC ENDIF ENDIF 90 CONTINUE PRINT *, ' ' PRINT *, 'PROGRAM ENDED; PRESS ' PAUSE STOP END C **************************************************** SUBROUTINE INIDAT C THIS SUBROUTINE USES THE DATA ARRAYS STARTN AND ENDN C TO CONSTRUCT AUXILIARY DATA ARRAYS FOUT, NXTOU, FIN, AND C NXTIN THAT ARE REQUIRED BY E-RELAX. IN THIS SUBROUTINE WE C ARBITRARILY ORDER THE ARCS LEAVING EACH NODE AND STORE C THIS INFORMATION IN FOUT AND NXTOU. SIMILARLY, WE ARBITRA- C RILLY ORDER THE ARCS ENTERING EACH NODE AND STORE THIS C INFORMATION IN FIN AND NXTIN. AT THE COMPLETION OF THE C CONSTRUCTION, WE HAVE THAT C C FOUT(I) = FIRST ARC LEAVING NODE I. C NXTOU(J) = NEXT ARC LEAVING THE HEAD NODE OF ARC J. C FIN(I) = FIRST ARC ENTERING NODE I. C NXTIN(J) = NEXT ARC ENTERING THE TAIL NODE OF ARC J. PARAMETER (MAXNODES=8000, MAXARCS=25000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,FACTOR,EPS COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER FINALIN(MAXNODES),FINALOU(MAXNODES) DO 20 NODE=1,N FIN(NODE)=0 FOUT(NODE)=0 FINALIN(NODE)=0 FINALOU(NODE)=0 20 CONTINUE DO 30 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) IF (FOUT(START).NE.0) THEN NXTOU(FINALOU(START))=ARC ELSE FOUT(START)=ARC END IF IF (FIN(END).NE.0) THEN NXTIN(FINALIN(END))=ARC ELSE FIN(END)=ARC END IF FINALOU(START)=ARC FINALIN(END)=ARC NXTIN(ARC)=0 NXTOU(ARC)=0 30 CONTINUE RETURN END C **************************************************************** C C SCALED E-RELAXATION C THIS VERSION USES A QUEUE TO SELECT NODES. C NODES JOIN THE QUEUE AT THE BOTTOM. C C **************************************************************** SUBROUTINE EPS_RELAX IMPLICIT NONE INTEGER N,NA,LARGE,FACTOR,EPS,ARC,ARCF,ARCB INTEGER NCYC,NITER,NODE,PRICE,PRICEINCR INTEGER I,J,K,LN,CAPOUT,CAPIN,PREVNODE,LASTQUEUE INTEGER SURPI,REJI,NEXT,START,END,PSTART,PEND,MAXPRICE INTEGER RESID,DEL,FLOW,NXTNODE,TOTALITER INTEGER REFPRICE,PRJ,XP,REFPRJ,PRICEJ,THRESHSURPLUS,SFACTOR INTEGER NXTQUEUE(1) INTEGER STARTN(1),ENDN(1),U(1),F(1),SURPLUS(1),FIN(1),FOUT(1) INTEGER NXTIN(1),NXTOU(1),P(1) INTEGER FPUSHF(1),NXTPUSHF(1),FPUSHB(1),NXTPUSHB(1) INTEGER COST(1) COMMON /SCALARS/ N,NA,LARGE,FACTOR,EPS, &THRESHSURPLUS,SFACTOR COMMON /STATS/ NCYC,NITER,TOTALITER COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /PRICES/ P COMMON /BLK3/ SURPLUS COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /BLK11/ FPUSHF COMMON /BLK12/ NXTPUSHF COMMON /BLK13/ FPUSHB COMMON /BLK14/ NXTPUSHB COMMON /COST/ COST COMMON /QUEUE/ NXTQUEUE C INITIALIZE COUNT OF NUMBER OF UP ITERATIONS PERFORMED NITER = 0 TOTALITER=0 NCYC=0 MAXPRICE =LARGE C REDUCE ARC CAPACITIES DO 40 NODE=1,N CAPOUT=0 ARC=FOUT(NODE) 41 IF (ARC.GT.0) THEN CAPOUT=MIN(LARGE,CAPOUT+U(ARC)) ARC=NXTOU(ARC) GO TO 41 END IF CAPOUT=MIN(LARGE,CAPOUT-SURPLUS(NODE)) IF (CAPOUT.LT.0) GOTO 400 CAPIN=0 ARC=FIN(NODE) 43 IF (ARC.GT.0) THEN IF (U(ARC) .GT. CAPOUT) THEN U(ARC)=CAPOUT ENDIF CAPIN=MIN(LARGE,CAPIN+U(ARC)) ARC=NXTIN(ARC) GO TO 43 END IF CAPIN=MIN(LARGE,CAPIN+SURPLUS(NODE)) IF (CAPIN.LT.0) GOTO 400 ARC=FOUT(NODE) 45 IF (ARC.GT.0) THEN IF (U(ARC) .GT. CAPIN) THEN U(ARC)=CAPIN ENDIF ARC=NXTOU(ARC) GO TO 45 END IF 40 CONTINUE C SET ARC FLOWS TO SATISFY E-CS DO 49 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) PSTART=P(START) PEND=P(END) IF (PSTART.GE.PEND+COST(ARC)+EPS) THEN SURPLUS(START)=SURPLUS(START)-U(ARC) SURPLUS(END)=SURPLUS(END)+U(ARC) F(ARC)=U(ARC) ELSE F(ARC)=0 END IF 49 CONTINUE C *********** START OF A NEW SCALING PHASE *********** 60 CONTINUE C ***** QUEUE INITIALIZATION ***** DO 82 NODE=1,N-1 NXTQUEUE(NODE)=NODE+1 82 CONTINUE NXTQUEUE(N)=1 I=1 PREVNODE=N LASTQUEUE=N C ************ START A NEW CYCLE OF UP ITERATIONS ************ 100 CONTINUE C TAKE UP NEXT NODE (NODE I) FOR ITERATION SURPI=SURPLUS(I) IF (SURPI.GT.THRESHSURPLUS) THEN C ARCF & ARCB ARE THE CURRENT VALUES OF THE STARTING ARCS OF THE C PUSH LISTS OF I TOTALITER=TOTALITER+1 ARCF = FPUSHF(I) ARCB = FPUSHB(I) PRICE = P(I) C PRINT*,'NODE ',I,' START PRICE = ',PRICE 115 IF ((ARCF .GT. 0) .OR. (ARCB .GT. 0)) THEN C START BY TRYING TO PUSH AWAY FLOW ON ARCS THAT WERE C ADMISSIBLE AT THE END OF THE LAST ITERATION (IF ANY) C AT THIS NODE. WE MUST CHECK THAT THEY ARE STILL C ADMISSIBLE. 120 IF ((SURPI .GT. 0) .AND. (ARCF .GT. 0)) THEN J = ENDN(ARCF) IF (PRICE-P(J)-COST(ARCF).EQ.EPS) THEN RESID = U(ARCF) - F(ARCF) IF (RESID.GT.0) THEN IF (SURPI .GE. RESID) THEN F(ARCF) = U(ARCF) SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + RESID ARCF = NXTPUSHF(ARCF) ELSE F(ARCF) = F(ARCF) + SURPI SURPLUS(J) = SURPLUS(J) + SURPI SURPI = 0 ENDIF IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCF = NXTPUSHF(ARCF) ENDIF ELSE ARCF = NXTPUSHF(ARCF) ENDIF GOTO 120 ENDIF 121 IF ((SURPI .GT. 0) .AND. (ARCB .GT. 0)) THEN J = STARTN(ARCB) IF (PRICE-P(J)+COST(ARCB).EQ.EPS) THEN RESID=F(ARCB) IF (RESID.GT.0) THEN IF (SURPI .GE. RESID) THEN SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + RESID F(ARCB) = 0 ARCB = NXTPUSHB(ARCB) ELSE F(ARCB) = F(ARCB) - SURPI SURPLUS(J) = SURPLUS(J) + SURPI SURPI = 0 ENDIF IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCB = NXTPUSHB(ARCB) ENDIF ELSE ARCB = NXTPUSHB(ARCB) ENDIF GOTO 121 ENDIF ENDIF C ***** END OF D-PUSHES; CHECK IF PRICE RISE IS NEEDED ***** IF ((ARCF .EQ. 0) .AND. (ARCB .EQ. 0)) THEN C ******************DO A PRICE RISE ********************* REFPRICE=PRICE PRICE = LARGE NITER=NITER+1 ARC = FOUT(I) 111 IF (ARC .GT. 0) THEN IF (F(ARC) .LT. U(ARC)) THEN XP = P(ENDN(ARC))+COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ARCF = ARC NXTPUSHF(ARC)=0 ELSE IF (XP.EQ.PRICE) THEN NXTPUSHF(ARC) = ARCF ARCF = ARC END IF ENDIF ENDIF ARC = NXTOU(ARC) GOTO 111 ENDIF ARC = FIN(I) 112 IF (ARC .GT. 0) THEN IF (F(ARC) .GT. 0) THEN XP = P(STARTN(ARC))-COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ARCB = ARC ARCF=0 NXTPUSHB(ARC)=0 ELSE IF (XP.EQ.PRICE) THEN NXTPUSHB(ARC) = ARCB ARCB = ARC END IF ENDIF ENDIF ARC = NXTIN(ARC) GOTO 112 ENDIF C IF PRICE=LARGE, THE PUSH LIST IS LIKELY EMPTY, SO SET PRICE TO AT LEAST THE C LEVEL OF THE BEST ARC IF (PRICE.EQ.LARGE) THEN ARCF=0 ARCB=0 ARC = FOUT(I) 211 IF (ARC.GT.0) THEN XP = P(ENDN(ARC))+COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ELSE IF (XP.GE.LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' PAUSE STOP END IF ENDIF ARC = NXTOU(ARC) GOTO 211 ENDIF ARC = FIN(I) 212 IF (ARC.GT.0) THEN XP = P(STARTN(ARC))-COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ELSE IF (XP.GE.LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' PAUSE STOP END IF ENDIF ARC = NXTIN(ARC) GOTO 212 ENDIF IF (SURPI.GT.0) GOTO 400 IF (PRICE.LT.INT(LARGE/10)) PRICE=INT(LARGE/10) END IF C SET PRICE OF NODE I PRICE=PRICE+EPS P(I) = PRICE FPUSHF(I) = ARCF FPUSHB(I) = ARCB X PRICEINCR=PRICE-REFPRICE X IF (PRICEINCR.LE.0) THEN X PRINT*,'NONPOSITIVE PRICE INCREMENT=',PRICEINCR X PAUSE X STOP X ENDIF C END OF PRICE RISE; IF THERE IS STILL A LOT OF SURPLUS, C TRY AND DO SOME PUSHING. IF (SURPI.GT.THRESHSURPLUS) GOTO 115 ELSE C IF NO PRICE RISE TOOK PLACE RESET THE START OF THE PUSH LISTS FPUSHF(I) = ARCF FPUSHB(I) = ARCB ENDIF C DO THE FINAL BOOKKEEPING OF THE ITERATION SURPLUS(I) = SURPI C IF THE PRICE IS TOO HIGH, THEN SOMETHING IS WRONG X IF (PRICE.GT.MAXPRICE) THEN X GO TO 400 X ENDIF ENDIF C ********** END OF UP ITERATION STARTING AT NODE I ********* C CHECK FOR THE END OF THE QUEUE. THIS STEP IS NOT NECESSARY FOR C THE ALGORITHM; IT IS INCLUDED FOR COLLECTING STATISTICS ABOUT C THE ALGORITHM'S BEHAVIOR, E.G. COUNTING THE NUMBER OF CYCLES X IF (I.EQ.LASTQUEUE) THEN X LASTQUEUE=PREVNODE X NCYC=NCYC+1 X PRINT*,'CYCLE # ',NCYC,' # OF PRICE RISES ',NITER X END IF C CHECK FOR TERMINATION OF SCALING PHASE. IF SCALING PHASE IS C NOT FINISHED, ADVANCE THE QUEUE AND RETURN TO TAKE ANOTHER NODE. NXTNODE=NXTQUEUE(I) IF (I.NE.NXTNODE) THEN NXTQUEUE(I)=0 NXTQUEUE(PREVNODE)=NXTNODE I=NXTNODE GO TO 100 END IF C ************** END OF SUBPROBLEM (SCALING PHASE) ************** X PRINT*,'END OF SCALING PHASE' C DO A DIAGNOSTIC CHECK X LN=0 X DO 500 NODE=1,N X IF (SURPLUS(NODE).NE.0) THEN X LN=LN+1 X ENDIF 500 CONTINUE X PRINT*,'NONZERO SURPLUS AT ',LN,' NODES ' X DO 600 ARC=1,NA X IF (F(ARC).GT.0) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).LT.-EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF X IF (F(ARC).LT.U(ARC)) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).GT.EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF 600 CONTINUE X PRINT*,'END OF CHECK OF OLD PHASE' C ****** IF EPSILON IS 1 TERMINATE; ELSE REDUCE EPSILON****** IF (EPS.EQ.1) THEN RETURN ELSE EPS=INT(EPS/FACTOR) IF (EPS.LT.1) EPS=1 END IF THRESHSURPLUS=INT(THRESHSURPLUS/SFACTOR) IF (EPS.EQ.1) THRESHSURPLUS=0 C RESET THE FLOWS & THE PUSH LISTS; FIND THE MINIMAL PRICE XP=LARGE DO 800 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 IF (P(NODE).LT.XP) XP=P(NODE) 800 CONTINUE C REDUCE ALL PRICES TO REDUCE DANGER OF OVERFLOW DEL=XP+INT(LARGE/10) IF (DEL.LT.0) DEL=0 DO 810 NODE=1,N P(NODE)=P(NODE)-DEL 810 CONTINUE X PRINT*,'MODIFYING ARC FLOWS TO SATISFY E-CS' DO 900 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) PSTART=P(START) PEND=P(END) IF (PSTART.GT.PEND+EPS+COST(ARC)) THEN RESID=U(ARC)-F(ARC) IF (RESID.GT.0) THEN SURPLUS(START)=SURPLUS(START)-RESID SURPLUS(END)=SURPLUS(END)+RESID F(ARC)=U(ARC) END IF ELSE IF (PSTART.LT.PEND-EPS+COST(ARC)) THEN FLOW=F(ARC) IF (FLOW.GT.0) THEN SURPLUS(START)=SURPLUS(START)+FLOW SURPLUS(END)=SURPLUS(END)-FLOW F(ARC)=0 END IF END IF END IF 900 CONTINUE C RETURN FOR ANOTHER PHASE X PRINT*,'CHECKING E-CS BEFORE STARTING NEW PHASE' X DO 700 ARC=1,NA X IF (F(ARC).GT.0) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).LT.-EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF X IF (F(ARC).LT.U(ARC)) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).GT.EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF 700 CONTINUE GO TO 60 400 PRINT*,'PROBLEM IS INFEASIBLE' PAUSE STOP END $\e$-RELAX-N This code is the same as $\e$-RELAX but allows relaxation iterations starting from nodes with negative as well as positive surplus. C ** SAMPLE CALLING PROGRAM FOR E-RELAX ** C C THIS PROGRAM WILL READ A PROBLEM FILE IN STANDARD FORMAT C AND SOLVE IT USING E-RELAX. C C THIS IS AN E-SCALED VERSION THAT USES A CYCLIC QUEUE C C ITERATIONS ARE DONE FROM BOTH POSITIVE AND NEGATIVE C SURPLUS NODES. AN ADAPTIVE MECHANISM IS USED TO DECIDE C WHETHER TO ITERATE FROM A NEGATIVE SURPLUS NODE C C *************************************************** C PARAMETER (MAXNODES=8000, MAXARCS=25000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,FACTOR,EPS, &THRESHSURPLUS,SFACTOR COMMON /STATS/ NCYC,NITER,TOTALITER COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /PRICES/ P COMMON /BLK3/ SURPLUS COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /BLK11/ FPUSHF COMMON /BLK12/ NXTPUSHF COMMON /BLK13/ FPUSHB COMMON /BLK14/ NXTPUSHB COMMON /COST/ COST COMMON /QUEUE/ NXTQUEUE COMMON /STATUS/ STATUS INTEGER NXTQUEUE(MAXNODES) INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER U(MAXARCS) INTEGER F(MAXARCS) INTEGER SURPLUS(MAXNODES) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER P(MAXNODES) INTEGER FPUSHF(MAXNODES) INTEGER NXTPUSHF(MAXARCS) INTEGER FPUSHB(MAXNODES) INTEGER NXTPUSHB(MAXARCS) INTEGER COST(MAXARCS) INTEGER STATUS(MAXNODES) REAL*8 TCOST,TT,TIMER PRINT*,'E-RELAXATION METHOD FOR MIN COST FLOW' PRINT*,'*************************************' PRINT *,'READING PROBLEM DATA' OPEN(13,FILE='FOR013.DAT',STATUS='OLD') REWIND(13) C READ NUMBER OF NODES AND ARCS READ(13,1010) N,NA C READ START, END, COST, AND CAPACITY OF EACH ARC C AND GENERATE MAX COST VALUE MAXCOST=0 DO 5 I=1,NA READ(13,1020) STARTN(I),ENDN(I),COST(I),U(I) COST(I)=COST(I)*(N+1) IF (ABS(COST(I)).GT.MAXCOST) MAXCOST=ABS(COST(I)) 5 CONTINUE C READ SUPPLY OF EACH NODE DO 8 I=1,N READ(13,1000) SURPLUS(I) 8 CONTINUE ENDFILE(13) REWIND(13) PRINT*,'END OF READING' 1000 FORMAT(1I8) 1010 FORMAT(2I8) 1020 FORMAT(4I8) C LARGE=1500000000 C INITIALIZE NODE PRICES, SURPLUSES, ETC DO 10 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 P(NODE)=0 STATUS(NODE)=0 10 CONTINUE C IN THE FOLLOWING STATEMENTS, THE E-SCALING PARAMETERS ARE SET. C THE FOLLOWING RANGES ARE RECOMMENDED: C STARTING EPSILON: BETWEEN MAXCOST/20 TO MAXCOST/2 C SCALING FACTOR: BETWEEN FOUR AND TEN 25 CONTINUE PRINT*,'ENTER STARTING EPSILON' PRINT*,'SHOULD BE BETWEEN 1 & ',MAXCOST READ*,EPS IF (EPS.LT.1) GO TO 25 PRINT*,'ENTER SCALING FACTOR' 30 CONTINUE READ*,FACTOR IF ((EPS.GT.1).AND.(FACTOR.LE.1)) THEN PRINT*,'ENTER SCALING FACTOR; SHOULD BE GREATER THAN 1' GO TO 30 END IF C MAXSURPLUS=-LARGE DO 920 NODE=1,N IF (SURPLUS(NODE).GT.MAXSURPLUS) THEN MAXSURPLUS=SURPLUS(NODE) END IF 920 CONTINUE C SET THE PARAMETER SSCALE TO 1 TO ALLOW SURPLUS SCALING C A HEURISTIC SCHEME IS USED TO SET THE SFACTOR AND C THRESHSURPLUS PARAMETERS THAT CONTROL SURPLUS SCALING SSCALE=1 IF (SSCALE.EQ.1) THEN SFACTOR=2+INT(FACTOR*MAXSURPLUS/MAXCOST) IF (SFACTOR.LT.4) SFACTOR=4 THRESHSURPLUS=INT(MAXSURPLUS/SFACTOR) IF (EPS.EQ.1) THRESHSURPLUS=0 ELSE THRESHSURPLUS=0 END IF PRINT*,'INITIALIZING DATA STRUCTURES' CALL INIDAT PRINT *,'**********************************' PRINT *,'CALLING E-RELAX FOR MCF (+ AND - SURPLUS NODES ITERATED)' TIMER = LONG(362) CALL EPS_RELAX_N TT =(LONG(362) - TIMER)/60 PRINT*, 'TOTAL TIME = ',TT, ' secs.' TCOST=0 DO 330 I=1,NA 330 TCOST=TCOST+F(I)*COST(I)/(N+1) WRITE(9,1100) TCOST 1100 FORMAT(' ','OPTIMAL COST =',F14.2) PRINT *,'**********************************' PRINT *,'# OF ITERATIONS = ',TOTALITER PRINT *,'# OF S. N. PRICE RISES = ',NITER PRINT *,'**********************************' C CHECK CORRECTNESS OF THE ANSWER IF (EPS.NE.1) THEN PRINT*,'* CAUTION * THE FINAL EPSILON IS EQUAL TO ',EPS END IF DO 80 NODE=1,N IF (SURPLUS(NODE).NE.0) THEN PRINT*,'NONZERO SURPLUS AT NODE ',NODE ENDIF 80 CONTINUE DO 90 ARC=1,NA IF (F(ARC).GT.0) THEN IF (P(STARTN(ARC))-P(ENDN(ARC)).LT.-EPS+COST(ARC)) THEN PRINT*,'E-CS VIOLATED AT ARC ',ARC ENDIF ENDIF IF (F(ARC).LT.U(ARC)) THEN IF (P(STARTN(ARC))-P(ENDN(ARC)).GT.EPS+COST(ARC)) THEN PRINT*,'E-CS VIOLATED AT ARC ',ARC ENDIF ENDIF 90 CONTINUE PRINT *, ' ' PRINT *, 'PROGRAM ENDED; PRESS ' PAUSE STOP END C **************************************************** SUBROUTINE INIDAT C THIS SUBROUTINE USES THE DATA ARRAYS STARTN AND ENDN C TO CONSTRUCT AUXILIARY DATA ARRAYS FOUT, NXTOU, FIN, AND C NXTIN THAT ARE REQUIRED BY E-RELAX-N. IN THIS SUBROUTINE WE C ARBITRARILY ORDER THE ARCS LEAVING EACH NODE AND STORE C THIS INFORMATION IN FOUT AND NXTOU. SIMILARLY, WE ARBITRA- C RILLY ORDER THE ARCS ENTERING EACH NODE AND STORE THIS C INFORMATION IN FIN AND NXTIN. AT THE COMPLETION OF THE C CONSTRUCTION, WE HAVE THAT C C FOUT(I) = FIRST ARC LEAVING NODE I. C NXTOU(J) = NEXT ARC LEAVING THE HEAD NODE OF ARC J. C FIN(I) = FIRST ARC ENTERING NODE I. C NXTIN(J) = NEXT ARC ENTERING THE TAIL NODE OF ARC J. PARAMETER (MAXNODES=8000, MAXARCS=25000) IMPLICIT INTEGER (A-Z) COMMON /SCALARS/ N,NA,LARGE,FACTOR,EPS COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU INTEGER STARTN(MAXARCS) INTEGER ENDN(MAXARCS) INTEGER FIN(MAXNODES) INTEGER FOUT(MAXNODES) INTEGER NXTIN(MAXARCS) INTEGER NXTOU(MAXARCS) INTEGER FINALIN(MAXNODES),FINALOU(MAXNODES) DO 20 NODE=1,N FIN(NODE)=0 FOUT(NODE)=0 FINALIN(NODE)=0 FINALOU(NODE)=0 20 CONTINUE DO 30 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) IF (FOUT(START).NE.0) THEN NXTOU(FINALOU(START))=ARC ELSE FOUT(START)=ARC END IF IF (FIN(END).NE.0) THEN NXTIN(FINALIN(END))=ARC ELSE FIN(END)=ARC END IF FINALOU(START)=ARC FINALIN(END)=ARC NXTIN(ARC)=0 NXTOU(ARC)=0 30 CONTINUE RETURN END C **************************************************************** C C SCALED E-RELAXATION C THIS VERSION USES A QUEUE TO SELECT NODES. C NODES JOIN THE QUEUE AT THE BOTTOM. C USES BOTH POSITIVE AND NEGATIVE SURPLUS ITERATIONS. C C **************************************************************** SUBROUTINE EPS_RELAX_N IMPLICIT NONE INTEGER N,NA,LARGE,FACTOR,EPS,ARC,ARCF,ARCB INTEGER NCYC,NITER,NODE,PRICE,PRICEINCR INTEGER I,J,K,LN,CAPOUT,CAPIN,PREVNODE,LASTQUEUE INTEGER SURPI,REJI,NEXT,START,END,PSTART,PEND,MAXPRICE INTEGER RESID,DEL,FLOW,NXTNODE,TOTALITER,PARSTAT INTEGER REFPRICE,PRJ,XP,REFPRJ,PRICEJ,THRESHSURPLUS,SFACTOR INTEGER NXTQUEUE(1) INTEGER STARTN(1),ENDN(1),U(1),F(1),SURPLUS(1),FIN(1),FOUT(1) INTEGER NXTIN(1),NXTOU(1),P(1) INTEGER FPUSHF(1),NXTPUSHF(1),FPUSHB(1),NXTPUSHB(1) INTEGER COST(1),STATUS(1) COMMON /SCALARS/ N,NA,LARGE,FACTOR,EPS, &THRESHSURPLUS,SFACTOR COMMON /STATS/ NCYC,NITER,TOTALITER COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /UBOUND/ U COMMON /FLOW/ F COMMON /PRICES/ P COMMON /BLK3/ SURPLUS COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU COMMON /BLK11/ FPUSHF COMMON /BLK12/ NXTPUSHF COMMON /BLK13/ FPUSHB COMMON /BLK14/ NXTPUSHB COMMON /COST/ COST COMMON /QUEUE/ NXTQUEUE COMMON /STATUS/ STATUS C C INITIALIZE COUNT OF NUMBER OF UP ITERATIONS PERFORMED C NITER = 0 TOTALITER=0 NCYC=0 MAXPRICE =LARGE C C REDUCE ARC CAPACITIES C DO 40 NODE=1,N CAPOUT=0 ARC=FOUT(NODE) 41 IF (ARC.GT.0) THEN CAPOUT=MIN(LARGE,CAPOUT+U(ARC)) ARC=NXTOU(ARC) GO TO 41 END IF CAPOUT=MIN(LARGE,CAPOUT-SURPLUS(NODE)) IF (CAPOUT.LT.0) GOTO 400 CAPIN=0 ARC=FIN(NODE) 43 IF (ARC.GT.0) THEN IF (U(ARC) .GT. CAPOUT) THEN U(ARC)=CAPOUT ENDIF CAPIN=MIN(LARGE,CAPIN+U(ARC)) ARC=NXTIN(ARC) GO TO 43 END IF CAPIN=MIN(LARGE,CAPIN+SURPLUS(NODE)) IF (CAPIN.LT.0) GOTO 400 ARC=FOUT(NODE) 45 IF (ARC.GT.0) THEN IF (U(ARC) .GT. CAPIN) THEN U(ARC)=CAPIN ENDIF ARC=NXTOU(ARC) GO TO 45 END IF 40 CONTINUE C C SET ARC FLOWS TO SATISFY E-CS C DO 49 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) PSTART=P(START) PEND=P(END) IF (PSTART.GE.PEND+COST(ARC)+EPS) THEN SURPLUS(START)=SURPLUS(START)-U(ARC) SURPLUS(END)=SURPLUS(END)+U(ARC) F(ARC)=U(ARC) ELSE F(ARC)=0 END IF 49 CONTINUE C *********** START OF A NEW SCALING PHASE *********** 60 CONTINUE C ***** QUEUE INITIALIZATION ***** DO 82 NODE=1,N-1 NXTQUEUE(NODE)=NODE+1 82 CONTINUE NXTQUEUE(N)=1 I=1 PREVNODE=N LASTQUEUE=N C SET THE STATUS PARAMETER FOR ALLOWING A NEG. SURPLUS ITERATION PARSTAT=0 C ************ START A NEW CYCLE OF UP ITERATIONS ************ 100 CONTINUE C TAKE UP NEXT NODE (NODE I) FOR ITERATION SURPI=SURPLUS(I) IF (SURPI.GT.THRESHSURPLUS) THEN C ARCF & ARCB ARE THE CURRENT VALUES OF THE STARTING ARCS OF THE C PUSH LISTS OF I TOTALITER=TOTALITER+1 ARCF = FPUSHF(I) ARCB = FPUSHB(I) STATUS(I)=STATUS(I)+1 PRICE = P(I) 115 IF ((ARCF .GT. 0) .OR. (ARCB .GT. 0)) THEN C START BY TRYING TO PUSH AWAY FLOW ON ARCS THAT WERE C ADMISSIBLE AT THE END OF THE LAST ITERATION (IF ANY) C AT THIS NODE. WE MUST CHECK THAT THEY ARE STILL C ADMISSIBLE. 120 IF ((SURPI .GT. 0) .AND. (ARCF .GT. 0)) THEN J = ENDN(ARCF) IF (PRICE-P(J)-COST(ARCF).EQ.EPS) THEN RESID = U(ARCF) - F(ARCF) IF (RESID.GT.0) THEN IF (SURPLUS(J).LT.0) PARSTAT=PARSTAT+1 IF (SURPI .GE. RESID) THEN F(ARCF) = U(ARCF) SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + RESID ARCF = NXTPUSHF(ARCF) ELSE F(ARCF) = F(ARCF) + SURPI SURPLUS(J) = SURPLUS(J) + SURPI SURPI = 0 ENDIF IF (ABS(SURPLUS(J)).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCF = NXTPUSHF(ARCF) ENDIF ELSE ARCF = NXTPUSHF(ARCF) ENDIF GOTO 120 ENDIF 121 IF ((SURPI .GT. 0) .AND. (ARCB .GT. 0)) THEN J = STARTN(ARCB) IF (PRICE-P(J)+COST(ARCB).EQ.EPS) THEN RESID=F(ARCB) IF (RESID.GT.0) THEN IF (SURPLUS(J).LT.0) PARSTAT=PARSTAT+1 IF (SURPI .GE. RESID) THEN SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + RESID F(ARCB) = 0 ARCB = NXTPUSHB(ARCB) ELSE F(ARCB) = F(ARCB) - SURPI SURPLUS(J) = SURPLUS(J) + SURPI SURPI = 0 ENDIF IF (ABS(SURPLUS(J)).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCB = NXTPUSHB(ARCB) ENDIF ELSE ARCB = NXTPUSHB(ARCB) ENDIF GOTO 121 ENDIF ENDIF C ***** END OF D-PUSHES; CHECK IF PRICE RISE IS NEEDED ***** IF ((ARCF .EQ. 0) .AND. (ARCB .EQ. 0)) THEN C ******************DO A PRICE RISE ********************* REFPRICE=PRICE PRICE = LARGE NITER=NITER+1 ARC = FOUT(I) 111 IF (ARC .GT. 0) THEN IF (F(ARC) .LT. U(ARC)) THEN XP = P(ENDN(ARC))+COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ARCF = ARC NXTPUSHF(ARC)=0 ELSE IF (XP.EQ.PRICE) THEN NXTPUSHF(ARC) = ARCF ARCF = ARC END IF ENDIF ENDIF ARC = NXTOU(ARC) GOTO 111 ENDIF ARC = FIN(I) 112 IF (ARC .GT. 0) THEN IF (F(ARC) .GT. 0) THEN XP = P(STARTN(ARC))-COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ARCB = ARC ARCF=0 NXTPUSHB(ARC)=0 ELSE IF (XP.EQ.PRICE) THEN NXTPUSHB(ARC) = ARCB ARCB = ARC END IF ENDIF ENDIF ARC = NXTIN(ARC) GOTO 112 ENDIF C IF PRICE=LARGE, THE PUSH LIST IS LIKELY EMPTY, SO SET PRICE TO AT LEAST THE C LEVEL OF THE BEST ARC IF (PRICE.EQ.LARGE) THEN ARCF=0 ARCB=0 ARC = FOUT(I) 211 IF (ARC.GT.0) THEN XP = P(ENDN(ARC))+COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ELSE PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' PAUSE STOP ENDIF ARC = NXTOU(ARC) GOTO 211 ENDIF ARC = FIN(I) 212 IF (ARC.GT.0) THEN XP = P(STARTN(ARC))-COST(ARC) IF (XP.LT.PRICE) THEN PRICE = XP ELSE PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' PAUSE STOP ENDIF ARC = NXTIN(ARC) GOTO 212 ENDIF IF (SURPI.GT.0) GOTO 400 IF (PRICE.LT.INT(LARGE/10)) PRICE=INT(LARGE/10) END IF C SET PRICE OF NODE I PRICE=PRICE+EPS P(I) = PRICE FPUSHF(I) = ARCF FPUSHB(I) = ARCB C END OF PRICE RISE; IF THERE IS STILL A LOT OF SURPLUS, C TRY AND DO SOME PUSHING. IF (SURPI.GT.THRESHSURPLUS) GOTO 115 ELSE C IF NO PRICE RISE TOOK PLACE RESET THE START OF THE PUSH LISTS FPUSHF(I) = ARCF FPUSHB(I) = ARCB ENDIF C DO THE FINAL BOOKKEEPING OF THE ITERATION SURPLUS(I) = SURPI C IF THE PRICE IS TOO HIGH, THEN SOMETHING IS WRONG X IF (PRICE.GT.MAXPRICE) THEN X GO TO 400 X ENDIF ENDIF C CHECK FOR A NEGATIVE SURPLUS ITERATION IF (SURPI.LT.(-THRESHSURPLUS)) THEN TOTALITER=TOTALITER+1 ARCF = FPUSHF(I) ARCB = FPUSHB(I) PRICE = P(I) 1115 IF ((ARCF .GT. 0) .OR. (ARCB .GT. 0)) THEN C START BY TRYING TO PUSH AWAY FLOW ON ARCS THAT WERE C ADMISSIBLE AT THE END OF THE LAST ITERATION (IF ANY) C AT THIS NODE. WE MUST CHECK THAT THEY ARE STILL C ADMISSIBLE. 1120 IF ((SURPI .LT. 0) .AND. (ARCF .GT. 0)) THEN J = ENDN(ARCF) IF (PRICE-P(J)-COST(ARCF).EQ.(-EPS)) THEN RESID = F(ARCF) IF (RESID.GT.0) THEN IF (SURPLUS(J).GT.0) PARSTAT=PARSTAT+1 IF ((-SURPI) .GE. RESID) THEN F(ARCF) = 0 SURPI = SURPI + RESID SURPLUS(J) = SURPLUS(J) - RESID ARCF = NXTPUSHF(ARCF) ELSE F(ARCF) = F(ARCF) + SURPI SURPLUS(J) = SURPLUS(J) + SURPI SURPI = 0 ENDIF IF (ABS(SURPLUS(J)).GT.THRESHSURPLUS) THEN IF ((NXTQUEUE(J).EQ.0).AND.(STATUS(J).LE.(PARSTAT/2))) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCF = NXTPUSHF(ARCF) ENDIF ELSE ARCF = NXTPUSHF(ARCF) ENDIF GOTO 1120 ENDIF 1121 IF ((SURPI .LT. 0) .AND. (ARCB .GT. 0)) THEN J = STARTN(ARCB) IF (PRICE-P(J)+COST(ARCB).EQ.(-EPS)) THEN RESID=U(ARCB)-F(ARCB) IF (RESID.GT.0) THEN IF (SURPLUS(J).GT.0) PARSTAT=PARSTAT+1 IF ((-SURPI) .GE. RESID) THEN SURPI = SURPI + RESID SURPLUS(J) = SURPLUS(J) - RESID F(ARCB) = U(ARCB) ARCB = NXTPUSHB(ARCB) ELSE F(ARCB) = F(ARCB) - SURPI SURPLUS(J) = SURPLUS(J) + SURPI SURPI = 0 ENDIF IF (ABS(SURPLUS(J)).GT.THRESHSURPLUS) THEN IF ((NXTQUEUE(J).EQ.0).AND.(STATUS(J).LE.(PARSTAT/2))) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCB = NXTPUSHB(ARCB) ENDIF ELSE ARCB = NXTPUSHB(ARCB) ENDIF GOTO 1121 ENDIF ENDIF C ***** END OF D-PUSHES; CHECK IF PRICE RISE IS NEEDED ***** IF ((ARCF .EQ. 0) .AND. (ARCB .EQ. 0)) THEN C ******************DO A PRICE RISE ********************* REFPRICE=PRICE PRICE = -LARGE NITER=NITER+1 ARC = FOUT(I) 1111 IF (ARC .GT. 0) THEN IF (F(ARC) .GT. 0) THEN XP = P(ENDN(ARC))+COST(ARC) IF (XP.GT.PRICE) THEN PRICE = XP ARCF = ARC NXTPUSHF(ARC)=0 ELSE IF (XP.EQ.PRICE) THEN NXTPUSHF(ARC) = ARCF ARCF = ARC END IF ENDIF ENDIF ARC = NXTOU(ARC) GOTO 1111 ENDIF ARC = FIN(I) 1112 IF (ARC .GT. 0) THEN IF (F(ARC) .LT. U(ARC)) THEN XP = P(STARTN(ARC))-COST(ARC) IF (XP.GT.PRICE) THEN PRICE = XP ARCB = ARC ARCF=0 NXTPUSHB(ARC)=0 ELSE IF (XP.EQ.PRICE) THEN NXTPUSHB(ARC) = ARCB ARCB = ARC END IF ENDIF ENDIF ARC = NXTIN(ARC) GOTO 1112 ENDIF C IF PRICE=-LARGE, THE PUSH LIST IS LIKELY EMPTY, SO SET PRICE TO AT MOST THE C LEVEL OF THE BEST ARC IF (PRICE.EQ.-LARGE) THEN ARCF=0 ARCB=0 ARC = FOUT(I) 1211 IF (ARC .GT. 0) THEN XP = P(ENDN(ARC))+COST(ARC) IF (XP.GT.PRICE) THEN PRICE = XP ELSE IF (XP.LE.-LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' PAUSE STOP END IF END IF ARC = NXTOU(ARC) GOTO 1211 END IF ARC = FIN(I) 1212 IF (ARC .GT. 0) THEN XP = P(STARTN(ARC))-COST(ARC) IF (XP.GT.PRICE) THEN PRICE = XP ELSE IF (XP.LE.-LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' PAUSE STOP END IF END IF ARC = NXTIN(ARC) GOTO 1212 ENDIF IF (SURPI.LT.0) GOTO 400 IF (PRICE.GT.-INT(LARGE/10)) PRICE=-INT(LARGE/10) END IF C SET PRICE OF NODE I PRICE=PRICE-EPS P(I) = PRICE FPUSHF(I) = ARCF FPUSHB(I) = ARCB C END OF PRICE RISE; IF THERE IS STILL A LOT OF SURPLUS, C TRY AND DO SOME PUSHING. IF (SURPI.LT.(-THRESHSURPLUS)) GOTO 1115 ENDIF C DO THE FINAL BOOKKEEPING OF THE ITERATION SURPLUS(I) = SURPI END IF C ********** END OF ITERATION STARTING AT NODE I ********* C CHECK FOR THE END OF THE QUEUE. THIS STEP IS NOT NECESSARY FOR C THE ALGORITHM; IT IS INCLUDED FOR COLLECTING STATISTICS ABOUT C THE ALGORITHM'S BEHAVIOR, E.G. COUNTING THE NUMBER OF CYCLES IF (I.EQ.LASTQUEUE) THEN LASTQUEUE=PREVNODE NCYC=NCYC+1 X PRINT*,'CYCLE # ',NCYC,' # OF PRICE RISES ',NITER END IF C CHECK FOR TERMINATION OF SCALING PHASE. IF SCALING PHASE IS C NOT FINISHED, ADVANCE THE QUEUE AND RETURN TO TAKE ANOTHER NODE. NXTNODE=NXTQUEUE(I) IF (I.NE.NXTNODE) THEN NXTQUEUE(I)=0 NXTQUEUE(PREVNODE)=NXTNODE I=NXTNODE GO TO 100 END IF C ************** END OF SUBPROBLEM (SCALING PHASE) ************** X PRINT*,'END OF SCALING PHASE' C DO A DIAGNOSTIC CHECK X LN=0 X DO 500 NODE=1,N X IF (SURPLUS(NODE).NE.0) THEN X LN=LN+1 X ENDIF 500 CONTINUE X PRINT*,'NONZERO SURPLUS AT ',LN,' NODES ' X DO 600 ARC=1,NA X IF (F(ARC).GT.0) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).LT.-EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF X IF (F(ARC).LT.U(ARC)) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).GT.EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF 600 CONTINUE X PRINT*,'END OF CHECK OF OLD PHASE' C ****** IF EPSILON IS 1 TERMINATE; ELSE REDUCE EPSILON****** IF (EPS.EQ.1) THEN RETURN ELSE EPS=INT(EPS/FACTOR) IF (EPS.LT.1) EPS=1 END IF THRESHSURPLUS=INT(THRESHSURPLUS/SFACTOR) IF (EPS.EQ.1) THRESHSURPLUS=0 C RESET THE FLOWS & THE PUSH LISTS DO 800 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 STATUS(NODE)=0 800 CONTINUE DO 900 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) PSTART=P(START) PEND=P(END) IF (PSTART.GT.PEND+EPS+COST(ARC)) THEN RESID=U(ARC)-F(ARC) IF (RESID.GT.0) THEN SURPLUS(START)=SURPLUS(START)-RESID SURPLUS(END)=SURPLUS(END)+RESID F(ARC)=U(ARC) END IF ELSE IF (PSTART.LT.PEND-EPS+COST(ARC)) THEN FLOW=F(ARC) IF (FLOW.GT.0) THEN SURPLUS(START)=SURPLUS(START)+FLOW SURPLUS(END)=SURPLUS(END)-FLOW F(ARC)=0 END IF END IF END IF 900 CONTINUE C RETURN FOR ANOTHER PHASE X PRINT*,'CHECKING E-CS BEFORE STARTING NEW PHASE' X DO 700 ARC=1,NA X IF (F(ARC).GT.0) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).LT.-EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF X IF (F(ARC).LT.U(ARC)) THEN X IF (P(STARTN(ARC))-P(ENDN(ARC)).GT.EPS+COST(ARC)) THEN X PRINT*,'E-CS VIOLATED AT ARC ',ARC X ENDIF X ENDIF 700 CONTINUE GO TO 60 400 PRINT*,'PROBLEM IS INFEASIBLE' PAUSE STOP END ##################################################### End 92001-2.for