Below is the syntax highlighted version of TwoPhaseSimplex.java
from §6.5 Reductions.
/****************************************************************************** * Compilation: javac TwoPhaseSimplex.java * Execution: java TwoPhaseSimplex * Dependencies: StdOut.java * * Given an m-by-n matrix A, an m-length vector b, and an * n-length vector c, solve the LP { max cx : Ax <= b, x >= 0 }. * Unlike Simplex.java, this version does not assume b >= 0, * so it needs to find a basic feasible solution in Phase I. * * Creates an (m+1)-by-(n+m+1) simplex tableaux with the * RHS in column m+n, the objective function in row m, and * slack variables in columns m through m+n-1. * ******************************************************************************/ public class TwoPhaseSimplex { private static final double EPSILON = 1.0E-8; private double[][] a; // tableaux // row m = objective function // row m+1 = artificial objective function // column n to n+m-1 = slack variables // column n+m to n+m+m-1 = artificial variables private int m; // number of constraints private int n; // number of original variables private int[] basis; // basis[i] = basic variable corresponding to row i // sets up the simplex tableaux public TwoPhaseSimplex(double[][] A, double[] b, double[] c) { m = b.length; n = c.length; a = new double[m+2][n+m+m+1]; for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) a[i][j] = A[i][j]; for (int i = 0; i < m; i++) a[i][n+i] = 1.0; for (int i = 0; i < m; i++) a[i][n+m+m] = b[i]; for (int j = 0; j < n; j++) a[m][j] = c[j]; // if negative RHS, multiply by -1 for (int i = 0; i < m; i++) { if (b[i] < 0) { a[i][n+m+m] = -b[i]; for (int j = 0; j <= n; j++) a[i][j] = -a[i][j]; a[i][n+i] = -1.0; } } // artificial variables form initial basis for (int i = 0; i < m; i++) a[i][n+m+i] = 1.0; for (int i = 0; i < m; i++) a[m+1][n+m+i] = -1.0; for (int i = 0; i < m; i++) pivot(i, n+m+i); basis = new int[m]; for (int i = 0; i < m; i++) basis[i] = n + m + i; // StdOut.println("before phase I"); // show(); phase1(); // StdOut.println("before phase II"); // show(); phase2(); // StdOut.println("after phase II"); // show(); // check optimality conditions assert check(A, b, c); } // run phase I simplex algorithm to find basic feasible solution private void phase1() { while (true) { // find entering column q int q = bland1(); if (q == -1) break; // optimal // find leaving row p int p = minRatioRule(q); assert p != -1 : "Entering column = " + q; // pivot pivot(p, q); // update basis basis[p] = q; // show(); } if (a[m+1][n+m+m] > EPSILON) throw new ArithmeticException("Linear program is infeasible"); } // run simplex algorithm starting from initial basic feasible solution private void phase2() { while (true) { // find entering column q int q = bland2(); if (q == -1) break; // optimal // find leaving row p int p = minRatioRule(q); if (p == -1) throw new ArithmeticException("Linear program is unbounded"); // pivot pivot(p, q); // update basis basis[p] = q; } } // lowest index of a non-basic column with a positive cost - using artificial objective function private int bland1() { for (int j = 0; j < n+m; j++) if (a[m+1][j] > EPSILON) return j; return -1; // optimal } // lowest index of a non-basic column with a positive cost private int bland2() { for (int j = 0; j < n+m; j++) if (a[m][j] > EPSILON) return j; return -1; // optimal } // find row p using min ratio rule (-1 if no such row) private int minRatioRule(int q) { int p = -1; for (int i = 0; i < m; i++) { // if (a[i][q] <= 0) continue; if (a[i][q] <= EPSILON) continue; else if (p == -1) p = i; else if ((a[i][n+m+m] / a[i][q]) < (a[p][n+m+m] / a[p][q])) p = i; } return p; } // pivot on entry (p, q) using Gauss-Jordan elimination private void pivot(int p, int q) { // everything but row p and column q for (int i = 0; i <= m+1; i++) for (int j = 0; j <= n+m+m; j++) if (i != p && j != q) a[i][j] -= a[p][j] * (a[i][q] / a[p][q]); // zero out column q for (int i = 0; i <= m+1; i++) if (i != p) a[i][q] = 0.0; // scale row p for (int j = 0; j <= n+m+m; j++) if (j != q) a[p][j] /= a[p][q]; a[p][q] = 1.0; } // return optimal objective value public double value() { return -a[m][n+m+m]; } // return primal solution vector public double[] primal() { double[] x = new double[n]; for (int i = 0; i < m; i++) if (basis[i] < n) x[basis[i]] = a[i][n+m+m]; return x; } // return dual solution vector public double[] dual() { double[] y = new double[m]; for (int i = 0; i < m; i++) { y[i] = -a[m][n+i]; if (y[i] == -0.0) y[i] = 0.0; } return y; } // is the solution primal feasible? private boolean isPrimalFeasible(double[][] A, double[] b) { double[] x = primal(); // check that x >= 0 for (int j = 0; j < x.length; j++) { if (x[j] < -EPSILON) { StdOut.println("x[" + j + "] = " + x[j] + " is negative"); return false; } } // check that Ax <= b for (int i = 0; i < m; i++) { double sum = 0.0; for (int j = 0; j < n; j++) { sum += A[i][j] * x[j]; } if (sum > b[i] + EPSILON) { StdOut.println("not primal feasible"); StdOut.println("b[" + i + "] = " + b[i] + ", A_i x = " + sum); return false; } } return true; } // is the solution dual feasible? private boolean isDualFeasible(double[][] A, double[] c) { double[] y = dual(); // check that y >= 0 for (int i = 0; i < y.length; i++) { if (y[i] < -EPSILON) { StdOut.println("y[" + i + "] = " + y[i] + " is negative"); return false; } } // check that yA >= c for (int j = 0; j < n; j++) { double sum = 0.0; for (int i = 0; i < m; i++) { sum += A[i][j] * y[i]; } if (sum < c[j] - EPSILON) { StdOut.println("not dual feasible"); StdOut.println("c[" + j + "] = " + c[j] + ", y A_j = " + sum); return false; } } return true; } // check that optimal value = cx = yb private boolean isOptimal(double[] b, double[] c) { double[] x = primal(); double[] y = dual(); double value = value(); // check that value = cx = yb double value1 = 0.0; for (int j = 0; j < x.length; j++) value1 += c[j] * x[j]; double value2 = 0.0; for (int i = 0; i < y.length; i++) value2 += y[i] * b[i]; if (Math.abs(value - value1) > EPSILON || Math.abs(value - value2) > EPSILON) { StdOut.println("value = " + value + ", cx = " + value1 + ", yb = " + value2); return false; } return true; } private boolean check(double[][]A, double[] b, double[] c) { return isPrimalFeasible(A, b) && isDualFeasible(A, c) && isOptimal(b, c); } // print tableaux public void show() { StdOut.println("m = " + m); StdOut.println("n = " + n); for (int i = 0; i <= m+1; i++) { for (int j = 0; j <= n+m+m; j++) { StdOut.printf("%7.2f ", a[i][j]); if (j == n+m-1 || j == n+m+m-1) StdOut.print(" |"); } StdOut.println(); } StdOut.print("basis = "); for (int i = 0; i < m; i++) StdOut.print(basis[i] + " "); StdOut.println(); StdOut.println(); } public static void test(double[][] A, double[] b, double[] c) { TwoPhaseSimplex lp; try { lp = new TwoPhaseSimplex(A, b, c); } catch (ArithmeticException e) { System.out.println(e); return; } StdOut.println("value = " + lp.value()); double[] x = lp.primal(); for (int i = 0; i < x.length; i++) StdOut.println("x[" + i + "] = " + x[i]); double[] y = lp.dual(); for (int j = 0; j < y.length; j++) StdOut.println("y[" + j + "] = " + y[j]); } // x0 = 12, x1 = 28, opt = 800 public static void test1() { double[] c = { 13.0, 23.0 }; double[] b = { 480.0, 160.0, 1190.0 }; double[][] A = { { 5.0, 15.0 }, { 4.0, 4.0 }, { 35.0, 20.0 }, }; test(A, b, c); } // dual of test1(): x0 = 12, x1 = 28, opt = 800 public static void test2() { double[] b = { -13.0, -23.0 }; double[] c = { -480.0, -160.0, -1190.0 }; double[][] A = { { -5.0, -4.0, -35.0 }, { -15.0, -4.0, -20.0 } }; test(A, b, c); } public static void test3() { double[][] A = { { -1, 1, 0 }, { 1, 4, 0 }, { 2, 1, 0 }, { 3, -4, 0 }, { 0, 0, 1 }, }; double[] c = { 1, 1, 1 }; double[] b = { 5, 45, 27, 24, 4 }; test(A, b, c); } // unbounded public static void test4() { double[] c = { 2.0, 3.0, -1.0, -12.0 }; double[] b = { 3.0, 2.0 }; double[][] A = { { -2.0, -9.0, 1.0, 9.0 }, { 1.0, 1.0, -1.0, -2.0 }, }; test(A, b, c); } // degenerate - cycles if you choose most positive objective function coefficient public static void test5() { double[] c = { 10.0, -57.0, -9.0, -24.0 }; double[] b = { 0.0, 0.0, 1.0 }; double[][] A = { { 0.5, -5.5, -2.5, 9.0 }, { 0.5, -1.5, -0.5, 1.0 }, { 1.0, 0.0, 0.0, 0.0 }, }; test(A, b, c); } // floating-point EPSILON needed in min-ratio test public static void test6() { double[] c = { -1, -1, -1, -1, -1, -1, -1, -1, -1 }; double[] b = { -0.9, 0.2, -0.2, -1.1, -0.7, -0.5, -0.1, -0.1, -1 }; double[][] A = { { -2, 1, 0, 0, 0, 0, 0, 0, 0 }, { 1, -2, -1, 0, 0, 0, 0, 0, 0 }, { 0, -1, -2, -1, 0, 0, 0, 0, 0 }, { 0, 0, -1, -2, -1, -1, 0, 0, 0 }, { 0, 0, 0, -1, -2, -1, 0, 0, 0 }, { 0, 0, 0, -1, -1, -2, -1, 0, 0 }, { 0, 0, 0, 0, 0, -1, -2, -1, 0 }, { 0, 0, 0, 0, 0, 0, -1, -2, -1 }, { 0, 0, 0, 0, 0, 0, 0, -1, -2 } }; test(A, b, c); } // test client public static void main(String[] args) { StdOut.println("----- test 1 --------------------"); test1(); StdOut.println(); StdOut.println("----- test 2 --------------------"); test2(); StdOut.println(); StdOut.println("----- test 3 --------------------"); test3(); StdOut.println(); StdOut.println("----- test 4 --------------------"); test4(); StdOut.println(); StdOut.println("----- test 5 --------------------"); test5(); StdOut.println(); StdOut.println("----- test 6 --------------------"); test6(); StdOut.println("----- test random ---------------"); int m = Integer.parseInt(args[0]); int n = Integer.parseInt(args[1]); double[] c = new double[n]; double[] b = new double[m]; double[][] A = new double[m][n]; for (int j = 0; j < n; j++) c[j] = StdRandom.uniformInt(1000); for (int i = 0; i < m; i++) b[i] = StdRandom.uniformInt(1000) - 200; for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) A[i][j] = StdRandom.uniformInt(100) - 20; test(A, b, c); } }