% FIGURE 4.6.1
      function [P,X,Y] = DLPRV(DOTA,B,C,N,M)
%
%  FUNCTION DLPRV 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)
%
%    THE LAST M COLUMNS OF A MUST CONTAIN AN IDENTITY MATRIX, AND
%    B(1),...,B(M) MUST BE NONNEGATIVE.
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    DOTA   - STRING CONTAINING THE NAME 
%             OF A USER-SUPPLIED FUNCTION:
%             DOTA(Z,J) SHOULD RETURN THE
%             DOT PRODUCT OF THE M-VECTOR Z
%             WITH COLUMN J OF A, FOR J <= N-M
%             (IDENTITY MATRIX ASSUMED IN LAST 
%             M COLUMNS). 
%
%    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 (N>M).
%
%    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.
%
%-----------------------------------------------------------------------
%                              EPS = MACHINE FLOATING POINT RELATIVE
%                                    PRECISION
%                              BIGNO = A VERY LARGE NUMBER
% *****************************
      EPS = eps;
      BIGNO = 1.e35;
% *****************************
      P = 0;
      X = zeros(N,1);
      Y = zeros(M,1);
%                              INITIALIZE Ab^(-1) TO IDENTITY
      for I=1:M
      for J=1:M
         ABINV(I,J) = 0.0;
         if (I == J)
            ABINV(I,J) = 1.0;
         end
      end
      end
%                              BASIS(1),...,BASIS(M) HOLD NUMBERS OF
%                              BASIS VARIABLES
      for I=1:M
%                              INITIAL BASIS = LAST M VARIABLES
         K = N-M+I;
         BASIS(I) = K;
%                              INITIALIZE Y TO Ab^(-T)*Cb = Cb
         Y(I) = C(K);
%                              INITIALIZE Xb TO Ab^(-1)*B = B
         XB(I) = B(I);
      end
%                              THRESH = SMALL NUMBER.  WE ASSUME SCALES
%                              OF A AND C ARE NOT *TOO* DIFFERENT
      THRESH = 0.0;
      for J=1:N
         THRESH = max(THRESH,abs(C(J)));
      end
      THRESH = 1000*EPS*THRESH;
%                              BEGINNING OF SIMPLEX STEP
      while (1 > 0)
%                              D^T = Y^T*A - C^T
         for J=1:N
            if (J <= N-M) 
               D(J) = feval(DOTA,Y,J) - C(J);
            else
               D(J) = Y(J-(N-M)) - C(J);
            end
         end
%                              FIND MOST NEGATIVE ENTRY IN D,
%                              IDENTIFYING PIVOT COLUMN JP
         CMIN = -THRESH;
         JP = 0;
         for J=1:N
            if (D(J) < CMIN)
               CMIN = D(J);
               JP = J;
            end
         end
%                              IF ALL ENTRIES NONNEGATIVE (ACTUALLY,
%                              IF GREATER THAN -THRESH) WE ARE THROUGH
         if (JP == 0)
            break
         end
%                              V = Ab^(-1)*Ajp (Ajp = COLUMN JP OF A)
         for I=1:M
            for J=1:M
               WK(J) = ABINV(I,J);
            end
            if (JP <= N-M)
               V(I) = feval(DOTA,WK,JP);
            else
               V(I) = WK(JP-(N-M));
            end
         end
%                              FIND SMALLEST POSITIVE RATIO
%                              Xb(I)/V(I), IDENTIFYING PIVOT ROW IP
         RATMIN = BIGNO;
         IP = 0;
         for I=1:M
            if (V(I) > THRESH)
               RATIO = XB(I)/V(I);
               if (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^(-1) = E^(-1)*Ab^(-1)
%                                     Xb = E^(-1)*Xb
         for J=1:M
            ABINV(IP,J) = ABINV(IP,J)/V(IP);
         end
         XB(IP) = XB(IP)/V(IP);
         for I=1:M
            if (I == IP)
               continue
            end
            for J=1:M
               ABINV(I,J) = ABINV(I,J) - V(I)*ABINV(IP,J);
            end
            XB(I) = XB(I) - V(I)*XB(IP);
         end
%                              CALCULATE Y = Ab^(-T)*Cb
         for I=1:M
            Y(I) = 0.0;
            for J=1:M
               K = BASIS(J);
               Y(I) = Y(I) + ABINV(J,I)*C(K);
            end
         end
      end
%                              END OF SIMPLEX STEP
%                              CALCULATE X
      for J=1:N
         X(J) = 0.0;
      end
      for I=1:M
         K = BASIS(I);
         X(K) = XB(I);
      end
%                              CALCULATE P
      P = 0.0;
      for I=1:N
         P = P + C(I)*X(I);
      end