ineqelim2.m

来自「线性M3网络及其示意图」· M 代码 · 共 54 行

function [Aout,bout]=IneqElim(A,b)
% A*x>b

n=size(A,1);

for i=1:2*n,   % 2n times?
    r=A(1,:);
    r=r';
    p=A(2:end,:);
    p=p';
    s=b(2:end);
    
    q=p\r;%  X = A\B is the solution to the equation AX = B
    
    if rank(p)~=rank([p,r])
        A(end+1,:)=A(1,:);
        b(end+1)=b(1);
        A(1,:)=[];
        b(1)=[];
        continue;
    end
        
    
    while (~isempty(find(q < 0)) & ~isempty(p))
        ql=find(q<0);
        p(:,ql)=[];
        s(ql)=[];
        
        if rank(p)~=rank([p,r])
            q=zeros(size(r,1));
            break;
        end        
    
        q=p\r;
    end
   
    if ~isempty(p)  &  b(1) <= q'*s  &  ~isequal(q,zeros(size(r,1))),
        A(1,:)=[];
        b(1)=[];
    else
        %m=max(max(A));
        %A(1,:)=A(1,:)*m;
        %b(1)=b(1)*m;
        A(end+1,:)=A(1,:);
        b(end+1)=b(1);
        A(1,:)=[];
        b(1)=[];
    end
end



Aout=A;  %gcd can't be used to double.
bout=b;