% FIGURE 4.6.2
      function [P,X,Y] = DLPRVS(A,IROW,JCOL,NZ,B,C,N,M)
%
%  FUNCTION DLPRVS USES THE REVISED SIMPLEX METHOD TO SOLVE THE PROBLEM
%
%             MAXIMIZE     P = C(1)*X(1) + ... + C(N)*X(N)
%
%    WITH X(1),...,X(N) NONNEGATIVE, AND
%
%           A(1,1)*X(1) + ... + A(1,N)*X(N)  = B(1)
%             .                   .             .
%             .                   .             .
%           A(M,1)*X(1) + ... + A(M,N)*X(N)  = B(M)
%
%    WHERE B(1),...,B(M) ARE ASSUMED TO BE NONNEGATIVE.
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - A(IZ) IS THE CONSTRAINT MATRIX 
%             ELEMENT IN ROW IROW(IZ), COLUMN 
%             JCOL(IZ), FOR IZ=1,...,NZ.
%
%    IROW   - (SEE A).
%
%    JCOL   - (SEE A).
%
%    NZ     - NUMBER OF NONZEROS IN A. 
%
%    B      - A VECTOR OF LENGTH M CONTAINING
%             THE RIGHT HAND SIDES OF THE
%             CONSTRAINTS.  THE COMPONENTS OF
%             B MUST ALL BE NONNEGATIVE.
%
%    C      - A VECTOR OF LENGTH N CONTAINING
%             THE COEFFICIENTS OF THE OBJECTIVE
%             FUNCTION.
%
%    N      - THE NUMBER OF UNKNOWNS.
%
%    M      - THE NUMBER OF CONSTRAINTS.
%
%    P      -                                   THE MAXIMUM OF THE
%                                               OBJECTIVE FUNCTION.
%
%    X      -                                   A VECTOR OF LENGTH N
%                                               WHICH CONTAINS THE LP
%                                               SOLUTION.
%
%    Y      -                                   A VECTOR OF LENGTH M
%                                               WHICH CONTAINS THE DUAL
%                                               SOLUTION.
%
%  NOTE: DLPRVS CALLS A SPARSE LINEAR SYSTEM SOLVER
%
%          X = SPARSOL(A,IROW,JCOL,NZ,B,N)
%
%  TO SOLVE THE N BY N SPARSE SYSTEM AX=B, WHERE A(IZ), IZ=1,...,NZ,
%  IS THE NONZERO ELEMENT OF A IN ROW IROW(IZ), COLUMN JCOL(IZ), B IS 
%  THE RIGHT HAND SIDE VECTOR AND X IS THE SOLUTION.
%
%  THE CALL TO SPARSOL SHOULD BE REPLACED BY A CALL TO A DIRECT LINEAR 
%  SYSTEM SOLVER (ITERATIVE SOLVERS NOT RECOMMENDED).
%
%-----------------------------------------------------------------------
%
%                              EPS = MACHINE FLOATING POINT RELATIVE
%                                    PRECISION
% *****************************
      EPS = eps;
% *****************************
      P = 0;
      X = zeros(N,1);
      Y = zeros(M,1);
%                              INITIALIZE Ab TO IDENTITY
      NZB = M;
      for I=1:M
         IROWB(I) = I;
         JCOLB(I) = I;
         AB(I) = 1.0;
      end
%                              OBJECTIVE FUNCTION COEFFICIENTS ARE
%                              CC(I,1) + CC(I,2)*ALPHA, WHERE "ALPHA" 
%                              IS TREATED AS INFINITY 
      for I=1:N+M
         CC(I,1) = 0.0;
         CC(I,2) = 0.0;
         if (I <= N) 
            CC(I,1) = C(I);
         else
            CC(I,2) = -1;
         end
      end
%                              BASIS(1),...,BASIS(M) HOLD NUMBERS OF
%                              BASIS VARIABLES.  INITIAL BASIS CONSISTS
%                              OF ARTIFICIAL VARIABLES ONLY
      for I=1:M
         K = N+I;
         BASIS(I) = K; 
%                              INITIALIZE Y TO Ab^(-T)*Cb = Cb
         YY(I,1) = CC(K,1);
         YY(I,2) = CC(K,2);
%                              INITIALIZE Xb TO Ab^(-1)*B = B
         XB(I) = B(I);
         if (B(I) < 0.0)
            disp ('***** ALL B(I) MUST BE NONNEGATIVE *****')
            return
         end
      end
%                              SIMPLEX METHOD CONSISTS OF TWO PHASES
      for IPHASE=1:2
         if (IPHASE == 1) 
%                              PHASE I:  ROW 2 OF D (WITH COEFFICIENTS OF
%                              ALPHA) SEARCHED FOR MOST NEGATIVE ENTRY
            MROW = 2;
            LIM = N+M;
         else
%                              PHASE II:  FIRST N ELEMENTS OF ROW 1 OF
%                              D SEARCHED FOR MOST NEGATIVE ENTRY
%                              (COEFFICIENTS OF ALPHA NONNEGATIVE NOW)
            MROW = 1;
            LIM = N;
%                              IF ANY ARTIFICIAL VARIABLES LEFT IN
%                              BASIS AT BEGINNING OF PHASE II, THERE
%                              IS NO FEASIBLE SOLUTION
            for I=1:M
               if (BASIS(I) > LIM) 
                  disp (' ***** NO FEASIBLE SOLUTION *****')
                  return
               end 
            end
         end 
%                              THRESH = SMALL NUMBER.  WE ASSUME SCALES
%                              OF A AND C ARE NOT *TOO* DIFFERENT 
         THRESH = 0.0;
         for J=1:LIM
            THRESH = max(THRESH,abs(CC(J,MROW)));
         end
         THRESH = 1000*EPS*THRESH;
%                              BEGINNING OF SIMPLEX STEP
         while (1 > 0)
%                              D^T = Y^T*A - C^T
            for IR=1:MROW
               for J=1:N+M
                  D(J,IR) = -CC(J,IR); 
               end
%                              LAST M COLUMNS OF A FORM IDENTITY MATRIX
               for J=1:M
                  D(N+J,IR) = D(N+J,IR) + YY(J,IR);
               end
%                              FIRST N COLUMNS STORED IN SPARSE A MATRIX
               for IZ=1:NZ
                  I = IROW(IZ);
                  J = JCOL(IZ);
                  D(J,IR) = D(J,IR) + A(IZ)*YY(I,IR);
               end
            end
%                              FIND MOST NEGATIVE ENTRY OF ROW MROW
%                              OF D, IDENTIFYING PIVOT COLUMN JP
            CMIN = -THRESH;
            JP = 0;
            for J=1:LIM
               if (D(J,MROW) < CMIN) 
                  CMIN = D(J,MROW);
                  JP = J;
               end
            end
%                              IF ALL ENTRIES NONNEGATIVE (ACTUALLY,
%                              IF GREATER THAN -THRESH) PHASE ENDS
            if (JP == 0)  
               break
            end
%                              COPY JP-TH COLUMN OF A ONTO VECTOR AP
            if (JP <= N) 
%                              JP-TH COLUMN IS PART OF SPARSE A MATRIX
               for I=1:M
                  AP(I) = 0;
               end
               for IZ=1:NZ
                  J = JCOL(IZ);
                  if (J == JP) 
                     I = IROW(IZ);
                     AP(I) = A(IZ);
                  end
               end
            else
%                              JP-TH COLUMN IS COLUMN OF FINAL IDENTITY
               for I=1:M
                  AP(I) = 0;
               end
               AP(JP-N) = 1;
            end
%                              SOLVE Ab*V = AP
            V = SPARSOL(AB,IROWB,JCOLB,NZB,AP,M);
%                              FIND SMALLEST POSITIVE RATIO
%                              Xb(I)/V(I), IDENTIFYING PIVOT ROW IP
            RATMIN = 0.0; 
            IP = 0;
            for I=1:M
               if (V(I) > THRESH) 
                  RATIO = XB(I)/V(I);
                  if (IP==0 || RATIO < RATMIN) 
                     RATMIN = RATIO;
                     IP = I;
                  end
               end
            end
%                              IF ALL RATIOS NONPOSITIVE, MAXIMUM
%                              IS UNBOUNDED
            if (IP == 0) 
               disp (' ***** UNBOUNDED MAXIMUM *****')
               return
            end 
%                              ADD X(JP) TO BASIS
            BASIS(IP) = JP;
%                              UPDATE Ab.  PUT NONZEROS OF JP-TH  
%                              COLUMN OF A (=AP) INTO COLUMN IP
%                              OF SPARSE MATRIX Ab
            NZBOLD = NZB;
            for I=1:M
               if (AP(I) ~= 0.0) 
                  NZB = NZB+1; 
                  IROWB(NZB) = I; 
                  JCOLB(NZB) = IP; 
                  AB(NZB) = AP(I);
               end
            end
            NZBNEW = NZB;
%                              REMOVE ELEMENTS OF OLD COLUMN IP 
            NZB = 0;
            for IZ=1:NZBNEW
               if (JCOLB(IZ) ~= IP || IZ > NZBOLD) 
                  NZB = NZB+1;
                  JCOLB(NZB) = JCOLB(IZ);
                  IROWB(NZB) = IROWB(IZ);
                  AB(NZB) = AB(IZ);
               end
            end
%                              SOLVE Ab*Xb = B
            XB = SPARSOL(AB,IROWB,JCOLB,NZB,B,M);
%                              SOLVE Ab^T*Y = Cb
            for IR=1:MROW
               for J=1:M
                  K = BASIS(J);
                  V(J) = CC(K,IR);
               end
               YY(1:M,IR) = SPARSOL(AB,JCOLB,IROWB,NZB,V,M);
            end
         end 
%                              END OF SIMPLEX STEP
      end
%                              END OF PHASE II; CALCULATE X
      for J=1:N
         X(J) = 0.0;
      end
      for I=1:M
         K = BASIS(I);
         X(K) = XB(I);
         Y(I) = YY(I,1);
      end
%                              CALCULATE P
      P = 0.0; 
      for I=1:N
         P = P + C(I)*X(I);
      end
