C/ FIGURE 4.6.1 SUBROUTINE DLPRV(A,IROW,JCOL,NZ,B,C,N,M,P,X,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DECLARATIONS FOR ARGUMENTS DOUBLE PRECISION A(NZ),B(M),C(N),P,X(N),Y(M) INTEGER N,M,NZ,IROW(NZ),JCOL(NZ) C DECLARATIONS FOR LOCAL VARIABLES DOUBLE PRECISION ABINV(M,M),CC(N+M,2),YY(M,2),D(N+M,2), & V(M),XB(M),AP(M) INTEGER BASIS(M) C C SUBROUTINE DLPRV USES THE REVISED SIMPLEX METHOD TO SOLVE THE PROBLEM C C MAXIMIZE P = C(1)*X(1) + ... + C(N)*X(N) C C WITH X(1),...,X(N) NONNEGATIVE, AND C C A(1,1)*X(1) + ... + A(1,N)*X(N) = B(1) C . . . C . . . C A(M,1)*X(1) + ... + A(M,N)*X(N) = B(M) C C WHERE B(1),...,B(M) ARE ASSUMED TO BE NONNEGATIVE. C C ARGUMENTS C C ON INPUT ON OUTPUT C -------- --------- C C A - A(IZ) IS THE CONSTRAINT MATRIX C ELEMENT IN ROW IROW(IZ), COLUMN C JCOL(IZ), FOR IZ=1,...,NZ. C C IROW - (SEE A). C C JCOL - (SEE A). C C NZ - NUMBER OF NONZEROS IN A. C C B - A VECTOR OF LENGTH M CONTAINING C THE RIGHT HAND SIDES OF THE C CONSTRAINTS. THE COMPONENTS OF C B MUST ALL BE NONNEGATIVE. C C C - A VECTOR OF LENGTH N CONTAINING C THE COEFFICIENTS OF THE OBJECTIVE C FUNCTION. C C N - THE NUMBER OF UNKNOWNS. C C M - THE NUMBER OF CONSTRAINTS. C C P - THE MAXIMUM OF THE C OBJECTIVE FUNCTION. C C X - A VECTOR OF LENGTH N C WHICH CONTAINS THE LP C SOLUTION. C C Y - A VECTOR OF LENGTH M C WHICH CONTAINS THE DUAL C SOLUTION. C C----------------------------------------------------------------------- C C EPS = MACHINE FLOATING POINT RELATIVE C PRECISION C ***************************** DATA EPS/2.D-16/ C ***************************** C INITIALIZE Ab**(-1) TO IDENTITY DO 5 I=1,M DO 5 J=1,M ABINV(I,J) = 0.0 IF (I.EQ.J) ABINV(I,J) = 1.0 5 CONTINUE C OBJECTIVE FUNCTION COEFFICIENTS ARE C CC(I,1) + CC(I,2)*ALPHA, WHERE "ALPHA" C IS TREATED AS INFINITY DO 10 I=1,N+M CC(I,1) = 0.0 CC(I,2) = 0.0 IF (I.LE.N) THEN CC(I,1) = C(I) ELSE CC(I,2) = -1 ENDIF 10 CONTINUE C BASIS(1),...,BASIS(M) HOLD NUMBERS OF C BASIS VARIABLES. INITIAL BASIS CONSISTS C OF ARTIFICIAL VARIABLES ONLY DO 15 I=1,M K = N+I BASIS(I) = K C INITIALIZE Y TO Ab**(-T)*Cb = Cb YY(I,1) = CC(K,1) YY(I,2) = CC(K,2) C INITIALIZE Xb TO Ab**(-1)*B = B XB(I) = B(I) IF (B(I).LT.0.0) THEN PRINT 12 12 FORMAT (' ***** ALL B(I) MUST BE NONNEGATIVE *****') RETURN ENDIF 15 CONTINUE C SIMPLEX METHOD CONSISTS OF TWO PHASES DO 130 IPHASE=1,2 IF (IPHASE.EQ.1) THEN C PHASE I: ROW 2 OF D (WITH COEFFICIENTS OF C ALPHA) SEARCHED FOR MOST NEGATIVE ENTRY MROW = 2 LIM = N+M ELSE C PHASE II: FIRST N ELEMENTS OF ROW 1 OF C D SEARCHED FOR MOST NEGATIVE ENTRY C (COEFFICIENTS OF ALPHA NONNEGATIVE NOW) MROW = 1 LIM = N C IF ANY ARTIFICIAL VARIABLES LEFT IN C BASIS AT BEGINNING OF PHASE II, THERE C IS NO FEASIBLE SOLUTION DO 25 I=1,M IF (BASIS(I).GT.LIM) THEN PRINT 20 20 FORMAT (' ***** NO FEASIBLE SOLUTION *****') RETURN ENDIF 25 CONTINUE ENDIF C THRESH = SMALL NUMBER. WE ASSUME SCALES C OF A AND C ARE NOT *TOO* DIFFERENT THRESH = 0.0 DO 30 J=1,LIM THRESH = MAX(THRESH,ABS(CC(J,MROW))) 30 CONTINUE THRESH = 1000*EPS*THRESH C BEGINNING OF SIMPLEX STEP 35 CONTINUE C D**T = Y**T*A - C**T DO 55 IR=1,MROW DO 40 J=1,N+M D(J,IR) = -CC(J,IR) 40 CONTINUE C LAST M COLUMNS OF A FORM IDENTITY MATRIX DO 45 J=1,M D(N+J,IR) = D(N+J,IR) + YY(J,IR) 45 CONTINUE C FIRST N COLUMNS STORED IN SPARSE A MATRIX DO 50 IZ=1,NZ I = IROW(IZ) J = JCOL(IZ) D(J,IR) = D(J,IR) + A(IZ)*YY(I,IR) 50 CONTINUE 55 CONTINUE C FIND MOST NEGATIVE ENTRY OF ROW MROW C OF D, IDENTIFYING PIVOT COLUMN JP CMIN = -THRESH JP = 0 DO 60 J=1,LIM IF (D(J,MROW).LT.CMIN) THEN CMIN = D(J,MROW) JP = J ENDIF 60 CONTINUE C IF ALL ENTRIES NONNEGATIVE (ACTUALLY, C IF GREATER THAN -THRESH) PHASE ENDS IF (JP.EQ.0) GO TO 130 C COPY JP-TH COLUMN OF A ONTO VECTOR AP IF (JP.LE.N) THEN C JP-TH COLUMN IS PART OF SPARSE A MATRIX DO 65 I=1,M AP(I) = 0 65 CONTINUE DO 70 IZ=1,NZ J = JCOL(IZ) IF (J.EQ.JP) THEN I = IROW(IZ) AP(I) = A(IZ) ENDIF 70 CONTINUE ELSE C JP-TH COLUMN IS COLUMN OF FINAL IDENTITY DO 75 I=1,M AP(I) = 0 75 CONTINUE AP(JP-N) = 1 ENDIF C V = Ab**(-1)*AP DO 85 I=1,M V(I) = 0.0 DO 80 J=1,M V(I) = V(I) + ABINV(I,J)*AP(J) 80 CONTINUE 85 CONTINUE C FIND SMALLEST POSITIVE RATIO C Xb(I)/V(I), IDENTIFYING PIVOT ROW IP RATMIN = 0.0 IP = 0 DO 90 I=1,M IF (V(I).GT.THRESH) THEN RATIO = XB(I)/V(I) IF (IP.EQ.0 .OR. RATIO.LT.RATMIN) THEN RATMIN = RATIO IP = I ENDIF ENDIF 90 CONTINUE C IF ALL RATIOS NONPOSITIVE, MAXIMUM C IS UNBOUNDED IF (IP.EQ.0) THEN PRINT 95 95 FORMAT (' ***** UNBOUNDED MAXIMUM *****') RETURN ENDIF C ADD X(JP) TO BASIS BASIS(IP) = JP C UPDATE Ab**(-1) = E**(-1)*Ab**(-1) C Xb = E**(-1)*Xb DO 100 J=1,M ABINV(IP,J) = ABINV(IP,J)/V(IP) 100 CONTINUE XB(IP) = XB(IP)/V(IP) DO 110 I=1,M IF (I.EQ.IP) GO TO 110 DO 105 J=1,M ABINV(I,J) = ABINV(I,J) - V(I)*ABINV(IP,J) 105 CONTINUE XB(I) = XB(I) - V(I)*XB(IP) 110 CONTINUE C CALCULATE Y = Ab**(-T)*Cb DO 125 IR=1,MROW DO 120 I=1,M YY(I,IR) = 0.0 DO 115 J=1,M K = BASIS(J) YY(I,IR) = YY(I,IR) + ABINV(J,I)*CC(K,IR) 115 CONTINUE 120 CONTINUE 125 CONTINUE GO TO 35 C END OF SIMPLEX STEP 130 CONTINUE C END OF PHASE II; CALCULATE X DO 135 J=1,N X(J) = 0.0 135 CONTINUE DO 140 I=1,M K = BASIS(I) X(K) = XB(I) Y(I) = YY(I,1) 140 CONTINUE C CALCULATE P P = 0.0 DO 145 I=1,N P = P + C(I)*X(I) 145 CONTINUE RETURN END