change Simplex

code/Simplex.cc

@@ -1,78 +1,59 @@
 // 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;
-
+// USAGE: create an LPSolver object with A, b, and c as
+// arguments. Then, call Solve(x).
 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) : 
+  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];
+        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 (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 (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;
@@ -80,10 +61,10 @@ struct LPSolver {
       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);
+        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();
@@ -92,32 +73,3 @@ struct LPSolver {
     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;
-}

main.tex

@@ -112,7 +112,7 @@ Temps de cuisson : $O(n)$
 {\scriptsize\lstinputlisting{code/FFT.cpp}} % OK
 
 \subsection{Simplex algorithm}
-{\scriptsize\lstinputlisting{code/Simplex.cc}}
+{\scriptsize\lstinputlisting{code/Simplex.cc}} % OK
 
 \section{Graphes}