单纯形

study from:

国家集训队2016论文

浅谈线性规划与对偶问题

(网络很多文章不适合初学者)

https://pan.baidu.com/s/1KtNmuqVSbQk08oqJiaQ4Ew

uoj179

http://uoj.ac/problem/179

Sol:

http://uoj.ac/submission/61872

https://github.com/sxysxy/math/blob/master/simplex/simplex.c (代码有点小问题,可以参考一下其注释)

其它题:

https://loj.ac/problem/149

https://loj.ac/problem/153

https://loj.ac/problem/154

注意怎么操作(哪里需要找到就退出,哪里需要找到最小/大值,使满足条件/更快到达结果)

  1 #include <cstdio>
  2 #include <algorithm>
  3 using namespace std;
  4 
  5 const int maxn=20+10;
  6 const int maxm=20+10;
  7 const double eps=1e-9;
  8 const double inf=1e18;
  9 double a[maxm][maxn],value[maxn+maxm];
 10 int base[maxm],unbase[maxn],n,m;
 11 
 12 int judge(double x)
 13 {
 14     if (x<-eps)
 15         return -1;
 16     else if (x>eps)
 17         return 1;
 18     else
 19         return 0;
 20 }
 21 
 22 void pivot(int x,int y)
 23 {
 24     int i,j;
 25     swap(base[x],unbase[y]);
 26     double k=a[x][y];
 27     a[x][y]=1.0;
 28     for (i=0;i<=n;i++)
 29         a[x][i]/=k;
 30     for (i=0;i<=m;i++)
 31         if (i!=x && judge(a[i][y]))
 32         {
 33             k=a[i][y];
 34             a[i][y]=0;
 35             a[i][0]-=(i?1:-1)*k*a[x][0];
 36             for (j=1;j<=n;j++)
 37                 a[i][j]-=k*a[x][j];
 38         }
 39 }
 40 
 41 bool init()
 42 {
 43     int i,x,y;
 44     for (i=1;i<=n;i++)
 45         unbase[i]=i;
 46     for (i=1;i<=m;i++)
 47         base[i]=n+i;
 48     double v;
 49     while (1)
 50     {
 51         x=y=0;
 52         v=-eps;
 53         for (i=1;i<=m;i++)
 54             if (a[i][0]<v)
 55             {
 56                 v=a[i][0];
 57                 x=i;
 58             }
 59         if (!x)
 60             return 1;
 61         for (i=1;i<=n;i++)
 62             if (a[x][i]<-eps)
 63             {
 64                 y=i;
 65                 break;
 66             }
 67         if (!y)
 68             return 0;
 69         pivot(x,y);
 70     }
 71     return 1;///
 72 }
 73 
 74 int work()
 75 {
 76     if (!init())
 77         return 0;
 78     int i,x,y;
 79     double v,t;
 80     while (1)
 81     {
 82         x=y=0;
 83         v=eps;
 84         for (i=1;i<=n;i++)
 85             if (a[0][i]>v)
 86             {
 87                 v=a[0][i];
 88                 y=i;
 89             }
 90         if (!y)
 91             return 1;
 92         v=inf;
 93         for (i=1;i<=m;i++)
 94         {
 95             t=a[i][0]/a[i][y];
 96             if (a[i][y]>eps && (!x || t<v))
 97             {
 98                 v=t;
 99                 x=i;
100             }
101         }
102         if (!x)
103             return -1;
104         pivot(x,y);
105     }
106     return 1;///
107 }
108 
109 int main()
110 {
111     int t,i,j;
112     scanf("%d%d%d",&n,&m,&t);
113     for (i=1;i<=n;i++)
114         scanf("%lf",&a[0][i]);
115     for (i=1;i<=m;i++)
116     {
117         for (j=1;j<=n;j++)
118             scanf("%lf",&a[i][j]);
119         scanf("%lf",&a[i][0]);
120     }
121     switch(work())
122     {
123         case -1:
124             printf("Unbounded");
125             break;
126         case 0:
127             printf("Infeasible");
128             break;
129         case 1:
130             printf("%.10f
",a[0][0]);
131             if (t==1)
132             {
133                 for (i=1;i<=m;i++)
134                     value[base[i]]=a[i][0];
135                 for (i=1;i<=n;i++)
136                     printf("%.10f%c",value[i],i==n?'
':' ');
137             }
138             break;
139     }
140     return 0;
141 }

原文地址:https://www.cnblogs.com/cmyg/p/10420915.html