Solve LP problem using Dual Simplex Method and Bland's Law.

dual_simplex.m

%% Dual Simplex Method
% /DAA/dual_simplex.m
% [Description] Solve LP problem using Dual Simplex Method and Bland's Law.
% [Author] chaonan99(R)
% [Last Modify] 2016/04/03
% [Contact Me] [email protected]

% [Note] The input should already have an initial basic solution for the
% original problem and a feasible solution for its dual problem.

%% code
function [x,z] = dual_simplex(A,b,c)
% A -- constraint matrix
% b -- constraint vector
% c -- checknumber
% x -- solution vector
% z -- optimal solution
format rat;
S = [A b'; c 0];
[m,n] = size(A);
% Check the validity of inputs
assert(n==length(c), 'Scale of A and c do not match!');
assert(m==length(b), 'Scale of A and b do not match!');
r = find(b<0);
len = length(r);
while len
    for k = 1:len
        d = find(A(r(k),:)<0, 1);
        % A(r,:)>=0, no feasible solution.
        if (isempty(d))
            error('No optimal solution!');
        end
    end
    rb = zeros(1,n);
    for p = 1:n
        if A(r(1),p)>=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) = c(p)/A(r(1),p);
        end
    end
    [~, col] = min(rb); % min method find the minimum index when there are more than one min value
    S(r(1),:) = S(r(1),:)./S(r(1),col);
    for i = 1:m+1
        if i ~= r(1)
            S(i,:) = S(i,:) - S(r(1),:)*S(i,col);
        end
    end
    A = S(1:m,1:n);
    b = S(1:m,n+1)';
    c = S(m+1,1:n);
    r = find(b<0);
    % break point
    len = length(r);
end
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+1, n+1);
end