123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123 |
- // Two-phase simplex algorithm for solving linear programs of the form
- //
- // maximize c^T x
- // subject to Ax <= b
- // x >= 0
- //
- // INPUT: A -- an m x n matrix
- // b -- an m-dimensional vector
- // c -- an n-dimensional vector
- // x -- a vector where the optimal solution will be stored
- //
- // OUTPUT: value of the optimal solution (infinity if unbounded
- // above, nan if infeasible)
- //
- // To use this code, create an LPSolver object with A, b, and c as
- // arguments. Then, call Solve(x).
- #include <iostream>
- #include <iomanip>
- #include <vector>
- #include <cmath>
- #include <limits>
- using namespace std;
- typedef long double DOUBLE;
- typedef vector<DOUBLE> VD;
- typedef vector<VD> VVD;
- typedef vector<int> VI;
- const DOUBLE EPS = 1e-9;
- struct LPSolver {
- int m, n;
- VI B, N;
- VVD D;
- LPSolver(const VVD &A, const VD &b, const VD &c) :
- m(b.size()), n(c.size()), N(n+1), B(m), D(m+2, VD(n+2)) {
- for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) D[i][j] = A[i][j];
- for (int i = 0; i < m; i++) { B[i] = n+i; D[i][n] = -1; D[i][n+1] = b[i]; }
- for (int j = 0; j < n; j++) { N[j] = j; D[m][j] = -c[j]; }
- N[n] = -1; D[m+1][n] = 1;
- }
-
- void Pivot(int r, int s) {
- for (int i = 0; i < m+2; i++) if (i != r)
- for (int j = 0; j < n+2; j++) if (j != s)
- D[i][j] -= D[r][j] * D[i][s] / D[r][s];
- for (int j = 0; j < n+2; j++) if (j != s) D[r][j] /= D[r][s];
- for (int i = 0; i < m+2; i++) if (i != r) D[i][s] /= -D[r][s];
- D[r][s] = 1.0 / D[r][s];
- swap(B[r], N[s]);
- }
- bool Simplex(int phase) {
- int x = phase == 1 ? m+1 : m;
- while (true) {
- int s = -1;
- for (int j = 0; j <= n; j++) {
- if (phase == 2 && N[j] == -1) continue;
- if (s == -1 || D[x][j] < D[x][s] || D[x][j] == D[x][s] && N[j] < N[s]) s = j;
- }
- if (D[x][s] >= -EPS) return true;
- int r = -1;
- for (int i = 0; i < m; i++) {
- if (D[i][s] <= 0) continue;
- if (r == -1 || D[i][n+1] / D[i][s] < D[r][n+1] / D[r][s] ||
- D[i][n+1] / D[i][s] == D[r][n+1] / D[r][s] && B[i] < B[r]) r = i;
- }
- if (r == -1) return false;
- Pivot(r, s);
- }
- }
- DOUBLE Solve(VD &x) {
- int r = 0;
- for (int i = 1; i < m; i++) if (D[i][n+1] < D[r][n+1]) r = i;
- if (D[r][n+1] <= -EPS) {
- Pivot(r, n);
- if (!Simplex(1) || D[m+1][n+1] < -EPS) return -numeric_limits<DOUBLE>::infinity();
- for (int i = 0; i < m; i++) if (B[i] == -1) {
- int s = -1;
- for (int j = 0; j <= n; j++)
- if (s == -1 || D[i][j] < D[i][s] || D[i][j] == D[i][s] && N[j] < N[s]) s = j;
- Pivot(i, s);
- }
- }
- if (!Simplex(2)) return numeric_limits<DOUBLE>::infinity();
- x = VD(n);
- for (int i = 0; i < m; i++) if (B[i] < n) x[B[i]] = D[i][n+1];
- return D[m][n+1];
- }
- };
- int main() {
- const int m = 4;
- const int n = 3;
- DOUBLE _A[m][n] = {
- { 6, -1, 0 },
- { -1, -5, 0 },
- { 1, 5, 1 },
- { -1, -5, -1 }
- };
- DOUBLE _b[m] = { 10, -4, 5, -5 };
- DOUBLE _c[n] = { 1, -1, 0 };
-
- VVD A(m);
- VD b(_b, _b + m);
- VD c(_c, _c + n);
- for (int i = 0; i < m; i++) A[i] = VD(_A[i], _A[i] + n);
- LPSolver solver(A, b, c);
- VD x;
- DOUBLE value = solver.Solve(x);
-
- cerr << "VALUE: "<< value << endl;
- cerr << "SOLUTION:";
- for (size_t i = 0; i < x.size(); i++) cerr << " " << x[i];
- cerr << endl;
- return 0;
- }
|