gt:         double[] b = { 4, 2, 36, 42 };
context: /******************************************************************************
   *  Compilation:  javac GaussianElimination.java
   *  Execution:    java GaussianElimination m n
   *  Dependencies: StdOut.java
   * 
   *  Gaussian elimination with partial pivoting for m-by-n system.
   *
   *  % java GaussianElimination m n
   *  -1.000000
   *  2.000000
   *  2.000000
   *
   *  3.000000
   *  -1.000000
   *  -2.000000
   * 
   *  System is infeasible
   *
   *  -6.250000
   *  -4.500000
   *  0.000000
   *  0.000000
   *  1.000000
   *
   *  System is infeasible
   *
   *  -1.375000
   *  1.625000
   *  0.000000
   *
   ******************************************************************************/
  
  package edu.princeton.cs.algs4;
  
  /**
   *  The {@code GaussianElimination} data type provides methods
   *  to solve a linear system of equations <em>Ax</em> = <em>b</em>,
   *  where <em>A</em> is an <em>m</em>-by-<em>n</em> matrix
   *  and <em>b</em> is a length <em>n</em> vector.
   *  <p>
   *  This is a bare-bones implementation that uses Gaussian elimination
   *  with partial pivoting.
   *  See <a href = "https://algs4.cs.princeton.edu/99scientific/GaussianEliminationLite.java.html">GaussianEliminationLite.java</a>
   *  for a stripped-down version that assumes the matrix <em>A</em> is square
   *  and nonsingular. See {@link GaussJordanElimination} for an alternate
   *  implementation that uses Gauss-Jordan elimination.
   *  For an industrial-strength numerical linear algebra library,
   *  see <a href = "http://math.nist.gov/javanumerics/jama/">JAMA</a>.
   *  <p>
   *  For additional documentation, see
   *  <a href="https://algs4.cs.princeton.edu/99scientific">Section 9.9</a>
   *  <i>Algorithms, 4th Edition</i> by Robert Sedgewick and Kevin Wayne.
   *
   *  @author Robert Sedgewick
   *  @author Kevin Wayne
   */
  public class GaussianElimination {
      private static final double EPSILON = 1e-8;
  
      private final int m;      // number of rows
      private final int n;      // number of columns
      private double[][] a;     // m-by-(n+1) augmented matrix
  
      /**
       * Solves the linear system of equations <em>Ax</em> = <em>b</em>,
       * where <em>A</em> is an <em>m</em>-by-<em>n</em> matrix and <em>b</em>
       * is a length <em>m</em> vector.
       *
       * @param  A the <em>m</em>-by-<em>n</em> constraint matrix
       * @param  b the length <em>m</em> right-hand-side vector
       * @throws IllegalArgumentException if the dimensions disagree, i.e.,
       *         the length of {@code b} does not equal {@code m}
       */
      public GaussianElimination(double[][] A, double[] b) {
          m = A.length;
          n = A[0].length;
  
          if (b.length != m) throw new IllegalArgumentException("Dimensions disagree");
  
          // build augmented matrix
          a = new double[m][n+1];
          for (int i = 0; i < m; i++)
              for (int j = 0; j < n; j++)
                  a[i][j] = A[i][j];
          for (int i = 0; i < m; i++)
              a[i][n] = b[i];
  
          forwardElimination();
  
          assert certifySolution(A, b);
      }
  
      // forward elimination
      private void forwardElimination() {
          for (int p = 0; p < Math.min(m, n); p++) {
  
              // find pivot row using partial pivoting
              int max = p;
              for (int i = p+1; i < m; i++) {
                  if (Math.abs(a[i][p]) > Math.abs(a[max][p])) {
                      max = i;
                  }
              }
  
              // swap
              swap(p, max);
  
              // singular or nearly singular
              if (Math.abs(a[p][p]) <= EPSILON) {
                  continue;
              }
  
              // pivot
              pivot(p);
          }
      }
  
      // swap row1 and row2
      private void swap(int row1, int row2) {
          double[] temp = a[row1];
          a[row1] = a[row2];
          a[row2] = temp;
      }
  
      // pivot on a[p][p]
      private void pivot(int p) {
          for (int i = p+1; i < m; i++) {
              double alpha = a[i][p] / a[p][p];
              for (int j = p; j <= n; j++) {
                  a[i][j] -= alpha * a[p][j];
              }
          }
      }
  
      /**
       * Returns a solution to the linear system of equations <em>Ax</em> = <em>b</em>.
       *      
       * @return a solution <em>x</em> to the linear system of equations
       *         <em>Ax</em> = <em>b</em>; {@code null} if no such solution
       */
      public double[] primal() {
          // back substitution
          double[] x = new double[n];
          for (int i = Math.min(n-1, m-1); i >= 0; i--) {
              double sum = 0.0;
              for (int j = i+1; j < n; j++) {
                  sum += a[i][j] * x[j];
              }
  
              if (Math.abs(a[i][i]) > EPSILON)
                  x[i] = (a[i][n] - sum) / a[i][i];
              else if (Math.abs(a[i][n] - sum) > EPSILON)
                  return null;
          }
  
          // redundant rows
          for (int i = n; i < m; i++) {
              double sum = 0.0;
              for (int j = 0; j < n; j++) {
                  sum += a[i][j] * x[j];
              }
              if (Math.abs(a[i][n] - sum) > EPSILON)
                  return null;
          }
          return x;
      }
  
      /**
       * Returns true if there exists a solution to the linear system of
       * equations <em>Ax</em> = <em>b</em>.
       *      
       * @return {@code true} if there exists a solution to the linear system
       *         of equations <em>Ax</em> = <em>b</em>; {@code false} otherwise
       */
      public boolean isFeasible() {
          return primal() != null;
      }
  
  
      // check that Ax = b
      private boolean certifySolution(double[][] A, double[] b) {
          if (!isFeasible()) return true;
          double[] x = primal();
          for (int i = 0; i < m; i++) {
              double sum = 0.0;
              for (int j = 0; j < n; j++) {
                  sum += A[i][j] * x[j];
              }
              if (Math.abs(sum - b[i]) > EPSILON) {
                  StdOut.println("not feasible");
                  StdOut.println("b[" + i + "] = " + b[i] + ", sum = " + sum);
                  return false;
              }
          }
          return true;
      }
  
  
      /**
       * Unit tests the {@code GaussianElimination} data type.
       */
      private static void test(String name, double[][] A, double[] b) {
          StdOut.println("----------------------------------------------------");
          StdOut.println(name);
          StdOut.println("----------------------------------------------------");
          GaussianElimination gaussian = new GaussianElimination(A, b);
          double[] x = gaussian.primal();
          if (gaussian.isFeasible()) {
              for (int i = 0; i < x.length; i++) {
                  StdOut.printf("%.6f\n", x[i]);
              }
          }
          else {
              StdOut.println("System is infeasible");
          }
          StdOut.println();
          StdOut.println();
      }
  
  
      // 3-by-3 nonsingular system
      private static void test1() {
          double[][] A = {
              { 0, 1,  1 },
              { 2, 4, -2 },
              { 0, 3, 15 }
          };
          double[] b = { 4, 2, 36 };
          test("test 1 (3-by-3 system, nonsingular)", A, b);
      }
  
      // 3-by-3 nonsingular system
      private static void test2() {
          double[][] A = {
              {  1, -3,   1 },
              {  2, -8,   8 },
              { -6,  3, -15 }
          };
          double[] b = { 4, -2, 9 };
          test("test 2 (3-by-3 system, nonsingular)", A, b);
      }
  
      // 5-by-5 singular: no solutions
      private static void test3() {
          double[][] A = {
              {  2, -3, -1,  2,  3 },
              {  4, -4, -1,  4, 11 },
              {  2, -5, -2,  2, -1 },
              {  0,  2,  1,  0,  4 },
              { -4,  6,  0,  0,  7 },
          };
          double[] b = { 4, 4, 9, -6, 5 };
          test("test 3 (5-by-5 system, no solutions)", A, b);
      }
  
      // 5-by-5 singular: infinitely many solutions
      private static void test4() {
          double[][] A = {
              {  2, -3, -1,  2,  3 },
              {  4, -4, -1,  4, 11 },
              {  2, -5, -2,  2, -1 },
              {  0,  2,  1,  0,  4 },
              { -4,  6,  0,  0,  7 },
          };
          double[] b = { 4, 4, 9, -5, 5 };
          test("test 4 (5-by-5 system, infinitely many solutions)", A, b);
      }
  
      // 3-by-3 singular: no solutions
      private static void test5() {
          double[][] A = {
              {  2, -1,  1 },
              {  3,  2, -4 },
              { -6,  3, -3 },
          };
          double[] b = { 1, 4, 2 };
          test("test 5 (3-by-3 system, no solutions)", A, b);
      }
  
      // 3-by-3 singular: infinitely many solutions
      private static void test6() {
          double[][] A = {
              {  1, -1,  2 },
              {  4,  4, -2 },
              { -2,  2, -4 },
          };
          double[] b = { -3, 1, 6 };
          test("test 6 (3-by-3 system, infinitely many solutions)", A, b);
      }
  
      // 4-by-3 full rank and feasible system
      private static void test7() {
          double[][] A = {
              { 0, 1,  1 },
              { 2, 4, -2 },
              { 0, 3, 15 },
              { 2, 8, 14 }
          };
  