单纯形法模板

54 阅读 0 评论 36 点赞

我是靠谱客的博主舒心鞋子，最近开发中收集的这篇文章主要介绍单纯形法模板，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

 1 #include <bits/stdc++.h>
 2 #define eps 1e-8
 3 #define INF 1e100
 4 using namespace std;
 5 long double a[105][105],x[205],ans,c[105];
 6 int n,m,xx,yy,t,B[105],N[105],d[105],tt;
 7 void pivot(int l,int e)
 8 {
 9 
swap(B[l],N[e]);
10
for(int i=0;i<=n;i++)if(i!=e)a[l][i]/=a[l][e];
11
a[l][e]=(long double)1/a[l][e];
12
for(int i=0;i<=m;i++)if(i!=l&&fabs(a[i][e])>eps)
13 
{
14
for(int j=0;j<=n;j++)if(j!=e)a[i][j]-=a[i][e]*a[l][j];
15
a[i][e]=-a[i][e]*a[l][e];
16 
}
17 }
18 int simplex()
19 {
20
int l,e;
21
while(1)
22 
{
23
e=n+1; tt=0;
24
for(int i=1;i<=n;i++)if(a[0][i]>eps)d[++tt]=i;
25
if(tt==0)break; e=d[rand()%tt+1]; long double temp=INF;
26
for(int i=1;i<=m;i++)if((a[i][e]>eps)and(a[i][0]/a[i][e]<temp)){ l=i; temp=a[i][0]/a[i][e]; }
27
if(temp==INF){ printf("Unboundedn"); return 0; }
28 
pivot(l,e);
29 
}
30
ans=-a[0][0]; return 1;
31 }
32 int main()
33 {
34
scanf("%d%d%dn",&n,&m,&t);
35
for(int i=1;i<=n;i++)scanf("%Lf",&c[i]); a[0][n+1]=-1;
36
for(int i=1;i<=n;i++)N[i]=i; N[n+1]=0;
37
for(int i=1;i<=m;i++)B[i]=n+i;
38
for(int i=1;i<=m;i++)
39 
{
40
for(int j=1;j<=n;j++)scanf("%Lf",&a[i][j]); a[i][n+1]=-1;
41
scanf("%Lf",&a[i][0]);
42 
}
43
n++; int mm=1; for(int i=2;i<=m;i++)if(a[i][0]<a[mm][0])mm=i;
44
if(a[mm][0]<-eps)pivot(mm,n);
45
if(!simplex())return 0;
46
if(ans<-eps){ printf("Infeasiblen"); return 0; }
47
xx=0; for(int i=1;i<=n;i++)if(N[i]==0){ xx=i; break; }
48
if(xx==0)
49 
{
50
yy=0;
51
for(int i=1;i<=m;i++)if(B[i]==0){ yy=i; break; }
52
for(xx=0;xx<=n;xx++)if(abs(a[yy][xx])>eps)break;
53
if(xx<=n)pivot(yy,xx);
54 
}
55
if(xx<=n)
56 
{
57
for(int i=xx;i<=n-1;i++)N[i]=N[i+1]; N[n]=0;
58
for(int i=1;i<=m;i++)for(int j=xx;j<=n-1;j++)a[i][j]=a[i][j+1];
59
n--; for(int i=0;i<=n;i++)a[0][i]=0;
60
for(int i=1;i<=n;i++)
61 
{
62
xx=0; for(int j=1;j<=n;j++)if(N[j]==i){ xx=j; break; }
63
if(xx>0)a[0][xx]+=c[i];else
64 
{
65
for(int j=1;j<=m;j++)if(B[j]==i){ xx=j; break; }
66
for(int j=0;j<=n;j++)a[0][j]=a[0][j]-c[i]*a[xx][j];
67 
}
68 
}
69
}else
70 
{
71
for(int i=yy;i<=m-1;i++)for(int j=0;j<=n;j++)a[i][j]=a[i+1][j];
72
m--;
73 
}
74
if(!simplex())return 0;
75
printf("%.10Lfn",ans);
76
if(t==1)
77 
{
78
for(int i=1;i<=m;i++)x[B[i]]=a[i][0];
79
for(int i=1;i<=n;i++)printf("%.10Lf ",x[i]);
80
printf("n");
81 
}
82
return 0;
83 }