001package algs91; // section 9.9 002import stdlib.*; 003/* *********************************************************************** 004 * Compilation: javac GaussJordanElimination.java 005 * Execution: java GaussJordanElimination N 006 * 007 * Finds a solutions to Ax = b using Gauss-Jordan elimination with partial 008 * pivoting. If no solution exists, find a solution to yA = 0, yb != 0, 009 * which serves as a certificate of infeasibility. 010 * 011 * % java GaussJordanElimination 012 * -1.000000 013 * 2.000000 014 * 2.000000 015 * 016 * 3.000000 017 * -1.000000 018 * -2.000000 019 * 020 * System is infeasible 021 * 022 * -6.250000 023 * -4.500000 024 * 0.000000 025 * 0.000000 026 * 1.000000 027 * 028 * System is infeasible 029 * 030 * -1.375000 031 * 1.625000 032 * 0.000000 033 * 034 * 035 *************************************************************************/ 036 037public class XGaussJordanElimination { 038 private static final double EPSILON = 1e-8; 039 040 private final int N; // N-by-N system 041 private final double[][] a; // N-by-N+1 augmented matrix 042 043 // Gauss-Jordan elimination with partial pivoting 044 public XGaussJordanElimination(double[][] A, double[] b) { 045 N = b.length; 046 047 // build augmented matrix 048 a = new double[N][N+N+1]; 049 for (int i = 0; i < N; i++) 050 for (int j = 0; j < N; j++) 051 a[i][j] = A[i][j]; 052 053 // only need if you want to find certificate of infeasibility (or compute inverse) 054 for (int i = 0; i < N; i++) 055 a[i][N+i] = 1.0; 056 057 for (int i = 0; i < N; i++) a[i][N+N] = b[i]; 058 059 solve(); 060 061 assert check(A, b); 062 } 063 064 private void solve() { 065 066 // Gauss-Jordan elimination 067 for (int p = 0; p < N; p++) { 068 // show(); 069 070 // find pivot row using partial pivoting 071 int max = p; 072 for (int i = p+1; i < N; i++) { 073 if (Math.abs(a[i][p]) > Math.abs(a[max][p])) { 074 max = i; 075 } 076 } 077 078 // exchange row p with row max 079 swap(p, max); 080 081 // singular or nearly singular 082 if (Math.abs(a[p][p]) <= EPSILON) { 083 continue; 084 // throw new Error("Matrix is singular or nearly singular"); 085 } 086 087 // pivot 088 pivot(p, p); 089 } 090 // show(); 091 } 092 093 // swap row1 and row2 094 private void swap(int row1, int row2) { 095 double[] temp = a[row1]; 096 a[row1] = a[row2]; 097 a[row2] = temp; 098 } 099 100 101 // pivot on entry (p, q) using Gauss-Jordan elimination 102 private void pivot(int p, int q) { 103 104 // everything but row p and column q 105 for (int i = 0; i < N; i++) { 106 double alpha = a[i][q] / a[p][q]; 107 for (int j = 0; j <= N+N; j++) { 108 if (i != p && j != q) a[i][j] -= alpha * a[p][j]; 109 } 110 } 111 112 // zero out column q 113 for (int i = 0; i < N; i++) 114 if (i != p) a[i][q] = 0.0; 115 116 // scale row p (ok to go from q+1 to N, but do this for consistency with simplex pivot) 117 for (int j = 0; j <= N+N; j++) 118 if (j != q) a[p][j] /= a[p][q]; 119 a[p][q] = 1.0; 120 } 121 122 // extract solution to Ax = b 123 public double[] primal() { 124 double[] x = new double[N]; 125 for (int i = 0; i < N; i++) { 126 if (Math.abs(a[i][i]) > EPSILON) 127 x[i] = a[i][N+N] / a[i][i]; 128 else if (Math.abs(a[i][N+N]) > EPSILON) 129 return null; 130 } 131 return x; 132 } 133 134 // extract solution to yA = 0, yb != 0 135 public double[] dual() { 136 double[] y = new double[N]; 137 for (int i = 0; i < N; i++) { 138 if ((Math.abs(a[i][i]) <= EPSILON) && (Math.abs(a[i][N+N]) > EPSILON)) { 139 for (int j = 0; j < N; j++) 140 y[j] = a[i][N+j]; 141 return y; 142 } 143 } 144 return null; 145 } 146 147 // does the system have a solution? 148 public boolean isFeasible() { 149 return primal() != null; 150 } 151 152 // print the tableaux 153 private void show() { 154 for (int i = 0; i < N; i++) { 155 for (int j = 0; j < N; j++) { 156 StdOut.format("%8.3f ", a[i][j]); 157 } 158 StdOut.format("| "); 159 for (int j = N; j < N+N; j++) { 160 StdOut.format("%8.3f ", a[i][j]); 161 } 162 StdOut.format("| %8.3f\n", a[i][N+N]); 163 } 164 StdOut.println(); 165 } 166 167 168 // check that Ax = b or yA = 0, yb != 0 169 private boolean check(double[][] A, double[] b) { 170 171 // check that Ax = b 172 if (isFeasible()) { 173 double[] x = primal(); 174 for (int i = 0; i < N; i++) { 175 double sum = 0.0; 176 for (int j = 0; j < N; j++) { 177 sum += A[i][j] * x[j]; 178 } 179 if (Math.abs(sum - b[i]) > EPSILON) { 180 StdOut.println("not feasible"); 181 StdOut.format("b[%d] = %8.3f, sum = %8.3f\n", i, b[i], sum); 182 return false; 183 } 184 } 185 return true; 186 } 187 188 // or that yA = 0, yb != 0 189 else { 190 double[] y = dual(); 191 for (int j = 0; j < N; j++) { 192 double sum = 0.0; 193 for (int i = 0; i < N; i++) { 194 sum += A[i][j] * y[i]; 195 } 196 if (Math.abs(sum) > EPSILON) { 197 StdOut.println("invalid certificate of infeasibility"); 198 StdOut.format("sum = %8.3f\n", sum); 199 return false; 200 } 201 } 202 double sum = 0.0; 203 for (int i = 0; i < N; i++) { 204 sum += y[i] * b[i]; 205 } 206 if (Math.abs(sum) < EPSILON) { 207 StdOut.println("invalid certificate of infeasibility"); 208 StdOut.format("yb = %8.3f\n", sum); 209 return false; 210 } 211 return true; 212 } 213 } 214 215 216 public static void test(double[][] A, double[] b) { 217 XGaussJordanElimination gaussian = new XGaussJordanElimination(A, b); 218 if (gaussian.isFeasible()) { 219 StdOut.println("Solution to Ax = b"); 220 double[] x = gaussian.primal(); 221 for (double element : x) { 222 StdOut.format("%10.6f\n", element); 223 } 224 } 225 else { 226 StdOut.println("Certificate of infeasibility"); 227 double[] y = gaussian.dual(); 228 for (double element : y) { 229 StdOut.format("%10.6f\n", element); 230 } 231 } 232 StdOut.println(); 233 } 234 235 236 // 3-by-3 nonsingular system 237 public static void test1() { 238 double[][] A = { 239 { 0, 1, 1 }, 240 { 2, 4, -2 }, 241 { 0, 3, 15 } 242 }; 243 double[] b = { 4, 2, 36 }; 244 test(A, b); 245 } 246 247 // 3-by-3 nonsingular system 248 public static void test2() { 249 double[][] A = { 250 { 1, -3, 1 }, 251 { 2, -8, 8 }, 252 { -6, 3, -15 } 253 }; 254 double[] b = { 4, -2, 9 }; 255 test(A, b); 256 } 257 258 // 5-by-5 singular: no solutions 259 // y = [ -1, 0, 1, 1, 0 ] 260 public static void test3() { 261 double[][] A = { 262 { 2, -3, -1, 2, 3 }, 263 { 4, -4, -1, 4, 11 }, 264 { 2, -5, -2, 2, -1 }, 265 { 0, 2, 1, 0, 4 }, 266 { -4, 6, 0, 0, 7 }, 267 }; 268 double[] b = { 4, 4, 9, -6, 5 }; 269 test(A, b); 270 } 271 272 // 5-by-5 singluar: infinitely many solutions 273 public static void test4() { 274 double[][] A = { 275 { 2, -3, -1, 2, 3 }, 276 { 4, -4, -1, 4, 11 }, 277 { 2, -5, -2, 2, -1 }, 278 { 0, 2, 1, 0, 4 }, 279 { -4, 6, 0, 0, 7 }, 280 }; 281 double[] b = { 4, 4, 9, -5, 5 }; 282 test(A, b); 283 } 284 285 // 3-by-3 singular: no solutions 286 // y = [ 1, 0, 1/3 ] 287 public static void test5() { 288 double[][] A = { 289 { 2, -1, 1 }, 290 { 3, 2, -4 }, 291 { -6, 3, -3 }, 292 }; 293 double[] b = { 1, 4, 2 }; 294 test(A, b); 295 } 296 297 // 3-by-3 singular: infinitely many solutions 298 public static void test6() { 299 double[][] A = { 300 { 1, -1, 2 }, 301 { 4, 4, -2 }, 302 { -2, 2, -4 }, 303 }; 304 double[] b = { -3, 1, 6 }; 305 test(A, b); 306 } 307 308 // sample client 309 public static void main(String[] args) { 310 311 try { test1(); } 312 catch (Exception e) { e.printStackTrace(); } 313 StdOut.println("--------------------------------"); 314 315 try { test2(); } 316 catch (Exception e) { e.printStackTrace(); } 317 StdOut.println("--------------------------------"); 318 319 try { test3(); } 320 catch (Exception e) { e.printStackTrace(); } 321 StdOut.println("--------------------------------"); 322 323 try { test4(); } 324 catch (Exception e) { e.printStackTrace(); } 325 StdOut.println("--------------------------------"); 326 327 try { test5(); } 328 catch (Exception e) { e.printStackTrace(); } 329 StdOut.println("--------------------------------"); 330 331 try { test6(); } 332 catch (Exception e) { e.printStackTrace(); } 333 StdOut.println("--------------------------------"); 334 335 // N-by-N random system (likely full rank) 336 int N = Integer.parseInt(args[0]); 337 double[][] A = new double[N][N]; 338 for (int i = 0; i < N; i++) 339 for (int j = 0; j < N; j++) 340 A[i][j] = StdRandom.uniform(1000); 341 double[] b = new double[N]; 342 for (int i = 0; i < N; i++) 343 b[i] = StdRandom.uniform(1000); 344 test(A, b); 345 346 StdOut.println("--------------------------------"); 347 348 // N-by-N random system (likely infeasible) 349 A = new double[N][N]; 350 for (int i = 0; i < N-1; i++) 351 for (int j = 0; j < N; j++) 352 A[i][j] = StdRandom.uniform(1000); 353 for (int i = 0; i < N-1; i++) { 354 double alpha = StdRandom.uniform(11) - 5.0; 355 for (int j = 0; j < N; j++) { 356 A[N-1][j] += alpha * A[i][j]; 357 } 358 } 359 b = new double[N]; 360 for (int i = 0; i < N; i++) 361 b[i] = StdRandom.uniform(1000); 362 test(A, b); 363 364 365 } 366 367}