Solve LP problem using Two-Stage Method and Bland's Law.

twostage.m

%% Two-Stage Method
% /DAA/twostage.m
% [Function] Solve LP problem using Two-Stage Method and Bland's Law.
% [Author] chaonan99(R)
% [Last modify] 2016/03/28
% [Contect me] [email protected]

% [Note] 
% 1. There are no examination on input size, which the user should ensure
% its correctness.
% 2. Vector b should not have negative value. You can multiply -1 to the
% line which has a negative b value.
% 3. The program calculate the MAX of z.

%% code
function [x,z] = twostage(A,b,c)
% A -- constraint matrix
% b -- constraint vector
% c -- checknumber
% x -- solution vector
% z -- optimal solution

[m,n] = size(A);
S = [A eye(m) b'; zeros(1,n) -ones(1,m) 0; c zeros(1,m+1)];
S(m+1,:) = S(m+1,:) + sum(S(1:m,:));

%% Stage 1
g = S(m+1,:);
r = find(g>0);
len = length(r);
while len
    rb = zeros(1,m);
    for p = 1:m
        if A(p,r(1))<=0 % In-base variable have the minimum index (Bland's law)
            rb(p) = Inf; % In order to use min to find out-base variable
        else
            rb(p) = b(p)/A(p,r(1));
        end
    end
    [~, row] = min(rb);
    S(row,:) = S(row,:)./S(row,r(1));
    for i = 1:m+2
        if i ~= row
            S(i,:) = S(i,:) - S(row,:)*S(i,r(1));
        end
    end
    A = S(1:m,1:n);
    b = S(1:m,n+m+1)';
    g = S(m+1,1:n);
    r = find(g>0);
    % break point
    len = length(r);
end

if S(m+1,m+n+1)>0
    error('No feasible solution!');
end

%% Stage 2
c = S(m+2,1:n);
r = find(c>0);
len = length(r);
while len
    for k = 1:len
        d = find(A(:,r(k))>0, 1);
        if (isempty(d))
            error('No optimal solution!');
        end
    end
    rb = zeros(1,m);
    for p = 1:m
        if A(p,r(1))<=0 % In-base variable have the minimum index (Bland's law)
            rb(p) = Inf; % In order to use min to find out-base variable
        else
            rb(p) = b(p)/A(p,r(1));
        end
    end
    [~, row] = min(rb); % min method find the minimum index when there are more than one min value
    S(row,:) = S(row,:)./S(row,r(1));
    for i = 1:m+2
        if i ~= row
            S(i,:) = S(i,:) - S(row,:)*S(i,r(1));
        end
    end
    A = S(1:m,1:n);
    b = S(1:m,n+m+1)';
    c = S(m+2,1:n);
    r = find(c>0);
    % break point
    len = length(r);
end

%% Outputing optimal solution and solution vector.
x = zeros(1,n);
for i = 1:n
    j = find(A(:,i)==1);
    k = find(A(:,i)==0);
    if length(j)==1 && length(k)==m-1
        x(i) = b(j);
    end
end
z = -S(m+2, n+m+1);
end

twostage_test.m

%% A test to twostage.m
% LP problem:
% min = 2*x1 + 3*x2 +5*x3 + 6*x4
% s.t. x1 + 2*x2 + 3*x3 + x4 >= 2
%    -2*x1 + x2 - x3 + 3*x4 <= -3
%        x_j >= 0 j = 1,2,3,4

% You should transform the problem into standard form (max value):
A = [1 2 3 1 -1 0
    2 -1 1 -3 0 -1];
b = [2 3];
c = [-2 -3 -5 -6 0 0];

[x,z] = twostage(A,b,c);

% Answer: 3.8 (opposite number of the result)