這裏有，裏面還有 說明 文檔。
以下是，采用大M法及字典序規則的 單純型法 算法 參考程序：
function [x,f]=Lps_Mlex(c,A,b,M,N,pre)
% [x,f]=Lps_Mlex(c,A,b,M,N,pre) 采用單純型法中的
% 大M法，並給字典序規則解下列線性規劃
% min f=ct*x s.t. Ax=b,x 所有分量 >=0
% M是壹個充分大的數，N是引進人工變量的個數，N應不超過
% （通常等於）約束等式的個數，pre 是精度
% 返回結果 x 是最優解，f 是最優解處的函數值
[m,n]=size(A);
if nargin<6,pre=0;end;if nargin<5,N=0;end
if N>m,error('N不能超過約束條件的個數m');
else,
A=[A,[zeros(N,m-N);eye(m-N)]];
c=[c(:)',zeros(1,m-N)];
A=[A,eye(m)];c=[c,M*ones(1,m)];
m1=n+m-N+1:n+2*m-N;B=A(:,m1);
x=zeros(n+2*m-N,1);x=x(:);
tb=find(b<0);b(tb)=-b(tb);A(tb,:)=-A(tb,:);
xB=B\b;x(m1)=xB;f=c*x;
criterion=c(m1)*(B\A)-c;[z1,z2]=max(criterion);
while z1>pre
az2=B\A(:,z2);
if z2<=pre,x=nan*ones(length(c));break
else,t1=find(az2>pre);p=[xB,B\eye(size(B))];pp=[];
for kk=1:length(t1);
pp(kk,:)=p(t1(kk),:)./az2(t1(kk));
end
tt1=min(xB(t1)./az2(t1));[tt0,tt2]=lex_min(pp);
t3=t1(tt2);B(:,t3)=A(:,z2);x(m1)=xB-tt1*az2;
m1(t3)=z2;x(z2)=tt1;f=c(m1)*xB;xB=x(m1);
criterion=c(m1)*(B\A)-c;[z1,z2]=max(criterion);
end,end,end;
if sum(x(n+m-N+1:n+2*-N))<=pre*m
x(m1)=xB;f=c(1:n)*x(1:n);x=x(1:n);
else,x=nan*ones(length(c));x=x(:);x=x(1:n);
end
% 例子
% f=[0;0;0;-3/4;20;-0.5;6];
% a=[0.25,-8,-1,9;0.5,-12,-0.5,3;0,0,1,0];
% a=[eye(3),a];b=[0;0;1];vlb=zeros(7,1);
% x=linprog(f,[],[],a,b,vlb); %用來做對比的
% xx=Lps_Mlex(f,a,b,100000,3);[x,xx]
function [y,k]=lex_min(x)
% [y,k]=lex_min(x) 按行求矩陣x 字典序最小行向量
% 返回值 y 是 矩陣 x 字典序最小行向量，k 是 y 在 x 中的行數
[mx,nx]=size(x);k=1;y=x(1,:);
if mx==1,k=1;y=x(1,:);
else,
[t1,t2]=min(x);t3=zeros(mx,nx);
for i=1:nx;t3(:,i)=t1(i)*ones(mx,1);end
t4=sum((t3~=x));tt=find(t4~=0);
k=t2(tt(1));y=x(k,:);
end
分別保存為 Lps_Mlex.m 和 lex_min.m 文件
這已經很簡單了