001package edu.pdx.cs410J.rmi;
002
003/**
004 * This class performs Gaussian Elimination to solve a system of
005 * linear equations of the form <code>Ax = b</code>.
006 */
007public class GaussianElimination {
008
009  /**
010   * Solve a system <code>Ax = b</code> using Gaussian elimination
011   *
012   * @throws IllegalArgumentException
013   *         <code>A</code> is not an <code>n x n</code> matrix,
014   *         <code>b</code> does not have length <code>n</code>
015   */
016  public static double[] solve(double[][] a, double[] b) {
017    // First validate that 
018
019    int rows = Matrix.countRows(a);
020    int columns = Matrix.countColumns(a);
021    if (rows != columns) {
022      String s = "Matrix a must have the same number of rows (" + rows
023        + ") and columns (" + columns + ")";
024      throw new IllegalArgumentException(s);
025    }
026
027    if (b.length != rows) {
028      String s = "Vector b must have " + rows + " rows (not " +
029      b.length + ")";
030      throw new IllegalArgumentException(s);
031    }
032
033    java.io.PrintWriter pw = new java.io.PrintWriter(System.out, true);
034    Matrix.print("Matrix A", a, pw);
035    Matrix.print("Vector b", b, pw);
036
037    int n = rows;
038    double[] x = new double[n];
039
040    // Transform "a" into a upper-diagnol matrix
041
042    // Divide each element in the row by a[0][0]
043    double d = a[0][0];
044    for(int h = 0; h < n; h++) {
045      a[0][h] = a[0][h] / d;
046    }
047    b[0] = b[0] / d;
048    
049    for (int k = 0; k < n; k++) {
050      for (int i = k + 1; i < n; i++) {
051        double multiplier = a[i][k] / a[k][k];
052        for (int j = k; j < n; j++) {
053          a[i][j] = a[i][j] - (multiplier * a[k][j]);
054        }
055        b[i] = b[i] - (multiplier * b[k]);
056      }
057    }
058
059    // Now solve for x from the bottom up
060    for (int i = n-1; i >= 0; i--) {
061      double v = b[i];
062      for (int j = i+1; j < n; j++) {
063        v -= (a[i][j] * x[j]);
064      }
065      x[i] = v / a[i][i];
066    }
067
068    Matrix.print("Vector x", x, pw);
069    
070    return x;
071  }
072
073}