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
