C....*...1.........2.........3.........4.........5.........6.........7.*.......8 C DLP 12/20/73 C C PURPOSE C SOLVE THE LINEAR PROGRAM: MAX C'X SUBJECT TO AX.LE.B AND X.GE.0. C C USAGE C CALL DLP(A,B,C,N,M,X,Y,D,EPS,IDX,IER) C C SUBROUTINES CALLED C DRAN C C ARGUMENTS C A - INPUT N BY M COEFFICIENT MATRIX STORED COLUMNWISE (STORAGE C MODE 0). DESTROYED ON RETURN. C REAL*8 C B - INPUT N BY 1 VECTOR CONTAINING THE RIGHT HAND SIDE OF THE C CONSTRAINTS AX.LE.B. DESTROYED ON RETURN. C REAL*8 C C - INPUT M BY 1 VECTOR CONTAINING THE COEFFICIENTS OF THE C OBJECTIVE FUNCTION C'X. DESTROYED ON RETURN. C REAL*8 C N - NUMBER OF INEQUALITY CONSTRAINTS. C INTEGER*4 C M - LENGTH OF THE VECTOR VALUED VARIABLE X. C INTEGER*4 C X - OUTPUT SOLUTION VECTOR OF LENGTH M IF THE PROGRAM IS C FEASIBLE AND BOUNDED. C REAL*8 C Y - OUTPUT SOLUTION VECTOR OF LENGTH N OF THE DUAL LINEAR C PROGRAM IF THE PROGRAM IS FEASIBLE AND BOUNDED. C REAL*8 C D - VALUE OF THE OBJECTIVE FUNCTION C'X PROVIDED THE PROGRAM IS C FEASIBLE AND BOUNDED. C REAL*8 C EPS - INPUT VALUE. VALUES LESS THAN TOL=EPS*MAX(A(I,J)) ARE C TAKEN TO BE ZERO. A REASONABLE CHOICE OF EPS IS 1.D-13. C REAL*8 C IDX - WORK VECTOR OF LENGTH M+N. C INTEGER*4 C IER - ERROR PARAMETER CODED AS FOLLOWS: C 0 - NO ERROR. PROGRAM IS FEASIBLE AND BOUNDED. C 1 - PROGRAM IS FEASIBLE BUT UNBOUNDED. C 2 - PROGRAM IS INFEASIBLE. C INTEGER*4 C C PROGRAMMER C DR. A. RONALD GALLANT C DEPARTMANT OF STATISTICS C NORTH CAROLINA STATE UNIVERSITY C RALEIGH, NORTH CAROLINA 27607 C C REFERENCE C G. OWEN. GAME THEORY. PHILADELPHIA: W.B. SAUNDERS, 1968. C C SUBROUTINE DLP(A,B,C,N,M,X,Y,D,EPS,IDX,IER) IMPLICIT REAL*8 (A-H,O-Z) save REAL*8 A(N,M),B(N),C(M),X(M),Y(N),D INTEGER*4 IDX(1) IER=3 IF(N.LE.0) RETURN IF(M.LE.0) RETURN CALL ZINDEX(1,I0,J0,X,Y,N,M,IDX,B,C) IX=191919 DO 10 I=1,N 10 B(I)=-B(I) D=0.D0 TOL=0.D0 DO 20 I=1,N DO 20 J=1,M TEST=DABS(A(I,J)) 20 IF(TEST.GT.TOL) TOL=TEST TOL=EPS*TOL C C OBTAIN A BASIC FEASIBLE POSITION (B.F.P) C 100 CONTINUE DO 110 K=1,N K0=N-(K-1) IF(B(K0).GT.TOL) GO TO 115 110 CONTINUE C FALLING THRU LOOP 110 IMPLYS B.F.P. OBTAINED. GO TO 200 115 DO 120 J=1,M L0=J IF(A(K0,J).LT.-TOL) GO TO 125 120 CONTINUE C FALLING THRU LOOP 120 IMPLYS B(K0) IS POSITIVE AND ALL ELEMENTS C OF ROW K0 OF A ARE POSTIVE. HENCE PRIMAL IS INFEASIBLE. C (NOTE LOOP 10) IER=2 RETURN 125 CALL ZPICK(1,X,M,IX,J0) DO 150 J=L0,M IF(A(K0,J).GE.-TOL) GO TO 150 CALL ZPICK(2,X,M,IX,J) TEST=B(K0)/A(K0,L0) DO 130 I=K0,N RATIO=0.D0 AIJ=A(I,J) IF(DABS(AIJ).GT.TOL) RATIO=B(I)/AIJ IF((RATIO.LT.-TOL).AND. (RATIO.GT.TEST)) TEST=RATIO Y(I)=RATIO 130 IF(RATIO.EQ.TEST) I0=I ISW=0 DO 140 I=K0,N 140 IF(DABS(Y(I)-TEST).LT.TOL) ISW=ISW+1 J0=J IF(ISW.EQ.1) GO TO 170 150 CONTINUE C FALLING THRU LOOP 150 IMPLYS THERE IS NO COLUMN WITH A(K0,J)>0 C AND A UNIQUE MAXIMUM OF THOSE B(I)/A(I,J)<0. (NOTE LOOP 10). C A COLUMN WITH A(K0,J)>0 IS CHOSEN RANDOMLY. CALL ZPICK(3,X,M,IX,J0) TEST=B(K0)/A(K0,J0) DO 160 I=K0,N RATIO=0.D0 AIJ0=A(I,J0) IF(DABS(AIJ0).GT.TOL) RATIO=B(I)/AIJ0 IF((RATIO.LT.-TOL).AND.(RATIO.GT.TEST)) TEST=RATIO 160 CONTINUE CALL ZPICK(1,Y,N,IX,I0) DO 165 I=K0,N RATIO=0.D0 AIJ0=A(I,J0) IF(DABS(AIJ0).GT.TOL) RATIO=B(I)/AIJ0 165 IF(DABS(RATIO-TEST).LT.TOL) CALL ZPICK(2,Y,N,IX,I) CALL ZPICK(3,Y,N,IX,I0) 170 CALL ZPIVOT(I0,J0,A,B,C,D,N,M) CALL ZINDEX(2,I0,J0,X,Y,N,M,IDX,B,C) GO TO 100 C C IMPROVE OBJECTIVE FUNCTION AND MAINTAIN A B.F.P. C 200 CONTINUE DO 210 J=1,M L0=J IF(C(L0).GE.TOL) GO TO 220 210 CONTINUE C FALLING THRU LOOP 210 IMPLYS SOLUTION OBTAINED GO TO 300 220 CALL ZPICK(1,X,M,IX,J0) DO 250 J=L0,M IF(C(J ).LT.TOL) GO TO 250 CALL ZPICK(2,X,M,IX,J) IER=1 TEST=-1.D62 DO 230 I=1,N RATIO=0.D0 AIJ=A(I,J) IF(DABS(AIJ).GT.TOL) RATIO=B(I)/AIJ IF(AIJ.GT.TOL) IER=0 IF((RATIO.LT.-TOL).AND.(RATIO.GT.TEST)) TEST=RATIO Y(I)=RATIO 230 IF(RATIO.EQ.TEST) I0=I C IF IER=1 THEN ALL A(I,J) IN COLUMN J ARE NEGATIVE, C(J) IS C POSITIVE, AND WE HAVE A B.F.P. SITUATION HENCE PRIMAL IS UNBOUNDED IF(IER.EQ.1) RETURN ISW=0 DO 240 I=1,N 240 IF(DABS(Y(I)-TEST).LT.TOL) ISW=ISW+1 J0=J IF(ISW.EQ.1) GO TO 270 250 CONTINUE C FALLING THRU LOOP 250 IMPLYS THERE IS NO COLUMN WITH C(J)>0 C AND A UNIQUE MAXIMUM OF THOSE B(I)/A(I,J)<0. (NOTE LOOP 10). C A COLUMN WITH C(J)>0 IS CHOSEN RANDOMLY. CALL ZPICK(3,X,M,IX,J0) TEST=-1.D62 DO 260 I=1,N RATIO=0.D0 AIJ0=A(I,J0) IF(DABS(AIJ0).GT.TOL) RATIO=B(I)/AIJ0 IF((RATIO.LT.-TOL).AND.(RATIO.GT.TEST)) TEST=RATIO 260 CONTINUE CALL ZPICK(1,Y,N,IX,I0) DO 265 I=1,N RATIO=0.D0 AIJ0=A(I,J0) IF(DABS(AIJ0).GT.TOL) RATIO=B(I)/AIJ0 265 IF(DABS(RATIO-TEST).LT.TOL) CALL ZPICK(2,Y,N,IX,I) CALL ZPICK(3,Y,N,IX,I0) 270 CALL ZPIVOT(I0,J0,A,B,C,D,N,M) CALL ZINDEX(2,I0,J0,X,Y,N,M,IDX,B,C) GO TO 100 300 CALL ZINDEX(3,I0,J0,X,Y,N,M,IDX,B,C) RETURN END SUBROUTINE ZPICK(ISW,INDEX,LNGTH,IX,K0) implicit real*8 (a-h,o-z) save INTEGER*4 INDEX(1) IF(ISW.EQ.1) GO TO 100 IF(ISW.EQ.2) GO TO 200 IF(ISW.EQ.3) GO TO 300 C C INITIALIZE INDEX C 100 DO 110 K=1,LNGTH 110 INDEX(K)=0 RETURN C C ENTER POSSIBLE PIVOT CHOICE C 200 INDEX(K0)=1 RETURN C C CHOOSE A PIVOT RANDOMLY C 300 CONTINUE K0=0 ISUM=0 DO 310 K=1,LNGTH II=INDEX(K) IF(II.EQ.1) K0=K 310 ISUM=ISUM+II IF(ISUM.EQ.1) RETURN c CALL RANDU(IX,IY,P) c IX=IY p=dran(ix) IP=P*ISUM ISUM=0 DO 320 K=1,LNGTH II=INDEX(K) IF(II.EQ.1) K0=K ISUM=ISUM+II 320 IF(ISUM.GT.IP) RETURN RETURN END SUBROUTINE ZPIVOT(I0,J0,A,B,C,D,N,M) IMPLICIT REAL*8 (A-H,O-Z) save REAL*8 A(N,M),B(N),C(M) C REMARK: A IS OWEN'S A', B IS HIS B', AND C HIS C'. PIVOT=A(I0,J0) DO 10 I=1,N DUMMY=A(I,J0)/PIVOT DO 10 J=1,M 10 IF((I.NE.I0).AND.(J.NE.J0)) A(I,J)=A(I,J)-A(I0,J)*DUMMY DO 20 I=1,N 20 IF(I.NE.I0) B(I)=B(I)-B(I0)*A(I,J0)/PIVOT DO 30 J=1,M 30 IF(J.NE.J0) C(J)=C(J)-A(I0,J)*C(J0)/PIVOT D=D-B(I0)*C(J0)/PIVOT DO 40 I=1,N 40 IF(I.NE.I0) A(I,J0)=-A(I,J0)/PIVOT DO 50 J=1,M 50 IF(J.NE.J0) A(I0,J)=A(I0,J)/PIVOT B(I0)=B(I0)/PIVOT C(J0)=-C(J0)/PIVOT A(I0,J0)=1.D0/PIVOT RETURN END SUBROUTINE ZINDEX(ISW,I0,J0,X,Y,N,M,IDX,B,C) IMPLICIT REAL*8 (A-H,O-Z) save REAL*8 X(1),Y(1),B(1),C(1) INTEGER*4 IDX(1) MPLUSN=M+N C C PERMUTATION VECTOR CORRESONDENCES C C PRIMAL: C IDX(I).LE.M IMPLYS X VARIABLE IS INDEXED C IDX(I).GT.M IMPLYS SLACK VARIABLE IS INDEXED C IDX(I) (I=1,...,M) REGERS TO A NON-BASIC VARIABLE (TOP BORDER) C IDX(I) (I=M+1,...,M+N) REFERS TO A BASIC VARIABLE (RIGHT BORDER) C C DUAL: C IDX(I).LE.M IMPLYS SLACK VARIABLE IS INDEXED C IDX(I).GT.M.IMPLYS Y VARIABLE IS INDEXED C IDX(I) (I=1,...,M) REFERS TO A BASIC VARIABLE (BOTTOM BORDER) C IDX(I) (I=M+1,...,M+N) REFERS TO A NON-BASIC VARIABLE (LEFT) C C IF(ISW.EQ.1) GO TO 100 IF(ISW.EQ.2) GO TO 200 IF(ISW.EQ.3) GO TO 300 C C INTIALIZE THE PERMUTATION VECTOR. C 100 DO 110 I=1,MPLUSN 110 IDX(I)=I RETURN C C BASIC VARIABLE IDX(M+I0) BECOMES NON-BASIC VARIABLE IDX(J0) C AND VICE-VERSA. C 200 ITEMP=IDX(M+I0) 3000 FORMAT('0',10I5) IDX(M+I0)=IDX(J0) IDX(J0)=ITEMP RETURN C C ASSIGN X AND Y THEIR VALUES C 300 DO 310 J=1,M IF(IDX(J).LE.M) X(IDX(J))=0.D0 310 IF(IDX(J).GT.M) Y(IDX(J)-M)=-C(J) DO 320 I=1,N II=M+I IF(IDX(II).LE.M) X(IDX(II))=-B(I) 320 IF(IDX(II).GT.M) Y(IDX(II)-M)=0.D0 RETURN END