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. C THE METHOD AND NUMERICAL C TEST RESULTS ARE DESCRIBED IN THE PAPERS: C D.P. Bertsekas and P. Tseng, C "An epsilon-Relaxation Method for Generalized Separable Convex C Cost Network Flow Problems", Math. Program., 88, 2000, 85-104. C 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 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. THE LINEAR COST COEFFICIENTS ARE STORED C IN THE ARRAY C A(I), I=1,...,NA, C WITH CDAT(I,1) = 0 RESERVED TO INDICATE THAT ARC I HAS C A LINEAR COST FUNCTION. C C VERSION 1: March 1, 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 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 enter now' C READ(*,*) SMAL C PRINT*, SMAL C end if PRINT*,'E-RELAXATION 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. ZERO).and.(smal.gt.1.0e-9)) then CDAT(I,1) = SMAL END IF 5 CONTINUE C READ SUPPLY OF EACH NODE (IS MAXSURPLUS USEFUL?) 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 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/5.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 = 5.0 WRITE(*,*) 'INITIAL 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 E-RELAX 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 SURPI,FLOW,RESID,UP,CAPOUT,CAPIN REAL*8 REF_THR,C,DC,DCINV,EPS2,GAIN LOGICAL*1 FLAG C ADDITIONAL DATA STRUCTURE USED TO COMPUTE DUAL COST REAL*8 PCOST,DCOST,RC,F1,DFCT(MAXNODES),DFCT1(MAXNODES) 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 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 NODE=1,N-1 NXTQUEUE(NODE)=NODE+1 82 CONTINUE NXTQUEUE(N)=1 I=1 PREVNODE=N LASTQUEUE=N FLAG = .FALSE. C ***** PHASE I: UP ITERATIONS ON NODES WITH POSITIVE SURPLUS **** 100 CONTINUE C TAKE UP NEXT NODE (NODE I) FOR ITERATION SURPI=SURPLUS(I) IF (SURPI .GT. THRESHSURPLUS) THEN C ARCF & ARCB ARE THE CURRENT VALUES OF THE STARTING ARCS OF THE C PUSH LISTS OF I TOTALITER=TOTALITER+1 ARCF = FPUSHF(I) ARCB = FPUSHB(I) QARCF = QFPUSHF(I) QARCB = QFPUSHB(I) PRICE = P(I) C PRINT*,'NODE',I,' PRICE=',PRICE,' surpi=',surpi C if ((i.eq.17) .and. (surpi .gt. 0)) then C print*,'0',arcf,arcb,qarcf,qarcb C end if 115 IF ((ARCF .GT. 0) .OR. (ARCB .GT. 0)) THEN C START BY TRYING TO PUSH AWAY FLOW ON ARCS THAT WERE C ADMISSIBLE AT THE END OF THE LAST ITERATION (IF ANY) C AT THIS NODE. WE MUST CHECK THAT THEY ARE STILL C ADMISSIBLE. 120 IF ((SURPI .GT. REF_THR) .AND. (ARCF .GT. 0)) THEN J = ENDN(ARCF) GAIN = G(ARCF) IF (PRICE-GAIN*P(J)-COST(ARCF).GT.EPS2) THEN C********TAKING CARE OF LINEAR FORWARD ARCS********************* RESID = U(ARCF) - F(ARCF) C if ((i.eq.17) .and. (surpi.gt.0.0)) then C print*,'0.5',resid,gain,price,p(j) C end if IF (RESID.GT.REF_THR) THEN IF (SURPI .GE. RESID) THEN F(ARCF) = U(ARCF) SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + GAIN*RESID ARCF = NXTPUSHF(ARCF) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(ARCF) = F(ARCF) + SURPI SURPLUS(J) = SURPLUS(J) + GAIN*SURPI SURPI = 0.0 c if ((i.eq.6) .and. (surpi .lt. 21.53)) then c print*,'j=',j,surpi,surplus(j),flag,prdcsr(j) c end if C C ***** IF FLOW PUSH IS NON-SATURATING, STORE ARCF AND CHECK FOR FLOW PUSH C ALONG CYCLE C IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=ARCF ELSE CALL CYCLE_PUSH1(I,J,ARCF) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCF = NXTPUSHF(ARCF) END IF ELSE ARCF = NXTPUSHF(ARCF) END IF GOTO 120 END IF 121 IF ((SURPI .GT. REF_THR) .AND. (ARCB .GT. 0)) THEN J = STARTN(ARCB) GAIN = G(ARCB) IF (GAIN*PRICE-P(J)+COST(ARCB).GT.EPS2) THEN C *********** TAKING CARE OF THE LINEAR COST REVERSE ARCS *************** RESID=F(ARCB) IF (RESID.GT.REF_THR) THEN IF (SURPI .GE. GAIN*RESID) THEN SURPI = SURPI - GAIN*RESID SURPLUS(J) = SURPLUS(J) + RESID F(ARCB) = 0 ARCB = NXTPUSHB(ARCB) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(ARCB) = F(ARCB) - SURPI/GAIN SURPLUS(J) = SURPLUS(J) + SURPI/GAIN SURPI = 0.0 C C ***** IF FLOW PUSH IS NON-SATURATING, STORE ARCB AND CHECK FOR FLOW PUSH C ALONG CYCLE C IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=-ARCB ELSE CALL CYCLE_PUSH1(I,J,-ARCB) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCB = NXTPUSHB(ARCB) END IF ELSE ARCB = NXTPUSHB(ARCB) END IF GOTO 121 END IF END IF 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. (surpi .gt. 0)) then C print*,'2',qarcf,qarcb,surpi,eps2 c end if IF ((QARCF+QARCB).GT.0) THEN 122 IF ((SURPI .GT. REF_THR) .AND. (QARCF .GT. 0)) THEN J = ENDN(QARCF) GAIN = G(QARCF) IF (PRICE-GAIN*P(J)-COST(QARCF).GT.EPS2) THEN C**************TAKING CARE OF NONLINEAR FORWARD ARCS********************* RESID=MIN(U(QARCF),DCINV(QARCF,PRICE-GAIN*P(J))) $ - F(QARCF) C if (PRICE-GAIN*P(J).gt.100) then C print*,qarcf,j,p(j),price C end if IF (RESID.GT.REF_THR) THEN IF (SURPI .GE. RESID) THEN F(QARCF) = F(QARCF)+RESID SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + GAIN*RESID COST(QARCF) = DC(QARCF,F(QARCF)) QARCF = QNXTPUSHF(QARCF) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(QARCF) = F(QARCF) + SURPI COST(QARCF) = DC(QARCF,F(QARCF)) SURPLUS(J) = SURPLUS(J) + GAIN*SURPI SURPI = 0.0 C C ***** IF FLOW PUSH IS NON-SATURATING, STORE QARCF AND CHECK FOR FLOW PUSH C ALONG CYCLE C IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=QARCF ELSE CALL CYCLE_PUSH1(I,J,QARCF) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF GOTO 122 END IF 123 IF ((SURPI .GT. REF_THR) .AND. (QARCB .GT. 0)) THEN J = STARTN(QARCB) GAIN = G(QARCB) IF (GAIN*PRICE-P(J)+COST(QARCB).GT.EPS2) THEN C ************* TAKING CARE OF ARCS WITH NON LINEAR COST ****************** RESID=F(QARCB)-MAX(0.0,DCINV(QARCB,P(J)-GAIN*PRICE)) C if (P(j)-GAIN*Price.gt.100) then C print*,qarcb,j,p(j),price C end if IF (RESID.GT.REF_THR) THEN IF (SURPI .GE. GAIN*RESID) THEN F(QARCB) = F(QARCB)-RESID SURPI = SURPI - GAIN*RESID SURPLUS(J) = SURPLUS(J) + RESID COST(QARCB) = DC(QARCB,F(QARCB)) QARCB = QNXTPUSHB(QARCB) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(QARCB) = F(QARCB) - SURPI/GAIN COST(QARCB) = DC(QARCB,F(QARCB)) SURPLUS(J) = SURPLUS(J) + SURPI/GAIN SURPI = 0.0 C C *** IF FLOW PUSH IS NON-SATURATING, STORE QARCB AND CHECK FOR FLOW PUSH C ALONG A CYCLE C IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=-QARCB ELSE CALL CYCLE_PUSH1(I,J,-QARCB) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF GOTO 123 END IF END IF 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. (surpi .gt. 0)) then C print*,'3',i,arcf,arcb,qarcf,qarcb,price c end if IF ((SURPI.GT.THRESHSURPLUS).AND. $ (ARCF+ARCB+QARCF+QARCB.EQ.0)) THEN IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF 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' C PAUSE 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' C PAUSE STOP END IF END IF ARC = NXTIN(ARC) GOTO 212 END IF c IF (SURPI.GT.REF_THR) print*,i,surpi IF (SURPI.GT.REF_THR) GOTO 400 IF (PRICE.LT.IDINT(LARGE/10)) PRICE=IDINT(LARGE/10) END IF C print*,'price rise at ',i,p(i),price 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.LE.0) THEN C PRINT*,'NONPOSITIVE PRICE INCREMENT=',PRICEINCR C PAUSE C STOP C ENDIF C DO A DIAGNOSTIC CHECK 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; IF THERE IS STILL A LOT OF SURPLUS, C TRY AND DO SOME PUSHING. IF (SURPI.GT.THRESHSURPLUS) GOTO 115 ELSE C IF NO PRICE RISE TOOK PLACE RESET THE START OF THE PUSH LISTS FPUSHF(I) = ARCF FPUSHB(I) = ARCB QFPUSHF(I) = QARCF QFPUSHB(I) = QARCB END IF C DO THE FINAL BOOKKEEPING OF THE ITERATION SURPLUS(I) = SURPI 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 ********** END OF UP ITERATION STARTING AT NODE I ********* C CHECK FOR THE END OF THE QUEUE. THIS STEP IS NOT NECESSARY FOR C THE ALGORITHM; IT IS INCLUDED FOR COLLECTING STATISTICS ABOUT C THE ALGORITHM'S BEHAVIOR, E.G. COUNTING THE NUMBER OF CYCLES 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 print*,' i ',i,' p(i)',p(i),surplus(i),flag C do 88 k=1,n C arc=fpushf(k) C89 if (arc.gt.0) then C if (endn(arc).eq.k) then C print*,k,arc,startn(arc),endn(arc) C stop C end if C arc=nxtpushf(arc) C end if C88 continue NXTNODE=NXTQUEUE(I) IF (I.NE.NXTNODE) THEN NXTQUEUE(I)=0 NXTQUEUE(PREVNODE)=NXTNODE C C IF LAST OPERATION IS A NONSATURATING FLOW PUSH TOWARDS J, ITERATE ON J NEXT C IF (FLAG) THEN K=0 240 IF (NXTNODE.NE.J) THEN PREVNODE=NXTNODE NXTNODE=NXTQUEUE(NXTNODE) K=K+1 GO TO 240 END IF I=J c print*,'k= ',k,' j=',j,' surp(j)=',surplus(j) END IF I=NXTNODE GO TO 100 END IF C ***** END OF PHASE I ***** C ***** PUSHLIST INITIALIZATION ***** DO 481 NODE=1,N FPUSHF(NODE)=0 FPUSHB(NODE)=0 QFPUSHF(NODE)=0 QFPUSHB(NODE)=0 PRDCSR(NODE)=0 481 CONTINUE C ***** QUEUE INITIALIZATION ***** DO 482 NODE=1,N-1 NXTQUEUE(NODE)=NODE+1 482 CONTINUE NXTQUEUE(N)=1 I=1 PREVNODE=N LASTQUEUE=N FLAG = .FALSE. C ***** PHASE II: DOWN ITERATIONS ON NODES WITH NEGATIVE SURPLUS ***** 500 CONTINUE C TAKE UP NEXT NODE (NODE I) FOR ITERATION SURPI=SURPLUS(I) IF (-SURPI .GT. THRESHSURPLUS) THEN C ARCF & ARCB ARE THE CURRENT VALUES OF THE STARTING ARCS OF THE C PUSH LISTS OF I TOTALITER=TOTALITER+1 ARCF = FPUSHF(I) ARCB = FPUSHB(I) QARCF = QFPUSHF(I) QARCB = QFPUSHB(I) PRICE = P(I) C PRINT*,'NODE',I,' PRICE=',PRICE,' surpi=',surpi c if ((i.eq.14) .and. (-surpi .gt. 0)) then c print*,'0',arcf,arcb,qarcf,qarcb c end if 515 IF ((ARCF .GT. 0) .OR. (ARCB .GT. 0)) THEN C START BY TRYING TO PUSH AWAY FLOW ON ARCS THAT WERE C ADMISSIBLE AT THE END OF THE LAST ITERATION (IF ANY) C AT THIS NODE. WE MUST CHECK THAT THEY ARE STILL C ADMISSIBLE. 520 IF ((-SURPI .GT. REF_THR) .AND. (ARCF .GT. 0)) THEN J = ENDN(ARCF) GAIN = G(ARCF) IF (-(PRICE-GAIN*P(J)-COST(ARCF)).GT.EPS2) THEN C******************TAKING CARE OF LINEAR FORWARD ARCS*************************** RESID = - F(ARCF) IF (-RESID.GT.REF_THR) THEN IF (-SURPI .GE. -RESID) THEN F(ARCF) = 0.0 SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + GAIN*RESID ARCF = NXTPUSHF(ARCF) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(ARCF) = F(ARCF) + SURPI SURPLUS(J) = SURPLUS(J) + GAIN*SURPI SURPI = 0.0 C C ***** IF FLOW PUSH IS NON-SATURATING, STORE ARCF AND CHECK FOR FLOW PUSH C ALONG CYCLE C IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=ARCF ELSE CALL CYCLE_PUSH2(I,J,ARCF) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCF = NXTPUSHF(ARCF) END IF ELSE ARCF = NXTPUSHF(ARCF) END IF GOTO 520 END IF 521 IF ((-SURPI .GT. REF_THR) .AND. (ARCB .GT. 0)) THEN J = STARTN(ARCB) GAIN = G(ARCB) IF (-(GAIN*PRICE-P(J)+COST(ARCB)).GT.EPS2) THEN C *********** TAKING CARE OF THE LINEAR COST REVERSE ARCS *************** RESID=F(ARCB)-U(ARCB) IF (-RESID.GT.REF_THR) THEN IF (-SURPI .GE. -GAIN*RESID) THEN SURPI = SURPI - GAIN*RESID SURPLUS(J) = SURPLUS(J) + RESID F(ARCB) = U(ARCB) ARCB = NXTPUSHB(ARCB) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(ARCB) = F(ARCB) - SURPI/GAIN SURPLUS(J) = SURPLUS(J) + SURPI/GAIN SURPI = 0.0 C C ***** IF FLOW PUSH IS NON-SATURATING, STORE ARCB AND CHECK FOR FLOW PUSH C ALONG CYCLE C IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=-ARCB ELSE CALL CYCLE_PUSH2(I,J,-ARCB) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE ARCB = NXTPUSHB(ARCB) END IF ELSE ARCB = NXTPUSHB(ARCB) END IF GOTO 521 END IF END IF 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)+SURPI) .GT. 1.d0 ) THEN c PRINT*,'SURPLUS ERR AT',I,-DFCT(I),SURPI c END IF c if ((i.eq.14) .and. (-surpi .gt. 0)) then C print*,'2',qarcf,qarcb,surpi,eps2 c end if IF ((QARCF+QARCB).GT.0) THEN 522 IF ((-SURPI .GT. REF_THR) .AND. (QARCF .GT. 0)) THEN J = ENDN(QARCF) GAIN = G(QARCF) IF (-(PRICE-GAIN*P(J)-COST(QARCF)).GT.EPS2) THEN C**************TAKING CARE OF NONLINEAR FORWARD ARCS******************************** RESID=MAX(0.0,DCINV(QARCF,PRICE-GAIN*P(J)))-F(QARCF) IF (-RESID.GT.REF_THR) THEN IF (-SURPI .GE. -RESID) THEN F(QARCF) = F(QARCF)+RESID SURPI = SURPI - RESID SURPLUS(J) = SURPLUS(J) + GAIN*RESID COST(QARCF) = DC(QARCF,F(QARCF)) QARCF = QNXTPUSHF(QARCF) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(QARCF) = F(QARCF) + SURPI COST(QARCF) = DC(QARCF,F(QARCF)) SURPLUS(J) = SURPLUS(J) + GAIN*SURPI SURPI = 0.0 C C ***** IF FLOW PUSH IS NON-SATURATING, STORE QARCF AND CHECK FOR FLOW PUSH C ALONG CYCLE C IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=QARCF ELSE CALL CYCLE_PUSH2(I,J,QARCF) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF ELSE QARCF = QNXTPUSHF(QARCF) END IF GOTO 522 END IF 523 IF ((-SURPI .GT. REF_THR) .AND. (QARCB .GT. 0)) THEN J = STARTN(QARCB) GAIN = G(QARCB) IF ((-(GAIN*PRICE-P(J)+COST(QARCB))-EPS2).GT.0.0) THEN C ************* TAKING CARE OF ARCS WITH NON LINEAR COST ****************** RESID=F(QARCB)-MIN(U(QARCB),DCINV(QARCB,P(J)-GAIN $ *PRICE)) c if (f(QARCb).gt.u(QARCb)-ref_thr) then c print*,'4',qarcb,f(QARCb),u(QARCb) c end if IF (-RESID.GT.REF_THR) THEN IF (-SURPI .GE. -GAIN*RESID) THEN F(QARCB) = F(QARCB)-RESID SURPI = SURPI - GAIN*RESID SURPLUS(J) = SURPLUS(J) + RESID COST(QARCB) = DC(QARCB,F(QARCB)) QARCB = QNXTPUSHB(QARCB) IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF ELSE F(QARCB) = F(QARCB) - SURPI/GAIN COST(QARCB) = DC(QARCB,F(QARCB)) SURPLUS(J) = SURPLUS(J) + SURPI/GAIN SURPI = 0.0 C C ***** IF FLOW PUSH IS NON-SATURATING, STORE QARCB AND CHECK FOR FLOW PUSH C ALONG A CYCLE C IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (FLAG) THEN IF (PRDCSR(J).EQ.0) THEN PRDCSR(J)=-QARCB ELSE CALL CYCLE_PUSH2(I,J,-QARCB) FLAG = .FALSE. END IF ELSE PRDCSR(J) = NA+1 FLAG = .TRUE. END IF ELSE IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF END IF END IF IF (-SURPLUS(J).GT.THRESHSURPLUS) THEN IF (NXTQUEUE(J).EQ.0) THEN NXTQUEUE(PREVNODE)=J NXTQUEUE(J)=I PREVNODE=J END IF END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF ELSE QARCB = QNXTPUSHB(QARCB) END IF GOTO 523 END IF END IF 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)+SURPI) .GT. 1.d0 ) THEN c PRINT*,'SURPLUS ERR AT',I,-DFCT(I),SURPI c stop c END IF C ***** END OF D-PUSHES; CHECK IF PRICE DROP IS NEEDED ***** c if ((i.eq.14) .and. (-surpi .gt. 0)) then C print*,'3',arcf,arcb,qarcf,qarcb,price,surpi c end if IF ((-SURPI.GT.THRESHSURPLUS).AND. $ (ARCF+ARCB+QARCF+QARCB.EQ.0)) THEN IF (FLAG) THEN FLAG = .FALSE. CALL CLR_PRDCSR(I) END IF 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) c if ((i.eq.12) .and. (arc.eq. 61)) then c print*,'3.5',xp,price,eps2,c(arc),g(arc) c end if C THE LINES BELOW DIFFERS FROM LAKIS' EVEN FOR ORDINARY NETWORKS 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' C PAUSE 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' C PAUSE STOP END IF END IF ARC = NXTIN(ARC) GOTO 612 END IF IF (-SURPI.GT.REF_THR) GOTO 400 IF (-PRICE.LT.IDINT(LARGE/10)) PRICE=-IDINT(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. IF (-SURPI.GT.THRESHSURPLUS) GOTO 515 ELSE C IF NO PRICE RISE TOOK PLACE RESET THE START OF THE PUSH LISTS FPUSHF(I) = ARCF FPUSHB(I) = ARCB QFPUSHF(I) = QARCF QFPUSHB(I) = QARCB END IF C DO THE FINAL BOOKKEEPING OF THE ITERATION SURPLUS(I) = SURPI 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 ********** END OF DOWN ITERATION STARTING AT NODE I ********* C CHECK FOR THE END OF THE QUEUE. THIS STEP IS NOT NECESSARY FOR C THE ALGORITHM; IT IS INCLUDED FOR COLLECTING STATISTICS ABOUT C THE ALGORITHM'S BEHAVIOR, E.G. COUNTING THE NUMBER OF CYCLES 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 write(*,444)(surplus(k),k=1,n) c write(*,444)(p(k),k=1,n) c444 format(10e8.1) 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 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 END IF c IF (-DFCT(K) .GT. THRESHSURPLUS + REF_THR) THEN c PRINT*,'SURPLUS TOO POSITIVE AT',K,-DFCT(K) c STOP c END IF c DEL=DEL+DABS(DFCT(K)) c635 CONTINUE c DO 639 K=1,N c IF (SURPLUS(K) .GT. THRESHSURPLUS + REF_THR) THEN c PRINT*,'SURPLUS TOO POSITIVE AT',K,SURPLUS(K) c STOP c END IF c639 CONTINUE C PRINT*,'i=',i,p(i),surplus(i), DEL NXTNODE=NXTQUEUE(I) IF (I.NE.NXTNODE) THEN NXTQUEUE(I)=0 NXTQUEUE(PREVNODE)=NXTNODE C C IF LAST OPERATION IS A NONSATURATING FLOW PUSH TOWARDS J, ITERATE ON J NEXT C IF (FLAG) THEN K=0 640 IF (NXTNODE.NE.J) THEN PREVNODE=NXTNODE NXTNODE=NXTQUEUE(NXTNODE) K=K+1 GO TO 640 END IF I=J END IF I=NXTNODE GO TO 500 END IF C ************** END OF SUBPROBLEM (SCALING PHASE) ************** c 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 ' 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,DELTA REAL*8 REF_THR print*,'do cycle push1 at',node,j,iarc C print*,' ' REF_THR = 1.0E-10 IF (IARC.GT.0) THEN GAIN1=G(IARC) DELTA=U(IARC)-F(IARC) ELSE GAIN1=1.D0/G(-IARC) DELTA=F(-IARC)/G(-IARC) END IF I=NODE C print*,i,iarc,delta,gain1 ARC = PRDCSR(I) 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) 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 DELTA = MIN(DELTA,RESID)*G(ARC) GAIN1 = GAIN1/G(ARC) END IF c print*,i,k,j,arc,delta,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 IF (GAIN1.LT. 1.D0 - REF_THR) THEN DELTA = MIN(DELTA,SURPLUS(J)/(1.D0-GAIN1)) END IF C print*,'push1',node,j,iarc,surplus(j) C C **** DO FLOW PUSH ALONG THE CYLE BY THE AMOUNT DELTA C SURPLUS(J)=SURPLUS(J) - DELTA*(1.D0-GAIN1) IF (IARC.GT.0) THEN GAIN1=GAIN1/G(IARC) F(IARC)=F(IARC)+DELTA*GAIN1 ELSE IARC=-IARC F(IARC)=F(IARC)-DELTA*GAIN1 GAIN1=GAIN1*G(IARC) END IF COST(IARC) = DC(IARC,F(IARC)) I=NODE ARC = PRDCSR(I) 20 IF (ARC.GT.0) THEN K = STARTN(ARC) GAIN1=GAIN1/G(ARC) F(ARC)=F(ARC)+DELTA*GAIN1 ELSE ARC = -ARC K = ENDN(ARC) F(ARC)=F(ARC)-DELTA*GAIN1 GAIN1=GAIN1*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 20 END IF C C ***** SET PREDECESSOR FOR ALL REMAINING NODES TO ZERO C I = K ARC = PRDCSR(I) 30 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 30 END IF PRDCSR(I)=0 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,DELTA REAL*8 REF_THR print*,'do cycle push2 at',node,j,iarc REF_THR = 1.0E-10 IF (IARC.GT.0) THEN GAIN1=G(IARC) DELTA=F(IARC) ELSE GAIN1=1.D0/G(-IARC) DELTA=(U(-IARC)-F(-IARC))/G(-IARC) END IF I=NODE ARC = PRDCSR(I) 10 IF (ARC.GT.0) THEN K = STARTN(ARC) C print*,CDAT(ARC,1) C print*,delta,arc,f(arc),u(arc),P(K)-G(ARC)*P(I) 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) GAIN1 = GAIN1*G(ARC) ELSE ARC = -ARC K = ENDN(ARC) C print*,CDAT(ARC,1) C print*,arc,f(arc),u(arc),P(i)-G(ARC)*P(k),a(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) GAIN1 = GAIN1/G(ARC) END IF 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 IF (GAIN1.LT. 1.D0 - REF_THR) THEN DELTA = MIN(DELTA,-SURPLUS(J)/(1.D0-GAIN1)) END IF c print*,'push2',node,j,iarc,surplus(j),delta,gain1 C C **** DO FLOW PUSH ALONG THE CYLE BY THE AMOUNT DELTA C SURPLUS(J)=SURPLUS(J) + DELTA*(1.D0-GAIN1) IF (IARC.GT.0) THEN GAIN1=GAIN1/G(IARC) F(IARC)=F(IARC)-DELTA*GAIN1 ELSE IARC=-IARC F(IARC)=F(IARC)+DELTA*GAIN1 GAIN1=GAIN1*G(IARC) END IF COST(IARC) = DC(IARC,F(IARC)) I=NODE ARC = PRDCSR(I) 20 IF (ARC.GT.0) THEN K = STARTN(ARC) GAIN1=GAIN1/G(ARC) F(ARC)=F(ARC)-DELTA*GAIN1 ELSE ARC = -ARC K = ENDN(ARC) F(ARC)=F(ARC)+DELTA*GAIN1 GAIN1=GAIN1*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 c print*,'*',i,k,arc,j,iarc C if (i.eq.0) stop IF (K.NE.J) THEN I = K ARC = PRDCSR(I) GO TO 20 END IF C C ***** SET PREDECESSOR FOR ALL REMAINING NODES TO ZERO C I = K ARC = PRDCSR(I) 30 IF (ARC.NE.NA+1) THEN IF (ARC.GT.0) THEN K=STARTN(ARC) ELSE K=ENDN(-ARC) END IF PRDCSR(I)=0 c print*,'&',i,k,arc,j,iarc C if (i.eq.0) stop I=K ARC=PRDCSR(I) GO TO 30 END IF PRDCSR(I)=0 c print*,'cycle push done!',surplus(j) END C C SUBROUTINE CLR_PRDCSR(NODE) IMPLICIT NONE INTEGER*4 NODE 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 LARGE,ZERO COMMON /SCALARS/ N,NA,LARGE,ZERO COMMON /BLK1/ STARTN COMMON /BLK2/ ENDN COMMON /BLK8/ PRDCSR INTEGER*4 I,K,ARC C print*,'clrprdcsr',node ARC=PRDCSR(NODE) I=NODE C C ***** SET PREDECESSOR FOR ALL NODES TO ZERO C 30 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 30 END IF PRDCSR(I)=0 c print*,'clr done!' C END