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

simplex.m

%% Simplex Method
% /DAA/simplex.m
% Solve LP problem using Simplex Method and Bland's Law.
% chaonan99(R)
% 2016/03/20

% [Note] The input should already have a initial basic feasible solution.
% [Thanks to] http://wenku.baidu.com/link?url=bUal8_VtNYHdIXF97JKo_r6vzYt3lGlSFO2W2Wz6wInLnw-3STZpbyxSoaOIu2XKE3z5O4D5xjErH5_tqXY-j3-DubAapkfCuyraVMxMNNO

%% code
function [x,z] = 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(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+1
        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+1)';
    c = S(m+1,1:n);
    r = find(c>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