线性规划模板题,贴个模板
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 |
// <simplex.cpp> - 05/28/16 08:01:08 // This file is created by XuYike's black technology automatically. // Copyright (C) 2015 ChangJun High School, Inc. // I don't know what this program is. #include <iostream> #include <vector> #include <algorithm> #include <cstring> #include <cstdio> #include <cmath> #include <ctime> #define double long double using namespace std; typedef long long lol; int gi(){ int res=0,fh=1;char ch=getchar(); while((ch>'9'||ch<'0')&&ch!='-')ch=getchar(); if(ch=='-')fh=-1,ch=getchar(); while(ch>='0'&&ch<='9')res=res*10+ch-'0',ch=getchar(); return fh*res; } const int MAXN=41; const int INF=1e9; const double eps=1e-12; double a[MAXN][MAXN]; int n,m,L[MAXN],E[MAXN],b[MAXN]; double val[MAXN]; void pivot(int l,int e){ swap(L[l],E[e]); double k=a[l][e];a[l][e]=1; for(int i=0;i<=n;i++)a[l][i]/=k; int len=0; for(int i=0;i<=n;i++)if(fabs(a[l][i])>eps)b[++len]=i; for(int i=0;i<=m;i++) if(i!=l&&fabs(a[i][e])>eps){ k=a[i][e];a[i][e]=0; for(int j=1;j<=len;j++)a[i][b[j]]-=k*a[l][b[j]]; } } double simplex(){ while(1){ int l,e; for(e=1;e<=n;e++)if(a[0][e]>eps)break;//选择第一个正的c[e],据说可以防止死循环 if(e==n+1)return -a[0][0];//若所有的c都为负数那么显然都取0为最大值 double tmp=INF; for(int i=1;i<=m;i++) if(a[i][e]>eps&&a[i][0]/a[i][e]<tmp)tmp=a[i][0]/a[i][e],l=i; //找到最小的b[i]/a[i][e],观察pivot函数可以发现这样做能够防止出现负的b if(tmp==INF)return INF;//找不到说明e可以无限增大增加答案而不破坏约束条件 pivot(l,e); } } bool init(){ for(int i=1;i<=n;i++)E[i]=i; while(1){ int l=0,e=0; for(int i=1;i<=m;i++)if(a[i][0]<-eps&&(!l||(rand()&1)))l=i; if(!l)return 1; for(int i=1;i<=n;i++)if(a[l][i]<-eps&&(!e||(rand()&1)))e=i; if(!e)return 0; pivot(l,e); } } int main(){ srand(time(NULL)); n=gi(),m=gi();int t=gi(); for(int i=1;i<=n;i++)scanf("%Lf",&a[0][i]); for(int i=1;i<=m;i++){ for(int j=1;j<=n;j++)scanf("%Lf",&a[i][j]); scanf("%Lf",&a[i][0]); } if(!init()){printf("Infeasible");return 0;} double ans=simplex(); if(ans==INF)printf("Unbounded"); else{ printf("%.8Lf\n",ans); if(!t)return 0; for(int i=1;i<=m;i++)val[L[i]]=a[i][0]; for(int i=1;i<=n;i++)printf("%.8Lf ",val[i]); } return 0; } |
膜拜会单纯形的神犇。
(顺便求友链)