% FIGURE 4.5.1
      function [P,X,Y] = DLPRG(A,B,C,N,M)
%
%  FUNCTION DLPRG USES THE 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      - THE M BY N CONSTRAINT COEFFICIENT
%             MATRIX.
%
%    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.
%
%-----------------------------------------------------------------------
%                              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);
%                              BASIS(1),...,BASIS(M) HOLD NUMBERS OF
%                              BASIS VARIABLES.  INITIAL BASIS CONSISTS
%                              OF ARTIFICIAL VARIABLES ONLY
      for I=1:M
         BASIS(I) = N+I;
      end
%                              INITIALIZE SIMPLEX TABLEAU
      for I=1:M+2
      for J=1:N+M+1
         TAB(I,J) = 0.0;
      end
      end
%                              LOAD A INTO UPPER LEFT HAND CORNER
%                              OF TABLEAU
      for I=1:M
      for J=1:N
         TAB(I,J) = A(I,J);
      end
      end
%                              LOAD M BY M IDENTITY TO RIGHT OF A
%                              AND LOAD B INTO LAST COLUMN
      for I=1:M
         TAB(I,N+I) = 1.0;
         TAB(I,N+M+1) = B(I);
      end
%                              ROW M+1 CONTAINS -C, INITIALLY
      for J=1:N
         TAB(M+1,J) = -C(J);
      end
%                              ROW M+2 CONTAINS COEFFICIENTS OF
%                              "ALPHA", WHICH IS TREATED AS +INFINITY
      for I=1:M
         TAB(M+2,N+I) = 1.0;
      end
%                              CLEAR "ALPHAS" IN LAST ROW
      for I=1:M
      for J=1:N+M+1
         TAB(M+2,J) = TAB(M+2,J) - TAB(I,J);
      end
      end
%                              SIMPLEX METHOD CONSISTS OF TWO PHASES
      for IPHASE=1:2
         if (IPHASE == 1)
%                              PHASE I:  ROW M+2 (WITH COEFFICIENTS OF
%                              ALPHA) SEARCHED FOR MOST NEGATIVE ENTRY
            MROW = M+2;
            LIM = N+M;
         else
%                              PHASE II:  FIRST N ELEMENTS OF ROW M+1
%                              SEARCHED FOR MOST NEGATIVE ENTRY
%                              (COEFFICIENTS OF ALPHA NONNEGATIVE NOW)
            MROW = M+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) > N)
                  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(TAB(MROW,J)));
         end
         THRESH = 1000*EPS*THRESH;
%                              BEGINNING OF SIMPLEX STEP
         while (1 > 0)
%                              FIND MOST NEGATIVE ENTRY IN ROW MROW,
%                              IDENTIFYING PIVOT COLUMN JP
            CMIN = -THRESH;
            JP = 0;
            for J=1:LIM
               if (TAB(MROW,J) < CMIN)
                  CMIN = TAB(MROW,J);
                  JP = J;
               end
            end
%                              IF ALL ENTRIES NONNEGATIVE (ACTUALLY,
%                              IF GREATER THAN -THRESH) PHASE ENDS
            if (JP == 0)
               break
            end
%                              FIND SMALLEST POSITIVE RATIO
%                              B(*)/TAB(*,JP), IDENTIFYING PIVOT
%                              ROW IP
            RATMIN = BIGNO;
            IP = 0;
            for I=1:M
               if (TAB(I,JP) > THRESH)
                  RATIO = TAB(I,N+M+1)/TAB(I,JP);
                  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;
%                              NORMALIZE PIVOT ROW TO MAKE TAB(IP,JP)=1
            AMULT = 1.0/TAB(IP,JP);
            for J=1:N+M+1
               TAB(IP,J) = AMULT*TAB(IP,J);
            end
%                              ADD MULTIPLES OF PIVOT ROW TO OTHER
%                              ROWS, TO KNOCK OUT OTHER ELEMENTS IN
%                              PIVOT COLUMN
            for I=1:MROW
               if (I == IP)
                  continue
               end
               AMULT = TAB(I,JP);
               for J=1:N+M+1
                  TAB(I,J) = TAB(I,J) - AMULT*TAB(IP,J);
               end
            end
         end
%                              END OF SIMPLEX STEP
      end
%                              END OF PHASE II; READ X,P,Y FROM
%                              FINAL TABLEAU
      for J=1:N
         X(J) = 0.0;
      end
      for I=1:M
         K = BASIS(I);
         X(K) = TAB(I,N+M+1);
      end
      P = TAB(M+1,N+M+1);
      for I=1:M
         Y(I) = TAB(M+1,N+I);
      end