C THIS FORTRAN CODE IS FOR SOLVING MINIMUM COST NETWORK FLOW C PROBLEM WITH POSITIVE ARC GAINS, SEPARABLE CONVEX DIFFERENTIABLE C ARC COSTS, AND BOUNDS. C C THE COST FUNCTION VALUE, ITS DERIVATIVE, AND INVERSE OF C ITS DERIVATIVE FOR EACH ARC I WITH FLOW F C ARE ASSUMED TO HAVE FORMULAS, SPECIFIED THROUGH THE C FUNCTIONS C(I,F), DC(I,F), DCINV(I,T), C AND PARAMETERIZED BY DATA VALUES C A(I) AND CDAT(I,J), J=1,...,NDAT, C STORED IN AN ASCII DATA FILE NL013.DAT. THESE ARE C SUPPLIED BY THE USER. C C THIS CODE IMPLEMENTS THE EPSILON-RELAXATION METHOD WITH C THE AUCTION/SEQUENTIAL-SHORTEST-PATH VARIATION C (I.E., FLOW PUSH ALONG PATHS OR CYCLES WITHOUT CHANGING THE SIGN C OF SURPLUSES; PRICE RISE AT A NODE). THE METHOD AND NUMERICAL C TEST RESULTS ARE DESCRIBED IN THE PAPER: C Francesca Guerriero and Paul Tseng, C "Implementation and Testing of Auction Methods for Solving C Separable Convex Cost Generalized Network Flow Problems", C J. Optim. Theory Appl., 115, 2002, 113-144. C C FOR ANY QUESTIONS, CONTACT PAUL TSENG (tseng@math.washington.edu). C THIS CODE IS PLACED IN THE PUBLIC DOMAIN TO BE USED FOR C RESEARCH PURPOSES. FOR ANY OTHER PURPOSES, PLEASE CONTACT C PAUL TSENG. C C THIS CODE DOES NOT HAVE A BUILT-IN SAFEGUARD TO ENSURE FINITE TERMINATION. C ONE SUCH SAFEGUARD WOULD BE TO COUNT THE NUMBER OF FLOW PUSHES ALONG C CYCLES WITH GAIN GREATER THAN 1. IF THIS NUMBER OF EXCEEDS A THRESHOLD, C THEN SWITCH TO THE TWO PHASE SCHEME: DO UP ITERATIONS ONLY UNTIL NO C NODE WITH SURPLUS EXCEEDING THRESHOLD REMAINS; THEN DO DOWN ITERATIONS C ONLY. HOWEVER, IN PRACTICE THIS NUMBER SEEMS TO REMAIN VERY SMALL C SO WE HAVE NOT BOTHERED TO BUILD THE SAFEGUARD. C C THIS VERSION TREATS NONLINEAR ARCS DIFFERENTLY. C RATHER THAN KEEPING ALL PUSH LIST ARCS IN THE SAME QUEUE C THIS VERSION MAINTAINS SEPERATE LISTS FOR NONLINEAR AND C LINEAR COST ARCS. FLOW IS FIRST PUSHED ON LINEAR ARCS C THEN ON NONLINEAR. C C VERSION 1: March 9, 2001. C VERSION 1.1: April 22, 2005. C * Array dimensions in subroutines are declared explicitly. C * A bug in initializing THRESHSURPLUS (which can be initialized to 0 C and thus cause infinite cycling when node supplies/demands are small) C is fixed. C *************************************************** C IMPLICIT INTEGER*4 (A-Z) PARAMETER (MAXNODES=5000, MAXARCS=20000, MAXNDAT=5) INTEGER*4 STARTN(MAXARCS) INTEGER*4 ENDN(MAXARCS) INTEGER*4 FIN(MAXNODES) INTEGER*4 FOUT(MAXNODES) INTEGER*4 NXTIN(MAXARCS) INTEGER*4 NXTOU(MAXARCS) INTEGER*4 FPUSHF(MAXNODES) INTEGER*4 NXTPUSHF(MAXARCS) INTEGER*4 FPUSHB(MAXNODES) INTEGER*4 NXTPUSHB(MAXARCS) INTEGER*4 QFPUSHF(MAXNODES),QNXTPUSHF(MAXARCS) INTEGER*4 QFPUSHB(MAXNODES),QNXTPUSHB(MAXARCS) INTEGER*4 NXTQUEUE(MAXNODES) INTEGER*4 PRDCSR(MAXNODES) REAL*8 FACTOR,EPS,THRESHSURPLUS,SFACTOR,MAXSURPLUS,LARGE,ZERO REAL*8 A(MAXARCS),CDAT(MAXARCS,MAXNDAT) REAL*8 U(MAXARCS) REAL*8 G(MAXARCS) REAL*8 COST(MAXARCS) REAL*8 F(MAXARCS) REAL*8 SURPLUS(MAXNODES) REAL*8 P(MAXNODES) COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /STATS/ NCYC,NITER,TOTALITER COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN 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 /QBU1/ QFPUSHF COMMON /QBU2/ QNXTPUSHF COMMON /QBU3/ QFPUSHB COMMON /QBU4/ QNXTPUSHB COMMON /QUEUE/ NXTQUEUE COMMON /BLK8/ PRDCSR COMMON /SC2/ FACTOR,EPS,THRESHSURPLUS,SFACTOR,MAXSURPLUS COMMON /LCOS/ A COMMON /QCO1/ CDAT COMMON /UBOUND/ U COMMON /QG/ G COMMON /COST/ COST COMMON /FLOW/ F COMMON /BLK3/ SURPLUS COMMON /PRICES/ P C EPS3 IS NEEDED ONLY IF DCINV IS LINKED FROM INEXACT INVERSE DCINV.F REAL*8 EPS3 COMMON /EPS3/EPS3 real*8 REF_THR,C,DC,DCINV REAL*8 MAXCOST,GAIN,TEMP REAL*8 TT,TIMER DOUBLE PRECISION TCOST,PRODUCT real*8 SMAL REAL*8 MAXVIOL real etime, tarray1(2) external etime LARGE=2000000000 ZERO=0.0 SMAL = 0.0 C Print*,'DO you want to have ill conditioning?(1)' C read(*,*),dummy C if (dummy .eq.1) then C print*, 'Will substitute zero quad. coeff with C $ the small coefficient you will enter now' C READ(*,*) SMAL C PRINT*, SMAL C end if PRINT*,'E-RELAXATION ASSP METHOD FOR MIN COST FLOW' PRINT*,'*************************************' PRINT *,'READING PROBLEM DATA' OPEN(13,FILE='NL013.DAT',STATUS='OLD') REWIND(13) READ(13,1010) N,NA,NDAT 1010 FORMAT(2I8,I3) PRINT*,' #NODES=',N,' #ARCS=',NA,' #COST PARAMS=',NDAT C READ START, END, COST PARAMETERS, CAPACITY, GAIN OF EACH ARC C AND GENERATE MAX COST VALUE DO 5 I=1,NA READ(13,1020) STARTN(I),ENDN(I),A(I),(CDAT(I,J),J=1,NDAT), $ U(I),G(I) 1020 FORMAT(2I8,8E15.9) if ((CDAT(I,1) .eq.0).and.(smal.gt.1.0e-9)) then CDAT(I,1) = SMAL END IF 5 CONTINUE C READ SUPPLY OF EACH NODE MAXSURPLUS=0 DO 8 I=1,N READ(13,1000) SURPLUS(I) MAXSURPLUS = MAX( MAXSURPLUS , DABS(SURPLUS(I)) ) 8 CONTINUE REWIND(13) CLOSE(13) PRINT*,'END OF READING' 1000 FORMAT(E15.9) MAXCOST=0.0 K=NA DO 9 I=1,K C THIS, USED IN PREVIOUS ASSP CODES, DIFFERS FROM ONE BELOW USED IN E-RELAX CODES. C TEMP = U(I) TEMP = MIN( MAXSURPLUS , U(I) ) C IF DERIVATIVE AT TEMP IS INFINITE, REDUCE TEMP IF (DC(I,TEMP).EQ.LARGE) THEN TEMP=TEMP/2. END IF MAXCOST=MAX( MAXCOST , DABS(DC(I,TEMP)) ) C CHECK FOR SELF-LOOP AND ADD AN ARC FOR EACH SUCH LOOP IF (STARTN(I).EQ.ENDN(I)) THEN N=N+1 SURPLUS(N)=0.0 ENDN(I)=N NA=NA+1 STARTN(NA)=N ENDN(NA)=STARTN(I) A(NA)=0.0 DO 11 J=1,NDAT CDAT(NA,J)=ZERO 11 CONTINUE U(NA)=U(I)*G(I) G(NA)=1.0 END IF 9 CONTINUE IF (NA.GT.K) PRINT*,'ADD ',NA-K,' ARCS FOR SELF-LOOPS' C INPUT NUMBER OF DIGITS OF ACCURACY DESIRED IN SOLUTION C PRINT*,' INPUT #DIGITS OF ACCURACY FOR SOLUTION' READ(*,*)ACCUR C INITIALIZE PRICES & PUSH LISTS DO 10 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 QFPUSHF(NODE)=0 QFPUSHB(NODE)=0 C P(NODE)=-(LARGE/10) P(NODE)=0.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 c WRITE(*,*)'ENTER STARTING EPSILON' c WRITE(*,*)'SHOULD BE BETWEEN 1 & ',MAXCOST c READ*,EPS c IF (EPS.LT.1) GO TO 25 EPS = MAXCOST/3.0 IF (EPS.LT.(1.0/(N+1.0))) EPS = 1.0/(N+1.0) C PRINT*,'ENTER MAXIMUM EPSILON (TYPICALLY A LARGE NO.)' C READ*,TEMP C EPS = MIN(EPS,TEMP) FACTOR = 3.0 WRITE(*,*) 'EPS, FACTOR:',EPS,FACTOR c WRITE(*,*)'ENTER SCALING FACTOR' c30 CONTINUE c READ*,FACTOR c IF ((EPS.GT.1).AND.(FACTOR.LE.1)) THEN c PRINT*,'ENTER SCALING FACTOR; SHOULD BE GREATER THAN 1' c GO TO 30 c END IF C MAXSURPLUS=1.0 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 REF_THR = 1.0E-10 SSCALE=1 IF (SSCALE.EQ.1) THEN SFACTOR=2.0+INT(FACTOR*MAXSURPLUS/MAXCOST) IF (SFACTOR.LT.4.0) SFACTOR=4.0 THRESHSURPLUS=MAX(INT(MAXSURPLUS/SFACTOR),REF_THR) C IF (EPS.LE.(1.0/(N+1.0))) THRESHSURPLUS=REF_THR ELSE THRESHSURPLUS=REF_THR END IF PRINT*,'INITIALIZING DATA STRUCTURES' CALL INIDAT PRINT *,'**********************************' PRINT *,'CALLING ASSP ALGORITHM FOR MCF WITH GAINS' TIMER = etime(tarray1) TIMER = TARRAY1(1) CALL EPS_RELAX(ACCUR) TT = etime(tarray1) TT = TARRAY1(1)-TIMER PRINT*, 'TOTAL TIME = ',TT, ' secs.' C open(17,FILE='FORT9',STATUS='unknown') C rewind(17) TCOST=0 DO 330 I=1,NA TCOST=TCOST+C(I,F(I)) 330 CONTINUE PRINT *,'**********************************' WRITE(*,*) 'OPTIMAL COST = ', TCOST PRINT *,'# OF ITERATIONS = ',TOTALITER PRINT *,'# OF S. N. PRICE RISES/DROPS = ',NITER PRINT *,'**********************************' C CHECK CORRECTNESS OF THE ANSWER DO 80 NODE=1,N IF (ABS(SURPLUS(NODE)).GT.MAXSURPLUS*(.1**ACCUR)) THEN PRINT*,'NONZERO SURPLUS AT NODE ',NODE, SURPLUS(NODE) ENDIF 80 CONTINUE NUM_POS_VI = 0 NUM_NEG_VI = 0 MAXVIOL = 0.0 DO 90 ARC=1,NA COST(ARC)=DC(ARC,F(ARC)) GAIN=G(ARC) IF (F(ARC).GT.REF_THR) THEN IF (P(STARTN(ARC))-GAIN*P(ENDN(ARC)).LT.-EPS+COST(ARC) $ -REF_THR) THEN PRINT*,'E-CS VIOLATED AT ARC ',ARC, $ P(STARTN(ARC))-GAIN*P(ENDN(ARC))+EPS-COST(ARC) NUM_POS_VI = NUM_POS_VI + 1 MAXVIOL = MAX(MAXVIOL, $ ABS(P(STARTN(ARC))-GAIN*P(ENDN(ARC))+EPS-COST(ARC))) ENDIF ENDIF IF (F(ARC).LT.(U(ARC)-REF_THR)) THEN IF (P(STARTN(ARC))-GAIN*P(ENDN(ARC)).GT.EPS+COST(ARC) $ +REF_THR) THEN c PRINT*,'E-CS VIOLATED AT ARC ',ARC, c $ P(STARTN(ARC))-GAIN*P(ENDN(ARC))-EPS-COST(ARC) NUM_NEG_VI = NUM_NEG_VI + 1 MAXVIOL = MAX(MAXVIOL, $ ABS(P(STARTN(ARC))-GAIN*P(ENDN(ARC))-EPS-COST(ARC))) ENDIF ENDIF 90 CONTINUE PRINT *, ' ' print*, 'NUMBER OF E-CS VIOLATIONS: ',NUM_POS_VI,' ',NUM_NEG_VI PRINT*, 'MAXIMUM VIOLATION: ', MAXVIOL C PRINT *, 'PROGRAM ENDED; PRESS ' C 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=5000, MAXARCS=20000) IMPLICIT INTEGER*4 (A-Z) INTEGER*4 STARTN(MAXARCS) INTEGER*4 ENDN(MAXARCS) INTEGER*4 FIN(MAXNODES) INTEGER*4 FOUT(MAXNODES) INTEGER*4 NXTIN(MAXARCS) INTEGER*4 NXTOU(MAXARCS) INTEGER*4 FINALIN(MAXNODES),FINALOU(MAXNODES) REAL*8 FACTOR,EPS,THRESHSURPLUS,SFACTOR,MAXSURPLUS,LARGE,ZERO COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /SC2/ FACTOR,EPS,THRESHSURPLUS,SFACTOR,MAXSURPLUS COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK4/ FIN COMMON /BLK5/ FOUT COMMON /BLK6/ NXTIN COMMON /BLK7/ NXTOU 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(ACCUR) IMPLICIT NONE INTEGER*4 ACCUR INTEGER*4 N,NA,MAXNODES,MAXARCS,MAXNDAT INTEGER*4 NCYC,NITER,TOTALITER PARAMETER (MAXNODES=5000, MAXARCS=20000, MAXNDAT=5) INTEGER*4 STARTN(MAXARCS),ENDN(MAXARCS),FIN(MAXNODES), $ FOUT(MAXNODES) INTEGER*4 NXTIN(MAXARCS),NXTOU(MAXARCS) INTEGER*4 FPUSHF(MAXNODES),NXTPUSHF(MAXARCS),FPUSHB(MAXNODES), $ NXTPUSHB(MAXARCS) INTEGER*4 QFPUSHF(MAXNODES),QNXTPUSHF(MAXARCS) INTEGER*4 QFPUSHB(MAXNODES),QNXTPUSHB(MAXARCS) INTEGER*4 NXTQUEUE(MAXNODES) INTEGER*4 PRDCSR(MAXNODES) REAL*8 FACTOR,EPS,THRESHSURPLUS,SFACTOR,MAXSURPLUS,LARGE,ZERO REAL*8 A(MAXARCS),CDAT(MAXARCS,MAXNDAT),U(MAXARCS),G(MAXARCS) REAL*8 COST(MAXARCS),F(MAXARCS),SURPLUS(MAXNODES),P(MAXNODES) COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /STATS/ NCYC,NITER,TOTALITER COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN 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 /QBU1/ QFPUSHF COMMON /QBU2/ QNXTPUSHF COMMON /QBU3/ QFPUSHB COMMON /QBU4/ QNXTPUSHB COMMON /QUEUE/ NXTQUEUE COMMON /BLK8/ PRDCSR COMMON /SC2/ FACTOR,EPS,THRESHSURPLUS,SFACTOR,MAXSURPLUS COMMON /LCOS/ A COMMON /QCO1/ CDAT COMMON /UBOUND/ U COMMON /QG/ G COMMON /COST/ COST COMMON /FLOW/ F COMMON /BLK3/ SURPLUS COMMON /PRICES/ P REAL*8 EPS3 COMMON /EPS3/EPS3 INTEGER*4 NODE,NXTNODE INTEGER*4 ARC,ARCF,ARCB INTEGER*4 QARCF,QARCB INTEGER*4 I,J,K,LN,PREVNODE,LASTQUEUE INTEGER*4 START,END REAL*8 REFPRICE,XP,PRICE,PRICEINCR REAL*8 PSTART,PEND,DEL,MAXPRICE REAL*8 FLOW,RESID,UP,CAPOUT,CAPIN REAL*8 REF_THR,C,DC,DCINV,EPS2,GAIN C ADDITIONAL DATA STRUCTURE USED TO COMPUTE DUAL COST REAL*8 PCOST,DCOST,RC,F1,DFCT(5000),DFCT1(5000) DO 1 I=1,N DFCT1(I)=-SURPLUS(I) PRDCSR(I)=0 1 CONTINUE C INITIALIZE COUNT OF NUMBER OF UP/DOWN ITERATIONS PERFORMED NITER = 0 TOTALITER=0 NCYC=0 MAXPRICE = LARGE REF_THR = 1.0E-10 C REDUCE ARC CAPACITIES DO 40 NODE=1,N CAPOUT=0.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.0) GOTO 400 CAPIN=0.0 ARC=FIN(NODE) 43 IF (ARC.GT.0) THEN IF (G(ARC)*U(ARC).GT.CAPOUT) THEN U(ARC)=CAPOUT/G(ARC) ENDIF CAPIN=MIN(LARGE,CAPIN+G(ARC)*U(ARC)) ARC=NXTIN(ARC) GO TO 43 END IF CAPIN=MIN(LARGE,CAPIN+SURPLUS(NODE)) IF (CAPIN.LT.0.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 EPS2 = EPS/2.0 EPS3=EPS2/2.0 DO 49 ARC=1,NA START=STARTN(ARC) END=ENDN(ARC) PSTART=P(START) PEND=P(END) GAIN=G(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN IF (PSTART.GE.GAIN*PEND+A(ARC)+EPS2) THEN UP=U(ARC) ELSE UP=0.0 END IF COST(ARC)=A(ARC) ELSE C THIS LINE IS NEEDED ONLY IF DCINV IS LINKED FROM DCINV.F F(ARC)=0.5*U(ARC) COST(ARC)=DC(ARC,F(ARC)) UP=MIN(U(ARC),DCINV(ARC,PSTART-GAIN*PEND)) UP=MAX(0.0,UP) COST(ARC)=DC(ARC,UP) END IF SURPLUS(START)=SURPLUS(START)-UP SURPLUS(END)=SURPLUS(END)+GAIN*UP F(ARC)=UP 49 CONTINUE C PRINT*, 'STARTING NEW SCALING PHASE' C *********** START OF A NEW SCALING PHASE *********** 60 CONTINUE EPS2 = EPS/2.0 EPS3=EPS2/2.0 DO 81 NODE=1,N PRDCSR(NODE)=0 81 CONTINUE C ***** QUEUE INITIALIZATION ***** DO 82 I=1,N-1 NXTQUEUE(I)=I+1 82 CONTINUE NXTQUEUE(N)=1 NODE=1 PREVNODE=N LASTQUEUE=N 100 CONTINUE C TAKE UP NEXT NODE FOR ITERATION IF (SURPLUS(NODE) .GT. THRESHSURPLUS) THEN TOTALITER=TOTALITER+1 PRDCSR(NODE) = NA+1 I = NODE C ***** UP ITERATION ON NODE I **** C ARCF & ARCB ARE THE CURRENT VALUES OF THE STARTING ARCS OF THE C PUSH LISTS OF I 110 ARCF = FPUSHF(I) ARCB = FPUSHB(I) QARCF = QFPUSHF(I) QARCB = QFPUSHB(I) PRICE = P(I) C PRINT*,'I=',I,' PRICE=',PRICE,' surpi=',surplus(i) c if ((i.eq.10) .and. (surplus(i) .gt. 0)) then C print*,'0',arcf,arcb,qarcf,qarcb c end if C ************* TAKING CARE OF LINEAR FORWARD ARCS*************************** 120 IF (ARCF .GT. 0) THEN J = ENDN(ARCF) GAIN = G(ARCF) IF (PRICE-GAIN*P(J)-COST(ARCF).GT.EPS2) THEN IF (U(ARCF)-F(ARCF).GT.REF_THR) THEN FPUSHF(I) = ARCF IF (-SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH1(J,ARCF) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J) = ARCF I = J GOTO 110 ELSE CALL CYCLE_PUSH1(I,J,ARCF) GOTO 620 END IF END IF ELSE ARCF = NXTPUSHF(ARCF) END IF ELSE ARCF = NXTPUSHF(ARCF) END IF GOTO 120 END IF FPUSHF(I) = ARCF C *********** TAKING CARE OF THE LINEAR COST REVERSE ARCS *************** 121 IF (ARCB .GT. 0) THEN J = STARTN(ARCB) GAIN = G(ARCB) IF (GAIN*PRICE-P(J)+COST(ARCB).GT.EPS2) THEN IF (F(ARCB).GT.REF_THR) THEN FPUSHB(I) = ARCB IF (-SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH1(J,-ARCB) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J) = -ARCB I = J GOTO 110 ELSE CALL CYCLE_PUSH1(I,J,-ARCB) GOTO 620 END IF END IF ELSE ARCB = NXTPUSHB(ARCB) END IF ELSE ARCB = NXTPUSHB(ARCB) END IF GOTO 121 END IF FPUSHB(I) = ARCB C DO A DIAGNOSTIC CHECK C DO 609 ARC=1,NA C IF (F(ARC).GT.REF_THR) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).LT.COST(ARC)-1.1*EPS) C $THEN C PRINT*,'1. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS,F(ARC) C stop C ENDIF C ENDIF C IF (F(ARC)+REF_THR.LT.U(ARC)) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).GT.COST(ARC)+1.1*EPS) C $THEN C PRINT*,'2. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS, C $ U(ARC)-F(ARC) C stop C ENDIF C ENDIF C IF ( (F(ARC).LT.-REF_THR) .OR. (F(ARC)-REF_THR.GT.U(ARC)) ) THEN C PRINT*,'FLOW OUTSIDE BDS AT',ARC,F(ARC),U(ARC) C STOP C END IF C IF ( DABS(DC(ARC,F(ARC))-COST(ARC)) .GT. REF_THR) THEN C PRINT*,'COST ERR AT',ARC,DC(ARC,F(ARC)),COST(ARC) C STOP C END IF C609 CONTINUE c if ((i.eq.13) .and. (surplus(i) .gt. 0)) then c print*,'2',qarcf,qarcb,surplus(i),eps2 c end if C**************TAKING CARE OF NONLINEAR FORWARD ARCS********************* 122 IF (QARCF .GT. 0) THEN J = ENDN(QARCF) GAIN = G(QARCF) IF (PRICE-GAIN*P(J)-COST(QARCF).GT.EPS2) THEN RESID=MIN(U(QARCF),DCINV(QARCF,PRICE-GAIN*P(J))) $ - F(QARCF) IF (RESID.GT.REF_THR) THEN QFPUSHF(I) = QARCF IF (-SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH1(J,QARCF) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J) = QARCF I = J GOTO 110 ELSE CALL CYCLE_PUSH1(I,J,QARCF) GOTO 620 END IF END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF GOTO 122 END IF QFPUSHF(I) = QARCF C ************* TAKING CARE OF NONLINEAR BACKWARD ARCS****************** 123 IF (QARCB .GT. 0) THEN J = STARTN(QARCB) GAIN = G(QARCB) IF (GAIN*PRICE-P(J)+COST(QARCB).GT.EPS2) THEN RESID=F(QARCB)-MAX(0.0,DCINV(QARCB,P(J)-GAIN*PRICE)) IF (RESID.GT.REF_THR) THEN QFPUSHB(I) = QARCB IF (-SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH1(J,-QARCB) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=-QARCB I = J GOTO 110 ELSE CALL CYCLE_PUSH1(I,J,-QARCB) GOTO 620 END IF END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF GOTO 123 END IF QFPUSHB(I) = QARCB C DO A DIAGNOSTIC CHECK C DO 608 ARC=1,NA C IF (F(ARC).GT.REF_THR) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).LT.COST(ARC)-1.1*EPS) C $THEN C PRINT*,'1. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS,F(ARC) C stop C ENDIF C ENDIF C IF (F(ARC)+REF_THR.LT.U(ARC)) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).GT.COST(ARC)+1.1*EPS) C $THEN C PRINT*,'2. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS, C $ U(ARC)-F(ARC) C stop C ENDIF C ENDIF C IF ( (F(ARC).LT.-REF_THR) .OR. (F(ARC)-REF_THR.GT.U(ARC)) ) THEN C PRINT*,'FLOW OUTSIDE BDS AT',ARC,F(ARC),U(ARC) C STOP C END IF C IF ( DABS(DC(ARC,F(ARC))-COST(ARC)) .GT. REF_THR) THEN C PRINT*,'COST ERR AT',ARC,DC(ARC,F(ARC)),COST(ARC) C STOP C END IF C608 CONTINUE C ***** END OF D-PUSHES; CHECK IF PRICE RISE IS NEEDED ***** c if ((i.eq.10) .and. (surplus(i) .gt. 0)) then C print*,'3',i,arcf,arcb,qarcf,qarcb,price c end if IF ((ARCF+ARCB+QARCF+QARCB) .EQ. 0) THEN C ******************DO A PRICE RISE ********************* REFPRICE=PRICE PRICE = LARGE NITER=NITER+1 ARC = FOUT(I) 111 IF (ARC .GT. 0) THEN UP = U(ARC) IF ((F(ARC)-UP).LT.(-REF_THR)) THEN XP = G(ARC)*P(ENDN(ARC))+COST(ARC) IF ((XP-PRICE).LT.-EPS) THEN PRICE = XP + EPS IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR FORWARD PUSH LIST*********************** QARCF = ARC QNXTPUSHF(ARC)=0 ARCF = 0 ELSE C*********** LINEAR FORWARD PUSH LIST*********************** ARCF = ARC NXTPUSHF(ARC)=0 QARCF = 0 END IF ELSE IF ((XP-PRICE).LT.-EPS2) THEN IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR FORWARD PUSH LIST*********************** IF (QARCF.EQ.0) THEN QARCF = ARC QNXTPUSHF(ARC)=0 ELSE QNXTPUSHF(ARC) = QARCF QARCF = ARC END IF ELSE C*********** LINEAR FORWARD PUSH LIST*********************** IF (ARCF.EQ.0) THEN ARCF = ARC NXTPUSHF(ARC)=0 ELSE NXTPUSHF(ARC) = ARCF ARCF = ARC END IF END IF END IF END IF END IF ARC = NXTOU(ARC) GOTO 111 END IF ARC = FIN(I) 112 IF (ARC .GT. 0) THEN IF (F(ARC) .GT. REF_THR) THEN XP = (P(STARTN(ARC))-COST(ARC))/G(ARC) IF ((XP-PRICE)*G(ARC) .LT. -EPS) THEN PRICE = XP + EPS/G(ARC) ARCF=0 QARCF = 0 IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR REVERSE PUSH LIST*********************** QARCB = ARC QNXTPUSHB(ARC)=0 ARCB = 0 ELSE C*********** LINEAR REVERSE PUSH LIST*********************** ARCB = ARC NXTPUSHB(ARC)=0 QARCB = 0 END IF ELSE IF ((XP-PRICE)*G(ARC) .LT. -EPS2) THEN IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR FORWARD PUSH LIST*********************** IF (QARCB.EQ.0) THEN QARCB = ARC QNXTPUSHB(ARC)=0 ELSE QNXTPUSHB(ARC) = QARCB QARCB = ARC END IF ELSE C*********** LINEAR FORWARD PUSH LIST*********************** IF (ARCB.EQ.0) THEN ARCB = ARC NXTPUSHB(ARC)=0 ELSE NXTPUSHB(ARC) = ARCB ARCB = ARC END IF END IF END IF END IF END IF ARC = NXTIN(ARC) GOTO 112 END IF 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.GE.LARGE) THEN ARCF=0 ARCB=0 QARCF=0 QARCB=0 ARC = FOUT(I) 211 IF (ARC.GT.0) THEN XP = G(ARC)*P(ENDN(ARC))+COST(ARC) IF ((XP-PRICE).LT. (-REF_THR)) THEN PRICE = XP ELSE IF (XP.GE.LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' STOP END IF END IF ARC = NXTOU(ARC) GOTO 211 END IF ARC = FIN(I) 212 IF (ARC.GT.0) THEN XP = P(STARTN(ARC))-COST(ARC) IF ((XP-G(ARC)*PRICE).LT.(-REF_THR)) THEN PRICE = XP/G(ARC) ELSE IF (XP.GE.LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' STOP END IF END IF ARC = NXTIN(ARC) GOTO 212 END IF IF (PRICE.LT.INT(LARGE/10)) PRICE=INT(LARGE/10) END IF C SET PRICE OF NODE I P(I) = PRICE FPUSHF(I) = ARCF FPUSHB(I) = ARCB QFPUSHF(I) = QARCF QFPUSHB(I) = QARCB C DO A DIAGNOSTIC CHECK C PRICEINCR=PRICE-REFPRICE C IF (PRICEINCR.LE.-REF_THR) THEN C PRINT*,'NONPOSITIVE PRICE INCREMENT=',PRICEINCR C PAUSE C STOP C ENDIF C DO 607 ARC=1,NA C IF (F(ARC).GT.REF_THR) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).LT.COST(ARC)-1.1*EPS) C $THEN C PRINT*,'1. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS,F(ARC) C stop C ENDIF C ENDIF C IF (F(ARC)+REF_THR.LT.U(ARC)) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).GT.COST(ARC)+1.1*EPS) C $THEN C PRINT*,'2. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS, C $ U(ARC)-F(ARC) C stop C ENDIF C ENDIF C IF ( (F(ARC).LT.-REF_THR) .OR. (F(ARC)-REF_THR.GT.U(ARC)) ) THEN C PRINT*,'FLOW OUTSIDE BDS AT',ARC,F(ARC),U(ARC) C STOP C END IF C IF ( DABS(DC(ARC,F(ARC))-COST(ARC)) .GT. REF_THR) THEN C PRINT*,'COST ERR AT',ARC,DC(ARC,F(ARC)),COST(ARC) C STOP C END IF C607 CONTINUE C END OF PRICE RISE. ARC = PRDCSR(I) IF (ARC.NE.NA+1) THEN PRDCSR(I) = 0 IF (ARC.GT.0) THEN I = STARTN(ARC) ELSE I = ENDN(-ARC) END IF END IF GOTO 110 END IF C IF THE PRICE IS TOO HIGH, THEN SOMETHING IS WRONG C IF (PRICE.GT.MAXPRICE) THEN C GO TO 400 C ENDIF ELSE IF (-SURPLUS(NODE) .GT. THRESHSURPLUS) THEN TOTALITER=TOTALITER+1 PRDCSR(NODE) = NA+1 I = NODE C ***** DOWN ITERATION ON NODE I ***** C ARCF & ARCB ARE THE CURRENT VALUES OF THE STARTING ARCS OF THE C PUSH LISTS OF I 510 ARCF = FPUSHF(I) ARCB = FPUSHB(I) QARCF = QFPUSHF(I) QARCB = QFPUSHB(I) PRICE = P(I) C PRINT*,'I=',I,' PRICE=',PRICE,' surpi=',surplus(i) c if ((i.eq.14) .and. (-surplus(i) .gt. 0)) then c print*,'0',arcf,arcb,qarcf,qarcb c end if C******************TAKING CARE OF LINEAR FORWARD ARCS*************************** 520 IF (ARCF .GT. 0) THEN J = ENDN(ARCF) GAIN = G(ARCF) IF (-(PRICE-GAIN*P(J)-COST(ARCF)).GT.EPS2) THEN IF (F(ARCF).GT.REF_THR) THEN FPUSHF(I) = ARCF IF (SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH2(J,ARCF) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=ARCF I = J GOTO 510 ELSE CALL CYCLE_PUSH2(I,J,ARCF) GOTO 620 END IF END IF ELSE ARCF = NXTPUSHF(ARCF) END IF ELSE ARCF = NXTPUSHF(ARCF) END IF GOTO 520 END IF FPUSHF(I) = ARCF C *********** TAKING CARE OF THE LINEAR COST REVERSE ARCS *************** 521 IF (ARCB .GT. 0) THEN J = STARTN(ARCB) GAIN = G(ARCB) IF (-(GAIN*PRICE-P(J)+COST(ARCB)).GT.EPS2) THEN IF (U(ARCB)-F(ARCB).GT.REF_THR) THEN FPUSHB(I) = ARCB IF (SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH2(J,-ARCB) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=-ARCB I = J GOTO 510 ELSE CALL CYCLE_PUSH2(I,J,-ARCB) GOTO 620 END IF END IF ELSE ARCB = NXTPUSHB(ARCB) END IF ELSE ARCB = NXTPUSHB(ARCB) END IF GOTO 521 END IF FPUSHB(I) = ARCB C DFCT(I)=DFCT1(I) C k=fout(i) C423 if (k.gt.0) then C dfct(i)=dfct(i)+f(k) C k=nxtou(k) C go to 423 C end if C k=fin(i) C424 if (k.gt.0) then C dfct(i)=dfct(i)-g(k)*f(k) C k=nxtin(k) C go to 424 C end if C IF ( DABS(DFCT(I)+SURPLUS(I)) .GT. 1.d0 ) THEN C PRINT*,'SURPLUS ERR AT',I,-DFCT(I),SURPLUS(I) C END IF c if ((i.eq.14) .and. (-surplus(i) .gt. 0)) then c print*,'2',qarcf,qarcb,surplus(i),eps2 c end if C **************TAKING CARE OF NONLINEAR FORWARD ARCS************************ 522 IF (QARCF .GT. 0) THEN J = ENDN(QARCF) GAIN = G(QARCF) IF (-(PRICE-GAIN*P(J)-COST(QARCF)).GT.EPS2) THEN RESID=F(QARCF)-MAX(0.0,DCINV(QARCF,PRICE-GAIN*P(J))) IF (RESID.GT.REF_THR) THEN QFPUSHF(I) = QARCF IF (SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH2(J,QARCF) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=QARCF I = J GOTO 510 ELSE CALL CYCLE_PUSH2(I,J,QARCF) GOTO 620 END IF END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF GOTO 522 END IF QFPUSHF(I) = QARCF C ************* TAKING CARE OF NONLINEAR BACKWARD ARCS****************** 523 IF (QARCB .GT. 0) THEN J = STARTN(QARCB) GAIN = G(QARCB) IF (-(GAIN*PRICE-P(J)+COST(QARCB)).GT.EPS2) THEN RESID=MIN(U(QARCB),DCINV(QARCB,P(J)-GAIN*PRICE)) $ - F(QARCB) IF (RESID.GT.REF_THR) THEN QFPUSHB(I) = QARCB IF (SURPLUS(J).GT.REF_THR) THEN CALL PATH_PUSH2(J,-QARCB) GOTO 620 ELSE IF (PRDCSR(J).EQ.0) THEN PRDCSR(J) = -QARCB I = J GOTO 510 ELSE CALL CYCLE_PUSH2(I,J,-QARCB) GOTO 620 END IF END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF GOTO 523 END IF QFPUSHB(I) = QARCB C DFCT(I)=DFCT1(I) C k=fout(i) C425 if (k.gt.0) then C dfct(i)=dfct(i)+f(k) C k=nxtou(k) C go to 425 C end if C k=fin(i) C426 if (k.gt.0) then C dfct(i)=dfct(i)-g(k)*f(k) C k=nxtin(k) C go to 426 C end if C IF ( DABS(DFCT(I)+SURPLUS(I)) .GT. 1.d0 ) THEN C PRINT*,'SURPLUS ERR AT',I,-DFCT(I),SURPLUS(I) C stop C END IF C ***** END OF D-PUSHES; CHECK IF PRICE DROP IS NEEDED ***** c if ((i.eq.14) .and. (-surplus(i) .gt. 0)) then c print*,'3',arcf,arcb,qarcf,qarcb,price,surplus(i) c end if IF ((ARCF+ARCB+QARCF+QARCB) .EQ. 0) THEN C ******************DO A PRICE DROP ********************* REFPRICE=PRICE PRICE = -LARGE NITER=NITER+1 ARC = FOUT(I) 511 IF (ARC .GT. 0) THEN IF (F(ARC).GT.REF_THR) THEN XP = G(ARC)*P(ENDN(ARC))+COST(ARC) IF (-(XP-PRICE).LT. -EPS) THEN PRICE = XP - EPS IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR FORWARD PUSH LIST*********************** QARCF = ARC QNXTPUSHF(ARC)=0 ARCF = 0 ELSE C*********** LINEAR FORWARD PUSH LIST*********************** ARCF = ARC NXTPUSHF(ARC)=0 QARCF = 0 END IF ELSE IF (-(XP-PRICE).LT.-EPS2) THEN IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR FORWARD PUSH LIST*********************** IF (QARCF.EQ.0) THEN QARCF = ARC QNXTPUSHF(ARC)=0 ELSE QNXTPUSHF(ARC) = QARCF QARCF = ARC END IF ELSE C*********** LINEAR FORWARD PUSH LIST*********************** IF (ARCF.EQ.0) THEN ARCF = ARC NXTPUSHF(ARC)=0 ELSE NXTPUSHF(ARC) = ARCF ARCF = ARC END IF END IF END IF END IF END IF ARC = NXTOU(ARC) GOTO 511 END IF ARC = FIN(I) 512 IF (ARC .GT. 0) THEN UP = U(ARC) IF (UP-F(ARC) .GT. REF_THR) THEN XP = (P(STARTN(ARC))-COST(ARC))/G(ARC) C THE LINES BELOW DIFFERS FROM LAKIS' EVEN FOR ORDINARY NETWORKS IF (-(XP-PRICE)*G(ARC).LT. -EPS) THEN PRICE = XP - EPS/G(ARC) ARCF=0 QARCF = 0 IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR REVERSE PUSH LIST*********************** QARCB = ARC QNXTPUSHB(ARC)=0 ARCB = 0 ELSE C*********** LINEAR REVERSE PUSH LIST*********************** ARCB = ARC NXTPUSHB(ARC)=0 QARCB = 0 END IF ELSE IF (-(XP-PRICE)*G(ARC).LT. -EPS2) THEN IF (CDAT(ARC,1).NE.ZERO) THEN C*********** NONLINEAR FORWARD PUSH LIST*********************** IF (QARCB.EQ.0) THEN QARCB = ARC QNXTPUSHB(ARC)=0 ELSE QNXTPUSHB(ARC) = QARCB QARCB = ARC END IF ELSE C*********** LINEAR FORWARD PUSH LIST*********************** IF (ARCB.EQ.0) THEN ARCB = ARC NXTPUSHB(ARC)=0 ELSE NXTPUSHB(ARC) = ARCB ARCB = ARC END IF END IF END IF END IF END IF ARC = NXTIN(ARC) GOTO 512 END IF 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.GE.LARGE) THEN ARCF=0 ARCB=0 QARCF=0 QARCB=0 ARC = FOUT(I) 611 IF (ARC.GT.0) THEN XP = G(ARC)*P(ENDN(ARC))+COST(ARC) IF (-(XP-PRICE).LT. (-REF_THR)) THEN PRICE = XP ELSE IF (-XP.GE.LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' STOP END IF END IF ARC = NXTOU(ARC) GOTO 611 ENDIF ARC = FIN(I) 612 IF (ARC.GT.0) THEN XP = P(STARTN(ARC))-COST(ARC) IF (-(XP-G(ARC)*PRICE).LT.(-REF_THR)) THEN PRICE = XP/G(ARC) ELSE IF (-XP.GE.LARGE) THEN PRINT*,'PRICE OF ',I,'HAS EXCEEDED THE INT. RANGE' STOP END IF END IF ARC = NXTIN(ARC) GOTO 612 END IF IF (-PRICE.LT.INT(LARGE/10)) PRICE=-INT(LARGE/10) END IF C SET PRICE OF NODE I P(I) = PRICE FPUSHF(I) = ARCF FPUSHB(I) = ARCB QFPUSHF(I) = QARCF QFPUSHB(I) = QARCB C PRICEINCR=PRICE-REFPRICE C IF (PRICEINCR.GE.REF_THR) THEN C PRINT*,'POSITIVE PRICE INCREMENT=',PRICEINCR,PRICE C STOP C ENDIF C C END OF PRICE RISE; IF THERE IS STILL A LOT OF SURPLUS, C TRY AND DO SOME PUSHING. c print*,'+',i,p(i),surplus(i),resid,flag,arcf,arcb,qarcf,qarcb ARC = PRDCSR(I) IF (ARC.NE.NA+1) THEN PRDCSR(I) = 0 IF (ARC.GT.0) THEN I = STARTN(ARC) ELSE I = ENDN(-ARC) END IF END IF GOTO 510 END IF C IF THE PRICE IS TOO HIGH, THEN SOMETHING IS WRONG C IF (-PRICE.GT.MAXPRICE) THEN C GO TO 400 C ENDIF END IF 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 C IF (I.EQ.LASTQUEUE) THEN C LASTQUEUE=PREVNODE C NCYC=NCYC+1 C PRINT*,'CYCLE # ',NCYC,' # OF PRICE RISES ',NITER C 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 C (UNLESS THE PREVIOUS OPERATION INVOLVES A NON-SATURAT. FLOW PUSH C ALONG AN ARC BETWEEN I AND J, IN WHICH CASE WE ITERATE ON J NEXT). C I = NODE C DFCT(I)=DFCT1(I) C k=fout(i) C421 if (k.gt.0) then C dfct(i)=dfct(i)+f(k) C k=nxtou(k) C go to 421 C end if C k=fin(i) C422 if (k.gt.0) then C dfct(i)=dfct(i)-g(k)*f(k) C k=nxtin(k) C go to 422 C end if C IF ( DABS(DFCT(I)+SURPLUS(I)) .GT. 1.d0 ) THEN C PRINT*,'SURPLUS ERR AT',I,-DFCT(I),SURPLUS(I) C stop C END IF C DO A DIAGNOSTIC CHECK C DO 600 ARC=1,NA C IF (F(ARC).GT.REF_THR) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).LT.COST(ARC)-1.1*EPS) C $THEN C PRINT*,'1. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS,F(ARC) C stop C ENDIF C ENDIF C IF (F(ARC)+REF_THR.LT.U(ARC)) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).GT.COST(ARC)+1.1*EPS) C $THEN C PRINT*,'2. E-CS VIOLATED AT ARC ',ARC, C $ COST(ARC)-P(STARTN(ARC))+G(ARC)*P(ENDN(ARC)),EPS, C $ U(ARC)-F(ARC) C stop C ENDIF C ENDIF C IF ( (F(ARC).LT.-REF_THR) .OR. (F(ARC)-REF_THR.GT.U(ARC)) ) THEN C PRINT*,'FLOW OUTSIDE BDS AT',ARC,F(ARC),U(ARC) C STOP C END IF C IF ( DABS(DC(ARC,F(ARC))-COST(ARC)) .GT. REF_THR) THEN C PRINT*,'COST ERR AT',ARC,DC(ARC,F(ARC)),COST(ARC) C STOP C END IF C600 CONTINUE C RETURN TO ITERATE ON THE NEXT NODE IN THE QUEUE 620 CONTINUE c write(*,444)(surplus(k),k=1,n) c write(*,444)(p(k),k=1,n) 444 format(10e8.1) C------------------------------------------------------------ C DO 641 K=1,N C641 DFCT(K)=DFCT1(K) C DO 631 K=1,NA C DFCT(STARTN(K))=DFCT(STARTN(K))+F(K) C DFCT(ENDN(K))=DFCT(ENDN(K))-G(K)*F(K) C631 CONTINUE C DEL=0.0 C DO 635 K=1,N C IF ( DABS(DFCT(K)+SURPLUS(K)) .GT. 0.1 ) THEN C PRINT*,'SURPLUS ERR AT',K,-DFCT(K),SURPLUS(K) C STOP C END IF C DEL=DEL+DABS(DFCT(K)) C635 CONTINUE C PRINT*,'NODE=',NODE,p(NODE),surplus(NODE), DEL C------------------------------------------------------------ NXTNODE=NXTQUEUE(NODE) IF (DABS(SURPLUS(NODE)).GT.THRESHSURPLUS) THEN GOTO 100 ELSE IF (NODE.NE.NXTNODE) THEN NXTQUEUE(NODE)=0 NXTQUEUE(PREVNODE)=NXTNODE NODE=NXTNODE GO TO 100 END IF END IF C ************** END OF SUBPROBLEM (SCALING PHASE) ************** PRINT*,'END OF SCALING PHASE',EPS,THRESHSURPLUS LN=0 DO 501 NODE=1,N IF (DABS(SURPLUS(NODE)).GT.THRESHSURPLUS) THEN LN=LN+1 ENDIF 501 CONTINUE IF (LN.GT.0) THEN PRINT*,'SURPLUS EXCEED. THRESHOLD AT ',LN,' NODES ' C STOP END IF C PRINT*,'END OF CHECK OF OLD PHASE' C ****** IF DUALITY GAP IS WITHIN .1**ACCUR, THEN STOP; C ELSE REDUCE EPSILON****** C PCOST=0. DO 330 K=1,NA 330 PCOST=PCOST+C(K,F(K)) PRINT*,'PRIMAL COST =',PCOST DCOST=0. DO 341 I=1,N 341 DFCT(I)=DFCT1(I) DO 331 K=1,NA IF (CDAT(K,1).EQ.ZERO) THEN RC=A(K)-P(STARTN(K))+G(K)*P(ENDN(K)) F1=0 IF (RC.LT.0) F1=U(K) ELSE F1=DMIN1(DCINV(K,P(STARTN(K))-G(K)*P(ENDN(K))),U(K)) IF (F1.LT.0) F1=0 END IF DFCT(STARTN(K))=DFCT(STARTN(K))+F1 DFCT(ENDN(K))=DFCT(ENDN(K))-G(K)*F1 DCOST=DCOST+C(K,F1) 331 CONTINUE DO 335 I=1,N 335 DCOST=DCOST-P(I)*DFCT(I) PRINT*,'DUAL COST =',DCOST c write(*,77)(p(i),i=1,n) 77 format(10e8.2) C IF ( (DABS(PCOST-DCOST).LE.DABS(PCOST)*(.1**ACCUR)) $ .AND. (THRESHSURPLUS.LE.MAXSURPLUS*(.1**ACCUR)) ) THEN RETURN ELSE EPS=(EPS/FACTOR) THRESHSURPLUS=(THRESHSURPLUS/SFACTOR) THRESHSURPLUS=MAX(THRESHSURPLUS,EPS*N/(NA+1.0)) IF ((THRESHSURPLUS.LE.REF_THR).OR.(EPS.LE.REF_THR)) THEN PRINT*,'EXIT: THRESHSURPLUS=',THRESHSURPLUS,' EPS=',EPS RETURN END IF END IF C IF (EPS.LT. 1.0/(n+1.0)) THRESHSURPLUS=REF_THR C RESET THE FLOWS & THE PUSH LISTS; FIND THE MINIMAL PRICE XP=LARGE DO 800 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 QFPUSHF(NODE)=0 QFPUSHB(NODE)=0 IF (P(NODE).LT.XP) XP=P(NODE) 800 CONTINUE C REDUCE ALL PRICES TO REDUCE DANGER OF OVERFLOW C DEL=XP+(LARGE/10) C IF (DEL.LT.0.0) DEL=0.0 C DO 810 NODE=1,N C P(NODE)=P(NODE)-DEL C810 CONTINUE C PRINT*,'MODIFYING PRICES & 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-G(ARC)*PEND-EPS-COST(ARC)).GE.0.0) THEN IF (CDAT(ARC,1).NE.ZERO) THEN UP = MIN(DCINV(ARC,pstart-g(arc)*pend),U(ARC)) COST(ARC) = DC(ARC,UP) ELSE UP = U(ARC) END IF RESID=UP-F(ARC) C IF (RESID.GT.REF_THR) THEN SURPLUS(START)=SURPLUS(START)-RESID SURPLUS(END)=SURPLUS(END)+G(ARC)*RESID F(ARC)=UP C END IF ELSE IF ((-PSTART+G(ARC)*PEND-EPS+COST(ARC)).GE.0.0) THEN IF (CDAT(ARC,1).NE.ZERO) THEN FLOW=DMAX1(DCINV(ARC,pstart-g(arc)*pend),0.D0) COST(ARC) = DC(ARC,FLOW) ELSE FLOW =0.0 END IF FLOW=F(ARC)-FLOW C IF (FLOW.GT.REF_THR) THEN SURPLUS(START)=SURPLUS(START)+FLOW SURPLUS(END)=SURPLUS(END)-G(ARC)*FLOW F(ARC)=F(ARC)-FLOW C END IF END IF END IF 900 CONTINUE C print*,'type 1 for next phase' C read(*,*)k C RETURN FOR ANOTHER PHASE C PRINT*,'CHECKING E-CS BEFORE STARTING NEW PHASE' C DO 700 ARC=1,NA C IF (F(ARC).GT.REF_THR) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).LT.-EPS+COST(ARC)) THEN C PRINT*,'E-CS VIOLATED AT ARC ',ARC C ENDIF C ENDIF C IF (F(ARC)+REF_THR.LT.U(ARC)) THEN C IF (P(STARTN(ARC))-G(ARC)*P(ENDN(ARC)).GT.EPS+COST(ARC)) THEN C PRINT*,'E-CS VIOLATED AT ARC ',ARC C ENDIF C ENDIF 700 CONTINUE GO TO 60 400 PRINT*,'PROBLEM IS INFEASIBLE' C PAUSE STOP END C C SUBROUTINE CYCLE_PUSH1(NODE,J,IARC) IMPLICIT NONE INTEGER*4 NODE,J,IARC INTEGER*4 N,NA,MAXNODES,MAXARCS,MAXNDAT PARAMETER (MAXNODES=5000, MAXARCS=20000, MAXNDAT=5) INTEGER*4 STARTN(MAXARCS),ENDN(MAXARCS) INTEGER*4 PRDCSR(MAXNODES) REAL*8 A(MAXARCS),CDAT(MAXARCS,MAXNDAT),U(MAXARCS),G(MAXARCS) REAL*8 COST(MAXARCS),P(MAXNODES),F(MAXARCS),SURPLUS(MAXNODES), $ LARGE,ZERO,DC,DCINV COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK8/ PRDCSR COMMON /LCOS/ A COMMON /QCO1/ CDAT COMMON /UBOUND/ U COMMON /QG/ G COMMON /COST/ COST COMMON /FLOW/ F COMMON /BLK3/ SURPLUS COMMON /PRICES/ P INTEGER*4 I,K,ARC REAL*8 RESID,GAIN1,DELTA1,GAIN2,DELTA2 REAL*8 REF_THR C print*,' ' C print*,'do cycle push1 at',node,j,iarc REF_THR = 1.0E-10 GAIN1=1.0 DELTA1=LARGE I = J ARC = IARC c print*,i,iarc,delta1,gain1 10 IF (ARC.GT.0) THEN K = STARTN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = U(ARC)-F(ARC) ELSE RESID = MIN(U(ARC),DCINV(ARC,P(K)-G(ARC)*P(I)))-F(ARC) END IF DELTA1 = MIN(DELTA1/G(ARC),RESID) GAIN1 = GAIN1*G(ARC) ELSE ARC = -ARC K = ENDN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = F(ARC) ELSE RESID = F(ARC)-MAX(0.0,DCINV(ARC,P(I)-G(ARC)*P(K))) END IF DELTA1 = MIN(DELTA1,RESID)*G(ARC) GAIN1 = GAIN1/G(ARC) END IF c print*,i,k,j,arc,delta1,gain1 C C ***** IF J IS NOT REACHED, MOVE TO THE PREDECESSOR NODE C IF (K.NE.J) THEN I = K ARC = PRDCSR(I) GO TO 10 END IF C print*,'out of cycle1',gain1,delta1 C IF (GAIN1.LT. 1.D0 - REF_THR) THEN DELTA2 = DELTA1*(1.D0-GAIN1) GAIN2 = 1.0 C C ***** COMPUTE MAXIMUM FLOW CHANGE ALONG PREDECESSOR ARCS C I = J ARC = PRDCSR(I) 20 IF (ARC.NE.NA+1) THEN IF (ARC.GT.0) THEN K=STARTN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = U(ARC)-F(ARC) ELSE RESID = MIN(U(ARC),DCINV(ARC,P(K)-G(ARC)*P(I)))-F(ARC) END IF DELTA2 = MIN(DELTA2/G(ARC),RESID) GAIN2 = GAIN2*G(ARC) ELSE ARC = -ARC K=ENDN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = F(ARC) ELSE RESID = F(ARC)-MAX(0.0,DCINV(ARC,P(I)-G(ARC)*P(K))) END IF DELTA2 = MIN(DELTA2,RESID)*G(ARC) GAIN2 = GAIN2/G(ARC) END IF I=K ARC=PRDCSR(I) GO TO 20 END IF DELTA2 = MIN(DELTA2,SURPLUS(I)) C print*,'NODE=',I,surplus(i),surplus(i)-delta2 SURPLUS(I) = SURPLUS(I)-DELTA2 DELTA2 = DELTA2*GAIN2 DELTA1 = DELTA2/(1.D0-GAIN1) C C ***** UPDATE FLOW & SET PREDECESSOR FOR ALL REMAINING NODES TO ZERO C I = J ARC = PRDCSR(I) 30 IF (ARC.NE.NA+1) THEN IF (ARC.GT.0) THEN K=STARTN(ARC) DELTA2 = DELTA2/G(ARC) F(ARC) = F(ARC) + DELTA2 ELSE ARC = -ARC K=ENDN(ARC) F(ARC) = F(ARC) - DELTA2 DELTA2 = DELTA2*G(ARC) END IF COST(ARC) = DC(ARC,F(ARC)) PRDCSR(I)=0 I=K ARC=PRDCSR(I) GO TO 30 END IF PRDCSR(I)=0 ELSE DELTA2 = 0.0 C C ***** SET PREDECESSOR FOR ALL REMAINING NODES TO ZERO C I = K ARC = PRDCSR(I) 40 IF (ARC.NE.NA+1) THEN IF (ARC.GT.0) THEN K=STARTN(ARC) ELSE K=ENDN(-ARC) END IF PRDCSR(I)=0 I=K ARC=PRDCSR(I) GO TO 40 END IF PRDCSR(I)=0 SURPLUS(J)=SURPLUS(J) - DELTA1*(1.D0-GAIN1) END IF C print*,'push1',node,j,iarc,surplus(j) C C **** DO FLOW PUSH ALONG THE CYLE BY THE AMOUNT DELTA1 C DELTA1=DELTA1*GAIN1 IF (IARC.GT.0) THEN DELTA1=DELTA1/G(IARC) F(IARC)=F(IARC)+DELTA1 ELSE IARC=-IARC F(IARC)=F(IARC)-DELTA1 DELTA1=DELTA1*G(IARC) END IF COST(IARC) = DC(IARC,F(IARC)) I=NODE ARC = PRDCSR(I) 50 IF (ARC.GT.0) THEN K = STARTN(ARC) DELTA1=DELTA1/G(ARC) F(ARC)=F(ARC)+DELTA1 ELSE ARC = -ARC K = ENDN(ARC) F(ARC)=F(ARC)-DELTA1 DELTA1=DELTA1*G(ARC) END IF COST(ARC) = DC(ARC,F(ARC)) C C ***** IF J IS NOT REACHED, MOVE TO THE PREDECESSOR NODE C PRDCSR(I) = 0 IF (K.NE.J) THEN I = K ARC = PRDCSR(I) GO TO 50 END IF C print*,'cycle push done!',surplus(j) C END C C SUBROUTINE CYCLE_PUSH2(NODE,J,IARC) IMPLICIT NONE INTEGER*4 NODE,J,IARC INTEGER*4 N,NA,MAXNODES,MAXARCS,MAXNDAT PARAMETER (MAXNODES=5000, MAXARCS=20000, MAXNDAT=5) INTEGER*4 STARTN(MAXARCS),ENDN(MAXARCS) INTEGER*4 PRDCSR(MAXNODES) REAL*8 A(MAXARCS),CDAT(MAXARCS,MAXNDAT),U(MAXARCS),G(MAXARCS) REAL*8 COST(MAXARCS),P(MAXNODES),F(MAXARCS),SURPLUS(MAXNODES), $LARGE,ZERO,DC,DCINV COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK8/ PRDCSR COMMON /LCOS/ A COMMON /QCO1/ CDAT COMMON /UBOUND/ U COMMON /QG/ G COMMON /COST/ COST COMMON /FLOW/ F COMMON /BLK3/ SURPLUS COMMON /PRICES/ P INTEGER*4 I,K,ARC REAL*8 RESID,GAIN1,DELTA1,GAIN2,DELTA2 REAL*8 REF_THR c print*,' ' C print*,'do cycle push2 at',node,j,iarc REF_THR = 1.0E-10 GAIN1=1.0 DELTA1=LARGE I = J ARC = IARC 10 IF (ARC.GT.0) THEN K = STARTN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = F(ARC) ELSE RESID = F(ARC)-MAX(0.0,DCINV(ARC,P(K)-G(ARC)*P(I))) END IF DELTA1 = MIN(DELTA1/G(ARC),RESID) GAIN1 = GAIN1*G(ARC) ELSE ARC = -ARC K = ENDN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = U(ARC)-F(ARC) ELSE RESID = MIN(U(ARC),DCINV(ARC,P(I)-G(ARC)*P(K)))-F(ARC) END IF DELTA1 = MIN(DELTA1,RESID)*G(ARC) GAIN1 = GAIN1/G(ARC) END IF c print*,i,k,arc,prdcsr(i),delta1,gain1 C C ***** IF J IS NOT REACHED, MOVE TO THE PREDECESSOR NODE C IF (K.NE.J) THEN I = K ARC = PRDCSR(I) GO TO 10 END IF C print*,'out of cycle2',gain1,delta1 IF (GAIN1.LT. 1.D0 - REF_THR) THEN DELTA2 = DELTA1*(1.D0-GAIN1) GAIN2 = 1.0 C C ***** COMPUTE MAXIMUM FLOW CHANGE ALONG PREDECESSOR ARCS C I = J ARC = PRDCSR(I) 20 IF (ARC.NE.NA+1) THEN IF (ARC.GT.0) THEN K=STARTN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = F(ARC) ELSE RESID = F(ARC)-MAX(0.0,DCINV(ARC,P(K)-G(ARC)*P(I))) END IF DELTA2 = MIN(DELTA2/G(ARC),RESID) GAIN2 = GAIN2*G(ARC) ELSE ARC = -ARC K=ENDN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = U(ARC)-F(ARC) ELSE RESID = MIN(U(ARC),DCINV(ARC,P(I)-G(ARC)*P(K)))-F(ARC) END IF DELTA2 = MIN(DELTA2,RESID)*G(ARC) GAIN2 = GAIN2/G(ARC) END IF c print*,i,k,arc,prdcsr(i),delta2,gain2 I=K ARC=PRDCSR(I) GO TO 20 END IF DELTA2 = MIN(DELTA2,-SURPLUS(I)) c print*,'NODE=',I,surplus(i),surplus(i)-delta2 SURPLUS(I) = SURPLUS(I)+DELTA2 DELTA2 = DELTA2*GAIN2 DELTA1 = DELTA2/(1.D0-GAIN1) C C ***** UPDATE FLOW & SET PREDECESSOR FOR ALL REMAINING NODES TO ZERO C I = J ARC = PRDCSR(I) 30 IF (ARC.NE.NA+1) THEN IF (ARC.GT.0) THEN K=STARTN(ARC) DELTA2 = DELTA2/G(ARC) c print*,'%arc2',arc,delta2,g(arc) F(ARC) = F(ARC) - DELTA2 ELSE ARC = -ARC K=ENDN(ARC) F(ARC) = F(ARC) + DELTA2 DELTA2 = DELTA2*G(ARC) END IF COST(ARC) = DC(ARC,F(ARC)) PRDCSR(I)=0 I=K ARC=PRDCSR(I) GO TO 30 END IF PRDCSR(I)=0 ELSE C C ***** SET PREDECESSOR FOR ALL REMAINING NODES TO ZERO C I = K ARC = PRDCSR(I) 40 IF (ARC.NE.NA+1) THEN IF (ARC.GT.0) THEN K=STARTN(ARC) ELSE K=ENDN(-ARC) END IF PRDCSR(I)=0 I=K ARC=PRDCSR(I) GO TO 40 END IF PRDCSR(I)=0 SURPLUS(J)=SURPLUS(J) + DELTA1*(1.D0-GAIN1) END IF c print*,'push2',node,j,iarc,surplus(j) C C **** DO FLOW PUSH ALONG THE CYLE BY THE AMOUNT DELTA1 C DELTA1=DELTA1*GAIN1 IF (IARC.GT.0) THEN DELTA1=DELTA1/G(IARC) c print*,'%iarc',iarc,delta1,g(iarc) F(IARC)=F(IARC)-DELTA1 ELSE IARC=-IARC F(IARC)=F(IARC)+DELTA1 DELTA1=DELTA1*G(IARC) END IF COST(IARC) = DC(IARC,F(IARC)) I=NODE ARC = PRDCSR(I) 50 IF (ARC.GT.0) THEN K = STARTN(ARC) DELTA1=DELTA1/G(ARC) c print*,'%arc1',arc,delta1,g(arc) F(ARC)=F(ARC)-DELTA1 ELSE ARC = -ARC K = ENDN(ARC) F(ARC)=F(ARC)+DELTA1 DELTA1=DELTA1*G(ARC) END IF COST(ARC) = DC(ARC,F(ARC)) C C ***** IF J IS NOT REACHED, MOVE TO THE PREDECESSOR NODE C PRDCSR(I) = 0 IF (K.NE.J) THEN I = K ARC = PRDCSR(I) GO TO 50 END IF C print*,'cycle push done!',surplus(j) C END C C THIS ROUTINE PERFORMS A FLOW PUSH ALONG A PATH TO J C SUBROUTINE PATH_PUSH1(J,IARC) IMPLICIT NONE INTEGER*4 J,IARC INTEGER*4 N,NA,MAXNODES,MAXARCS,MAXNDAT PARAMETER (MAXNODES=5000, MAXARCS=20000, MAXNDAT=5) INTEGER*4 STARTN(MAXARCS),ENDN(MAXARCS) INTEGER*4 PRDCSR(MAXNODES) REAL*8 A(MAXARCS),CDAT(MAXARCS,MAXNDAT),U(MAXARCS),G(MAXARCS) REAL*8 COST(MAXARCS),P(MAXNODES),F(MAXARCS),SURPLUS(MAXNODES), $LARGE,ZERO,DC,DCINV COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK8/ PRDCSR COMMON /LCOS/ A COMMON /QCO1/ CDAT COMMON /UBOUND/ U COMMON /QG/ G COMMON /COST/ COST COMMON /FLOW/ F COMMON /BLK3/ SURPLUS COMMON /PRICES/ P INTEGER*4 I,K,ARC REAL*8 RESID,GAIN,DELTA REAL*8 REF_THR C print*,'do PATH push1 TO',j,iarc REF_THR = 1.0E-10 GAIN=1.0 DELTA=-SURPLUS(J) I = J ARC = IARC 10 IF (ARC.GT.0) THEN K = STARTN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = U(ARC)-F(ARC) ELSE RESID = MIN(U(ARC),DCINV(ARC,P(K)-G(ARC)*P(I)))-F(ARC) END IF DELTA = MIN(DELTA/G(ARC),RESID) GAIN = GAIN*G(ARC) ELSE ARC = -ARC K = ENDN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = F(ARC) ELSE RESID = F(ARC)-MAX(0.0,DCINV(ARC,P(I)-G(ARC)*P(K))) END IF DELTA = MIN(DELTA,RESID)*G(ARC) GAIN = GAIN/G(ARC) END IF C print*,i,k,j,arc,delta,gain C C ***** IF J IS NOT REACHED, MOVE TO THE PREDECESSOR NODE C I = K ARC = PRDCSR(I) IF (ARC.NE.NA+1) GO TO 10 DELTA = MIN(DELTA,SURPLUS(I)) C print*,'push1',I,j,iarc,surplus(I) C C **** DO FLOW PUSH ALONG THE PATH BY THE AMOUNT DELTA C SURPLUS(I)=SURPLUS(I) - DELTA SURPLUS(J)=SURPLUS(J) + DELTA*GAIN I = J ARC = IARC 20 IF (ARC.GT.0) THEN K = STARTN(ARC) GAIN=GAIN/G(ARC) F(ARC)=F(ARC)+DELTA*GAIN ELSE ARC = -ARC K = ENDN(ARC) F(ARC)=F(ARC)-DELTA*GAIN GAIN=GAIN*G(ARC) END IF COST(ARC) = DC(ARC,F(ARC)) PRDCSR(I) = 0 I = K ARC = PRDCSR(I) IF (ARC .NE. NA+1) GO TO 20 PRDCSR(I)=0 C print*,'PATH push done!' END C C SUBROUTINE PATH_PUSH2(J,IARC) IMPLICIT NONE INTEGER*4 J,IARC INTEGER*4 N,NA,MAXNODES,MAXARCS,MAXNDAT PARAMETER (MAXNODES=5000, MAXARCS=20000, MAXNDAT=5) INTEGER*4 STARTN(MAXARCS),ENDN(MAXARCS) INTEGER*4 PRDCSR(MAXNODES) REAL*8 A(MAXARCS),CDAT(MAXARCS,MAXNDAT),U(MAXARCS),G(MAXARCS) REAL*8 COST(MAXARCS),P(MAXNODES),F(MAXARCS),SURPLUS(MAXNODES), $LARGE,ZERO,DC,DCINV COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK8/ PRDCSR COMMON /LCOS/ A COMMON /QCO1/ CDAT COMMON /UBOUND/ U COMMON /QG/ G COMMON /COST/ COST COMMON /FLOW/ F COMMON /BLK3/ SURPLUS COMMON /PRICES/ P INTEGER*4 I,K,ARC REAL*8 RESID,GAIN,DELTA REAL*8 REF_THR C print*,'do PATH push2 TO',j,iarc REF_THR = 1.0E-10 GAIN=1.0 DELTA=SURPLUS(J) I = J ARC = IARC 10 IF (ARC.GT.0) THEN K = STARTN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = F(ARC) ELSE RESID = F(ARC)-MAX(0.0,DCINV(ARC,P(K)-G(ARC)*P(I))) END IF DELTA = MIN(DELTA/G(ARC),RESID) GAIN = GAIN*G(ARC) ELSE ARC = -ARC K = ENDN(ARC) IF (CDAT(ARC,1).EQ.ZERO) THEN RESID = U(ARC)-F(ARC) ELSE RESID = MIN(U(ARC),DCINV(ARC,P(I)-G(ARC)*P(K)))-F(ARC) END IF DELTA = MIN(DELTA,RESID)*G(ARC) GAIN = GAIN/G(ARC) END IF C print*,i,k,j,arc,delta,gain C C ***** IF J IS NOT REACHED, MOVE TO THE PREDECESSOR NODE C I = K ARC = PRDCSR(I) IF (ARC.NE.NA+1) GO TO 10 DELTA = MIN(DELTA,-SURPLUS(I)) C print*,'push2',I,j,iarc,surplus(I) C C **** DO FLOW PUSH ALONG THE PATH BY THE AMOUNT DELTA C SURPLUS(I)=SURPLUS(I) + DELTA SURPLUS(J)=SURPLUS(J) - DELTA*GAIN I = J ARC = IARC 20 IF (ARC.GT.0) THEN K = STARTN(ARC) GAIN=GAIN/G(ARC) F(ARC)=F(ARC)-DELTA*GAIN ELSE ARC = -ARC K = ENDN(ARC) F(ARC)=F(ARC)+DELTA*GAIN GAIN=GAIN*G(ARC) END IF COST(ARC) = DC(ARC,F(ARC)) PRDCSR(I) = 0 I = K ARC = PRDCSR(I) IF (ARC .NE. NA+1) GO TO 20 PRDCSR(I)=0 C print*,'PATH push done!' END