线形规划的单纯形法子程序，给你一堆，自己慢慢看。来自于《Delphi常用数值算法》，网上有
该书光盘的下载，不记得网址了。
Procedure SIMPLX(var A:matrx2;
M, N, MP, NP, M1, M2, M3:integer;
var ICASE:integer;var IZROV, IPOSV:array of integer);
Label 1,2,3,4;
const
EPS =0.000001;
var
L1,L2,L3:array[0..100] of integer;
K,I,NL1,NL2,IR,KP,M12,IP,IQ,KH:integer;
Q1,BMAX:real;
begin
If M <> M1 + M2 + M3 then
begin
ShowMessage(' Bad input constraint counts');
Exit;
end;
NL1:=N;
For K:=1 To N do
begin
L1[K]:=K;
IZROV[K]:=K;
end;
NL2:=M;
For I:=1 To M do
begin
If A[I + 1, 1] < 0 then
begin
ShowMessage('Bad input tableau.');
Exit;
end;
L2:=I;
IPOSV:=N + I;
end;
For I:=1 To M2 do
L3:=1;
IR:=0;
If M2 + M3 = 0 then
GoTo 3;
IR:=1;
For K:=1 To N + 1 do
begin
Q1:=0;
For I:=M1 + 1 To M do
Q1:=Q1 + A[I + 1, K];
A[M + 2, K]:=-Q1;
end;
repeat
SIMP1(A, MP, NP, M + 1, L1, NL1, 0, KP, BMAX);
If (BMAX <= EPS) And (A[M + 2, 1] < -EPS) then
begin
ICASE:=-1;
Exit;
end
else
If (BMAX <= EPS) And (A[M + 2, 1] <= EPS) then
begin
M12:=M1 + M2 + 1;
If M12 <= M then
begin
For IP:=M12 To M do
begin
If IPOSV[IP] = IP + N then
begin
SIMP1(A, MP, NP, IP, L1, NL1, 1, KP, BMAX);
If BMAX > 0 then
GoTo 1;
end;
end;
end;
IR:=0;
M12:=M12 - 1;
If M1 + 1 > M12 then
Exit;
For I:=M1 + 1 To M12 do
begin
If L3[I - M1] = 1 then
begin
For K:=1 To N + 1 do
A[I + 1, K]:=-A[I + 1, K];
end;
end;
goto 3;
end;
SIMP2(A, M, N, MP, NP, L2, NL2, IP, KP, Q1);
If IP = 0 then
begin
ICASE:=-1;
Exit;
end;
1: SIMP3(A, MP, NP, M + 1, N, IP, KP);
If IPOSV[IP] >= N + M1 + M2 + 1 then
begin
For K:=1 To NL1 do
If L1[K] = KP then
goto 4;
4: NL1:=NL1 - 1;
For IQ:=K To NL1 do
L1[IQ]:=L1[IQ + 1];
end
else
begin
If IPOSV[IP] < N + M1 + 1 then
GoTo 2;
KH:=IPOSV[IP] - M1 - N;
If L3[KH] = 0 then
GoTo 2;
L3[KH]:=0;
end;
A[M + 2, KP + 1]:=A[M + 2, KP + 1] + 1;
For I:=1 To M + 2 do
A[I, KP + 1]:=-A[I, KP + 1];
2: IQ:=IZROV[KP];
IZROV[KP]:=IPOSV[IP];
IPOSV[IP]:=IQ;
until IR <> 0;
3: SIMP1(A, MP, NP, 0, L1, NL1, 0, KP, BMAX);
If BMAX <= 0 then
begin
ICASE:=0;
Exit;
end;
SIMP2(A, M, N, MP, NP, L2, NL2, IP, KP, Q1);
If IP = 0 then
begin
ICASE:=1;
Exit;
end;
SIMP3(A, MP, NP, M, N, IP, KP);
GoTo 2;
end;
Procedure SIMP1(var A:matrx2;
MP, NP, MM:integer;
LL:array of integer;
NLL, IABF:integer;var KP:integer;
var BMAX:real);
var
K:integer;
TEST:real;
begin
KP:=LL[1];
BMAX:=A[MM + 1, KP + 1];
For K:=2 To NLL do
begin
If IABF = 0 then
TEST:=A[MM + 1, LL[K] + 1] - BMAX
else
TEST:=Abs(A[MM + 1, LL[K] + 1]) - Abs(BMAX);
If TEST > 0 then
begin
BMAX:=A[MM + 1, LL[K] + 1];
KP:=LL[K];
end;
end;
end;
Procedure SIMP2(var A:matrx2;
M, N, MP, NP:integer;
L2:array of integer;
NL2:integer;var IP, KP:integer;
Q1:real);
label 1,2;
const
EPS = 0.000001;
var
I,J,II,K:integer;
FLAG,Q,QP,Q0:real;
begin
IP:=0;
FLAG:=0;
For I:=1 To NL2 do
begin
If A[L2 + 1, KP + 1] < -EPS then
FLAG:=1;
If FLAG = 1 then
goto 1;
end;
1: If FLAG = 0 then
Exit;
Q1:=-A[L2 + 1, 1] / A[L2 + 1, KP + 1];
IP:=L2;
For I:=I + 1 To NL2 do
begin
II:=L2;
If A[II + 1, KP + 1] < -EPS then
begin
Q:=-A[II + 1, 1] / A[II + 1, KP + 1];
If Q < Q1 then
begin
IP:=II;
Q1:=Q;
end
else
If Q = Q1 then
begin
For K:=1 To N do
begin
QP:=-A[IP + 1, K + 1] / A[IP + 1, KP + 1];
Q0:=-A[II + 1, K + 1] / A[II + 1, KP + 1];
If Q0 <> QP then
goto 2;
end;
2: If Q0 < QP then
IP:=II;
end;
end;
end;
end;
Procedure SIMP3(var A:matrx2;
MP, NP, I1, K1:integer;var IP, KP:integer);
var
PIV:real;
II,KK:integer;
begin
PIV:=1 / A[IP + 1, KP + 1];
For II:=1 To I1 + 1 do
begin
If II - 1 <> IP then
begin
A[II, KP + 1]:=A[II, KP + 1] * PIV;
For KK:=1 To K1 + 1 do
begin
If KK - 1 <> KP then
A[II, KK]:=A[II, KK] - A[IP + 1, KK] * A[II, KP + 1];
end;
end;
end;
For KK:=1 To K1 + 1 do
If KK - 1 <> KP then
A[IP + 1, KK]:=-A[IP + 1, KK] * PIV;
A[IP + 1, KP + 1]:=PIV;
end;
