001package algs9; // 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}