001package algs91; // section 9.9
002import stdlib.*;
003/* ***********************************************************************
004 *  Compilation:  javac GaussianElimination.java
005 *  Execution:    java GaussianElimination
006 *
007 *  Gaussian elimination with partial pivoting.
008 *
009 *  % java GaussianElimination
010 *  -1.0
011 *  2.0
012 *  2.0
013 *
014 *************************************************************************/
015
016public class GaussianElimination {
017        private static final double EPSILON = 1e-10;
018
019        // Gaussian elimination with partial pivoting
020        public static double[] lsolve(double[][] A, double[] b) {
021                int N  = b.length;
022
023                for (int p = 0; p < N; p++) {
024
025                        // find pivot row and swap
026                        int max = p;
027                        for (int i = p + 1; i < N; i++) {
028                                if (Math.abs(A[i][p]) > Math.abs(A[max][p])) {
029                                        max = i;
030                                }
031                        }
032                        double[] temp = A[p]; A[p] = A[max]; A[max] = temp;
033                        double   t    = b[p]; b[p] = b[max]; b[max] = t;
034
035                        // singular or nearly singular
036                        if (Math.abs(A[p][p]) <= EPSILON) {
037                                throw new Error("Matrix is singular or nearly singular");
038                        }
039
040                        // pivot within A and b
041                        for (int i = p + 1; i < N; i++) {
042                                double alpha = A[i][p] / A[p][p];
043                                b[i] -= alpha * b[p];
044                                for (int j = p; j < N; j++) {
045                                        A[i][j] -= alpha * A[p][j];
046                                }
047                        }
048                }
049
050                // back substitution
051                double[] x = new double[N];
052                for (int i = N - 1; i >= 0; i--) {
053                        double sum = 0.0;
054                        for (int j = i + 1; j < N; j++) {
055                                sum += A[i][j] * x[j];
056                        }
057                        x[i] = (b[i] - sum) / A[i][i];
058                }
059                return x;
060        }
061
062
063        // sample client
064        public static void main(String[] args) {
065                int N = 3;
066                double[][] A = {
067                                { 0, 1,  1 },
068                                { 2, 4, -2 },
069                                { 0, 3, 15 }
070                };
071                double[] b = { 4, 2, 36 };
072                double[] x = lsolve(A, b);
073
074
075                // print results
076                for (int i = 0; i < N; i++) {
077                        StdOut.println(x[i]);
078                }
079
080        }
081
082}