function [x,T,B,flg] = LPsolve(A,b,c)
%
% This is a combined Phase I/Phase II calculation for solving LP.
% We assume non-degeneracy in the entire computation.
% Written by Ming Gu for Math 170
% March 2007
%
% 
% Prepare for initial Phase I computation.
%
flg       = 0;
[m,n]     = size(A);
mask      = find(b<0);
A(mask,:) = - A(mask,:);
b(mask)   = - b(mask);
B         = (n+1:n+m)';
x         = [zeros(n,1);b];
c1        = [zeros(n,1);ones(m,1)];
y         = ones(m,1);
T         = [b eye(m); c1'*x,y'];
A1        = [A eye(m)];
%
% Starting Simplex Method
%
f = ['Starting Initial Phase I Simplex Iteration... '];
disp(f);
[ITER1,B,T,flg1] = ILPRevisedSimplex(A1,B,c1,T);
f = ['Phase I took ',num2str(ITER1),' Iterations. '];
disp(f);
if (flg1 ~= 0)
    f = ['Initial Phase I Failed, flg =  ',num2str(flg1)];
    disp(f);
    flg = flg1;
    return;
end
if (max(B)>n)
    f = ['Initial Phase I Detected Primal Infeasibility, max(B) = ',num2str(max(B))];
    disp(f);
    flg = 3;
    return;
end
x   = zeros(n,1);
x(B)= T(1:m,1);
y   = T(1:m,2:m+1)'*c(B);
T(end,:) = [c'*x,y'];
f = ['Starting Phase II Iteration... '];
disp(f);
A1      = A;
b1      = b;
c1      = c;
xtol    = 1e-12;
bsave   = b;
x       = zeros(size(A1,2),1);
x(B)    = A1(:,B)\bsave;
T       =[x(B),inv(A1(:,B));x'*c1,c1(B)'/A1(:,B)];
[ITER,B,T,flg2] = ILPRevisedSimplex(A1,B,c1,T);
f = [' Phase II took ',num2str(ITER),' Iterations. '];
disp(f);
if (flg2 ~= 0)
    f = ['Phase II Failed, flg =  ',num2str(flg2)];
    disp(f);
    flg = flg2;
    return;
end
y = T(end,2:end)';
x = zeros(size(A1,2),1);
x(B) = T(1:end-1,1);