clear;clc;close all
% 两阶段法
disp('在运行程序之前，请确认已经将LP化为标准型');
tic;
A=[1 1 1 1 0 0 0;-2 1 -1 0 -1 1 0;0 3 1 0 0 0 1];
c0=[-3 0 1 0 0 -1 -1];
b=[4;1;9];
% A=[2 2 1 0 0;4 0 0 1 0;0 5 0 0 1];
% c=[2 3 0 0 0];
% b=[12;16;15];
% A=input('请输入系数矩阵A=');
% c=input('请输入目标函数的系数c=');
% b=input('请输入右端向量b=');
[m,n]=size(A);
% C_B B b x
bigm=zeros(m,m+n);
% 第一阶段
c=[0 0 0 0 0 -1 -1];
% 加入松弛变量以后，直接就会形成单位阵
% 第一步：形成初始的单纯形表
% bigm(:,2)=n-m+1:n;
total_m=nchoosek(1:n,m);
for i=1:size(total_m,1)
    temp_m=A(:,total_m(i,:));
    if all(temp_m-eye(m)==0)
        bigm(:,2)=total_m(i,:);
        break
    end
end
bigm(:,1)=c(bigm(:,2));
bigm(:,3)=b;
bigm(:,4:end)=A;
sigma=c-bigm(:,1)'*A;
% 第二步：进行最优性检验
while any(sigma>0)
    % 第三步：从一个基可行解转换到另一个目标函数值更大的基可行解，列出新的单纯形表
    % （1）确定换入基变量
    [dmax,zyl]=max(sigma);
    % （2）确定换出基变量
    tempv=bigm(:,3)./bigm(:,zyl+3);
    zyh0=find(tempv==min(tempv(tempv>0)));
    zyh=zyh0(end);
    % （3）用换入基变量代替换出基变量
    bigm(zyh,2)=zyl;
    bigm(zyh,1)=c(zyl);
    bigm(zyh,3:end)=bigm(zyh,3:end)/bigm(zyh,zyl+3);
    otherh=setdiff(find(bigm(:,zyl+3)),zyh);
    bigm(otherh,3:end)=bigm(otherh,3:end)-bsxfun(@times,bigm(zyh,3:end),bigm(otherh,zyl+3));
    sigma=c-bigm(:,1)'*bigm(:,4:end);
end
% 第二阶段
c=[-3 0 1 0 0];
s_temp=setdiff(1:n,1:length(c));
bigm(:,1)=c(bigm(:,2));
bigm(:,s_temp+3)=[];
sigma=c-bigm(:,1)'*bigm(:,4:end);
% 第二步：进行最优性检验
while any(sigma>0)
    % 第三步：从一个基可行解转换到另一个目标函数值更大的基可行解，列出新的单纯形表
    % （1）确定换入基变量
    [dmax,zyl]=max(sigma);
    % （2）确定换出基变量
    tempv=bigm(:,3)./bigm(:,zyl+3);
    zyh0=find(tempv==min(tempv(tempv>0)));
    zyh=zyh0(end);
    % （3）用换入基变量代替换出基变量
    bigm(zyh,2)=zyl;
    bigm(zyh,1)=c(zyl);
    bigm(zyh,3:end)=bigm(zyh,3:end)/bigm(zyh,zyl+3);
    otherh=setdiff(find(bigm(:,zyl+3)),zyh);
    bigm(otherh,3:end)=bigm(otherh,3:end)-bsxfun(@times,bigm(zyh,3:end),bigm(otherh,zyl+3));
    sigma=c-bigm(:,1)'*bigm(:,4:end);
end
n=length(c);
zyj=zeros(2,n);
zyj(1,:)=1:n;
zyj(2,bigm(:,2))=bigm(:,3);
zyz=c*zyj(2,:)';
toc;
fprintf('x%d=%f\n',zyj);
fprintf('Z=%f\n',zyz);