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