function [xopt,zmin,is_bounded,has_soln,basis] = linprog(c,A,b,D,d,E,e)
% usage: [xopt,zmin,is_bounded,sol,basis] = linprog(c,A,b,D,d,E,e)
% description: This program finds a solution to the
% standard linear programming problem:
%   minimize    z = c'*x
%   subject to  A*x <= b, D*x >= d, E*x = e, x >= 0
% using the two phase method, where the simplex
% method is used at each stage. Here x is the vector of
% decision variables. The return variables are:
% xopt: an optimal solution of decision variables
% zmin: the optimal value, 
% is_bounded: 1 if the solution is bounded; 0 otherwise,
% has_soln: 1 if the problem is solvable; 0 otherwise,
% basis:  values of slack/surplus variables.
% Arguments c,A,b are mandatory, pairs L,d and E,e optional.

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

% local variables:
% j,k: index variables

% a little bit of bulletproofing
if (nargin < 3)
  error('linprog: need at least three arguments');
end
[mm,nn] = size(A);
if (nn ~= 0)
  if (length(c) ~= nn | size(c,2) ~= 1)
    error('linprog: matrix A must have same column number as c');
  end
  if (length(b) ~= mm | size(b,2) > 1)
    error('linprog: vector b must have same row number as A');
  end
else
  nn = length(c);
end
if (exist('D','var'))
  if  (~exist('d','var'))
    error('linprog: need matching arguments D,d')
  else
    if (size(D,2) ~= nn & size(D,2) ~= 0)
      error('linprog: matrix D must have same column number as c');
    end
    if  (length(d) ~= size(D,1) | size(d,2) > 1)
      error('linprog: vector d must have same row number as D')
    end
  end
end
if (exist('E','var'))
  if  (~exist('e','var'))
    error('lingprog: need matching arguments E,e')
  else
    if (size(E,2) ~= nn)
      error('linprog: E must have same column number as c');
    end
    if  (length(e) ~= size(E,1) | size(e,2) > 1)
      error('linprog: column vector e must have length equal to row number of E')
    end
 end
end
% Form augmented system
if (exist('D','var'))
  A = [A;-D];
  b = [b;-d];
end
if (exist('E','var'))
  A = [A;E;-E];
  b = [b;e;-e];
end
A = [A,eye(size(A,1))];
c = [c;zeros(size(A,1),1)];
[xopt,zmin,is_bounded,has_soln,basis] = lp(c,A,b);
if (length(xopt)>0)
  basis = xopt(nn+1:length(xopt));
  xopt = xopt(1:nn);
else
  basis = [];
  xopt = [];
end