Documentation

%Attempts to solve the linear program:
%
%   maximize y = c*x
%   subject to: A*x <= b,
%   where x >= 0 and b >= 0
%
%   using the dictionary form of the simplex algorithm.
%
%   c must be a double, real row vector.
%   b must be a double, real column vector.
%   A must be a double, real, non-empty matrix with
%     size(A) = [length(b),length(c)].
%
%   Compare to
%   x = linprog(-c(:),A,b,[],[],zeros(length(c),1))
%   y = c*x
%   
%   Copyright 2011-2018 The MathWorks, Inc.

function [x,y] = lpsolve(c,A,b)
%#codegen
narginchk(3,3);
m = int32(size(A,1));
n = int32(size(A,2));
if ~validInputs(c,A,b)
    % Something is wrong with the inputs. Return NaNs.
    x = nan(n,1);
    y = nan;
    return
end
cn = c;
nIdx = 1:n; % indices of non-basic variables
bIdx = n + (1:m); % indices of basic variables
% Perform iterations.
j = findPivotColumn(cn,nIdx);
unbounded = false;
while j > 0 && ~unbounded
    [i,r] = findPivotRow(j,A,b);
    if i == 0 || isinf(r)
        unbounded = true;
    else
        b(i) = r;
        % Take account of the variables entering and leaving the basis.
        entering = nIdx(j);
        leaving = bIdx(i);
        bIdx(i) = entering;
        nIdx(j) = leaving;
        [cn,A,b] = pivot(cn,A,b,i,j);
        j = findPivotColumn(cn,nIdx);
    end
end
% Construct the solution.
if unbounded
    x = nan(n,1);
    y = inf;
else
    x = zeros(n,1);
    y = 0;
    for k = 1:m
        j = bIdx(k);
        if j <= n
            x(j) = b(k);
            y = y + c(j)*x(j);
        end
    end
end

%-------------------------------------------------------------------------------

function p = validInputs(c,A,b)
coder.inline('never');
% All inputs must be double.
p = isa(c,'double') && isa(A,'double') && isa(b,'double');
% All inputs must be real.
p = p && isreal(c) && isreal(A) && isreal(b);
% A must be a non-empty matrix.
p = p && ismatrix(A) && ~isempty(A);
% c must be a row-vector with length = size(A,2).
p = p && isrow(c) && size(c,2) == size(A,2);
% b must be a column-vector with length = size(A,1).
p = p && iscolumn(b) && size(b,1) == size(A,1);
% The elements of b must be non-negative.
p = p && all(b(:) >= 0);

%-------------------------------------------------------------------------------

function [c,A,b] = pivot(c,A,b,i,j)
% Pivot on element i,j.
coder.inline('never');
m = int32(size(A,1));
aij = A(i,j);
A(i,j) = 1;
A(i,:) = A(i,:)/aij;
for k = 1:m
    if k ~= i
        akj = A(k,j);
        A(k,j) = 0;
        A(k,:) = A(k,:) - akj*A(i,:);
        b(k) = b(k) - akj*b(i);
    end
end
cj = c(j);
c(j) = 0;
c = c - cj*A(i,:);

%-------------------------------------------------------------------------------

function j = findPivotColumn(c,idx)
% Find a suitable variable to enter the basis using the smallest
% subscript rule to prevent cycling.
coder.inline('never');
j = int32(0);
nc = int32(numel(c));
lastidx = int32(0);
for k = 1:nc
    if c(k) > 0 && (lastidx < 1 || idx(k) < lastidx)
        j = k;
        lastidx = idx(k);
    end
end

%-------------------------------------------------------------------------------

function [i,r] = findPivotRow(j,A,b)
% Find the pivot row in column j of A and also return r = b(i)/A(i,j).
coder.inline('never');
r = inf;
i = int32(0);
nb = int32(numel(b));
for k = 1:nb
    if A(k,j) > 0
        rtmp = b(k)/A(k,j);
        if rtmp < r
            r = rtmp;
            i = k;
        end
    end
end

%-------------------------------------------------------------------------------

% Test function for lpsolve
%
%   Copyright 2011-2018 The MathWorks, Inc.

function lpsolve_test()
c = [2 3 1 1];
A = [2 3 1 -1;1 0 2 1;0 2 1 1];
b = [27;9;18];
[x,y] = lpsolve(c,A,b);
x % display x
y % display y