C/ FIGURE 4.6.2
      SUBROUTINE DLPRVS(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 AB(NZ),CC(N+M,2),YY(M,2),D(N+M,2),
     & V(M),XB(M),AP(M)
      INTEGER BASIS(M),IROWB(NZ),JCOLB(NZ)
C
C  SUBROUTINE DLPRVS 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  NOTE: DLPRVS CALLS A SPARSE LINEAR SYSTEM SOLVER
C
C          SUBROUTINE SPARSOL(A,IROW,JCOL,NZ,X,B,N)
C          DOUBLE PRECISION A(NZ),X(N),B(N)
C          INTEGER IROW(NZ),JCOL(NZ)
C
C  TO SOLVE THE N BY N SPARSE SYSTEM AX=B, WHERE A(IZ), IZ=1,...,NZ,
C  IS THE NONZERO ELEMENT OF A IN ROW IROW(IZ), COLUMN JCOL(IZ), B IS 
C  THE RIGHT HAND SIDE VECTOR AND X IS THE SOLUTION.
C
C  THE CALL TO SPARSOL SHOULD BE REPLACED BY A CALL TO A DIRECT LINEAR 
C  SYSTEM SOLVER (ITERATIVE SOLVERS NOT RECOMMENDED).  
C          
C-----------------------------------------------------------------------
C
C                              EPS = MACHINE FLOATING POINT RELATIVE
C                                    PRECISION
C *****************************
      DATA EPS/2.D-16/
C *****************************
C                              INITIALIZE Ab TO IDENTITY
      NZB = M
      DO 5 I=1,M
         IROWB(I) = I
         JCOLB(I) = I
         AB(I) = 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                              SOLVE Ab*V = AP
            CALL SPARSOL(AB,IROWB,JCOLB,NZB,V,AP,M)
C                              FIND SMALLEST POSITIVE RATIO
C                              Xb(I)/V(I), IDENTIFYING PIVOT ROW IP
            RATMIN = 0.0 
            IP = 0
            DO 80 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
   80       CONTINUE
C                              IF ALL RATIOS NONPOSITIVE, MAXIMUM
C                              IS UNBOUNDED
            IF (IP.EQ.0) THEN
               PRINT 85 
   85          FORMAT (' ***** UNBOUNDED MAXIMUM *****')
               RETURN
            ENDIF
C                              ADD X(JP) TO BASIS
            BASIS(IP) = JP
C                              UPDATE Ab.  PUT NONZEROS OF JP-TH  
C                              COLUMN OF A (=AP) INTO COLUMN IP
C                              OF SPARSE MATRIX Ab
            NZBOLD = NZB
            DO 90 I=1,M
               IF (AP(I).NE.0.0) THEN
                  NZB = NZB+1 
                  IROWB(NZB) = I 
                  JCOLB(NZB) = IP 
                  AB(NZB) = AP(I)
               ENDIF
   90       CONTINUE
            NZBNEW = NZB
C                              REMOVE ELEMENTS OF OLD COLUMN IP 
            NZB = 0
            DO 95 IZ=1,NZBNEW
               IF (JCOLB(IZ).NE.IP .OR. IZ.GT.NZBOLD) THEN
                  NZB = NZB+1
                  JCOLB(NZB) = JCOLB(IZ)
                  IROWB(NZB) = IROWB(IZ)
                  AB(NZB) = AB(IZ)
               ENDIF
   95       CONTINUE
C                              SOLVE Ab*Xb = B
            CALL SPARSOL(AB,IROWB,JCOLB,NZB,XB,B,M)
C                              SOLVE Ab**T*Y = Cb
            DO 105 IR=1,MROW
               DO 100 J=1,M
                  K = BASIS(J)
                  V(J) = CC(K,IR)
  100          CONTINUE 
               CALL SPARSOL(AB,JCOLB,IROWB,NZB,YY(1,IR),V,M)
  105       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
