C C SAMPLE CALLING PROGRAM FOR RELAX-IV C C--------------------------------------------------------------- C C PURPOSE - THIS PROGRAM READS IN DATA FOR A LINEAR COST C ORDINARY NETWORK FLOW PROBLEM FROM THE FILE `RELAX4.INP', C CALLS THE ROUTINE INIDAT TO CONSTRUCT LINKED LIST FOR THE PROBLEM, C AND THEN CALLS THE ROUTINE RELAX4 TO SOLVE THE PROBLEM. C C--------------------------------------------------------------- C PROGRAM MAIN IMPLICIT INTEGER (A-Z) C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C C INPUT PARAMETERS C C N = NUMBER OF NODES C NA = NUMBER OF ARCS C LARGE = A VERY LARGE INTEGER TO REPRESENT INFINITY C STARTN(J) = STARTING NODE FOR ARC J, J = 1,...,NA C ENDN(J) = ENDING NODE FOR ARC J, J = 1,...,NA C C(J) = COST OF ARC J, J = 1,...,NA C INTEGER STARTN(MAXNA),ENDN(MAXNA),C(MAXNA) COMMON /INPUT/N,NA,LARGE COMMON /ARRAYS/STARTN/ARRAYE/ENDN/ARRAYC/C C C UPDATED PARAMETERS C C U(J) = CAPACITY OF ARC J ON INPUT AND RESIDUAL CAPACITY C ON OUTPUT, J = 1,...,NA C B(I) = DEMAND AT NODE I ON INPUT AND ZERO ON OUTPUT, C I = 1,...,N C INTEGER U(MAXNA),B(MAXNN) COMMON /ARRAYU/U/ARRAYB/B C C OUTPUT PARAMETERS C C X(J) = FLOW ON ARC J, J = 1,...,NA C RC(J) = REDUCED COST OF ARC J, J = 1,...,NA C NMULTINODE = NUMBER OF MULTINODE RELAXATION ITERATIONS IN RELAX4 C ITER = NUMBER OF RELAXATION ITERATIONS IN RELAX4 C NUM_AUGM = NUMBER OF FLOW AUGMENTATION STEPS IN RELAX4 C NUM_ASCNT = NUMBER OF MULTINODE ASCENT STEPS IN RELAX4 C NSP = NUMBER OF AUCTION/SHORTEST PATH ITERATIONS C TCOST = COST OF FLOW C INTEGER X(MAXNA),RC(MAXNA) REAL*8 TCOST COMMON /ARRAYX/X/ARRAYRC/RC COMMON /OUTPUT/NMULTINODE,ITER,NUM_AUGM,NUM_ASCNT,NSP C C WORKING PARAMETERS C INTEGER I1(MAXNN),I2(MAXNN),I3(MAXNN),I4(MAXNA) INTEGER I5(MAXNN),I6(MAXNA),I7(MAXNA) INTEGER TFSTOU(MAXNN),TNXTOU(MAXNA) INTEGER TFSTIN(MAXNN),TNXTIN(MAXNA) INTEGER I14(MAXNN),I15(MAXNN),I16(MAXNN),I17(MAXNN) LOGICAL*1 SCAN(MAXNN),MARK(MAXNN) LOGICAL*1 REPEAT COMMON /BLK1/I1/BLK2/I2/BLK3/I3/BLK4/I4 COMMON /BLK5/I5/BLK6/I6/BLK7/I7 COMMON /BLK8/SCAN/BLK9/MARK COMMON /BLK10/TFSTOU/BLK11/TNXTOU/BLK12/TFSTIN/BLK13/TNXTIN COMMON /BLK14/I14/BLK15/I15/BLK16/I16/BLK17/I17 COMMON /CR/CRASH COMMON /BLKR/REPEAT C C OPTIONAL WORKING PARAMETERS (FOR SENSITIVITY ANALYSIS ONLY) C INTEGER CAP(MAXNA) COMMON /BLKCAP/CAP C C DECLARE TIMING VARIABLES FOR UNIX SYSTEM C REAL*4 TIME0,TIME1,TIME2 COMMON /T/TIME0,TIME1 C C--------------------------------------------------------------- C C READ PROBLEM DATA FROM FILE RELAX4.INP C PRINT*,'READ PROBLEM DATA FROM RELAX4.INP' OPEN(13,FILE='RELAX4.INP',STATUS='OLD') REWIND(13) C C READ NUMBER OF NODES AND ARCS C READ(13,1000) N,NA C C READ START NODE, END NODE, COST, AND CAPACITY OF EACH ARC C DO 20 I=1,NA READ(13,1000) STARTN(I),ENDN(I),C(I),U(I) 20 CONTINUE C C READ SUPPLY OF EACH NODE; CONVERT IT TO DEMAND C DO 30 I=1,N READ(13,1000) B(I) B(I)=-B(I) 30 CONTINUE C 1000 FORMAT(4I8) REWIND(13) CLOSE(13) C PRINT*,'END OF READING' PRINT*,'NUMBER OF NODES =',N,', NUMBER OF ARCS =',NA C C SET LARGE TO A LARGE INTEGER FOR YOUR MACHINE C LARGE=500000000 DANGER_THRESH=LARGE/10 C C CHECK DATA IS WITHIN RANGE OF MACHINE C FLAG1=0 FLAG2=0 FLAG3=0 DO 40 I=1,NA IF (ABS(C(I)).GT.LARGE) FLAG1=1 IF (U(I).GT.LARGE) FLAG2=1 IF (ABS(C(I)).GT.DANGER_THRESH) FLAG3=1 40 CONTINUE IF (FLAG1.EQ.1) THEN PRINT*,'SOME COSTS EXCEED THE ALLOWABLE RANGE' PRINT*,'PROGRAM CANNOT RUN; PRESS TO EXIT' PAUSE STOP END IF IF (FLAG2.EQ.1) THEN PRINT*,'SOME ARC CAPACITIES EXCEED THE ALLOWABLE RANGE' PRINT*,'PROGRAM CANNOT RUN; PRESS TO EXIT' PAUSE STOP END IF IF (FLAG3.EQ.1) THEN PRINT*,'SOME COSTS ARE DANGEROUSLY LARGE' PRINT*,'PROGRAM MAY NOT RUN CORRECTLY' END IF C C--------------------------------------------------------------- C C CONSTRUCT LINKED LISTS FOR THE PROBLEM C PRINT*,'CONSTRUCT LINKED LISTS FOR THE PROBLEM' CALL INIDAT C C--------------------------------------------------------------- C C INITIALIZE DUAL PRICES C (DEFAULT: ALL DUAL PRICES = 0, SO REDUCED COST IS SET C EQUAL TO COST) C IN THE COURSE OF THE ALGORITHM THE REDUCED COST OF AN C ARC IS ADJUSTED TO C RD(ARC)=C(ARC)+IMPLICIT PRICE OF ENDNODE(ARC) - C IMPLICIT PRICE OF STARTNODE(ARC) C THE FOLLOWING COMPLEMENTARY SLACKNESS CONDITION C IS ENFORCED BY RELAX IN ITS INITIALIZATION PHASE AND C IS MAINTAINED THROUGHOUT ITS COURSE C IF RD(ARC).GT.0 THEN FLOW OF ARC MUST BE EQUAL TO 0 C IF RD(ARC).LT.0 THEN FLOW OF ARC MUST BE EQUAL TO ARC CAPACITY C DO 60 I=1,NA 60 RC(I)=C(I) C C NOTE REGARDING REOPTIMIZATION: C THE LOGICAL VARIABLE REPEAT IS SET TO .FALSE. IF THE PROBLEM C IS SOLVED FOR THE FIRST TIME C C REPEAT IS SET TO .TRUE. IF THE PROBLEM IS REOPTIMIZED C FROM A HOT START. IN THIS CASE THE USER MUST ENSURE THAT C THE ARC FLOWS ARE RESET SO THAT C THE COMPLEMENTARY SLACKNESS CONDITION (SEE ABOVE) IS C SATISFIED UPON CALLING RELAX AND C THE DFCT ARRAY PROPERLY CORRESPONDS TO THE RESET ARC FLOWS C C HERE WE SPECIFY THAT WE ARE SOLVING THE PROBLEM FROM SCRATCH C REPEAT=.FALSE. C C STORE CAPACITY OF ARCS IN CAP C (OPTIONAL IF SENSITIVITY ANALYSIS WILL NOT BE ACTIVATED) C DO 70 I=1,NA 70 CAP(I)=U(I) C C--------------------------------------------------------------- C C SET CRASH EQUAL TO 1 TO ACTIVATE AN AUCTION/SHORTEST PATH SUBROUTINE FOR C GETTING THE INITIAL PRICE-FLOW PAIR. THIS IS RECOMMENDED FOR DIFFICULT C PROBLEMS WHERE THE DEFAULT INITIALIZATION YIELDS C LONG SOLUTION TIMES. C PRINT*,'ENTER THE INITIALIZATION DESIRED' PRINT*,' <0> FOR THE DEFAULT INITIALIZATION' PRINT*,' <1> FOR AUCTION INITIALIZATION' READ*,CRASH C C CALL RELAX4 TO SOLVE THE PROBLEM C 75 CONTINUE PRINT*,'CALLING RELAX4 TO SOLVE THE PROBLEM' PRINT*,'***********************************' C C INITIALIZE SYSTEM TIMER C C TIME0 = LONG(362)/60.0 TIME0 = SECNDS(0.0) C CALL RELAX4 C C CALL SYSTEM TIMER TO DISPLAY EXECUTION TIME FOR RELAX4 C C TIME2 = LONG(362)/60.0 - TIME0 TIME2 = SECNDS(TIME0) PRINT*,'TOTAL SOLUTION TIME =',TIME2,' SECS.' PRINT*,'TIME IN INITIALIZATION =',TIME1,' SECS.' C C--------------------------------------------------------------- C C CHECK CORRECTNESS OF OUTPUT PARAMETERS C DO 80 NODE=1,N IF (B(NODE).NE.0) THEN PRINT*,'NONZERO SURPLUS AT NODE ',NODE END IF 80 CONTINUE DO 90 ARC=1,NA IF (X(ARC).GT.0) THEN IF (RC(ARC).GT.0) THEN PRINT*,'COMPLEMENTARY SLACKNESS VIOLATED AT ARC ',ARC ENDIF ENDIF IF (U(ARC).GT.0) THEN IF (RC(ARC).LT.0) THEN PRINT*,'COMPLEMENTARY SLACKNESS VIOLATED AT ARC ',ARC ENDIF ENDIF 90 CONTINUE C C COMPUTE AND DISPLAY COST OF FLOWS (IN DOUBLE PRECISION) C TCOST=FLOAT(0) DO 100 I=1,NA TCOST=TCOST + DFLOAT(X(I)*C(I)) 100 CONTINUE PRINT*,'OPTIMAL COST = ',TCOST C C DISPLAY RELAX4 STATISTICS C IF (CRASH.EQ.1) THEN PRINT*,'NUMBER OF AUCTION/SHORTEST PATH ITERATIONS =',NSP END IF PRINT*,'NUMBER OF ITERATIONS = ',ITER PRINT*,'NUMBER OF MULTINODE ITERATIONS = ',NMULTINODE PRINT*,'NUMBER OF MULTINODE ASCENT STEPS = ',NUM_ASCNT PRINT*,'NUMBER OF REGULAR AUGMENTATIONS = ',NUM_AUGM PRINT*,'***********************************' C C TO ACTIVATE SENSITIVITY ANALYSIS, INSERT THE FOLLOWING C THREE LINES HERE. C CALL SENSTV REPEAT=.TRUE. GO TO 75 C END C C SUBROUTINE INIDAT IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C PURPOSE - THIS ROUTINE CONSTRUCTS TWO LINKED LISTS FOR C THE NETWORK TOPOLOGY: ONE LIST (GIVEN BY FOU, NXTOU) FOR C THE OUTGOING ARCS OF NODES AND ONE LIST (GIVEN BY FIN, C NXTIN) FOR THE INCOMING ARCS OF NODES. THESE TWO LISTS C ARE REQUIRED BY RELAX4. C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C C INPUT PARAMETERS C C N = NUMBER OF NODES C NA = NUMBER OF ARCS C STARTN(J) = STARTING NODE FOR ARC J, J = 1,...,NA C ENDN(J) = ENDING NODE FOR ARC J, J = 1,...,NA C INTEGER STARTN(MAXNA),ENDN(MAXNA) COMMON /ARRAYS/STARTN/ARRAYE/ENDN COMMON /INPUT/N,NA,LARGE C C OUTPUT PARAMETERS C C FOU(I) = FIRST ARC OUT OF NODE I, I = 1,...,N C NXTOU(J) = NEXT ARC OUT OF THE STARTING NODE OF ARC J, C J = 1,...,NA C FIN(I) = FIRST ARC INTO NODE I, I = 1,...,N C NXTIN(J) = NEXT ARC INTO THE ENDING NODE OF ARC J, C J = 1,...,NA C INTEGER FOU(MAXNN),NXTOU(MAXNA),FIN(MAXNN),NXTIN(MAXNA) COMMON /BLK3/FOU/BLK4/NXTOU/BLK5/FIN/BLK6/NXTIN C C WORKING PARAMETERS C INTEGER TEMPIN(MAXNN),TEMPOU(MAXNN) COMMON /BLK1/TEMPIN/BLK2/TEMPOU C C--------------------------------------------------------------- C DO 10 I=1,N FIN(I)=0 FOU(I)=0 TEMPIN(I)=0 10 TEMPOU(I)=0 DO 20 I=1,NA NXTIN(I)=0 NXTOU(I)=0 I1=STARTN(I) I2=ENDN(I) IF (FOU(I1).NE.0) THEN NXTOU(TEMPOU(I1))=I ELSE FOU(I1)=I END IF TEMPOU(I1)=I IF (FIN(I2).NE.0) THEN NXTIN(TEMPIN(I2))=I ELSE FIN(I2)=I END IF 20 TEMPIN(I2)=I RETURN END C C SUBROUTINE RELAX4 IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C RELAX-IV (VERSION OF OCTOBER 1994) C C RELEASE NOTE - THIS VERSION OF RELAXATION CODE HAS OPTION FOR C A SPECIAL CRASH PROCEDURE FOR C THE INITIAL PRICE-FLOW PAIR. THIS IS RECOMMENDED FOR DIFFICULT C PROBLEMS WHERE THE DEFAULT INITIALIZATION C RESULTS IN LONG RUNNING TIMES. C CRASH =1 CORRESPONDS TO AN AUCTION/SHORTEST PATH METHOD C C THESE INITIALIZATIONS ARE RECOMMENDED IN THE ABSENCE OF ANY C PRIOR INFORMATION ON A FAVORABLE INITIAL FLOW-PRICE VECTOR PAIR C THAT SATISFIES COMPLEMENTARY SLACKNESS C C THE RELAXATION PORTION OF THE CODE DIFFERS FROM THE CODE RELAXT-III C AND OTHER EARLIER RELAXATION CODES IN THAT IT MAINTAINS C THE SET OF NODES WITH NONZERO DEFICIT IN A FIFO QUEUE. C LIKE ITS PREDECESSOR RELAXT-III, THIS CODE MAINTAINS A LINKED LIST C OF BALANCED (I.E., OF ZERO REDUCED COST) ARCS SO TO REDUCE C THE WORK IN LABELING AND SCANNING. C UNLIKE RELAXT-III, IT DOES NOT USE SELECTIVELY C SHORTEST PATH ITERATIONS FOR INITIALIZATION. C C--------------------------------------------------------------- C C PURPOSE - THIS ROUTINE IMPLEMENTS THE RELAXATION METHOD C OF BERTSEKAS AND TSENG (SEE [1], [2]) FOR LINEAR C COST ORDINARY NETWORK FLOW PROBLEMS. C C [1] BERTSEKAS, D. P., "A UNIFIED FRAMEWORK FOR PRIMAL-DUAL METHODS ..." C MATHEMATICAL PROGRAMMING, VOL. 32, 1985, PP. 125-145. C [2] BERTSEKAS, D. P., AND TSENG, P., "RELAXATION METHODS FOR C MINIMUM COST ..." OPERATIONS RESEARCH, VOL. 26, 1988, PP. 93-114. C C THE RELAXATION METHOD IS ALSO DESCRIBED IN THE BOOKS: C C [3] BERTSEKAS, D. P., "LINEAR NETWORK OPTIMIZATION: ALGORITHMS AND CODES" C MIT PRESS, 1991. C [4] BERTSEKAS, D. P. AND TSITSIKLIS, J. N., "PARALLEL AND DISTRIBUTED C COMPUTATION: NUMERICAL METHODS", PRENTICE-HALL, 1989. C [5] BERTSEKAS, D. P., "NETWORK OPTIMIZATION: C CONTINUOUS AND DISCRETE MODELS", ATHENA SCIENTIFIC, 1998. C C C C--------------------------------------------------------------- C C SOURCE - THIS CODE WAS WRITTEN BY DIMITRI P. BERTSEKAS C AND PAUL TSENG, WITH A CONTRIBUTION BY JONATHAN ECKSTEIN C IN THE PHASE II INITIALIZATION. THE ROUTINE AUCTION WAS WRITTEN C BY DIMITRI P. BERTSEKAS AND IS BASED ON THE METHOD DESCRIBED IN C THE PAPER: C C [6] BERTSEKAS, D. P., "AN AUCTION/SEQUENTIAL SHORTEST PATH ALGORITHM C FOR THE MINIMUM COST FLOW PROBLEM", LIDS REPORT P-2146, MIT, NOV. 1992. C C FOR INQUIRIES ABOUT THE CODE, PLEASE CONTACT: C C DIMITRI P. BERTSEKAS C LABORATORY FOR INFORMATION AND DECISION SYSTEMS C MASSACHUSETTS INSTITUTE OF TECHNOLOGY C CAMBRIDGE, MA 02139 C (617) 253-7267, DIMITRIB@MIT.EDU C C--------------------------------------------------------------- C C USER GUIDELINES - C C THIS ROUTINE IS IN THE PUBLIC DOMAIN TO BE USED ONLY FOR RESEARCH C PURPOSES. IT CANNOT BE USED AS PART OF A COMMERCIAL PRODUCT, OR C TO SATISFY IN ANY PART COMMERCIAL DELIVERY REQUIREMENTS TO C GOVERNMENT OR INDUSTRY, WITHOUT PRIOR AGREEMENT WITH THE AUTHORS. C USERS ARE REQUESTED TO ACKNOWLEDGE THE AUTHORSHIP OF THE CODE, C AND THE RELAXATION METHOD. C C NO MODIFICATION SHOULD BE MADE TO THIS CODE OTHER C THAN THE MINIMAL NECESSARY C TO MAKE IT COMPATIBLE WITH THE FORTRAN COMPILERS OF SPECIFIC C MACHINES. C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C C INPUT PARAMETERS (SEE NOTES 1, 2, 4) C C N = NUMBER OF NODES C NA = NUMBER OF ARCS C LARGE = A VERY LARGE INTEGER TO REPRESENT INFINITY C (SEE NOTE 3) C REPEAT = .TRUE. IF INITIALIZATION IS TO BE SKIPPED C (.FALSE. OTHERWISE) C CRASH = 0 IF DEFAULT INITIALIZATION IS USED C 1 IF AUCTION INITIALIZATION IS USED C STARTN(J) = STARTING NODE FOR ARC J, J = 1,...,NA C ENDN(J) = ENDING NODE FOR ARC J, J = 1,...,NA C FOU(I) = FIRST ARC OUT OF NODE I, I = 1,...,N C NXTOU(J) = NEXT ARC OUT OF THE STARTING NODE OF ARC J, C J = 1,...,NA C FIN(I) = FIRST ARC INTO NODE I, I = 1,...,N C NXTIN(J) = NEXT ARC INTO THE ENDING NODE OF ARC J, C J = 1,...,NA C INTEGER STARTN(MAXNA),ENDN(MAXNA) INTEGER FOU(MAXNN),NXTOU(MAXNA),FIN(MAXNN),NXTIN(MAXNA) LOGICAL*1 REPEAT COMMON /INPUT/N,NA,LARGE COMMON /ARRAYS/STARTN/ARRAYE/ENDN COMMON /BLK3/FOU/BLK4/NXTOU/BLK5/FIN/BLK6/NXTIN COMMON /BLKR/REPEAT COMMON /CR/CRASH C C UPDATED PARAMETERS (SEE NOTES 1, 3, 4) C C RC(J) = REDUCED COST OF ARC J, J = 1,...,NA C U(J) = CAPACITY OF ARC J ON INPUT C AND (CAPACITY OF ARC J) - X(J) ON OUTPUT, C J = 1,...,NA C DFCT(I) = DEMAND AT NODE I ON INPUT C AND ZERO ON OUTPUT, I = 1,...,N C INTEGER RC(MAXNA),U(MAXNA),DFCT(MAXNN) COMMON /ARRAYRC/RC/ARRAYU/U/ARRAYB/DFCT C C OUTPUT PARAMETERS (SEE NOTES 1, 3, 4) C C X(J) = FLOW ON ARC J, J = 1,...,NA C NMULTINODE = NUMBER OF MULTINODE RELAXATION ITERATIONS IN RELAX4 C ITER = NUMBER OF RELAXATION ITERATIONS IN RELAX4 C NUM_AUGM = NUMBER OF FLOW AUGMENTATION STEPS IN RELAX4 C NUM_ASCNT = NUMBER OF MULTINODE ASCENT STEPS IN RELAX4 C NSP = NUMBER OF AUCTION/SHORTEST PATH ITERATIONS C INTEGER X(MAXNA) COMMON /ARRAYX/X COMMON/OUTPUT/NMULTINODE,ITER,NUM_AUGM,NUM_ASCNT,NSP C C WORKING PARAMETERS (SEE NOTES 1, 4, 5) C INTEGER LABEL(MAXNN),PRDCSR(MAXNN),SAVE(MAXNA) INTEGER TFSTOU(MAXNN),TNXTOU(MAXNA),TFSTIN(MAXNN),TNXTIN(MAXNA) INTEGER DDPOS(MAXNN),DDNEG(MAXNN),NXTQUEUE(MAXNN) INTEGER I15(MAXNN),I16(MAXNN),I17(MAXNN) LOGICAL*1 SCAN(MAXNN),MARK(MAXNN) LOGICAL*1 FEASBL,QUIT,SWITCH,POSIT,PCHANGE COMMON /BLK1/LABEL/BLK2/PRDCSR/BLK7/SAVE COMMON /BLK10/TFSTOU/BLK11/TNXTOU/BLK12/TFSTIN/BLK13/TNXTIN COMMON /BLK14/NXTQUEUE/BLK15/I15/BLK16/I16/BLK17/I17 COMMON /BLK8/SCAN/BLK9/MARK EQUIVALENCE(DDPOS,TFSTOU),(DDNEG,TFSTIN) C C TIMING PARAMETERS C REAL*4 TIME0,TIME1 COMMON /T/TIME0,TIME1 C C NOTE 1 - C TO RUN IN LIMITED MEMORY SYSTEMS, DECLARE THE ARRAYS C STARTN, ENDN, NXTIN, NXTOU, FIN, FOU, LABEL, C PRDCSR, SAVE, TFSTOU, TNXTOU, TFSTIN, TNXTIN, C DDPOS,DDNEG,NXTQUEUE AS INTEGER*2 INSTEAD. C C NOTE 2 - C THIS ROUTINE MAKES NO EFFORT TO INITIALIZE WITH A FAVORABLE X C FROM AMONGST THOSE FLOW VECTORS THAT SATISFY COMPLEMENTARY SLACKNESS C WITH THE INITIAL REDUCED COST VECTOR RC. C IF A FAVORABLE X IS KNOWN, THEN IT CAN BE PASSED, TOGETHER C WITH THE CORRESPONDING ARRAYS U AND DFCT, TO THIS ROUTINE C DIRECTLY. THIS, HOWEVER, REQUIRES THAT THE CAPACITY C TIGHTENING PORTION AND THE FLOW INITIALIZATION PORTION C OF THIS ROUTINE (UP TO LINE LABELED 90) BE SKIPPED. C C NOTE 3 - C ALL PROBLEM DATA SHOULD BE LESS THAN LARGE IN MAGNITUDE, C AND LARGE SHOULD BE LESS THAN, SAY, 1/4 THE LARGEST INTEGER*4 C OF THE MACHINE USED. THIS WILL GUARD PRIMARILY AGAINST C OVERFLOW IN UNCAPACITATED PROBLEMS WHERE THE ARC CAPACITIES C ARE TAKEN FINITE BUT VERY LARGE. NOTE, HOWEVER, THAT AS IN C ALL CODES OPERATING WITH INTEGERS, OVERFLOW MAY OCCUR IF SOME C OF THE PROBLEM DATA TAKES VERY LARGE VALUES. C C NOTE 4 - C EACH COMMON BLOCK CONTAINS JUST ONE ARRAY, SO THE ARRAYS IN RELAX4 C CAN BE DIMENSIONED TO 1 AND TAKE THEIR DIMENSION FROM THE C MAIN CALLING ROUTINE. WITH THIS TRICK, RELAX4 NEED NOT BE RECOMPILED C IF THE ARRAY DIMENSIONS IN THE CALLING ROUTINE CHANGE. C IF YOUR FORTRAN COMPILER DOES NOT SUPPORT THIS FEATURE, THEN C CHANGE THE DIMENSION OF ALL THE ARRAYS TO BE THE SAME AS THE ONES C DECLARED IN THE MAIN CALLING PROGRAM. C C NOTE 5 - C DDPOS AND DDNEG ARE ARRAYS THAT GIVE THE DIRECTIONAL DERIVATIVES FOR C ALL POSITIVE AND NEGATIVE SINGLE-NODE PRICE CHANGES. THESE ARE USED C ONLY IN PHASE II OF THE INITIALIZATION PROCEDURE, BEFORE THE C LINKED LIST OF BALANCED ARCS COMES TO PLAY. THEREFORE, TO REDUCE C STORAGE, THEY ARE EQUIVALENCE TO TFSTOU AND TFSTIN, C WHICH ARE OF THE SAME SIZE (NUMBER OF NODES) AND ARE USED C ONLY AFTER THE TREE COMES INTO USE. C C--------------------------------------------------------------- C C INITIALIZATION PHASE I C C IN THIS PHASE, WE REDUCE THE ARC CAPACITIES BY AS MUCH AS C POSSIBLE WITHOUT CHANGING THE PROBLEM; C THEN WE SET THE INITIAL FLOW ARRAY X, TOGETHER WITH C THE CORRESPONDING ARRAYS U AND DFCT. C C THIS PHASE AND PHASE II (FROM HERE UP TO LINE LABELED 90) C CAN BE SKIPPED (BY SETTING REPEAT TO .TRUE.) IF THE CALLING PROGRAM C PLACES IN COMMON USER-CHOSEN VALUES FOR THE ARC FLOWS, THE RESIDUAL ARC C CAPACITIES, AND THE NODAL DEFICITS. WHEN THIS IS DONE, C IT IS CRITICAL THAT THE FLOW AND THE REDUCED COST FOR EACH ARC C SATISFY COMPLEMENTARY SLACKNESS C AND THE DFCT ARRAY PROPERLY CORRESPOND TO THE INITIAL ARC/FLOWS. C IF (REPEAT) GO TO 90 C DO 10 NODE=1,N NODE_DEF=DFCT(NODE) DDPOS(NODE)=NODE_DEF DDNEG(NODE)=-NODE_DEF MAXCAP=0 SCAPOU=0 ARC=FOU(NODE) 11 IF (ARC.GT.0) THEN IF (SCAPOU.LE.LARGE-U(ARC)) THEN SCAPOU=SCAPOU+U(ARC) ELSE GO TO 10 END IF ARC=NXTOU(ARC) GO TO 11 END IF IF (SCAPOU.LE.LARGE-NODE_DEF) THEN CAPOUT=SCAPOU+NODE_DEF ELSE GO TO 10 END IF IF (CAPOUT.LT.0) THEN C C PROBLEM IS INFEASIBLE - EXIT C PRINT*,'EXIT DURING CAPACITY ADJUSTMENT' PRINT*,'EXOGENOUS FLOW INTO NODE',NODE, $ ' EXCEEDS OUT CAPACITY' CALL PRINTFLOWS(NODE) GO TO 4400 END IF C SCAPIN=0 ARC=FIN(NODE) 12 IF (ARC.GT.0) THEN U(ARC)=MIN0(U(ARC),CAPOUT) IF (MAXCAP.LT.U(ARC)) MAXCAP=U(ARC) IF (SCAPIN.LE.LARGE-U(ARC)) THEN SCAPIN=SCAPIN+U(ARC) ELSE GO TO 10 END IF ARC=NXTIN(ARC) GO TO 12 END IF IF (SCAPIN.LE.LARGE+NODE_DEF) THEN CAPIN=SCAPIN-NODE_DEF ELSE GO TO 10 END IF IF (CAPIN.LT.0) THEN C C PROBLEM IS INFEASIBLE - EXIT C PRINT*,'EXIT DURING CAPACITY ADJUSTMENT' PRINT*,'EXOGENOUS FLOW OUT OF NODE',NODE, $ ' EXCEEDS IN CAPACITY' CALL PRINTFLOWS(NODE) GO TO 4400 END IF C ARC=FOU(NODE) 15 IF (ARC.GT.0) THEN U(ARC)=MIN0(U(ARC),CAPIN) ARC=NXTOU(ARC) GO TO 15 END IF 10 CONTINUE C C--------------------------------------------------------------- C C INITIALIZATION PHASE II C C IN THIS PHASE, WE INITIALIZE THE PRICES AND FLOWS BY EITHER CALLING C THE ROUTINE AUCTION OR BY PERFORMING ONLY SINGLE NODE (COORDINATE) C RELAXATION ITERATIONS. C IF (CRASH.EQ.1) THEN NSP=0 CALL AUCTION GO TO 70 END IF C C INITIALIZE THE ARC FLOWS TO SATISFY COMPLEMENTARY SLACKNESS WITH THE C PRICES. U(ARC) IS THE RESIDUAL CAPACITY OF ARC, AND X(ARC) IS THE FLOW. C THESE TWO ALWAYS ADD UP TO THE TOTAL CAPACITY FOR ARC. C ALSO COMPUTE THE DIRECTIONAL DERIVATIVES FOR EACH COORDINATE C AND COMPUTE THE ACTUAL DEFICITS. C DO 20 ARC=1,NA X(ARC) = 0 IF (RC(ARC).LE. 0) THEN T = U(ARC) T1 = STARTN(ARC) T2 = ENDN(ARC) DDPOS(T1) = DDPOS(T1) + T DDNEG(T2) = DDNEG(T2) + T IF (RC(ARC).LT. 0) THEN X(ARC) = T U(ARC) = 0 DFCT(T1) = DFCT(T1) + T DFCT(T2) = DFCT(T2) - T DDNEG(T1) = DDNEG(T1) - T DDPOS(T2) = DDPOS(T2) - T END IF END IF 20 CONTINUE C C MAKE 2 OR 3 PASSES THROUGH ALL NODES, PERFORMING ONLY C SINGLE NODE RELAXATION ITERATIONS. THE NUMBER OF C PASSES DEPENDS ON THE DENSITY OF THE NETWORK C IF (NA.GT.N*10) THEN NUMPASSES=2 ELSE NUMPASSES=3 END IF C DO 30 PASSES = 1,NUMPASSES DO 40 NODE=1,N IF (DFCT(NODE).EQ. 0) GO TO 40 IF (DDPOS(NODE).LE. 0) THEN C C COMPUTE DELPRC, THE STEPSIZE TO THE NEXT BREAKPOINT C IN THE DUAL COST AS THE PRICE OF NODE IS INCREASED. C [SINCE THE REDUCED COST OF ALL OUTGOING (RESP., C INCOMING) ARCS WILL DECREASE (RESP., INCREASE) AS C THE PRICE OF NODE IS INCREASED, THE NEXT BREAKPOINT IS C THE MINIMUM OF THE POSITIVE REDUCED COST ON OUTGOING C ARCS AND OF THE NEGATIVE REDUCED COST ON INCOMING ARCS.] C DELPRC = LARGE ARC = FOU(NODE) 51 IF (ARC.GT.0) THEN TRC = RC(ARC) IF ((TRC.GT. 0).AND.(TRC.LT.DELPRC)) THEN DELPRC = TRC END IF ARC = NXTOU(ARC) GOTO 51 END IF ARC = FIN(NODE) 52 IF (ARC.GT.0) THEN TRC = RC(ARC) IF ((TRC.LT.0).AND.(TRC.GT.-DELPRC)) THEN DELPRC = -TRC END IF ARC = NXTIN(ARC) GOTO 52 END IF C C IF NO BREAKPOINT IS LEFT AND DUAL ASCENT IS STILL C POSSIBLE, THE PROBLEM IS INFEASIBLE. C IF (DELPRC.GE.LARGE) THEN IF (DDPOS(NODE).EQ.0) GOTO 40 GOTO 4400 END IF C C DELPRC IS THE STEPSIZE TO NEXT BREAKPOINT. INCREASE C PRICE OF NODE BY DELPRC AND COMPUTE THE STEPSIZE TO C THE NEXT BREAKPOINT IN THE DUAL COST. C 53 NXTBRK = LARGE C C LOOK AT ALL ARCS OUT OF NODE. C ARC = FOU(NODE) 54 IF (ARC.GT.0) THEN TRC = RC(ARC) IF (TRC .EQ. 0) THEN T1 = ENDN(ARC) T = U(ARC) IF (T.GT.0) THEN DFCT(NODE) = DFCT(NODE) + T DFCT(T1) = DFCT(T1) - T X(ARC) = T U(ARC) = 0 ELSE T = X(ARC) END IF DDNEG(NODE) = DDNEG(NODE) - T DDPOS(T1) = DDPOS(T1) - T END IF C C DECREASE THE REDUCED COST ON ALL OUTGOING ARCS. C TRC = TRC - DELPRC IF ((TRC.GT.0).AND.(TRC.LT.NXTBRK)) THEN NXTBRK = TRC ELSE IF (TRC.EQ.0) THEN C C ARC GOES FROM INACTIVE TO BALANCED. UPDATE THE C RATE OF DUAL ASCENT AT NODE AND AT ITS NEIGHBOR. C DDPOS(NODE) = DDPOS(NODE) + U(ARC) DDNEG(ENDN(ARC)) = DDNEG(ENDN(ARC)) + U(ARC) END IF RC(ARC) = TRC ARC = NXTOU(ARC) GOTO 54 END IF C C LOOK AT ALL ARCS INTO NODE. C ARC = FIN(NODE) 55 IF (ARC.GT.0) THEN TRC = RC(ARC) IF (TRC.EQ.0) THEN T1 = STARTN(ARC) T = X(ARC) IF (T.GT.0) THEN DFCT(NODE) = DFCT(NODE) + T DFCT(T1) = DFCT(T1) - T U(ARC) = T X(ARC) = 0 ELSE T = U(ARC) END IF DDPOS(T1) = DDPOS(T1) - T DDNEG(NODE) = DDNEG(NODE) - T END IF C C INCREASE THE REDUCED COST ON ALL INCOMING ARCS. C TRC = TRC + DELPRC IF ((TRC.LT.0).AND.(TRC.GT.-NXTBRK)) THEN NXTBRK = -TRC ELSE IF (TRC.EQ.0) THEN C C ARC GOES FROM ACTIVE TO BALANCED. UPDATE THE C RATE OF DUAL ASCENT AT NODE AND AT ITS NEIGHBOR. C DDNEG(STARTN(ARC)) = DDNEG(STARTN(ARC)) + X(ARC) DDPOS(NODE) = DDPOS(NODE) + X(ARC) END IF RC(ARC) = TRC ARC = NXTIN(ARC) GOTO 55 END IF C C IF PRICE OF NODE CAN BE INCREASED FURTHER WITHOUT DECREASING C THE DUAL COST (EVEN IF THE DUAL COST DOESN'T INCREASE), C RETURN TO INCREASE THE PRICE FURTHER. C IF ((DDPOS(NODE).LE.0).AND.(NXTBRK.LT.LARGE)) THEN DELPRC = NXTBRK GOTO 53 END IF ELSE IF (DDNEG(NODE).LE.0) THEN C C COMPUTE DELPRC, THE STEPSIZE TO THE NEXT BREAKPOINT C IN THE DUAL COST AS THE PRICE OF NODE IS DECREASED. C [SINCE THE REDUCED COST OF ALL OUTGOING (RESP., C INCOMING) ARCS WILL INCREASE (RESP., DECREASE) AS C THE PRICE OF NODE IS DECREASED, THE NEXT BREAKPOINT IS C THE MINIMUM OF THE NEGATIVE REDUCED COST ON OUTGOING C ARCS AND OF THE POSITIVE REDUCED COST ON INCOMING ARCS.] C DELPRC = LARGE ARC = FOU(NODE) 61 IF (ARC.GT.0) THEN TRC = RC(ARC) IF ((TRC.LT.0).AND.(TRC.GT.-DELPRC)) THEN DELPRC = -TRC ENDIF ARC = NXTOU(ARC) GOTO 61 END IF ARC = FIN(NODE) 62 IF (ARC.GT.0) THEN TRC = RC(ARC) IF ((TRC.GT.0).AND.(TRC.LT.DELPRC)) THEN DELPRC = TRC END IF ARC = NXTIN(ARC) GOTO 62 END IF C C IF NO BREAKPOINT IS LEFT AND DUAL ASCENT IS STILL C POSSIBLE, THE PROBLEM IS INFEASIBLE. C IF (DELPRC.EQ.LARGE) THEN IF (DDNEG(NODE).EQ.0) GOTO 40 GOTO 4400 END IF C C DELPRC IS THE STEPSIZE TO NEXT BREAKPOINT. DECREASE C PRICE OF NODE BY DELPRC AND COMPUTE THE STEPSIZE TO C THE NEXT BREAKPOINT IN THE DUAL COST. C 63 NXTBRK = LARGE C C LOOK AT ALL ARCS OUT OF NODE. C ARC = FOU(NODE) 64 IF (ARC.GT.0) THEN TRC = RC(ARC) IF (TRC.EQ.0) THEN T1 = ENDN(ARC) T = X(ARC) IF (T.GT.0) THEN DFCT(NODE) = DFCT(NODE) - T DFCT(T1) = DFCT(T1) + T U(ARC) = T X(ARC) = 0 ELSE T = U(ARC) END IF DDPOS(NODE) = DDPOS(NODE) - T DDNEG(T1) = DDNEG(T1) - T END IF C C INCREASE THE REDUCED COST ON ALL OUTGOING ARCS. C TRC = TRC + DELPRC IF ((TRC.LT.0).AND.(TRC.GT.-NXTBRK)) THEN NXTBRK = -TRC ELSE IF (TRC.EQ.0) THEN C C ARC GOES FROM ACTIVE TO BALANCED. UPDATE THE C RATE OF DUAL ASCENT AT NODE AND AT ITS NEIGHBOR. C DDNEG(NODE) = DDNEG(NODE) + X(ARC) DDPOS(ENDN(ARC)) = DDPOS(ENDN(ARC)) + X(ARC) END IF RC(ARC) = TRC ARC = NXTOU(ARC) GOTO 64 END IF C C LOOK AT ALL ARCS INTO NODE. C ARC = FIN(NODE) 65 IF (ARC.GT.0) THEN TRC = RC(ARC) IF (TRC.EQ.0) THEN T1 = STARTN(ARC) T = U(ARC) IF (T.GT.0) THEN DFCT(NODE) = DFCT(NODE) - T DFCT(T1) = DFCT(T1) + T X(ARC) = T U(ARC) = 0 ELSE T = X(ARC) END IF DDNEG(T1) = DDNEG(T1) - T DDPOS(NODE) = DDPOS(NODE) - T END IF C C DECREASE THE REDUCED COST ON ALL INCOMING ARCS. C TRC = TRC - DELPRC IF ((TRC.GT.0).AND.(TRC.LT.NXTBRK)) THEN NXTBRK = TRC ELSE IF (TRC.EQ.0) THEN C C ARC GOES FROM INACTIVE TO BALANCED. UPDATE THE C RATE OF DUAL ASCENT AT NODE AND AT ITS NEIGHBOR. C DDPOS(STARTN(ARC)) = DDPOS(STARTN(ARC)) + U(ARC) DDNEG(NODE) = DDNEG(NODE) + U(ARC) END IF RC(ARC) = TRC ARC = NXTIN(ARC) GOTO 65 END IF C C IF PRICE OF NODE CAN BE DECREASED FURTHER WITHOUT DECREASING C THE DUAL COST (EVEN IF THE DUAL COST DOESN'T INCREASE), C RETURN TO DECREASE THE PRICE FURTHER. C IF ((DDNEG(NODE).LE.0).AND.(NXTBRK.LT.LARGE)) THEN DELPRC = NXTBRK GOTO 63 END IF END IF 40 CONTINUE 30 CONTINUE C C 70 CONTINUE C C READ TIME FOR INITIALIZATION C C TIME1 = LONG(362)/60.0 - TIME0 TIME1 = SECNDS(TIME0) C C--------------------------------------------------------------- C C INITIALIZE TREE DATA STRUCTURE. DO 80 I=1,N TFSTOU(I)=0 TFSTIN(I)=0 80 CONTINUE DO 81 I=1,NA TNXTIN(I)=-1 TNXTOU(I)=-1 IF (RC(I).EQ.0) THEN TNXTOU(I)=TFSTOU(STARTN(I)) TFSTOU(STARTN(I))=I TNXTIN(I)=TFSTIN(ENDN(I)) TFSTIN(ENDN(I))=I END IF 81 CONTINUE C C INITIALIZE OTHER VARIABLES. C 90 FEASBL=.TRUE. ITER=0 NMULTINODE=0 NUM_AUGM=0 NUM_ASCNT=0 NUM_PASSES=0 NUMNZ=N NUMNZ_NEW=0 SWITCH=.FALSE. DO 91 I=1,N MARK(I)=.FALSE. SCAN(I)=.FALSE. 91 CONTINUE NLABEL=0 C C RELAX4 USES AN ADAPTIVE STRATEGY TO DECIDE WHETHER TO C CONTINUE THE SCANNING PROCESS AFTER A MULTINODE PRICE CHANGE. C THE THRESHOLD PARAMETER TP AND TS THAT CONTROL C THIS STRATEGY ARE SET IN THE NEXT TWO LINES. C TP=10 TS=INT(N/15) C C INITIALIZE THE QUEUE OF NODES WITH NONZERO DEFICIT C DO 92 NODE=1,N-1 NXTQUEUE(NODE)=NODE+1 92 CONTINUE NXTQUEUE(N)=1 NODE=N LASTQUEUE=N C C--------------------------------------------------------------- C C START THE RELAXATION ALGORITHM. C 100 CONTINUE C C CODE FOR ADVANCING THE QUEUE OF NONZERO DEFICIT NODES C PREVNODE=NODE NODE=NXTQUEUE(NODE) DEFCIT=DFCT(NODE) IF (NODE.EQ.LASTQUEUE) THEN NUMNZ=NUMNZ_NEW NUMNZ_NEW=0 LASTQUEUE=PREVNODE NUM_PASSES=NUM_PASSES+1 END IF C C CODE FOR DELETING A NODE FROM THE QUEUE C IF (DEFCIT.EQ.0) THEN NXTNODE=NXTQUEUE(NODE) IF (NODE.EQ.NXTNODE) THEN RETURN ELSE NXTQUEUE(PREVNODE)=NXTNODE NXTQUEUE(NODE)=0 NODE=NXTNODE GO TO 100 END IF ELSE POSIT = (DEFCIT.GT.0) END IF C ITER=ITER+1 NUMNZ_NEW=NUMNZ_NEW+1 C IF (POSIT) THEN C C ATTEMPT A SINGLE NODE ITERATION FROM NODE WITH POSITIVE DEFICIT C PCHANGE = .FALSE. INDEF=DEFCIT DELX=0 NB=0 C C CHECK OUTGOING (PROBABLY) BALANCED ARCS FROM NODE. C ARC=TFSTOU(NODE) 4500 IF (ARC.GT.0) THEN IF ((RC(ARC).EQ.0).AND.(X(ARC).GT.0)) THEN DELX = DELX + X(ARC) NB = NB + 1 SAVE(NB) = ARC ENDIF ARC = TNXTOU(ARC) GOTO 4500 END IF C C CHECK INCOMING ARCS. C ARC = TFSTIN(NODE) 4501 IF (ARC.GT.0) THEN IF ((RC(ARC).EQ.0).AND.(U(ARC).GT.0)) THEN DELX = DELX + U(ARC) NB = NB + 1 SAVE(NB) = -ARC ENDIF ARC = TNXTIN(ARC) GOTO 4501 END IF C C END OF INITIAL NODE SCAN. C 4018 CONTINUE C C IF NO PRICE CHANGE IS POSSIBLE, EXIT. C IF (DELX.GT.DEFCIT) THEN QUIT = (DEFCIT .LT. INDEF) GO TO 4016 END IF C C RELAX4 SEARCHES ALONG THE ASCENT DIRECTION FOR THE C BEST PRICE BY CHECKING THE SLOPE OF THE DUAL COST C AT SUCCESSIVE BREAK POINTS. FIRST, WE C COMPUTE THE DISTANCE TO THE NEXT BREAK POINT. C DELPRC = LARGE ARC = FOU(NODE) 4502 IF (ARC .GT. 0) THEN RDCOST = RC(ARC) IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) THEN DELPRC = -RDCOST END IF ARC = NXTOU(ARC) GOTO 4502 END IF ARC = FIN(NODE) 4503 IF (ARC .GT. 0) THEN RDCOST = RC(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) THEN DELPRC = RDCOST END IF ARC = NXTIN(ARC) GOTO 4503 END IF C C CHECK IF PROBLEM IS INFEASIBLE. C IF ((DELX.LT.DEFCIT).AND.(DELPRC.EQ.LARGE)) THEN C C THE DUAL COST CAN BE DECREASED WITHOUT BOUND. C GO TO 4400 END IF C C SKIP FLOW ADJUSTEMT IF THERE IS NO FLOW TO MODIFY. C IF (DELX.EQ.0) GO TO 4014 C C ADJUST THE FLOW ON THE BALANCED ARCS INCIDENT TO NODE TO C MAINTAIN COMPLEMENTARY SLACKNESS AFTER THE PRICE CHANGE. C DO 4013 J=1,NB ARC=SAVE(J) IF (ARC.GT.0) THEN NODE2=ENDN(ARC) T1=X(ARC) DFCT(NODE2)=DFCT(NODE2)+T1 IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF U(ARC)=U(ARC)+T1 X(ARC)=0 ELSE NARC=-ARC NODE2=STARTN(NARC) T1=U(NARC) DFCT(NODE2)=DFCT(NODE2)+T1 IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF X(NARC)=X(NARC)+T1 U(NARC)=0 END IF 4013 CONTINUE DEFCIT=DEFCIT-DELX 4014 IF (DELPRC.EQ.LARGE) THEN QUIT=.TRUE. GO TO 4019 END IF C C NODE CORRESPONDS TO A DUAL ASCENT DIRECTION. DECREASE C THE PRICE OF NODE BY DELPRC AND COMPUTE THE STEPSIZE TO THE C NEXT BREAKPOINT IN THE DUAL COST. C NB=0 PCHANGE = .TRUE. DP=DELPRC DELPRC=LARGE DELX=0 ARC=FOU(NODE) 4504 IF (ARC.GT.0) THEN RDCOST=RC(ARC)+DP RC(ARC)=RDCOST IF (RDCOST.EQ.0) THEN NB=NB+1 SAVE(NB)=ARC DELX=DELX+X(ARC) END IF IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) DELPRC=-RDCOST ARC=NXTOU(ARC) GOTO 4504 END IF ARC=FIN(NODE) 4505 IF (ARC.GT.0) THEN RDCOST=RC(ARC)-DP RC(ARC)=RDCOST IF (RDCOST.EQ.0) THEN NB=NB+1 SAVE(NB)=-ARC DELX=DELX+U(ARC) END IF IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) DELPRC=RDCOST ARC=NXTIN(ARC) GOTO 4505 END IF C C RETURN TO CHECK IF ANOTHER PRICE CHANGE IS POSSIBLE. C GO TO 4018 C C PERFORM FLOW AUGMENTATION AT NODE. C 4016 DO 4011 J=1,NB ARC=SAVE(J) IF (ARC.GT.0) THEN C C ARC IS AN OUTGOING ARC FROM NODE. C NODE2=ENDN(ARC) T1=DFCT(NODE2) IF (T1.LT.0) THEN C C DECREASE THE TOTAL DEFICIT BY DECREASING FLOW OF ARC. C QUIT=.TRUE. T2=X(ARC) DX=MIN0(DEFCIT,-T1,T2) DEFCIT=DEFCIT-DX DFCT(NODE2)=T1+DX IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF X(ARC)=T2-DX U(ARC)=U(ARC)+DX IF (DEFCIT.EQ.0) GO TO 4019 END IF ELSE C C -ARC IS AN INCOMING ARC TO NODE. C NARC=-ARC NODE2=STARTN(NARC) T1=DFCT(NODE2) IF (T1.LT.0) THEN C C DECREASE THE TOTAL DEFICIT BY INCREASING FLOW OF -ARC. C QUIT=.TRUE. T2=U(NARC) DX=MIN0(DEFCIT,-T1,T2) DEFCIT=DEFCIT-DX DFCT(NODE2)=T1+DX IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF X(NARC)=X(NARC)+DX U(NARC)=T2-DX IF (DEFCIT.EQ.0) GO TO 4019 END IF END IF 4011 CONTINUE 4019 DFCT(NODE)=DEFCIT C C RECONSTRUCT THE LINKED LIST OF BALANCE ARCS INCIDENT TO THIS NODE. C FOR EACH ADJACENT NODE, WE ADD ANY NEWLY BLANCED ARCS C TO THE LIST, BUT DO NOT BOTHER REMOVING FORMERLY BALANCED ONES C (THEY WILL BE REMOVED THE NEXT TIME EACH ADJACENT NODE IS SCANNED). C IF (PCHANGE) THEN ARC = TFSTOU(NODE) TFSTOU(NODE) = 0 4506 IF (ARC .GT. 0) THEN NXTARC = TNXTOU(ARC) TNXTOU(ARC) = -1 ARC = NXTARC GOTO 4506 END IF ARC = TFSTIN(NODE) TFSTIN(NODE) = 0 4507 IF (ARC .GT. 0) THEN NXTARC = TNXTIN(ARC) TNXTIN(ARC) = -1 ARC = NXTARC GOTO 4507 END IF C C NOW ADD THE CURRENTLY BALANCED ARCS TO THE LIST FOR THIS NODE C (WHICH IS NOW EMPTY), AND THE APPROPRIATE ADJACENT ONES. C DO 4508 J=1,NB ARC = SAVE(J) IF (ARC.LE.0) ARC=-ARC IF (TNXTOU(ARC) .LT. 0) THEN TNXTOU(ARC) = TFSTOU(STARTN(ARC)) TFSTOU(STARTN(ARC)) = ARC END IF IF (TNXTIN(ARC) .LT. 0) THEN TNXTIN(ARC) = TFSTIN(ENDN(ARC)) TFSTIN(ENDN(ARC)) = ARC END IF 4508 CONTINUE END IF C C END OF SINGLE NODE ITERATION FOR POSITIVE DEFICIT NODE. C ELSE C C ATTEMPT A SINGLE NODE ITERATION FROM NODE WITH NEGATIVE DEFICIT C PCHANGE = .FALSE. DEFCIT=-DEFCIT INDEF=DEFCIT DELX=0 NB=0 C ARC = TFSTIN(NODE) 4509 IF (ARC .GT. 0) THEN IF ((RC(ARC) .EQ. 0) .AND. (X(ARC) .GT. 0)) THEN DELX = DELX + X(ARC) NB = NB + 1 SAVE(NB) = ARC END IF ARC = TNXTIN(ARC) GOTO 4509 END IF ARC=TFSTOU(NODE) 4510 IF (ARC .GT. 0) THEN IF ((RC(ARC) .EQ. 0) .AND. (U(ARC) .GT. 0)) THEN DELX = DELX + U(ARC) NB = NB + 1 SAVE(NB) = -ARC END IF ARC = TNXTOU(ARC) GOTO 4510 END IF C 4028 CONTINUE IF (DELX.GE.DEFCIT) THEN QUIT = (DEFCIT .LT. INDEF) GO TO 4026 END IF C C COMPUTE DISTANCE TO NEXT BREAKPOINT. C DELPRC = LARGE ARC = FIN(NODE) 4511 IF (ARC .GT. 0) THEN RDCOST = RC(ARC) IF ((RDCOST .LT. 0) .AND. (RDCOST.GT.-DELPRC)) THEN DELPRC = -RDCOST END IF ARC = NXTIN(ARC) GOTO 4511 END IF ARC = FOU(NODE) 4512 IF (ARC .GT. 0) THEN RDCOST = RC(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) THEN DELPRC = RDCOST END IF ARC = NXTOU(ARC) GOTO 4512 END IF C C CHECK IF PROBLEM IS INFEASIBLE. C IF ((DELX.LT.DEFCIT).AND.(DELPRC.EQ.LARGE)) THEN GO TO 4400 END IF IF (DELX.EQ.0) GO TO 4024 C C FLOW AUGMENTATION IS POSSIBLE. C DO 4023 J=1,NB ARC=SAVE(J) IF (ARC.GT.0) THEN NODE2=STARTN(ARC) T1=X(ARC) DFCT(NODE2)=DFCT(NODE2)-T1 IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF U(ARC)=U(ARC)+T1 X(ARC)=0 ELSE NARC=-ARC NODE2=ENDN(NARC) T1=U(NARC) DFCT(NODE2)=DFCT(NODE2)-T1 IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF X(NARC)=X(NARC)+T1 U(NARC)=0 END IF 4023 CONTINUE DEFCIT=DEFCIT-DELX 4024 IF (DELPRC.EQ.LARGE) THEN QUIT=.TRUE. GO TO 4029 END IF C C PRICE INCREASE AT NODE IS POSSIBLE. C NB=0 PCHANGE = .TRUE. DP=DELPRC DELPRC=LARGE DELX=0 ARC=FIN(NODE) 4513 IF (ARC.GT.0) THEN RDCOST=RC(ARC)+DP RC(ARC)=RDCOST IF (RDCOST.EQ.0) THEN NB=NB+1 SAVE(NB)=ARC DELX=DELX+X(ARC) END IF IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) DELPRC=-RDCOST ARC=NXTIN(ARC) GOTO 4513 END IF ARC=FOU(NODE) 4514 IF (ARC.GT.0) THEN RDCOST=RC(ARC)-DP RC(ARC)=RDCOST IF (RDCOST.EQ.0) THEN NB=NB+1 SAVE(NB)=-ARC DELX=DELX+U(ARC) END IF IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) DELPRC=RDCOST ARC=NXTOU(ARC) GOTO 4514 END IF GO TO 4028 C C PERFORM FLOW AUGMENTATION AT NODE. C 4026 DO 4021 J=1,NB ARC=SAVE(J) IF (ARC.GT.0) THEN C C ARC IS AN INCOMING ARC TO NODE. C NODE2=STARTN(ARC) T1=DFCT(NODE2) IF (T1.GT.0) THEN QUIT=.TRUE. T2=X(ARC) DX=MIN0(DEFCIT,T1,T2) DEFCIT=DEFCIT-DX DFCT(NODE2)=T1-DX IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF X(ARC)=T2-DX U(ARC)=U(ARC)+DX IF (DEFCIT.EQ.0) GO TO 4029 END IF ELSE C C -ARC IS AN OUTGOING ARC FROM NODE. C NARC=-ARC NODE2=ENDN(NARC) T1=DFCT(NODE2) IF (T1.GT.0) THEN QUIT=.TRUE. T2=U(NARC) DX=MIN0(DEFCIT,T1,T2) DEFCIT=DEFCIT-DX DFCT(NODE2)=T1-DX IF (NXTQUEUE(NODE2).EQ.0) THEN NXTQUEUE(PREVNODE)=NODE2 NXTQUEUE(NODE2)=NODE PREVNODE=NODE2 END IF X(NARC)=X(NARC)+DX U(NARC)=T2-DX IF (DEFCIT.EQ.0) GO TO 4029 END IF END IF 4021 CONTINUE 4029 DFCT(NODE)=-DEFCIT C C RECONSTRUCT THE LIST OF BALANCED ARCS INCIDENT TO NODE. C IF (PCHANGE) THEN ARC = TFSTOU(NODE) TFSTOU(NODE) = 0 4515 IF (ARC .GT. 0) THEN NXTARC = TNXTOU(ARC) TNXTOU(ARC) = -1 ARC = NXTARC GOTO 4515 END IF ARC = TFSTIN(NODE) TFSTIN(NODE) = 0 4516 IF (ARC .GT. 0) THEN NXTARC = TNXTIN(ARC) TNXTIN(ARC) = -1 ARC = NXTARC GOTO 4516 END IF C C NOW ADD THE CURRENTLY BALANCED ARCS TO THE LIST FOR THIS NODE C (WHICH IS NOW EMPTY), AND THE APPROPRIATE ADJACENT ONES. C DO 4517 J=1,NB ARC = SAVE(J) IF (ARC.LE.0) ARC=-ARC IF (TNXTOU(ARC) .LT. 0) THEN TNXTOU(ARC) = TFSTOU(STARTN(ARC)) TFSTOU(STARTN(ARC)) = ARC END IF IF (TNXTIN(ARC) .LT. 0) THEN TNXTIN(ARC) = TFSTIN(ENDN(ARC)) TFSTIN(ENDN(ARC)) = ARC END IF 4517 CONTINUE END IF C C END OF SINGLE NODE ITERATION FOR A NEGATIVE DEFICIT NODE. C END IF C IF (QUIT.OR.(NUM_PASSES.LE.3)) GO TO 100 C C DO A MULTINODE ITERATION FROM NODE. C NMULTINODE=NMULTINODE+1 C C IF NUMBER OF NONZERO DEFICIT NODES IS SMALL, CONTINUE C LABELING UNTIL A FLOW AUGMENTATION IS DONE. C SWITCH = (NUMNZ.LT.TP) C C UNMARK NODES LABELED EARLIER. C DO 4090 J=1,NLABEL NODE2=LABEL(J) MARK(NODE2)=.FALSE. SCAN(NODE2)=.FALSE. 4090 CONTINUE C C INITIALIZE LABELING. C NLABEL=1 LABEL(1)=NODE MARK(NODE)=.TRUE. PRDCSR(NODE)=0 C C SCAN STARTING NODE. C SCAN(NODE)=.TRUE. NSCAN=1 DM=DFCT(NODE) DELX=0 DO 4095 J=1,NB ARC=SAVE(J) IF (ARC.GT.0) THEN IF (POSIT) THEN NODE2=ENDN(ARC) ELSE NODE2=STARTN(ARC) END IF IF (.NOT.MARK(NODE2)) THEN NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 PRDCSR(NODE2)=ARC MARK(NODE2)=.TRUE. DELX=DELX+X(ARC) END IF ELSE NARC=-ARC IF (POSIT) THEN NODE2=STARTN(NARC) ELSE NODE2=ENDN(NARC) END IF IF (.NOT.MARK(NODE2)) THEN NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 PRDCSR(NODE2)=ARC MARK(NODE2)=.TRUE. DELX=DELX+U(NARC) END IF END IF 4095 CONTINUE C C START SCANNING A LABELED BUT UNSCANNED NODE. C 4120 NSCAN=NSCAN+1 C C CHECK TO SEE IF SWITCH NEEDS TO BE SET TO TRUE SO TO C CONTINUE SCANNING EVEN AFTER A PRICE CHANGE. C SWITCH = SWITCH .OR. $((NSCAN .GT. TS).AND.(NUMNZ.LT.TS)) C C SCANNING WILL CONTINUE UNTIL EITHER AN OVERESTIMATE OF THE RESIDUAL C CAPACITY ACROSS THE CUT CORRESPONDING TO THE SCANNED SET OF NODES (CALLED C DELX) EXCEEDS THE ABSOLUTE VALUE OF THE TOTAL DEFICIT OF THE SCANNED C NODES (CALLED DM), OR ELSE AN AUGMENTING PATH IS FOUND. ARCS THAT ARE C IN THE TREE BUT ARE NOT BALANCED ARE REMOVED AS PART OF THE SCANNING C PROCESS. C I=LABEL(NSCAN) SCAN(I)=.TRUE. NAUGNOD=0 IF (POSIT) THEN C C SCANNING NODE I IN CASE OF POSITIVE DEFICIT. C PRVARC=0 ARC = TFSTOU(I) 4518 IF (ARC.GT.0) THEN C C ARC IS AN OUTGOING ARC FROM NODE. C IF (RC(ARC) .EQ. 0) THEN IF (X(ARC) .GT. 0) THEN NODE2=ENDN(ARC) IF (.NOT. MARK(NODE2)) THEN C C NODE2 IS NOT LABELED, SO ADD NODE2 TO THE LABELED SET. C PRDCSR(NODE2)=ARC IF (DFCT(NODE2).LT.0) THEN NAUGNOD=NAUGNOD+1 SAVE(NAUGNOD)=NODE2 END IF NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. DELX=DELX+X(ARC) END IF END IF PRVARC = ARC ARC = TNXTOU(ARC) ELSE TMPARC = ARC ARC = TNXTOU(ARC) TNXTOU(TMPARC) = -1 IF (PRVARC .EQ. 0) THEN TFSTOU(I) = ARC ELSE TNXTOU(PRVARC) = ARC END IF END IF GOTO 4518 END IF PRVARC = 0 ARC=TFSTIN(I) 4519 IF (ARC.GT.0) THEN C C ARC IS AN INCOMING ARC INTO NODE. C IF (RC(ARC) .EQ. 0) THEN IF (U(ARC) .GT. 0) THEN NODE2=STARTN(ARC) IF (.NOT. MARK(NODE2)) THEN C C NODE2 IS NOT LABELED, SO ADD NODE2 TO THE LABELED SET. C PRDCSR(NODE2)=-ARC IF (DFCT(NODE2).LT.0) THEN NAUGNOD=NAUGNOD+1 SAVE(NAUGNOD)=NODE2 END IF NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. DELX=DELX+U(ARC) END IF END IF PRVARC = ARC ARC = TNXTIN(ARC) ELSE TMPARC = ARC ARC = TNXTIN(ARC) TNXTIN(TMPARC) = -1 IF (PRVARC .EQ. 0) THEN TFSTIN(I) = ARC ELSE TNXTIN(PRVARC) = ARC END IF END IF GOTO 4519 END IF C C CORRECT THE RESIDUAL CAPACITY OF THE SCANNED NODE CUT. C ARC=PRDCSR(I) IF (ARC.GT.0) THEN DELX=DELX-X(ARC) ELSE DELX=DELX-U(-ARC) END IF C C END OF SCANNING OF NODE I FOR POSITIVE DEFICIT CASE. C ELSE C C SCANNING NODE I FOR NEGATIVE DEFICIT CASE. C PRVARC = 0 ARC=TFSTIN(I) 4520 IF (ARC.GT.0) THEN IF (RC(ARC) .EQ. 0) THEN IF (X(ARC) .GT. 0) THEN NODE2=STARTN(ARC) IF (.NOT. MARK(NODE2)) THEN PRDCSR(NODE2)=ARC IF (DFCT(NODE2).GT.0) THEN NAUGNOD=NAUGNOD+1 SAVE(NAUGNOD)=NODE2 END IF NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. DELX=DELX+X(ARC) END IF END IF PRVARC = ARC ARC = TNXTIN(ARC) ELSE TMPARC = ARC ARC = TNXTIN(ARC) TNXTIN(TMPARC) = -1 IF (PRVARC .EQ. 0) THEN TFSTIN(I) = ARC ELSE TNXTIN(PRVARC) = ARC END IF END IF GOTO 4520 END IF C PRVARC = 0 ARC = TFSTOU(I) 4521 IF (ARC.GT.0) THEN IF (RC(ARC) .EQ. 0) THEN IF (U(ARC) .GT. 0) THEN NODE2=ENDN(ARC) IF (.NOT. MARK(NODE2)) THEN PRDCSR(NODE2)=-ARC IF (DFCT(NODE2).GT.0) THEN NAUGNOD=NAUGNOD+1 SAVE(NAUGNOD)=NODE2 END IF NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. DELX=DELX+U(ARC) END IF END IF PRVARC = ARC ARC = TNXTOU(ARC) ELSE TMPARC = ARC ARC = TNXTOU(ARC) TNXTOU(TMPARC) = -1 IF (PRVARC .EQ. 0) THEN TFSTOU(I) = ARC ELSE TNXTOU(PRVARC) = ARC END IF END IF GOTO 4521 END IF C ARC=PRDCSR(I) IF (ARC.GT.0) THEN DELX=DELX-X(ARC) ELSE DELX=DELX-U(-ARC) END IF END IF C C ADD DEFICIT OF NODE SCANNED TO DM. C DM=DM+DFCT(I) C C CHECK IF THE SET OF SCANNED NODES CORRESPOND C TO A DUAL ASCENT DIRECTION; IF YES, PERFORM A C PRICE ADJUSTMENT STEP, OTHERWISE CONTINUE LABELING. C IF (NSCAN.LT.NLABEL) THEN IF (SWITCH) GO TO 4210 IF ((DELX.GE.DM).AND.(DELX.GE.-DM)) GO TO 4210 END IF C C TRY A PRICE CHANGE. C [NOTE THAT SINCE DELX-ABS(DM) IS AN OVERESTIMATE OF ASCENT SLOPE, WE C MAY OCCASIONALLY TRY A DIRECTION THAT IS NOT AN ASCENT DIRECTION. C IN THIS CASE, THE ASCNT ROUTINES RETURN WITH QUIT=.FALSE., C SO WE CONTINUE LABELING NODES. C IF (POSIT) THEN CALL ASCNT1(DM,DELX,NLABEL,FEASBL, $ SWITCH,NSCAN,NODE,PREVNODE) NUM_ASCNT=NUM_ASCNT+1 ELSE CALL ASCNT2(DM,DELX,NLABEL,FEASBL, $ SWITCH,NSCAN,NODE,PREVNODE) NUM_ASCNT=NUM_ASCNT+1 END IF IF (.NOT.FEASBL) GO TO 4400 IF (.NOT.SWITCH) GO TO 100 C C STORE THOSE NEWLY LABELED NODES TO WHICH FLOW AUGMENTATION IS POSSIBLE. C NAUGNOD=0 DO 530 J=NSCAN+1,NLABEL NODE2=LABEL(J) IF (POSIT.AND.(DFCT(NODE2).LT.0)) THEN NAUGNOD=NAUGNOD+1 SAVE(NAUGNOD)=NODE2 ELSE IF ((.NOT.POSIT).AND.(DFCT(NODE2).GT.0)) THEN NAUGNOD=NAUGNOD+1 SAVE(NAUGNOD)=NODE2 END IF 530 CONTINUE C C CHECK IF FLOW AUGMENTATION IS POSSIBLE. C IF NOT, RETURN TO SCAN ANOTHER NODE. C 4210 CONTINUE C IF (NAUGNOD.EQ.0) GO TO 4120 C DO 4096 J=1,NAUGNOD NUM_AUGM=NUM_AUGM+1 AUGNOD=SAVE(J) IF (POSIT) THEN C C DO THE AUGMENTATION FROM NODE WITH POSITIVE DEFICIT. C DX=-DFCT(AUGNOD) IB=AUGNOD 1500 IF (IB.NE.NODE) THEN ARC=PRDCSR(IB) IF (ARC.GT.0) THEN DX=MIN0(DX,X(ARC)) IB=STARTN(ARC) ELSE DX=MIN0(DX,U(-ARC)) IB=ENDN(-ARC) END IF GOTO 1500 END IF DX=MIN0(DX,DFCT(NODE)) IF (DX .GT. 0) THEN C C INCREASE (DECREASE) THE FLOW OF ALL FORWARD (BACKWARD) C ARCS IN THE FLOW AUGMENTING PATH. ADJUST NODE DEFICIT ACCORDINGLY. C IF (NXTQUEUE(AUGNOD).EQ.0) THEN NXTQUEUE(PREVNODE)=AUGNOD NXTQUEUE(AUGNOD)=NODE PREVNODE=AUGNOD END IF DFCT(AUGNOD)=DFCT(AUGNOD)+DX DFCT(NODE)=DFCT(NODE)-DX IB=AUGNOD 1501 IF (IB.NE.NODE) THEN ARC=PRDCSR(IB) IF (ARC.GT.0) THEN X(ARC)=X(ARC)-DX U(ARC)=U(ARC)+DX IB=STARTN(ARC) ELSE NARC=-ARC X(NARC)=X(NARC)+DX U(NARC)=U(NARC)-DX IB=ENDN(NARC) END IF GOTO 1501 END IF END IF ELSE C C DO THE AUGMENTATION FROM NODE WITH NEGATIVE DEFICIT. C DX=DFCT(AUGNOD) IB=AUGNOD 1502 IF (IB.NE.NODE) THEN ARC=PRDCSR(IB) IF (ARC.GT.0) THEN DX=MIN0(DX,X(ARC)) IB=ENDN(ARC) ELSE DX=MIN0(DX,U(-ARC)) IB=STARTN(-ARC) END IF GOTO 1502 END IF DX=MIN0(DX,-DFCT(NODE)) IF (DX .GT. 0) THEN C C UPDATE THE FLOW AND DEFICITS. C IF (NXTQUEUE(AUGNOD).EQ.0) THEN NXTQUEUE(PREVNODE)=AUGNOD NXTQUEUE(AUGNOD)=NODE PREVNODE=AUGNOD END IF DFCT(AUGNOD)=DFCT(AUGNOD)-DX DFCT(NODE)=DFCT(NODE)+DX IB=AUGNOD 1503 IF (IB.NE.NODE) THEN ARC=PRDCSR(IB) IF (ARC.GT.0) THEN X(ARC)=X(ARC)-DX U(ARC)=U(ARC)+DX IB=ENDN(ARC) ELSE NARC=-ARC X(NARC)=X(NARC)+DX U(NARC)=U(NARC)-DX IB=STARTN(NARC) END IF GOTO 1503 END IF END IF END IF IF (DFCT(NODE).EQ.0) GO TO 100 IF (DFCT(AUGNOD).NE.0) SWITCH=.FALSE. 4096 CONTINUE C C IF NODE STILL HAS NONZERO DEFICIT AND ALL NEWLY C LABELED NODES HAVE SAME SIGN FOR THEIR DEFICIT AS C NODE, WE CAN CONTINUE LABELING. IN THIS CASE, CONTINUE C LABELING ONLY WHEN FLOW AUGMENTATION IS DONE C RELATIVELY INFREQUENTLY. C IF (SWITCH.AND.(ITER.GT.8*NUM_AUGM)) GO TO 4120 C C RETURN TO DO ANOTHER RELAXATION ITERATION. C GO TO 100 C C PROBLEM IS FOUND TO BE INFEASIBLE C 4400 PRINT*,' PROBLEM IS FOUND TO BE INFEASIBLE.' PRINT*, 'PROGRAM ENDED; PRESS TO EXIT' PAUSE STOP C END C C SUBROUTINE AUCTION IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C PURPOSE - THIS SUBROUTINE USES A VERSION OF THE AUCTION C ALGORITHM FOR MIN COST NETWORK FLOW TO COMPUTE A C GOOD INITIAL FLOW AND PRICES FOR THE PROBLEM. C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C C INPUT PARAMETERS C C N = NUMBER OF NODES C NA = NUMBER OF ARCS C LARGE = A VERY LARGE INTEGER TO REPRESENT INFINITY C (SEE NOTE 3) C STARTN(I) = STARTING NODE FOR THE I-TH ARC, I = 1,...,NA C ENDN(I) = ENDING NODE FOR THE I-TH ARC, I = 1,...,NA C FOU(I) = FIRST ARC LEAVING I-TH NODE, I = 1,...,N C NXTOU(I) = NEXT ARC LEAVING THE STARTING NODE OF J-TH ARC, C I = 1,...,NA C FIN(I) = FIRST ARC ENTERING I-TH NODE, I = 1,...,N C NXTIN(I) = NEXT ARC ENTERING THE ENDING NODE OF J-TH ARC, C I = 1,...,NA C INTEGER STARTN(MAXNA),ENDN(MAXNA) INTEGER FOU(MAXNN),NXTOU(MAXNA),FIN(MAXNN),NXTIN(MAXNA) COMMON /INPUT/N,NA,LARGE COMMON /ARRAYS/STARTN/ARRAYE/ENDN COMMON /BLK3/FOU/BLK4/NXTOU/BLK5/FIN/BLK6/NXTIN COMMON /CR/CRASH C C UPDATED PARAMETERS C C RC(J) = REDUCED COST OF ARC J, J = 1,...,NA C U(J) = RESIDUAL CAPACITY OF ARC J, C J = 1,...,NA C X(J) = FLOW ON ARC J, J = 1,...,NA C DFCT(I) = DEFICIT AT NODE I, I = 1,...,N C INTEGER RC(MAXNA),U(MAXNA),X(MAXNA),DFCT(MAXNN) COMMON /ARRAYRC/RC/ARRAYU/U/ARRAYX/X/ARRAYB/DFCT C C OUTPUT PARAMETERS C COMMON /OUTPUT/NMULTINODE,ITER,NUM_AUGM,NUM_ASCNT,NSP C C WORKING PARAMETERS C INTEGER P(MAXNN),PRDCSR(MAXNN),SAVE(MAXNA) INTEGER FPUSHF(MAXNN),NXTPUSHF(MAXNA) INTEGER FPUSHB(MAXNN),NXTPUSHB(MAXNA) INTEGER NXTQUEUE(MAXNN),EXTEND_ARC(MAXNN) INTEGER SB_LEVEL(MAXNN),SB_ARC(MAXNN) LOGICAL*1 PATH_ID(MAXNN) COMMON /BLK1/P/BLK2/PRDCSR/BLK7/SAVE COMMON /BLK10/FPUSHF/BLK11/NXTPUSHF/BLK12/FPUSHB/BLK13/NXTPUSHB COMMON /BLK14/NXTQUEUE/BLK15/EXTEND_ARC COMMON /BLK16/SB_LEVEL/BLK17/SB_ARC COMMON /BLK9/PATH_ID C C START INITIALIZATION USING AUCTION C NAUG=0 PASS=0 THRESH_DFCT=0 C C FACTOR DETERMINES BY HOW MUCH EPSILON IS REDUCED AT EACH MINIMIZATION C FACTOR=3 C C NUM_PASSES DETERMINES HOW MANY AUCTION SCALING PHASES ARE PERFORMED C NUM_PASSES=1 C SET ARC FLOWS TO SATISFY CS AND CALCULATE MAXCOST AND MINCOST MAXCOST=-LARGE MINCOST=LARGE DO 49 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) RDCOST=RC(ARC) IF (MAXCOST.LT.RDCOST) MAXCOST=RDCOST IF (MINCOST.GT.RDCOST) MINCOST=RDCOST IF (RDCOST.LT.0) THEN DFCT(START)=DFCT(START)+U(ARC) DFCT(END)=DFCT(END)-U(ARC) X(ARC)=U(ARC) U(ARC)=0 ELSE X(ARC)=0 END IF 49 CONTINUE C C SET INITIAL EPSILON C IF ((MAXCOST-MINCOST).GE.8) THEN EPS=INT((MAXCOST-MINCOST)/8) ELSE EPS=1 END IF C C SET INITIAL PRICES TO ZERO C DO 48 NODE=1,N P(NODE)=0 48 CONTINUE C C INITIALIZATION USING AUCTION/SHORTEST PATHS. C START OF THE FIRST SCALING PHASE. C 100 CONTINUE PASS=PASS+1 IF ((PASS.EQ.NUM_PASSES).OR.(EPS.EQ.1)) CRASH=0 NOLIST=0 C C CONSTRUCT LIST OF POSITIVE SURPLUS NODES AND QUEUE OF NEGATIVE SURPLUS C NODES C DO 110 NODE=1,N PRDCSR(NODE)=0 PATH_ID(NODE)=.FALSE. EXTEND_ARC(NODE)=0 SB_LEVEL(NODE)=-LARGE NXTQUEUE(NODE)=NODE+1 IF (DFCT(NODE).GT.0) THEN NOLIST=NOLIST+1 SAVE(NOLIST)=NODE END IF 110 CONTINUE C NXTQUEUE(N)=1 ROOT=1 PREVNODE=N LASTQUEUE=N C C INITIALIZATION WITH DOWN ITERATIONS FOR NEGATIVE SURPLUS NODES C DO 150 I=1,NOLIST NODE=SAVE(I) NSP=NSP+1 C C BUILD THE LIST OF ARCS W/ ROOM FOR PUSHING FLOW C AND FIND PROPER PRICE FOR DOWN ITERATION C BSTLEVEL=-LARGE FPUSHF(NODE)=0 ARC=FOU(NODE) 152 IF (ARC.GT.0) THEN IF (U(ARC).GT.0) THEN IF (FPUSHF(NODE).EQ.0) THEN FPUSHF(NODE)=ARC NXTPUSHF(ARC)=0 LAST=ARC ELSE NXTPUSHF(LAST)=ARC NXTPUSHF(ARC)=0 LAST=ARC END IF END IF IF (X(ARC).GT.0) THEN NEW_LEVEL = P(ENDN(ARC)) + RC(ARC) IF (NEW_LEVEL.GT.BSTLEVEL) THEN BSTLEVEL=NEW_LEVEL EXTARC=ARC END IF END IF ARC=NXTOU(ARC) GO TO 152 END IF C FPUSHB(NODE)=0 ARC=FIN(NODE) 154 IF (ARC.GT.0) THEN IF (X(ARC).GT.0) THEN IF (FPUSHB(NODE).EQ.0) THEN FPUSHB(NODE)=ARC NXTPUSHB(ARC)=0 LAST=ARC ELSE NXTPUSHB(LAST)=ARC NXTPUSHB(ARC)=0 LAST=ARC END IF END IF IF (U(ARC).GT.0) THEN NEW_LEVEL = P(STARTN(ARC)) - RC(ARC) IF (NEW_LEVEL.GT.BSTLEVEL) THEN BSTLEVEL=NEW_LEVEL EXTARC=-ARC END IF END IF ARC=NXTIN(ARC) GO TO 154 END IF EXTEND_ARC(NODE)=EXTARC P(NODE)=BSTLEVEL-EPS 150 CONTINUE C C START THE AUGMENTATION CYCLES OF THE NEW SCALING PHASE. C 200 CONTINUE IF (DFCT(ROOT).GE.THRESH_DFCT) GOTO 3000 TERM=ROOT PATH_ID(ROOT)=.TRUE. C C MAIN FORWARD ALGORITHM WITH ROOT AS ORIGIN. C 500 CONTINUE C START OF A NEW FORWARD ITERATION C PTERM=P(TERM) EXTARC=EXTEND_ARC(TERM) IF (EXTARC.EQ.0) THEN C C BUILD THE LIST OF ARCS W/ ROOM FOR PUSHING FLOW C FPUSHF(TERM)=0 ARC=FOU(TERM) 502 IF (ARC.GT.0) THEN IF (U(ARC).GT.0) THEN IF (FPUSHF(TERM).EQ.0) THEN FPUSHF(TERM)=ARC NXTPUSHF(ARC)=0 LAST=ARC ELSE NXTPUSHF(LAST)=ARC NXTPUSHF(ARC)=0 LAST=ARC END IF END IF ARC=NXTOU(ARC) GO TO 502 END IF C FPUSHB(TERM)=0 ARC=FIN(TERM) 504 IF (ARC.GT.0) THEN IF (X(ARC).GT.0) THEN IF (FPUSHB(TERM).EQ.0) THEN FPUSHB(TERM)=ARC NXTPUSHB(ARC)=0 LAST=ARC ELSE NXTPUSHB(LAST)=ARC NXTPUSHB(ARC)=0 LAST=ARC END IF END IF ARC=NXTIN(ARC) GO TO 504 END IF GO TO 600 END IF C C SPECULATIVE PATH EXTENSION ATTEMPT C NOTE: ARC>0 MEANS THAT ARC IS ORIENTED FROM THE ROOT TO THE DESTINATIONS C ARC<0 MEANS THAT ARC IS ORIENTED FROM THE DESTINATIONS TO THE ROOT C EXTARC=0 OR PRDARC=0, MEANS THE EXTENSION ARC OR THE PREDECESSOR ARC, C RESPECTIVELY, HAS NOT BEEN ESTABLISHED C 510 CONTINUE IF (EXTARC.GT.0) THEN IF (U(EXTARC).EQ.0) THEN SECLEVEL=SB_LEVEL(TERM) GO TO 580 END IF END=ENDN(EXTARC) BSTLEVEL=P(END)+RC(EXTARC) IF (PTERM.GE.BSTLEVEL) THEN IF (PATH_ID(END)) GOTO 1200 TERM=END PRDCSR(TERM)=EXTARC PATH_ID(TERM)=.TRUE. C C IF NEGATIVE SURPLUS NODE IS FOUND, DO AN AUGMENTATION C IF (DFCT(TERM).GT.0) GOTO 2000 C C RETURN FOR ANOTHER ITERATION C GO TO 500 END IF ELSE EXTARC=-EXTARC IF (X(EXTARC).EQ.0) THEN SECLEVEL=SB_LEVEL(TERM) GO TO 580 END IF START=STARTN(EXTARC) BSTLEVEL=P(START)-RC(EXTARC) IF (PTERM.GE.BSTLEVEL) THEN IF (PATH_ID(START)) GOTO 1200 TERM=START PRDCSR(TERM)=-EXTARC PATH_ID(TERM)=.TRUE. C C IF NEGATIVE SURPLUS NODE IS FOUND, DO AN AUGMENTATION C IF (DFCT(TERM).GT.0) GOTO 2000 C C RETURN FOR ANOTHER ITERATION C GO TO 500 END IF END IF C C SECOND BEST LOGIC TEST APPLIED TO SAVE A FULL NODE SCAN C IF OLD BEST LEVEL CONTINUES TO BE BEST GO FOR ANOTHER CONTRACTION C 550 SECLEVEL=SB_LEVEL(TERM) IF (BSTLEVEL.LE.SECLEVEL) GOTO 800 C C IF SECOND BEST CAN BE USED DO EITHER A CONTRACTION C OR START OVER WITH A SPECULATIVE EXTENSION C 580 IF (SECLEVEL.GT.-LARGE) THEN EXTARC=SB_ARC(TERM) IF (EXTARC.GT.0) THEN IF (U(EXTARC).EQ.0) GOTO 600 BSTLEVEL=P(ENDN(EXTARC))+RC(EXTARC) ELSE IF (X(-EXTARC).EQ.0) GOTO 600 BSTLEVEL=P(STARTN(-EXTARC))-RC(-EXTARC) END IF IF (BSTLEVEL.EQ.SECLEVEL) THEN SB_LEVEL(TERM)=-LARGE EXTEND_ARC(TERM)=EXTARC GOTO 800 END IF END IF C C EXTENSION/CONTRACTION ATTEMPT WAS UNSUCCESSFUL, SO SCAN TERMINAL NODE C 600 CONTINUE NSP=NSP+1 BSTLEVEL=LARGE SECLEVEL=LARGE ARC=FPUSHF(TERM) 700 IF (ARC.GT.0) THEN NEW_LEVEL = P(ENDN(ARC)) + RC(ARC) IF (NEW_LEVEL.LT.SECLEVEL) THEN IF (NEW_LEVEL.LT.BSTLEVEL) THEN SECLEVEL=BSTLEVEL BSTLEVEL=NEW_LEVEL SECARC=EXTARC EXTARC=ARC ELSE SECLEVEL=NEW_LEVEL SECARC=ARC END IF END IF ARC=NXTPUSHF(ARC) GOTO 700 END IF ARC=FPUSHB(TERM) 710 IF (ARC.GT.0) THEN NEW_LEVEL = P(STARTN(ARC)) - RC(ARC) IF (NEW_LEVEL.LT.SECLEVEL) THEN IF (NEW_LEVEL.LT.BSTLEVEL) THEN SECLEVEL=BSTLEVEL BSTLEVEL=NEW_LEVEL SECARC=EXTARC EXTARC=-ARC ELSE SECLEVEL=NEW_LEVEL SECARC=-ARC END IF END IF ARC=NXTPUSHB(ARC) GOTO 710 END IF SB_LEVEL(TERM)=SECLEVEL SB_ARC(TERM)=SECARC EXTEND_ARC(TERM)=EXTARC C C END OF NODE SCAN. C IF THE TERMINAL NODE IS THE ROOT, ADJUST ITS PRICE AND CHANGE ROOT C 800 IF (TERM.EQ.ROOT) THEN P(TERM)=BSTLEVEL+EPS IF (PTERM.GE.LARGE) THEN PRINT*,'NO PATH TO THE DESTINATION' PRINT*,' PROBLEM IS FOUND TO BE INFEASIBLE.' PRINT*, 'PROGRAM ENDED; PRESS TO EXIT' PAUSE STOP END IF PATH_ID(ROOT)=.FALSE. PREVNODE=ROOT ROOT=NXTQUEUE(ROOT) GO TO 200 END IF C C CHECK WHETHER EXTENSION OR CONTRACTION C PRD=PRDCSR(TERM) IF (PRD.GT.0) THEN PR_TERM=STARTN(PRD) PREVLEVEL=P(PR_TERM)-RC(PRD) ELSE PR_TERM=ENDN(-PRD) PREVLEVEL=P(PR_TERM)+RC(-PRD) END IF C IF (PREVLEVEL.GT.BSTLEVEL) THEN C C PATH EXTENSION C IF (PREVLEVEL.GE.BSTLEVEL+EPS) THEN P(TERM)=BSTLEVEL+EPS ELSE P(TERM)=PREVLEVEL END IF IF (EXTARC.GT.0) THEN END=ENDN(EXTARC) IF (PATH_ID(END)) GOTO 1200 TERM=END ELSE START=STARTN(-EXTARC) IF (PATH_ID(START)) GOTO 1200 TERM=START END IF PRDCSR(TERM)=EXTARC PATH_ID(TERM)=.TRUE. C C IF NEGATIVE SURPLUS NODE IS FOUND, DO AN AUGMENTATION C IF (DFCT(TERM).GT.0) GOTO 2000 C C RETURN FOR ANOTHER ITERATION C GO TO 500 ELSE C C PATH CONTRACTION. C P(TERM)=BSTLEVEL+EPS PATH_ID(TERM)=.FALSE. TERM=PR_TERM IF (PR_TERM.NE.ROOT) THEN IF (BSTLEVEL.LE.PTERM+EPS) THEN GOTO 2000 END IF END IF PTERM=P(TERM) EXTARC=PRD IF (PRD.GT.0) THEN BSTLEVEL=BSTLEVEL+EPS+RC(PRD) ELSE BSTLEVEL=BSTLEVEL+EPS-RC(-PRD) END IF C C DO A SECOND BEST TEST AND IF THAT FAILS, DO A FULL NODE SCAN C GOTO 550 END IF C C A CYCLE IS ABOUT TO FORM; DO A RETREAT SEQUENCE. C 1200 CONTINUE C NODE=TERM 1600 IF (NODE.NE.ROOT) THEN PATH_ID(NODE)=.FALSE. PRD=PRDCSR(NODE) IF (PRD.GT.0) THEN PR_TERM=STARTN(PRD) IF (P(PR_TERM).EQ.P(NODE)+RC(PRD)+EPS) THEN NODE=PR_TERM GOTO 1600 END IF ELSE PR_TERM=ENDN(-PRD) IF (P(PR_TERM).EQ.P(NODE)-RC(-PRD)+EPS) THEN NODE=PR_TERM GOTO 1600 END IF END IF C C DO A FULL SCAN AND PRICE RISE AT PR_TERM C NSP=NSP+1 BSTLEVEL=LARGE SECLEVEL=LARGE ARC=FPUSHF(PR_TERM) 1700 IF (ARC.GT.0) THEN NEW_LEVEL = P(ENDN(ARC)) + RC(ARC) IF (NEW_LEVEL.LT.SECLEVEL) THEN IF (NEW_LEVEL.LT.BSTLEVEL) THEN SECLEVEL=BSTLEVEL BSTLEVEL=NEW_LEVEL SECARC=EXTARC EXTARC=ARC ELSE SECLEVEL=NEW_LEVEL SECARC=ARC END IF END IF ARC=NXTPUSHF(ARC) GOTO 1700 END IF C ARC=FPUSHB(PR_TERM) 1710 IF (ARC.GT.0) THEN NEW_LEVEL = P(STARTN(ARC)) - RC(ARC) IF (NEW_LEVEL.LT.SECLEVEL) THEN IF (NEW_LEVEL.LT.BSTLEVEL) THEN SECLEVEL=BSTLEVEL BSTLEVEL=NEW_LEVEL SECARC=EXTARC EXTARC=-ARC ELSE SECLEVEL=NEW_LEVEL SECARC=-ARC END IF END IF ARC=NXTPUSHB(ARC) GOTO 1710 END IF SB_LEVEL(PR_TERM)=SECLEVEL SB_ARC(PR_TERM)=SECARC EXTEND_ARC(PR_TERM)=EXTARC P(PR_TERM)=BSTLEVEL+EPS IF (PR_TERM.EQ.ROOT) THEN PREVNODE=ROOT PATH_ID(ROOT)=.FALSE. ROOT=NXTQUEUE(ROOT) GOTO 200 END IF PATH_ID(PR_TERM)=.FALSE. PRD=PRDCSR(PR_TERM) IF (PRD.GT.0) THEN TERM=STARTN(PRD) ELSE TERM=ENDN(-PRD) END IF IF (TERM.EQ.ROOT) THEN PREVNODE=ROOT PATH_ID(ROOT)=.FALSE. ROOT=NXTQUEUE(ROOT) GOTO 200 ELSE GOTO 2000 END IF END IF C C END OF AUCTION/SHORTEST PATH ROUTINE. C DO AUGMENTATION FROM ROOT AND CORRECT THE PUSH LISTS C 2000 CONTINUE INCR=-DFCT(ROOT) NODE = ROOT 2050 EXTARC=EXTEND_ARC(NODE) PATH_ID(NODE)=.FALSE. IF (EXTARC.GT.0) THEN NODE=ENDN(EXTARC) IF (INCR.GT.U(EXTARC)) INCR=U(EXTARC) ELSE NODE=STARTN(-EXTARC) IF (INCR.GT.X(-EXTARC)) INCR=X(-EXTARC) END IF IF (NODE.NE.TERM) GOTO 2050 PATH_ID(TERM)=.FALSE. IF (DFCT(TERM).GT.0) THEN IF (INCR.GT.DFCT(TERM)) INCR=DFCT(TERM) END IF C NODE = ROOT 2100 EXTARC=EXTEND_ARC(NODE) IF (EXTARC.GT.0) THEN END=ENDN(EXTARC) C C ADD ARC TO THE REDUCED GRAPH C IF (X(EXTARC).EQ.0) THEN NXTPUSHB(EXTARC)=FPUSHB(END) FPUSHB(END)=EXTARC NEW_LEVEL=P(NODE)-RC(EXTARC) IF (SB_LEVEL(END).GT.NEW_LEVEL) THEN SB_LEVEL(END)=NEW_LEVEL SB_ARC(END)=-EXTARC END IF END IF X(EXTARC)=X(EXTARC)+INCR U(EXTARC)=U(EXTARC)-INCR C C REMOVE ARC FROM THE REDUCED GRAPH C IF (U(EXTARC).EQ.0) THEN NAS=NAS+1 ARC=FPUSHF(NODE) IF (ARC.EQ.EXTARC) THEN FPUSHF(NODE)=NXTPUSHF(ARC) ELSE PREVARC=ARC ARC=NXTPUSHF(ARC) 2200 IF (ARC.GT.0) THEN IF (ARC.EQ.EXTARC) THEN NXTPUSHF(PREVARC)=NXTPUSHF(ARC) GO TO 2250 END IF PREVARC=ARC ARC=NXTPUSHF(ARC) GOTO 2200 END IF END IF END IF 2250 NODE=END ELSE EXTARC=-EXTARC START=STARTN(EXTARC) C C ADD ARC TO THE REDUCED GRAPH C IF (U(EXTARC).EQ.0) THEN NXTPUSHF(EXTARC)=FPUSHF(START) FPUSHF(START)=EXTARC NEW_LEVEL=P(NODE)+RC(EXTARC) IF (SB_LEVEL(START).GT.NEW_LEVEL) THEN SB_LEVEL(START)=NEW_LEVEL SB_ARC(START)=EXTARC END IF END IF U(EXTARC)=U(EXTARC)+INCR X(EXTARC)=X(EXTARC)-INCR C C REMOVE ARC FROM THE REDUCED GRAPH C IF (X(EXTARC).EQ.0) THEN NAS=NAS+1 ARC=FPUSHB(NODE) IF (ARC.EQ.EXTARC) THEN FPUSHB(NODE)=NXTPUSHB(ARC) ELSE PREVARC=ARC ARC=NXTPUSHB(ARC) 2300 IF (ARC.GT.0) THEN IF (ARC.EQ.EXTARC) THEN NXTPUSHB(PREVARC)=NXTPUSHB(ARC) GO TO 2350 END IF PREVARC=ARC ARC=NXTPUSHB(ARC) GOTO 2300 END IF END IF END IF 2350 NODE=START END IF IF (NODE.NE.TERM) GOTO 2100 DFCT(TERM)=DFCT(TERM)-INCR DFCT(ROOT)=DFCT(ROOT)+INCR C C INSERT TERM IN THE QUEUE IF IT HAS A LARGE ENOUGH SURPLUS C IF (DFCT(TERM).LT.THRESH_DFCT) THEN IF (NXTQUEUE(TERM).EQ.0) THEN NXTNODE=NXTQUEUE(ROOT) IF ((P(TERM).GE.P(NXTNODE)).AND.(ROOT.NE.NXTNODE)) THEN NXTQUEUE(ROOT)=TERM NXTQUEUE(TERM)=NXTNODE ELSE NXTQUEUE(PREVNODE)=TERM NXTQUEUE(TERM)=ROOT PREVNODE=TERM END IF END IF END IF C C IF ROOT HAS A LARGE ENOUGH SURPLUS, KEEP IT C IN THE QUEUE AND RETURN FOR ANOTHER ITERATION C IF (DFCT(ROOT).LT.THRESH_DFCT) THEN PREVNODE=ROOT ROOT=NXTQUEUE(ROOT) GO TO 200 END IF C C END OF AUGMENTATION CYCLE C 3000 CONTINUE C C CHECK FOR TERMINATION OF SCALING PHASE. IF SCALING PHASE IS C NOT FINISHED, ADVANCE THE QUEUE AND RETURN TO TAKE ANOTHER NODE. C NXTNODE=NXTQUEUE(ROOT) IF (ROOT.NE.NXTNODE) THEN NXTQUEUE(ROOT)=0 NXTQUEUE(PREVNODE)=NXTNODE ROOT=NXTNODE GO TO 200 END IF C C END OF SUBPROBLEM (SCALING PHASE). C 3600 CONTINUE C C REDUCE EPSILON. C EPS=INT(EPS/FACTOR) IF (EPS.LT.1) EPS=1 THRESH_DFCT=INT(THRESH_DFCT/FACTOR) IF (EPS.EQ.1) THRESH_DFCT=0 C C IF ANOTHER AUCTION SCALING PHASE REMAINS, RESET THE FLOWS & THE PUSH LISTS C ELSE RESET ARC FLOWS TO SATISFY CS AND COMPUTE REDUCED COSTS C IF (CRASH.EQ.1) THEN DO 3800 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) PSTART=P(START) PEND=P(END) IF (PSTART.GT.PEND+EPS+RC(ARC)) THEN RESID=U(ARC) IF (RESID.GT.0) THEN DFCT(START)=DFCT(START)+RESID DFCT(END)=DFCT(END)-RESID X(ARC)=X(ARC)+RESID U(ARC)=0 END IF ELSE IF (PSTART.LT.PEND-EPS+RC(ARC)) THEN FLOW=X(ARC) IF (FLOW.GT.0) THEN DFCT(START)=DFCT(START)-FLOW DFCT(END)=DFCT(END)+FLOW X(ARC)=0 U(ARC)=U(ARC)+FLOW END IF END IF END IF 3800 CONTINUE C C RETURN FOR ANOTHER PHASE C 3850 CONTINUE GOTO 100 ELSE CRASH=1 DO 3900 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) RED_COST=RC(ARC)+P(END)-P(START) IF (RED_COST.LT.0) THEN RESID=U(ARC) IF (RESID.GT.0) THEN DFCT(START)=DFCT(START)+RESID DFCT(END)=DFCT(END)-RESID X(ARC)=X(ARC)+RESID U(ARC)=0 END IF ELSE IF (RED_COST.GT.0) THEN FLOW=X(ARC) IF (FLOW.GT.0) THEN DFCT(START)=DFCT(START)-FLOW DFCT(END)=DFCT(END)+FLOW X(ARC)=0 U(ARC)=U(ARC)+FLOW END IF END IF END IF RC(ARC)=RED_COST 3900 CONTINUE END IF RETURN END C C SUBROUTINE PRINTFLOWS(NODE) IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C PURPOSE - THIS ROUTINE PRINTS THE DEFICIT AND THE FLOWS C OF ARCS INCIDENT TO NODE. IT IS USED FOR DIAGNOSTIC C PURPOSES IN CASE OF AN INFEASIBLE PROBLEM HERE. C IT CAN BE USED ALSO FOR MORE GENERAL DIAGNOSTIC C PURPOSES. C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C COMMON/ARRAYS/STARTN/ARRAYE/ENDN/ARRAYU/U/ARRAYX/X $/ARRAYB/DFCT/BLK3/FOU/BLK4/NXTOU/BLK5/FIN/BLK6/NXTIN C INTEGER STARTN(MAXNA),ENDN(MAXNA),U(MAXNA),X(MAXNA),DFCT(MAXNN) INTEGER FOU(MAXNN),NXTOU(MAXNA) INTEGER FIN(MAXNN),NXTIN(MAXNA) C C--------------------------------------------------------------- C PRINT*,'DEFICIT (I.E., NET FLOW OUT) OF NODE =',DFCT(NODE) PRINT*,'FLOWS AND CAPACITIES OF INCIDENT ARCS OF NODE',NODE C C CHECK ALL ARCS LEAVING NODE C IF (FOU(NODE).EQ.0) THEN PRINT*,'NO OUTGOING ARCS' ELSE ARC=FOU(NODE) 5 IF (ARC.GT.0) THEN PRINT*,'ARC',ARC,' BETWEEN NODES',NODE,ENDN(ARC) PRINT*,'FLOW =',X(ARC) PRINT*,'RESIDUAL CAPACITY =',U(ARC) ARC=NXTOU(ARC) GO TO 5 END IF END IF C C CHECK ALL ARCS INCOMING TO NODE C IF (FIN(NODE).EQ.0) THEN PRINT*,'NO INCOMING ARCS' ELSE ARC=FIN(NODE) 10 IF (ARC.GT.0) THEN PRINT*,'ARC',ARC,' BETWEEN NODES',STARTN(ARC),NODE PRINT*,'FLOW =',X(ARC) PRINT*,'RESIDUAL CAPACITY =',U(ARC) ARC=NXTIN(ARC) GO TO 10 END IF END IF C RETURN END C C SUBROUTINE ASCNT1(DM,DELX,NLABEL,FEASBL,SWITCH, $NSCAN,CURNODE,PREVNODE) IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C PURPOSE - THIS SUBROUTINE PERFORMS THE MULTI-NODE PRICE C ADJUSTMENT STEP FOR THE CASE WHERE THE SCANNED NODES C HAVE POSITIVE DEFICIT. IT FIRST CHECKS IF DECREASING C THE PRICE OF THE SCANNED NODES INCREASES THE DUAL COST. C IF YES, THEN IT DECREASES THE PRICE OF ALL SCANNED NODES. C THERE ARE TWO POSSIBILITIES FOR PRICE DECREASE: C IF SWITCH=.TRUE., THEN THE SET OF SCANNED NODES C CORRESPONDS TO AN ELEMENTARY DIRECTION OF MAXIMAL C RATE OF ASCENT, IN WHICH CASE THE PRICE OF ALL SCANNED C NODES ARE DECREASED UNTIL THE NEXT BREAKPOINT IN THE C DUAL COST IS ENCOUNTERED. AT THIS POINT, SOME ARC C BECOMES BALANCED AND MORE NODE(S) ARE ADDED TO THE C LABELED SET AND THE SUBROUTINE IS EXITED. C IF SWITCH=.FALSE., THEN THE PRICE OF ALL SCANNED NODES C ARE DECREASED UNTIL THE RATE OF ASCENT BECOMES C NEGATIVE (THIS CORRESPONDS TO THE PRICE ADJUSTMENT C STEP IN WHICH BOTH THE LINE SEARCH AND THE DEGENERATE C ASCENT ITERATION ARE IMPLEMENTED). C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C C INPUT PARAMETERS C C DM = TOTAL DEFICIT OF SCANNED NODES C SWITCH = .TRUE. IF LABELING IS TO CONTINUE AFTER PRICE CHANGE C NSCAN = NUMBER OF SCANNED NODES C CURNODE = MOST RECENTLY SCANNED NODE C N = NUMBER OF NODES C NA = NUMBER OF ARCS C LARGE = A VERY LARGE INTEGER TO REPRESENT INFINITY C (SEE NOTE 3) C STARTN(I) = STARTING NODE FOR THE I-TH ARC, I = 1,...,NA C ENDN(I) = ENDING NODE FOR THE I-TH ARC, I = 1,...,NA C FOU(I) = FIRST ARC LEAVING I-TH NODE, I = 1,...,N C NXTOU(I) = NEXT ARC LEAVING THE STARTING NODE OF J-TH ARC, C I = 1,...,NA C FIN(I) = FIRST ARC ENTERING I-TH NODE, I = 1,...,N C NXTIN(I) = NEXT ARC ENTERING THE ENDING NODE OF J-TH ARC, C I = 1,...,NA C INTEGER STARTN(MAXNA),ENDN(MAXNA) INTEGER FOU(MAXNN),NXTOU(MAXNA),FIN(MAXNN),NXTIN(MAXNA) COMMON /INPUT/N,NA,LARGE COMMON /ARRAYS/STARTN/ARRAYE/ENDN COMMON /BLK3/FOU/BLK4/NXTOU/BLK5/FIN/BLK6/NXTIN C C UPDATED PARAMETERS C C DELX = A LOWER ESTIMATE OF THE TOTAL FLOW ON BALANCED ARCS C IN THE SCANNED-NODES CUT C NLABEL = NUMBER OF LABELED NODES C FEASBL = .FALSE. IF PROBLEM IS FOUND TO BE INFEASIBLE C PREVNODE = THE NODE BEFORE CURNODE IN QUEUE C RC(J) = REDUCED COST OF ARC J, J = 1,...,NA C U(J) = RESIDUAL CAPACITY OF ARC J, C J = 1,...,NA C X(J) = FLOW ON ARC J, J = 1,...,NA C DFCT(I) = DEFICIT AT NODE I, I = 1,...,N C LABEL(K) = K-TH NODE LABELED, K = 1,NLABEL C PRDCSR(I) = PREDECESSOR OF NODE I IN TREE OF LABELED NODES C (O IF I IS UNLABELED), I = 1,...,N C TFSTOU(I) = FIRST BALANCED ARC OUT OF NODE I, I = 1,...,N C TNXTOU(J) = NEXT BALANCED ARC OUT OF THE STARTING NODE OF ARC J, C J = 1,...,NA C TFSTIN(I) = FIRST BALANCED ARC INTO NODE I, I = 1,...,N C TNXTIN(J) = NEXT BALANCED ARC INTO THE ENDING NODE OF ARC J, C J = 1,...,NA C NXTQUEUE(I) = NODE FOLLOWING NODE I IN THE FIFO QUEUE C (0 IF NODE IS NOT IN THE QUEUE), I = 1,...,N C SCAN(I) = .TRUE. IF NODE I IS SCANNED, I = 1,...,N C MARK(I) = .TRUE. IF NODE I IS LABELED, I = 1,...,N C INTEGER RC(MAXNA),U(MAXNA),X(MAXNA),DFCT(MAXNN) INTEGER LABEL(MAXNN),PRDCSR(MAXNN) INTEGER TFSTOU(MAXNN),TNXTOU(MAXNA),TFSTIN(MAXNN),TNXTIN(MAXNA) INTEGER NXTQUEUE(MAXNN) LOGICAL*1 SCAN(MAXNN),MARK(MAXNN) COMMON /ARRAYRC/RC/ARRAYU/U/ARRAYX/X/ARRAYB/DFCT COMMON /BLK1/LABEL/BLK2/PRDCSR COMMON /BLK10/TFSTOU/BLK11/TNXTOU/BLK12/TFSTIN/BLK13/TNXTIN COMMON /BLK14/NXTQUEUE COMMON /BLK8/SCAN/BLK9/MARK C C WORKING PARAMETERS C INTEGER SAVE(MAXNA) LOGICAL*1 SWITCH,FEASBL COMMON /BLK7/SAVE C C STORE THE ARCS BETWEEN THE SET OF SCANNED NODES AND C ITS COMPLEMENT IN SAVE AND COMPUTE DELPRC, THE STEPSIZE C TO THE NEXT BREAKPOINT IN THE DUAL COST IN THE DIRECTION C OF DECREASING PRICES OF THE SCANNED NODES. C [THE ARCS ARE STORED INTO SAVE BY LOOKING AT THE ARCS C INCIDENT TO EITHER THE SET OF SCANNED NODES OR ITS C COMPLEMENT, DEPENDING ON WHETHER NSCAN>N/2 OR NOT. C THIS IMPROVES THE EFFICIENCY OF STORING.] C DELPRC=LARGE DLX=0 NSAVE=0 IF (NSCAN.LE.N/2) THEN DO 1 I=1,NSCAN NODE=LABEL(I) ARC=FOU(NODE) 500 IF (ARC.GT.0) THEN C C ARC POINTS FROM SCANNED NODE TO AN UNSCANNED NODE. C NODE2=ENDN(ARC) IF (.NOT.SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE2).NE.ARC)) $ DLX=DLX+X(ARC) IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) $ DELPRC=-RDCOST END IF ARC=NXTOU(ARC) GOTO 500 END IF ARC=FIN(NODE) 501 IF (ARC.GT.0) THEN C C ARC POINTS FROM UNSCANNED NODE TO SCANNED NODE. C NODE2=STARTN(ARC) IF (.NOT.SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=-ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE2).NE.-ARC)) $ DLX=DLX+U(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) $ DELPRC=RDCOST END IF ARC=NXTIN(ARC) GOTO 501 END IF 1 CONTINUE ELSE DO 2 NODE=1,N IF (SCAN(NODE)) GO TO 2 ARC=FIN(NODE) 502 IF (ARC.GT.0) THEN NODE2=STARTN(ARC) IF (SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE).NE.ARC)) $ DLX=DLX+X(ARC) IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) $ DELPRC=-RDCOST END IF ARC=NXTIN(ARC) GOTO 502 END IF ARC=FOU(NODE) 503 IF (ARC.GT.0) THEN NODE2=ENDN(ARC) IF (SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=-ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE).NE.-ARC)) $ DLX=DLX+U(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) $ DELPRC=RDCOST END IF ARC=NXTOU(ARC) GOTO 503 END IF 2 CONTINUE END IF C C CHECK IF THE SET OF SCANNED NODES TRULY CORRESPONDS C TO A DUAL ASCENT DIRECTION. [HERE DELX+DLX IS THE EXACT C SUM OF THE FLOW ON ARCS FROM THE SCANNED SET TO THE C UNSCANNED SET PLUS THE (CAPACITY - FLOW) ON ARCS FROM C THE UNSCANNED SET TO THE SCANNED SET.] C IF THIS WERE NOT THE CASE, SET SWITCH TO .TRUE. C AND EXIT SUBROUTINE. C IF ((.NOT.SWITCH).AND.(DELX+DLX.GE.DM)) THEN SWITCH=.TRUE. RETURN END IF DELX=DELX+DLX C C CHECK THAT THE PROBLEM IS FEASIBLE. C 4 IF (DELPRC.EQ.LARGE) THEN C C WE CAN INCREASE THE DUAL COST WITHOUT BOUND, SO C THE PRIMAL PROBLEM IS INFEASIBLE. C FEASBL=.FALSE. RETURN END IF C C DECREASE THE PRICES OF THE SCANNED NODES, ADD MORE C NODES TO THE LABELED SET AND CHECK IF A NEWLY LABELED NODE C HAS NEGATIVE DEFICIT. C IF (SWITCH) THEN DO 7 I=1,NSAVE ARC=SAVE(I) IF (ARC.GT.0) THEN RC(ARC)=RC(ARC)+DELPRC IF (RC(ARC).EQ.0) THEN NODE2=ENDN(ARC) IF (TNXTOU(ARC) .LT. 0) THEN TNXTOU(ARC) = TFSTOU(STARTN(ARC)) TFSTOU(STARTN(ARC)) = ARC END IF IF (TNXTIN(ARC) .LT. 0) THEN TNXTIN(ARC) = TFSTIN(NODE2) TFSTIN(NODE2) = ARC END IF IF (.NOT.MARK(NODE2)) THEN PRDCSR(NODE2)=ARC NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. END IF END IF ELSE ARC=-ARC RC(ARC)=RC(ARC)-DELPRC IF (RC(ARC).EQ.0) THEN NODE2=STARTN(ARC) IF (TNXTOU(ARC) .LT. 0) THEN TNXTOU(ARC) = TFSTOU(NODE2) TFSTOU(NODE2) = ARC END IF IF (TNXTIN(ARC) .LT. 0) THEN TNXTIN(ARC) = TFSTIN(ENDN(ARC)) TFSTIN(ENDN(ARC)) = ARC END IF IF (.NOT.MARK(NODE2)) THEN PRDCSR(NODE2)=-ARC NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. END IF END IF END IF 7 CONTINUE RETURN ELSE C C DECREASE THE PRICES OF THE SCANNED NODES BY DELPRC. C ADJUST FLOW TO MAINTAIN COMPLEMENTARY SLACKNESS WITH C THE PRICES. C NB = 0 DO 6 I=1,NSAVE ARC=SAVE(I) IF (ARC.GT.0) THEN T1=RC(ARC) IF (T1.EQ.0) THEN T2=X(ARC) T3=STARTN(ARC) DFCT(T3)=DFCT(T3)-T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF T3=ENDN(ARC) DFCT(T3)=DFCT(T3)+T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF U(ARC)=U(ARC)+T2 X(ARC)=0 END IF RC(ARC)=T1+DELPRC IF (RC(ARC).EQ.0) THEN DELX=DELX+X(ARC) NB = NB + 1 PRDCSR(NB) = ARC END IF ELSE ARC=-ARC T1=RC(ARC) IF (T1.EQ.0) THEN T2=U(ARC) T3=STARTN(ARC) DFCT(T3)=DFCT(T3)+T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF T3=ENDN(ARC) DFCT(T3)=DFCT(T3)-T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF X(ARC)=X(ARC)+T2 U(ARC)=0 END IF RC(ARC)=T1-DELPRC IF (RC(ARC).EQ.0) THEN DELX=DELX+U(ARC) NB = NB + 1 PRDCSR(NB) = ARC END IF END IF 6 CONTINUE END IF C IF (DELX.LE.DM) THEN C C THE SET OF SCANNED NODES STILL CORRESPONDS TO A C DUAL (POSSIBLY DEGENERATE) ASCENT DIRECTON. COMPUTE C THE STEPSIZE DELPRC TO THE NEXT BREAKPOINT IN THE C DUAL COST. C DELPRC=LARGE DO 10 I=1,NSAVE ARC=SAVE(I) IF (ARC.GT.0) THEN RDCOST=RC(ARC) IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) DELPRC=-RDCOST ELSE ARC=-ARC RDCOST=RC(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) DELPRC=RDCOST END IF 10 CONTINUE IF ((DELPRC.NE.LARGE).OR.(DELX.LT.DM)) GO TO 4 END IF C C ADD NEW BALANCED ARCS TO THE SUPERSET OF BALANCED ARCS. C DO 9 I=1,NB ARC=PRDCSR(I) IF (TNXTIN(ARC).EQ.-1) THEN J=ENDN(ARC) TNXTIN(ARC)=TFSTIN(J) TFSTIN(J)=ARC END IF IF (TNXTOU(ARC).EQ.-1) THEN J=STARTN(ARC) TNXTOU(ARC)=TFSTOU(J) TFSTOU(J)=ARC END IF 9 CONTINUE RETURN END C C SUBROUTINE ASCNT2(DM,DELX,NLABEL,FEASBL,SWITCH, $NSCAN,CURNODE,PREVNODE) IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C PURPOSE - THIS ROUTINE IS ANALOGOUS TO ASCNT BUT FOR C THE CASE WHERE THE SCANNED NODES HAVE NEGATIVE DEFICIT. C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C INTEGER STARTN(MAXNA),ENDN(MAXNA) INTEGER FOU(MAXNN),NXTOU(MAXNA),FIN(MAXNN),NXTIN(MAXNA) COMMON /INPUT/N,NA,LARGE COMMON /ARRAYS/STARTN/ARRAYE/ENDN COMMON /BLK3/FOU/BLK4/NXTOU/BLK5/FIN/BLK6/NXTIN INTEGER RC(MAXNA),U(MAXNA),X(MAXNA),DFCT(MAXNN) INTEGER LABEL(MAXNN),PRDCSR(MAXNN) INTEGER TFSTOU(MAXNN),TNXTOU(MAXNA),TFSTIN(MAXNN),TNXTIN(MAXNA) INTEGER NXTQUEUE(MAXNN) LOGICAL*1 SCAN(MAXNN),MARK(MAXNN) COMMON /ARRAYRC/RC/ARRAYU/U/ARRAYX/X/ARRAYB/DFCT COMMON /BLK1/LABEL/BLK2/PRDCSR COMMON /BLK10/TFSTOU/BLK11/TNXTOU/BLK12/TFSTIN/BLK13/TNXTIN COMMON /BLK14/NXTQUEUE COMMON /BLK8/SCAN/BLK9/MARK INTEGER SAVE(MAXNA) LOGICAL*1 SWITCH,FEASBL COMMON /BLK7/SAVE C C STORE THE ARCS BETWEEN THE SET OF SCANNED NODES AND C ITS COMPLEMENT IN SAVE AND COMPUTE DELPRC, THE STEPSIZE C TO THE NEXT BREAKPOINT IN THE DUAL COST IN THE DIRECTION C OF INCREASING PRICES OF THE SCANNED NODES. C DELPRC=LARGE DLX=0 NSAVE=0 IF (NSCAN.LE.N/2) THEN DO 1 I=1,NSCAN NODE=LABEL(I) ARC=FIN(NODE) 500 IF (ARC.GT.0) THEN NODE2=STARTN(ARC) IF (.NOT.SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE2).NE.ARC)) $ DLX=DLX+X(ARC) IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) $ DELPRC=-RDCOST END IF ARC=NXTIN(ARC) GOTO 500 END IF ARC=FOU(NODE) 501 IF (ARC.GT.0) THEN NODE2=ENDN(ARC) IF (.NOT.SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=-ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE2).NE.-ARC)) $ DLX=DLX+U(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) $ DELPRC=RDCOST END IF ARC=NXTOU(ARC) GOTO 501 END IF 1 CONTINUE ELSE DO 2 NODE=1,N IF (SCAN(NODE)) GO TO 2 ARC=FOU(NODE) 502 IF (ARC.GT.0) THEN NODE2=ENDN(ARC) IF (SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE).NE.ARC)) $ DLX=DLX+X(ARC) IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) $ DELPRC=-RDCOST END IF ARC=NXTOU(ARC) GOTO 502 END IF ARC=FIN(NODE) 503 IF (ARC.GT.0) THEN NODE2=STARTN(ARC) IF (SCAN(NODE2)) THEN NSAVE=NSAVE+1 SAVE(NSAVE)=-ARC RDCOST=RC(ARC) IF ((RDCOST.EQ.0).AND.(PRDCSR(NODE).NE.-ARC)) $ DLX=DLX+U(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) $ DELPRC=RDCOST END IF ARC=NXTIN(ARC) GOTO 503 END IF 2 CONTINUE END IF C IF ((.NOT.SWITCH).AND.(DELX+DLX.GE.-DM)) THEN SWITCH=.TRUE. RETURN END IF DELX=DELX+DLX C C CHECK THAT THE PROBLEM IS FEASIBLE. C 4 IF (DELPRC.EQ.LARGE) THEN FEASBL=.FALSE. RETURN END IF C C INCREASE THE PRICES OF THE SCANNED NODES, ADD MORE C NODES TO THE LABELED SET AND CHECK IF A NEWLY LABELED NODE C HAS POSITIVE DEFICIT. C IF (SWITCH) THEN DO 7 I=1,NSAVE ARC=SAVE(I) IF (ARC.GT.0) THEN RC(ARC)=RC(ARC)+DELPRC IF (RC(ARC).EQ.0) THEN NODE2=STARTN(ARC) IF (TNXTOU(ARC) .LT. 0) THEN TNXTOU(ARC) = TFSTOU(NODE2) TFSTOU(NODE2) = ARC END IF IF (TNXTIN(ARC) .LT. 0) THEN TNXTIN(ARC) = TFSTIN(ENDN(ARC)) TFSTIN(ENDN(ARC)) = ARC END IF IF (.NOT.MARK(NODE2)) THEN PRDCSR(NODE2)=ARC NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. END IF END IF ELSE ARC=-ARC RC(ARC)=RC(ARC)-DELPRC IF (RC(ARC).EQ.0) THEN NODE2=ENDN(ARC) IF (TNXTOU(ARC) .LT. 0) THEN TNXTOU(ARC) = TFSTOU(STARTN(ARC)) TFSTOU(STARTN(ARC)) = ARC END IF IF (TNXTIN(ARC) .LT. 0) THEN TNXTIN(ARC) = TFSTIN(NODE2) TFSTIN(NODE2) = ARC END IF IF (.NOT.MARK(NODE2)) THEN PRDCSR(NODE2)=-ARC NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 MARK(NODE2)=.TRUE. END IF END IF END IF 7 CONTINUE RETURN ELSE NB = 0 DO 6 I=1,NSAVE ARC=SAVE(I) IF (ARC.GT.0) THEN T1=RC(ARC) IF (T1.EQ.0) THEN T2=X(ARC) T3=STARTN(ARC) DFCT(T3)=DFCT(T3)-T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF T3=ENDN(ARC) DFCT(T3)=DFCT(T3)+T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF U(ARC)=U(ARC)+T2 X(ARC)=0 END IF RC(ARC)=T1+DELPRC IF (RC(ARC).EQ.0) THEN DELX=DELX+X(ARC) NB = NB + 1 PRDCSR(NB) = ARC END IF ELSE ARC=-ARC T1=RC(ARC) IF (T1.EQ.0) THEN T2=U(ARC) T3=STARTN(ARC) DFCT(T3)=DFCT(T3)+T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF T3=ENDN(ARC) DFCT(T3)=DFCT(T3)-T2 IF (NXTQUEUE(T3).EQ.0) THEN NXTQUEUE(PREVNODE)=T3 NXTQUEUE(T3)=CURNODE PREVNODE=T3 END IF X(ARC)=X(ARC)+T2 U(ARC)=0 END IF RC(ARC)=T1-DELPRC IF (RC(ARC).EQ.0) THEN DELX=DELX+U(ARC) NB = NB + 1 PRDCSR(NB) = ARC END IF END IF 6 CONTINUE END IF C IF (DELX.LE.-DM) THEN DELPRC=LARGE DO 10 I=1,NSAVE ARC=SAVE(I) IF (ARC.GT.0) THEN RDCOST=RC(ARC) IF ((RDCOST.LT.0).AND.(RDCOST.GT.-DELPRC)) DELPRC=-RDCOST ELSE ARC=-ARC RDCOST=RC(ARC) IF ((RDCOST.GT.0).AND.(RDCOST.LT.DELPRC)) DELPRC=RDCOST END IF 10 CONTINUE IF ((DELPRC.NE.LARGE).OR.(DELX.LT.-DM)) GO TO 4 END IF C C ADD NEW BALANCED ARCS TO THE SUPERSET OF BALANCED ARCS. C DO 9 I=1,NB ARC=PRDCSR(I) IF (TNXTIN(ARC).EQ.-1) THEN J=ENDN(ARC) TNXTIN(ARC)=TFSTIN(J) TFSTIN(J)=ARC END IF IF (TNXTOU(ARC).EQ.-1) THEN J=STARTN(ARC) TNXTOU(ARC)=TFSTOU(J) TFSTOU(J)=ARC END IF 9 CONTINUE RETURN END C C SUBROUTINE SENSTV IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C PURPOSE - THIS SUBROUTINE ALLOWS THE USER TO INTERACTIVELY C EITHER CHANGE NODE SUPPLY, OR CHANGE FLOW UPPER BOUND C OF AN EXISTING ARC, OR CHANGE COST OF AN EXISTING ARC, C OR DELETE AN EXISTING ARC, OR ADD AN ARC. C [NOTE: IF THE OPERATING SYSTEM RESETS ALL LOCAL VARIABLES C TO SOME DEFAULT VALUE EACH TIME THIS SUBROUTINE IS CALLED, C THEN THE USER MUST MAKE THE FOLLOWING CURRENTLY LOCAL C VARIABLES: C DELARC, DARC, DU, ADDARC, AARC C GLOBAL (BY EITHER PUTTING THEM IN A COMMON BLOCK OR C PASSING THEM THROUGH THE CALLING PARAMETER).] C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C C INPUT PARAMETERS C C N = NUMBER OF NODES C NA = NUMBER OF ARCS C LARGE = A VERY LARGE INTEGER TO REPRESENT INFINITY C STARTN(J) = STARTING NODE FOR ARC J, J = 1,...,NA C ENDN(J) = ENDING NODE FOR ARC J, J = 1,...,NA C FOU(I) = FIRST ARC OUT OF NODE I, I = 1,...,N C NXTOU(J) = NEXT ARC OUT OF THE STARTING NODE OF ARC J, C J = 1,...,NA C FIN(I) = FIRST ARC INTO NODE I, I = 1,...,N C NXTIN(J) = NEXT ARC INTO THE ENDING NODE OF ARC J, C J = 1,...,NA C REPEAT = .TRUE. IF CAP(J)=X(J)+U(J), J=1,...,NA, ON INPUT C (.FALSE. OTHERWISE) C INTEGER STARTN(MAXNA),ENDN(MAXNA) INTEGER FOU(MAXNN),NXTOU(MAXNA),FIN(MAXNN),NXTIN(MAXNA) LOGICAL*1 REPEAT COMMON /INPUT/N,NA,LARGE COMMON /ARRAYS/STARTN/ARRAYE/ENDN COMMON /BLK3/FOU/BLK4/NXTOU/BLK5/FIN/BLK6/NXTIN COMMON /BLKR/REPEAT C C UPDATED PARAMETERS C C C(J) = COST OF ARC J, J = 1,...,NA C CAP(J) = CAPACITY OF ARC J, J = 1,...,NA C RC(J) = REDUCED COST OF ARC J, J = 1,...,NA C X(J) = FLOW ON ARC J, J = 1,...,NA C U(J) = CAP(J) - X(J) ON OUTPUT, J = 1,...,NA C DFCT(I) = DEMAND AT NODE I ON INPUT C AND ZERO ON OUTPUT, I = 1,...,N C TFSTOU(I) = FIRST BALANCED ARC OUT OF NODE I, I = 1,...,N C TNXTOU(J) = NEXT BALANCED ARC OUT OF THE STARTING NODE OF ARC J, C J = 1,...,NA C TFSTIN(I) = FIRST BALANCED ARC INTO NODE I, I = 1,...,N C TNXTIN(J) = NEXT BALANCED ARC INTO THE ENDING NODE OF ARC J, C J = 1,...,NA C INTEGER C(MAXNA),CAP(MAXNA) INTEGER RC(MAXNA),X(MAXNA),U(MAXNA),DFCT(MAXNN) INTEGER TFSTOU(MAXNN),TNXTOU(MAXNA),TFSTIN(MAXNN),TNXTIN(MAXNA) COMMON /ARRAYC/C/BLKCAP/CAP COMMON /ARRAYRC/RC/ARRAYX/X/ARRAYU/U/ARRAYB/DFCT COMMON/BLK10/TFSTOU/BLK11/TNXTOU/BLK12/TFSTIN/BLK13/TNXTIN C C WORKING PARAMETERS C INTEGER LABEL(MAXNN),PRICE(MAXNN),SAVE(MAXNA) LOGICAL*1 ADDARC,DELARC,MARK(MAXNN) COMMON/BLK1/LABEL/BLK2/PRICE/BLK7/SAVE COMMON/BLK9/MARK C IF (.NOT.REPEAT) THEN C C RESTORE THE ARC CAPACITY TO THAT OF THE ORIGINAL PROBLEM C [RECALL THAT, IN SOLVING THE ORIGINAL PROBLEM, RELAX4 C MAY HAVE DECREASED THE ARC CAPACITIES DURING C INITIALIZATION PHASE I.] THEN UPDATE FLOW AND THE C NODE DEFICITS TO MATCH THIS "NEW" CAPACITY. C DO 10 I=1,NA IF (RC(I).LT.0) THEN DFCT(STARTN(I))=DFCT(STARTN(I))+CAP(I)-X(I) DFCT(ENDN(I))=DFCT(ENDN(I))-CAP(I)+X(I) X(I)=CAP(I) ELSE U(I)=CAP(I)-X(I) END IF 10 CONTINUE REPEAT=.TRUE. END IF 20 WRITE(6,30) WRITE(6,40) WRITE(6,50) WRITE(6,60) WRITE(6,70) WRITE(6,80) IF (ADDARC) WRITE(6,90) AARC IF (DELARC) WRITE(6,100) DARC WRITE(6,105) 30 FORMAT(' INPUT 0 TO SOLVE THE MODIFIED PROBLEM') 40 FORMAT(' 1 TO CHANGE NODE FLOW SUPPLY') 50 FORMAT(' 2 TO CHANGE ARC FLOW UPPER BOUND') 60 FORMAT(' 3 TO CHANGE ARC COST') 70 FORMAT(' 4 TO DELETE AN ARC') 80 FORMAT(' 5 TO ADD AN ARC') 90 FORMAT(' 6 TO DELETE LAST ARC',I8,' ADDED') 100 FORMAT(' 7 TO RESTORE LAST ARC',I8,' DELETED') 105 FORMAT(' 8 TO QUIT PROGRAM') READ(5,*)SEL IF (SEL.EQ.8) THEN STOP ELSE IF (SEL.EQ.0) THEN RETURN ELSE IF (SEL.EQ.1) THEN C C CHANGE NODE FLOW SUPPLY. C 110 WRITE(6,120) 120 FORMAT(' INPUT NODE # WHERE FLOW SUPPLY IS INCREASED') READ(5,*)NODE IF ((NODE.LE.0).OR.(NODE.GT.N)) GO TO 110 WRITE(6,130) 130 FORMAT(' INPUT AMOUNT OF INCREASE (<0 VALUE ALLOWED)') READ(5,*)DELSUP DFCT(NODE)=DFCT(NODE)-DELSUP 140 WRITE(6,150) 150 FORMAT(' INPUT NODE NO. WHERE FLOW SUPPLY IS DECREASED') READ(5,*)NODE IF ((NODE.LE.0).OR.(NODE.GT.N)) GO TO 140 DFCT(NODE)=DFCT(NODE)+DELSUP ELSE IF (SEL.EQ.2) THEN C C CHANGE ARC FLOW CAPACITY. C [NOTE: U IS THE RESIDUAL CAPACITY, I.E., U = CAPACITY - X.] C 160 WRITE(6,170) 170 FORMAT(' INPUT ARC NO. AND THE INCREASE IN UPPER BOUND') READ(5,*)ARC,DELUB IF ((ARC.LE.0).OR.(ARC.GT.NA)) GO TO 160 IF (RC(ARC).LT.0) THEN C C ARC IS ACTIVE, SO MAINTAIN FLOW AT (NEW) CAPACITY. C DFCT(STARTN(ARC))=DFCT(STARTN(ARC))+DELUB DFCT(ENDN(ARC))=DFCT(ENDN(ARC))-DELUB X(ARC)=X(ARC)+DELUB IF (X(ARC).LT.0) WRITE(6,180) ELSE IF (RC(ARC).EQ.0) THEN IF (U(ARC).GE.-DELUB) THEN U(ARC)=U(ARC)+DELUB ELSE C C NEW CAPACITY IS LESS THAN CURRENT FLOW, SO DECREASE C FLOW TO NEW CAPACITY. C DEL=-DELUB-U(ARC) DFCT(STARTN(ARC))=DFCT(STARTN(ARC))-DEL DFCT(ENDN(ARC))=DFCT(ENDN(ARC))+DEL X(ARC)=X(ARC)-DEL IF (X(ARC).LT.0) WRITE(6,180) U(ARC)=0 END IF ELSE U(ARC)=U(ARC)+DELUB IF (U(ARC).LT.0) WRITE(6,180) 180 FORMAT(' FLOW UPPER BOUND IS NOW < 0') END IF ELSE IF (SEL.EQ.3) THEN C C CHANGE ARC COST. C 190 WRITE(6,200) 200 FORMAT(' INPUT ARC NO. & INCREASE IN COST') READ(5,*)ARC,DELC IF ((ARC.LE.0).OR.(ARC.GT.NA)) GO TO 190 IF ((RC(ARC).GE.0).AND.(RC(ARC)+DELC.LT.0)) THEN C C ARC BECOMES ACTIVE, SO INCREASE FLOW TO CAPACITY. C DFCT(STARTN(ARC))=DFCT(STARTN(ARC))+U(ARC) DFCT(ENDN(ARC))=DFCT(ENDN(ARC))-U(ARC) X(ARC)=U(ARC)+X(ARC) U(ARC)=0 ELSE IF ((RC(ARC).LE.0).AND.(RC(ARC)+DELC.GT.0))THEN C C ARC BECOMES INACTIVE, SO DECREASE FLOW TO ZERO. C DFCT(STARTN(ARC))=DFCT(STARTN(ARC))-X(ARC) DFCT(ENDN(ARC))=DFCT(ENDN(ARC))+X(ARC) U(ARC)=U(ARC)+X(ARC) X(ARC)=0 END IF RC(ARC)=RC(ARC)+DELC C(ARC)=C(ARC)+DELC C C IF ARC BECOMES BALANCED, CHECK TO ADD ARC TO TFSTOU, TFSTIN,.... C IF ((RC(ARC).EQ.0).AND.(DELC.NE.0)) THEN CALL ADDTR(ARC) END IF C ELSE IF ((SEL.EQ.4).OR.(SEL.EQ.6)) THEN C C DELETE AN ARC. C IF (SEL.EQ.6) THEN IF (.NOT.ADDARC) GO TO 20 ADDARC=.FALSE. ARC=AARC ELSE 210 WRITE(6,220) 220 FORMAT(' INPUT ARC NO. FOR DELETION') READ(5,*)ARC IF ((ARC.LE.0).OR.(ARC.GT.NA)) GO TO 210 DELARC=.TRUE. DARC=ARC DU=U(ARC)+X(ARC) END IF C C REMOVE ARC FROM THE ARRAY FIN, FOU, NXTIN, NXTOU. C ARC1=FOU(STARTN(ARC)) IF (ARC1.EQ.ARC) THEN FOU(STARTN(ARC))=NXTOU(ARC1) ELSE 230 ARC2=NXTOU(ARC1) IF (ARC2.EQ.ARC) THEN NXTOU(ARC1)=NXTOU(ARC2) GO TO 240 END IF ARC1=ARC2 IF (NXTOU(ARC1).GT.0) GO TO 230 END IF 240 ARC1=FIN(ENDN(ARC)) IF (ARC1.EQ.ARC) THEN FIN(ENDN(ARC))=NXTIN(ARC1) ELSE 250 ARC2=NXTIN(ARC1) IF (ARC2.EQ.ARC) THEN NXTIN(ARC1)=NXTIN(ARC2) GO TO 260 END IF ARC1=ARC2 IF (NXTIN(ARC1).GT.0) GO TO 250 END IF C C REMOVE ARC FROM THE ARRAY TFSTIN, TFSTOU, TNXTIN, TNXTOU. C 260 ARC1=TFSTOU(STARTN(ARC)) IF (ARC1.EQ.0) GO TO 262 IF (ARC1.EQ.ARC) THEN TFSTOU(STARTN(ARC))=TNXTOU(ARC1) ELSE 261 ARC2=TNXTOU(ARC1) IF (ARC2.EQ.ARC) THEN TNXTOU(ARC1)=TNXTOU(ARC2) GO TO 262 END IF ARC1=ARC2 IF (TNXTOU(ARC1).GT.0) GO TO 261 END IF 262 ARC1=TFSTIN(ENDN(ARC)) IF (ARC1.EQ.0) GO TO 264 IF (ARC1.EQ.ARC) THEN TFSTIN(ENDN(ARC))=TNXTIN(ARC1) ELSE 263 ARC2=TNXTIN(ARC1) IF (ARC2.EQ.ARC) THEN TNXTIN(ARC1)=TNXTIN(ARC2) GO TO 264 END IF ARC1=ARC2 IF (TNXTIN(ARC1).GT.0) GO TO 263 END IF 264 TNXTOU(ARC) = -1 TNXTIN(ARC) = -1 C C REMOVE FLOW OF ARC FROM NETWORK BY SETTING ITS FLOW C AND CAPACITY TO 0. C DFCT(STARTN(ARC))=DFCT(STARTN(ARC))-X(ARC) DFCT(ENDN(ARC))=DFCT(ENDN(ARC))+X(ARC) X(ARC)=0 U(ARC)=0 ELSE IF ((SEL.EQ.5).OR.(SEL.EQ.7)) THEN IF (SEL.EQ.7) THEN IF (.NOT.DELARC) GO TO 20 IARC=DARC IH=STARTN(IARC) IT=ENDN(IARC) DELARC=.FALSE. IU=DU ELSE 270 WRITE(6,280)NA+1 280 FORMAT(' INPUT HEAD & TAIL NODES OF NEW ARC',I8) READ(5,*)IH,IT IF ((IH.LE.0).OR.(IH.GT.N).OR.(IT.LE.0).OR.(IT.GT.N))GO TO 270 290 WRITE(6,300) 300 FORMAT(' INPUT COST & FLOW UPPER BD') READ(5,*)IC,IU IF (IU.LT.0) GO TO 290 ADDARC=.TRUE. AARC=NA+1 NA=NA+1 C(NA)=IC STARTN(NA)=IH ENDN(NA)=IT IARC=NA END IF C C DETERMINE THE DUAL PRICES AT IH AND IT: C FIRST SET THE PRICE AT IH TO ZERO AND THEN CONSTRUCT THE C PRICE AT OTHER NODES BY BREADTH-FIRST SEARCH AND USING C THE FACT THAT C RC(ARC) = C(ARC) - PRICE(STARTN(ARC)) + PRICE(ENDN(ARC)). C THE NETWORK (GIVEN BY FOU, NXTOU, FIN, NXTIN) NEED NOT BE CONNECTED. C NSCAN=0 NLABEL=1 LABEL(1)=IH PRICE(IH)=0 DO 310 I=1,N 310 MARK(I)=.FALSE. MARK(IH)=.TRUE. 320 IF (NLABEL.GT.NSCAN) THEN NSCAN=NSCAN+1 NODE=LABEL(NSCAN) ARC=FOU(NODE) 330 IF (ARC.GT.0) THEN NODE2=ENDN(ARC) IF (.NOT.MARK(NODE2)) THEN MARK(NODE2)=.TRUE. PRICE(NODE2)=RC(ARC)-C(ARC)+PRICE(NODE) IF (NODE2.EQ.IT) GO TO 370 NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 END IF ARC=NXTOU(ARC) GO TO 330 END IF ARC=FIN(NODE) 340 IF (ARC.GT.0) THEN NODE2=STARTN(ARC) IF (.NOT.MARK(NODE2)) THEN MARK(NODE2)=.TRUE. PRICE(NODE2)=C(ARC)-RC(ARC)+PRICE(NODE) IF (NODE2.EQ.IT) GO TO 370 NLABEL=NLABEL+1 LABEL(NLABEL)=NODE2 END IF ARC=NXTIN(ARC) GO TO 340 END IF GO TO 320 END IF PRICE(IT)=-C(IARC) C C COMPUTE REDUCED COST OF THE NEW ARC AND UPDATE FLOW AND DEFICIT C ACCORDINGLY. C 370 RC(IARC)=C(IARC)+PRICE(IT) IF (RC(IARC).LT.0) THEN DFCT(IH)=DFCT(IH)+IU DFCT(IT)=DFCT(IT)-IU X(IARC)=IU U(IARC)=0 ELSE X(IARC)=0 U(IARC)=IU END IF NXTOU(IARC)=FOU(IH) FOU(IH)=IARC NXTIN(IARC)=FIN(IT) FIN(IT)=IARC IF (RC(IARC).EQ.0) THEN TNXTOU(IARC)=TFSTOU(IH) TFSTOU(IH)=IARC TNXTIN(IARC)=TFSTIN(IT) TFSTIN(IT)=IARC END IF END IF GO TO 20 END C C SUBROUTINE ADDTR(ARC) IMPLICIT INTEGER (A-Z) C C--------------------------------------------------------------- C C PURPOSE - THIS SUBROUTINE CHECKS IF ARC IS IN THE ARRAYS TFSTOU, TNXTOU, C TFSTIN, TNXTIN AND, IF NOT, ADDS ARC TO THEM. C C--------------------------------------------------------------- C C MAXNN = DIMENSION OF NODE-LENGTH ARRAYS C MAXNA = DIMENSION OF ARC-LENGTH ARRAYS C PARAMETER (MAXNN=10000, MAXNA=70000) C COMMON/ARRAYS/STARTN/ARRAYE/ENDN COMMON/BLK10/TFSTOU/BLK11/TNXTOU/BLK12/TFSTIN/BLK13/TNXTIN INTEGER STARTN(MAXNA),ENDN(MAXNA) INTEGER TFSTOU(MAXNN),TNXTOU(MAXNA),TFSTIN(MAXNN),TNXTIN(MAXNA) C NODE=STARTN(ARC) ARC1=TFSTOU(NODE) 10 IF (ARC1.GT.0) THEN IF (ARC1.EQ.ARC) GO TO 20 ARC1=TNXTOU(ARC1) GO TO 10 END IF TNXTOU(ARC)=TFSTOU(NODE) TFSTOU(NODE)=ARC 20 NODE=ENDN(ARC) ARC1=TFSTIN(NODE) 30 IF (ARC1.GT.0) THEN IF (ARC1.EQ.ARC) RETURN ARC1=TNXTIN(ARC1) GO TO 30 END IF TNXTIN(ARC)=TFSTIN(NODE) TFSTIN(NODE)=ARC RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SAMPLE INPUT FILE RELAX4.INP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 20 80 1 10 38 115 1 11 17 120 8 18 82 115 8 1 56 197 8 11 48 126 10 8 69 115 10 16 39 115 10 15 78 179 10 14 82 131 10 12 100 146 10 9 74 101 2 15 27 166 2 1 75 195 2 16 57 129 2 11 38 181 2 8 67 185 2 14 39 162 12 18 1 166 12 14 99 126 12 10 85 155 12 2 15 157 12 11 59 148 15 12 36 166 15 16 18 166 15 13 100 114 15 4 1 130 15 7 96 100 15 18 30 155 3 6 71 100 3 5 12 153 6 7 68 100 6 17 89 100 6 12 13 163 7 9 93 100 7 13 26 176 7 16 24 148 7 20 45 169 9 20 2 100 9 16 58 100 9 1 77 180 9 19 13 145 9 6 63 104 4 13 81 100 4 1 61 160 4 12 4 133 4 3 23 136 4 14 77 188 13 19 95 100 13 17 43 100 13 11 82 104 5 11 31 119 5 8 18 103 5 7 84 173 5 20 95 119 5 6 27 175 11 14 62 119 11 15 15 194 11 10 83 185 11 4 52 184 11 7 94 165 14 19 87 119 14 16 92 119 14 12 52 128 16 14 80 145 16 12 89 148 16 17 83 173 17 19 33 189 17 4 20 132 17 14 15 136 18 17 92 114 18 2 3 157 19 16 77 146 19 14 48 166 20 9 69 160 20 13 9 154 20 18 32 157 20 17 80 165 20 16 76 185 20 8 54 117 20 7 3 131 115 166 5 95 119 0 0 0 0 0 0 0 0 0 0 -183 -73 -161 -80 -3