单纯形法模板

Posted 2020-10-19 G . H . O . S . T . R . E . A

 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("Unbounded\n"); return 0; }
28         pivot(l,e);
29     }
30     ans=-a[0][0]; return 1;
31 }
32 int main()
33 {
34     scanf("%d%d%d\n",&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("Infeasible\n"); 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("%.10Lf\n",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 }

View Code