function [xopt,zmin,is_bounded,has_soln,basis] = lp(c,A,b)
% Usage: [xopt,zmin,is_bounded,has_soln,basis] = lp(c,A,b)
% Description: This program finds a solution of the
% standard linear programming problem:
%   minimize    z = c'x
%   subject to  Ax = b, x >= 0
% using the two phase method, where the simplex
% method is used at each stage.  Output:
%   xopt: an optimal solution.
%   zmin: the optimal value.
%   is_bounded: 1 if the solution is bounded; 0 if unbounded.
%   has_soln: 1 if the problem is solvable; 0 if unsolvable.
% basis: indices of the basis of the solution.

% Contact: Tom Shores, tshores@math.unl.edu

[m,n] = size(A);
% Phase one
A = ((1-2*(b<0))*ones(1,n)).*A;
b = (1-2*(b<0)).*b;
if (m == 1)
  d = -A;
else
  d = -sum(A);
end
w0 = sum(b);
H = [A b;c' 0;d -w0]; % The initial simplex table of phase one
epsilon = norm(H,inf)*10*eps; % relatively near zero
index = 1:n;
basis = n+1:n+m;
[n1,n2] = size(H);
has_soln = 0;
while has_soln == 0
  [fm,jp] = min(H(n1,1:n2-1));
  if fm >= 0
    is_bounded = 1;    % bounded solution
    has_soln = 1;
  else
    [hm,ip] = max(H(1:n1-2,jp));
    if hm <= 0
      is_bounded = 0;  % unbounded solution
      has_soln = 1;
    else
      for i = 1:n1-2
        if H(i,jp) > 0
          h1(i) = H(i,n2)/H(i,jp);
        else
          h1(i) = Inf;
        end
      end
      [minh1,ip] = min(h1);
      basis(ip) = index(jp);
      piv = H(ip,jp);
      if (abs(piv) <= epsilon)
        has_soln = 0; xopt = []; zmin = []; is_bounded = [];
        error('lp: this problem is singular.');
      else
        H(ip,:) = H(ip,:)/piv;
        for i = 1:n1
          if i ~= ip
            H(i,:) = H(i,:) - H(i,jp)*H(ip,:);
          end
        end
      end
    end
  end
end
if (H(m+2,n+1) < -epsilon)
  has_soln = 0; xopt = []; zmin = []; is_bounded = [];
  error('lp: this problem is unsolvable.');
else
  has_soln = 1;
  j = 0;
  for i = 1:n
    j = j+1;
    if H(m+2,j) > 1e-10
      H = [H(:,1:j-1) H(:,j+1:size(H,2))];
      index = [index(1:j-1) index(j+1:length(index))];
      j = j-1;
    end
  end
  H(m+2,:) = [];
  % Phase 2
  if length(index) > 0
    [n1,n2] = size(H);
    has_soln = 0;
    while has_soln == 0
      [fm,jp] = min(H(n1,1:n2-1));
      if (fm >= 0)
        is_bounded = 1;    % bounded solution
        has_soln = 1;
      else
        [hm,ip] = max(H(1:n1-1,jp));
        if hm <= 0
          is_bounded = 0;  % unbounded solution
          has_soln = 1;
        else
          for i = 1:n1-1
            if H(i,jp) > 0
              h1(i) = H(i,n2)/H(i,jp);
            else 
              h1(i) = Inf;
            end
          end
          [minh1,ip] = min(h1);
          basis(ip) = index(jp);
          piv = H(ip,jp);
          if piv == 0
            has_soln = 0; xopt = []; zmin = []; is_bounded = [];
            error('lp: this problem is singular.');
          else
            H(ip,:) = H(ip,:)/piv;
            for i = 1:n1
              if i ~= ip
                H(i,:) = H(i,:) - H(i,jp)*H(ip,:);
              end
            end
          end
        end
      end
    end
    if (is_bounded == 1)
      xopt = zeros(n+m,1);
      [n1,n2] = size(H);
      for i = 1:m
        xopt(basis(i)) = H(i,n2);
      end
      xopt = [xopt(1:n,1)];
      zmin = -H(n1,n2);
    else
      xopt = [];
      zmin = -Inf;
    end
  else
    xopt = zeros(n+m,1);
    zmin = 0;
  end
end