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
