Merge branch 'master' of https://source.jalview.org/git/jalviewjs.git
[jalviewjs.git] / src / javajs / util / Matrix.java
index b188b7b..cea5a04 100644 (file)
-package javajs.util;\r
-\r
-/**\r
- * \r
- * streamlined and refined for Jmol by Bob Hanson\r
- * \r
- * from http://math.nist.gov/javanumerics/jama/\r
- * \r
- * Jama = Java Matrix class.\r
- * \r
- * @author The MathWorks, Inc. and the National Institute of Standards and\r
- *         Technology.\r
- * @version 5 August 1998\r
- */\r
-\r
-public class Matrix implements Cloneable {\r
-\r
-  public double[][] a;\r
-  protected int m, n;\r
-\r
-  /**\r
-   * Construct a matrix quickly without checking arguments.\r
-   * \r
-   * @param a\r
-   *        Two-dimensional array of doubles or null\r
-   * @param m\r
-   *        Number of rows.\r
-   * @param n\r
-   *        Number of colums.\r
-   */\r
-\r
-  public Matrix(double[][] a, int m, int n) {\r
-    this.a = (a == null ? new double[m][n] : a);\r
-    this.m = m;\r
-    this.n = n;\r
-  }\r
-\r
-  /**\r
-   * Get row dimension.\r
-   * \r
-   * @return m, the number of rows.\r
-   */\r
-\r
-  public int getRowDimension() {\r
-    return m;\r
-  }\r
-\r
-  /**\r
-   * Get column dimension.\r
-   * \r
-   * @return n, the number of columns.\r
-   */\r
-\r
-  public int getColumnDimension() {\r
-    return n;\r
-  }\r
-\r
-  /**\r
-   * Access the internal two-dimensional array.\r
-   * \r
-   * @return Pointer to the two-dimensional array of matrix elements.\r
-   */\r
-\r
-  public double[][] getArray() {\r
-    return a;\r
-  }\r
-\r
-  /**\r
-   * Copy the internal two-dimensional array.\r
-   * \r
-   * @return Two-dimensional array copy of matrix elements.\r
-   */\r
-\r
-  public double[][] getArrayCopy() {\r
-    double[][] x = new double[m][n];\r
-    for (int i = m; --i >= 0;)\r
-      for (int j = n; --j >= 0;)\r
-        x[i][j] = a[i][j];\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * Make a deep copy of a matrix\r
-   * \r
-   * @return copy\r
-   */\r
-\r
-  public Matrix copy() {\r
-    Matrix x = new Matrix(null, m, n);\r
-    double[][] c = x.a;\r
-    for (int i = m; --i >= 0;)\r
-      for (int j = n; --j >= 0;)\r
-        c[i][j] = a[i][j];\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * Clone the Matrix object.\r
-   */\r
-\r
-  @Override\r
-  public Object clone() {\r
-    return copy();\r
-  }\r
-\r
-  /**\r
-   * Get a submatrix.\r
-   * \r
-   * @param i0\r
-   *        Initial row index\r
-   * @param j0\r
-   *        Initial column index\r
-   * @param nrows\r
-   *        Number of rows\r
-   * @param ncols\r
-   *        Number of columns\r
-   * @return submatrix\r
-   * \r
-   */\r
-\r
-  public Matrix getSubmatrix(int i0, int j0, int nrows, int ncols) {\r
-    Matrix x = new Matrix(null, nrows, ncols);\r
-    double[][] xa = x.a;\r
-    for (int i = nrows; --i >= 0;)\r
-      for (int j = ncols; --j >= 0;)\r
-        xa[i][j] = a[i0 + i][j0 + j];\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * Get a submatrix for a give number of columns and selected row set.\r
-   * \r
-   * @param r\r
-   *        Array of row indices.\r
-   * @param n\r
-   *        number of rows \r
-   * @return submatrix\r
-   */\r
-\r
-  public Matrix getMatrixSelected(int[] r, int n) {\r
-    Matrix x = new Matrix(null, r.length, n);\r
-    double[][] xa = x.a;\r
-    for (int i = r.length; --i >= 0;) {\r
-      double[] b = a[r[i]];\r
-      for (int j = n; --j >= 0;)\r
-        xa[i][j] = b[j];\r
-    }\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * Matrix transpose.\r
-   * \r
-   * @return A'\r
-   */\r
-\r
-  public Matrix transpose() {\r
-    Matrix x = new Matrix(null, n, m);\r
-    double[][] c = x.a;\r
-    for (int i = m; --i >= 0;)\r
-      for (int j = n; --j >= 0;)\r
-        c[j][i] = a[i][j];\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * add two matrices\r
-   * @param b\r
-   * @return new Matrix this + b\r
-   */\r
-  public Matrix add(Matrix b) {\r
-    return scaleAdd(b, 1);\r
-  }\r
-\r
-  /**\r
-   * subtract two matrices\r
-   * @param b\r
-   * @return new Matrix this - b\r
-   */\r
-  public Matrix sub(Matrix b) {\r
-    return scaleAdd(b, -1);\r
-  }\r
-  \r
-  /**\r
-   * X = A + B*scale\r
-   * @param b \r
-   * @param scale \r
-   * @return X\r
-   * \r
-   */\r
-  public Matrix scaleAdd(Matrix b, double scale) {\r
-    Matrix x = new Matrix(null, m, n);\r
-    double[][] xa = x.a;\r
-    double[][] ba = b.a;\r
-    for (int i = m; --i >= 0;)\r
-      for (int j = n; --j >= 0;)\r
-        xa[i][j] = ba[i][j] * scale + a[i][j];\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * Linear algebraic matrix multiplication, A * B\r
-   * \r
-   * @param b\r
-   *        another matrix\r
-   * @return Matrix product, A * B or null for wrong dimension\r
-   */\r
-\r
-  public Matrix mul(Matrix b) {\r
-    if (b.m != n)\r
-      return null;\r
-    Matrix x = new Matrix(null, m, b.n);\r
-    double[][] xa = x.a;\r
-    double[][] ba = b.a;\r
-    for (int j = b.n; --j >= 0;)\r
-      for (int i = m; --i >= 0;) {\r
-        double[] arowi = a[i];\r
-        double s = 0;\r
-        for (int k = n; --k >= 0;)\r
-          s += arowi[k] * ba[k][j];\r
-        xa[i][j] = s;\r
-      }\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * Matrix inverse or pseudoinverse\r
-   * \r
-   * @return inverse (m == n) or pseudoinverse (m != n)\r
-   */\r
-\r
-  public Matrix inverse() {\r
-    return new LUDecomp(m, n).solve(identity(m, m), n);\r
-  }\r
-\r
-  /**\r
-   * Matrix trace.\r
-   * \r
-   * @return sum of the diagonal elements.\r
-   */\r
-\r
-  public double trace() {\r
-    double t = 0;\r
-    for (int i = Math.min(m, n); --i >= 0;)\r
-      t += a[i][i];\r
-    return t;\r
-  }\r
-\r
-  /**\r
-   * Generate identity matrix\r
-   * \r
-   * @param m\r
-   *        Number of rows.\r
-   * @param n\r
-   *        Number of columns.\r
-   * @return An m-by-n matrix with ones on the diagonal and zeros elsewhere.\r
-   */\r
-\r
-  public static Matrix identity(int m, int n) {\r
-    Matrix x = new Matrix(null, m, n);\r
-    double[][] xa = x.a;\r
-    for (int i = Math.min(m, n); --i >= 0;)\r
-      xa[i][i] = 1;\r
-    return x;\r
-  }\r
-\r
-  /**\r
-   * similarly to M3/M4 standard rotation/translation matrix\r
-   * we set a rotationTranslation matrix to be:\r
-   * \r
-   * [   nxn rot    nx1 trans\r
-   * \r
-   *     1xn  0     1x1 1      ]\r
-   * \r
-   * \r
-   * @return rotation matrix\r
-   */\r
-  public Matrix getRotation() {\r
-    return getSubmatrix(0, 0, m - 1, n - 1);\r
-  }\r
-\r
-  public Matrix getTranslation() {\r
-    return getSubmatrix(0, n - 1, m - 1, 1);\r
-  }\r
-\r
-  public static Matrix newT(T3 r, boolean asColumn) {\r
-    return (asColumn ? new Matrix(new double[][] { new double[] { r.x },\r
-        new double[] { r.y }, new double[] { r.z } }, 3, 1) : new Matrix(\r
-        new double[][] { new double[] { r.x, r.y, r.z } }, 1, 3));\r
-  }\r
-\r
-  @Override\r
-  public String toString() {\r
-    String s = "[\n";\r
-    for (int i = 0; i < m; i++) {\r
-      s += "  [";\r
-      for (int j = 0; j < n; j++)\r
-        s += " " + a[i][j];\r
-      s += "]\n";\r
-    }\r
-    s += "]";\r
-    return s;\r
-  }\r
-\r
-  /**\r
-   * \r
-   * Edited down by Bob Hanson for minimum needed by Jmol -- just constructor\r
-   * and solve\r
-   * \r
-   * LU Decomposition.\r
-   * <P>\r
-   * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit\r
-   * lower triangular matrix L, an n-by-n upper triangular matrix U, and a\r
-   * permutation vector piv of length m so that A(piv,:) = L*U. If m < n, then L\r
-   * is m-by-m and U is m-by-n.\r
-   * <P>\r
-   * The LU decompostion with pivoting always exists, even if the matrix is\r
-   * singular, so the constructor will never fail. The primary use of the LU\r
-   * decomposition is in the solution of square systems of simultaneous linear\r
-   * equations. This will fail if isNonsingular() returns false.\r
-   */\r
-\r
-  private class LUDecomp {\r
-\r
-    /* ------------------------\r
-       Class variables\r
-     * ------------------------ */\r
-\r
-    /**\r
-     * Array for internal storage of decomposition.\r
-     * \r
-     */\r
-    private double[][] LU;\r
-\r
-    /**\r
-     * Internal storage of pivot vector.\r
-     * \r
-     */\r
-    private int[] piv;\r
-\r
-    private int pivsign;\r
-\r
-    /* ------------------------\r
-       Constructor\r
-     * ------------------------ */\r
-\r
-    /**\r
-     * LU Decomposition Structure to access L, U and piv.\r
-     * @param m \r
-     * @param n \r
-     * \r
-     */\r
-\r
-    protected LUDecomp(int m, int n) {\r
-      \r
-      // Use a "left-looking", dot-product, Crout/Doolittle algorithm.\r
-\r
-      LU = getArrayCopy();\r
-      piv = new int[m];\r
-      for (int i = m; --i >= 0;)\r
-        piv[i] = i;\r
-      pivsign = 1;\r
-      double[] LUrowi;\r
-      double[] LUcolj = new double[m];\r
-\r
-      // Outer loop.\r
-\r
-      for (int j = 0; j < n; j++) {\r
-\r
-        // Make a copy of the j-th column to localize references.\r
-\r
-        for (int i = m; --i >= 0;)\r
-          LUcolj[i] = LU[i][j];\r
-\r
-        // Apply previous transformations.\r
-\r
-        for (int i = m; --i >= 0;) {\r
-          LUrowi = LU[i];\r
-\r
-          // Most of the time is spent in the following dot product.\r
-\r
-          int kmax = Math.min(i, j);\r
-          double s = 0.0;\r
-          for (int k = kmax; --k >= 0;)\r
-            s += LUrowi[k] * LUcolj[k];\r
-\r
-          LUrowi[j] = LUcolj[i] -= s;\r
-        }\r
-\r
-        // Find pivot and exchange if necessary.\r
-\r
-        int p = j;\r
-        for (int i = m; --i > j;)\r
-          if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p]))\r
-            p = i;\r
-        if (p != j) {\r
-          for (int k = n; --k >= 0;) {\r
-            double t = LU[p][k];\r
-            LU[p][k] = LU[j][k];\r
-            LU[j][k] = t;\r
-          }\r
-          int k = piv[p];\r
-          piv[p] = piv[j];\r
-          piv[j] = k;\r
-          pivsign = -pivsign;\r
-        }\r
-\r
-        // Compute multipliers.\r
-\r
-        if (j < m & LU[j][j] != 0.0)\r
-          for (int i = m; --i > j;)\r
-            LU[i][j] /= LU[j][j];\r
-      }\r
-    }\r
-\r
-    /* ------------------------\r
-       default Methods\r
-     * ------------------------ */\r
-\r
-    /**\r
-     * Solve A*X = B\r
-     * \r
-     * @param b\r
-     *        A Matrix with as many rows as A and any number of columns.\r
-     * @param n \r
-     * @return X so that L*U*X = B(piv,:) or null for wrong size or singular matrix\r
-     */\r
-\r
-    protected Matrix solve(Matrix b, int n) {\r
-      for (int j = 0; j < n; j++)\r
-        if (LU[j][j] == 0)\r
-          return null; // matrix is singular\r
-\r
-      // Copy right hand side with pivoting\r
-      int nx = b.n;\r
-      Matrix x = b.getMatrixSelected(piv, nx);\r
-      double[][] a = x.a;\r
-\r
-      // Solve L*Y = B(piv,:)\r
-      for (int k = 0; k < n; k++)\r
-        for (int i = k + 1; i < n; i++)\r
-          for (int j = 0; j < nx; j++)\r
-            a[i][j] -= a[k][j] * LU[i][k];\r
-\r
-      // Solve U*X = Y;\r
-      for (int k = n; --k >= 0;) {\r
-        for (int j = nx; --j >= 0;)\r
-          a[k][j] /= LU[k][k];\r
-        for (int i = k; --i >= 0;)\r
-          for (int j = nx; --j >= 0;)\r
-            a[i][j] -= a[k][j] * LU[i][k];\r
-      }\r
-      return x;\r
-    }\r
-  }\r
-\r
-}\r
+package javajs.util;
+
+/**
+ * 
+ * streamlined and refined for Jmol by Bob Hanson
+ * 
+ * from http://math.nist.gov/javanumerics/jama/
+ * 
+ * Jama = Java Matrix class.
+ * 
+ * @author The MathWorks, Inc. and the National Institute of Standards and
+ *         Technology.
+ * @version 5 August 1998
+ */
+
+public class Matrix implements Cloneable {
+
+  public double[][] a;
+  protected int m, n;
+
+  /**
+   * Construct a matrix quickly without checking arguments.
+   * 
+   * @param a
+   *        Two-dimensional array of doubles or null
+   * @param m
+   *        Number of rows.
+   * @param n
+   *        Number of colums.
+   */
+
+  public Matrix(double[][] a, int m, int n) {
+    this.a = (a == null ? new double[m][n] : a);
+    this.m = m;
+    this.n = n;
+  }
+
+  /**
+   * Get row dimension.
+   * 
+   * @return m, the number of rows.
+   */
+
+  public int getRowDimension() {
+    return m;
+  }
+
+  /**
+   * Get column dimension.
+   * 
+   * @return n, the number of columns.
+   */
+
+  public int getColumnDimension() {
+    return n;
+  }
+
+  /**
+   * Access the internal two-dimensional array.
+   * 
+   * @return Pointer to the two-dimensional array of matrix elements.
+   */
+
+  public double[][] getArray() {
+    return a;
+  }
+
+  /**
+   * Copy the internal two-dimensional array.
+   * 
+   * @return Two-dimensional array copy of matrix elements.
+   */
+
+  public double[][] getArrayCopy() {
+    double[][] x = new double[m][n];
+    for (int i = m; --i >= 0;)
+      for (int j = n; --j >= 0;)
+        x[i][j] = a[i][j];
+    return x;
+  }
+
+  /**
+   * Make a deep copy of a matrix
+   * 
+   * @return copy
+   */
+
+  public Matrix copy() {
+    Matrix x = new Matrix(null, m, n);
+    double[][] c = x.a;
+    for (int i = m; --i >= 0;)
+      for (int j = n; --j >= 0;)
+        c[i][j] = a[i][j];
+    return x;
+  }
+
+  /**
+   * Clone the Matrix object.
+   */
+
+  @Override
+  public Object clone() {
+    return copy();
+  }
+
+  /**
+   * Get a submatrix.
+   * 
+   * @param i0
+   *        Initial row index
+   * @param j0
+   *        Initial column index
+   * @param nrows
+   *        Number of rows
+   * @param ncols
+   *        Number of columns
+   * @return submatrix
+   * 
+   */
+
+  public Matrix getSubmatrix(int i0, int j0, int nrows, int ncols) {
+    Matrix x = new Matrix(null, nrows, ncols);
+    double[][] xa = x.a;
+    for (int i = nrows; --i >= 0;)
+      for (int j = ncols; --j >= 0;)
+        xa[i][j] = a[i0 + i][j0 + j];
+    return x;
+  }
+
+  /**
+   * Get a submatrix for a give number of columns and selected row set.
+   * 
+   * @param r
+   *        Array of row indices.
+   * @param n
+   *        number of rows 
+   * @return submatrix
+   */
+
+  public Matrix getMatrixSelected(int[] r, int n) {
+    Matrix x = new Matrix(null, r.length, n);
+    double[][] xa = x.a;
+    for (int i = r.length; --i >= 0;) {
+      double[] b = a[r[i]];
+      for (int j = n; --j >= 0;)
+        xa[i][j] = b[j];
+    }
+    return x;
+  }
+
+  /**
+   * Matrix transpose.
+   * 
+   * @return A'
+   */
+
+  public Matrix transpose() {
+    Matrix x = new Matrix(null, n, m);
+    double[][] c = x.a;
+    for (int i = m; --i >= 0;)
+      for (int j = n; --j >= 0;)
+        c[j][i] = a[i][j];
+    return x;
+  }
+
+  /**
+   * add two matrices
+   * @param b
+   * @return new Matrix this + b
+   */
+  public Matrix add(Matrix b) {
+    return scaleAdd(b, 1);
+  }
+
+  /**
+   * subtract two matrices
+   * @param b
+   * @return new Matrix this - b
+   */
+  public Matrix sub(Matrix b) {
+    return scaleAdd(b, -1);
+  }
+  
+  /**
+   * X = A + B*scale
+   * @param b 
+   * @param scale 
+   * @return X
+   * 
+   */
+  public Matrix scaleAdd(Matrix b, double scale) {
+    Matrix x = new Matrix(null, m, n);
+    double[][] xa = x.a;
+    double[][] ba = b.a;
+    for (int i = m; --i >= 0;)
+      for (int j = n; --j >= 0;)
+        xa[i][j] = ba[i][j] * scale + a[i][j];
+    return x;
+  }
+
+  /**
+   * Linear algebraic matrix multiplication, A * B
+   * 
+   * @param b
+   *        another matrix
+   * @return Matrix product, A * B or null for wrong dimension
+   */
+
+  public Matrix mul(Matrix b) {
+    if (b.m != n)
+      return null;
+    Matrix x = new Matrix(null, m, b.n);
+    double[][] xa = x.a;
+    double[][] ba = b.a;
+    for (int j = b.n; --j >= 0;)
+      for (int i = m; --i >= 0;) {
+        double[] arowi = a[i];
+        double s = 0;
+        for (int k = n; --k >= 0;)
+          s += arowi[k] * ba[k][j];
+        xa[i][j] = s;
+      }
+    return x;
+  }
+
+  /**
+   * Matrix inverse or pseudoinverse
+   * 
+   * @return inverse (m == n) or pseudoinverse (m != n)
+   */
+
+  public Matrix inverse() {
+    return new LUDecomp(m, n).solve(identity(m, m), n);
+  }
+
+  /**
+   * Matrix trace.
+   * 
+   * @return sum of the diagonal elements.
+   */
+
+  public double trace() {
+    double t = 0;
+    for (int i = Math.min(m, n); --i >= 0;)
+      t += a[i][i];
+    return t;
+  }
+
+  /**
+   * Generate identity matrix
+   * 
+   * @param m
+   *        Number of rows.
+   * @param n
+   *        Number of columns.
+   * @return An m-by-n matrix with ones on the diagonal and zeros elsewhere.
+   */
+
+  public static Matrix identity(int m, int n) {
+    Matrix x = new Matrix(null, m, n);
+    double[][] xa = x.a;
+    for (int i = Math.min(m, n); --i >= 0;)
+      xa[i][i] = 1;
+    return x;
+  }
+
+  /**
+   * similarly to M3/M4 standard rotation/translation matrix
+   * we set a rotationTranslation matrix to be:
+   * 
+   * [   nxn rot    nx1 trans
+   * 
+   *     1xn  0     1x1 1      ]
+   * 
+   * 
+   * @return rotation matrix
+   */
+  public Matrix getRotation() {
+    return getSubmatrix(0, 0, m - 1, n - 1);
+  }
+
+  public Matrix getTranslation() {
+    return getSubmatrix(0, n - 1, m - 1, 1);
+  }
+
+  public static Matrix newT(T3 r, boolean asColumn) {
+    return (asColumn ? new Matrix(new double[][] { new double[] { r.x },
+        new double[] { r.y }, new double[] { r.z } }, 3, 1) : new Matrix(
+        new double[][] { new double[] { r.x, r.y, r.z } }, 1, 3));
+  }
+
+  @Override
+  public String toString() {
+    String s = "[\n";
+    for (int i = 0; i < m; i++) {
+      s += "  [";
+      for (int j = 0; j < n; j++)
+        s += " " + a[i][j];
+      s += "]\n";
+    }
+    s += "]";
+    return s;
+  }
+
+  /**
+   * 
+   * Edited down by Bob Hanson for minimum needed by Jmol -- just constructor
+   * and solve
+   * 
+   * LU Decomposition.
+   * <P>
+   * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit
+   * lower triangular matrix L, an n-by-n upper triangular matrix U, and a
+   * permutation vector piv of length m so that A(piv,:) = L*U. If m < n, then L
+   * is m-by-m and U is m-by-n.
+   * <P>
+   * The LU decompostion with pivoting always exists, even if the matrix is
+   * singular, so the constructor will never fail. The primary use of the LU
+   * decomposition is in the solution of square systems of simultaneous linear
+   * equations. This will fail if isNonsingular() returns false.
+   */
+
+  private class LUDecomp {
+
+    /* ------------------------
+       Class variables
+     * ------------------------ */
+
+    /**
+     * Array for internal storage of decomposition.
+     * 
+     */
+    private double[][] LU;
+
+    /**
+     * Internal storage of pivot vector.
+     * 
+     */
+    private int[] piv;
+
+    private int pivsign;
+
+    /* ------------------------
+       Constructor
+     * ------------------------ */
+
+    /**
+     * LU Decomposition Structure to access L, U and piv.
+     * @param m 
+     * @param n 
+     * 
+     */
+
+    protected LUDecomp(int m, int n) {
+      
+      // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+
+      LU = getArrayCopy();
+      piv = new int[m];
+      for (int i = m; --i >= 0;)
+        piv[i] = i;
+      pivsign = 1;
+      double[] LUrowi;
+      double[] LUcolj = new double[m];
+
+      // Outer loop.
+
+      for (int j = 0; j < n; j++) {
+
+        // Make a copy of the j-th column to localize references.
+
+        for (int i = m; --i >= 0;)
+          LUcolj[i] = LU[i][j];
+
+        // Apply previous transformations.
+
+        for (int i = m; --i >= 0;) {
+          LUrowi = LU[i];
+
+          // Most of the time is spent in the following dot product.
+
+          int kmax = Math.min(i, j);
+          double s = 0.0;
+          for (int k = kmax; --k >= 0;)
+            s += LUrowi[k] * LUcolj[k];
+
+          LUrowi[j] = LUcolj[i] -= s;
+        }
+
+        // Find pivot and exchange if necessary.
+
+        int p = j;
+        for (int i = m; --i > j;)
+          if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p]))
+            p = i;
+        if (p != j) {
+          for (int k = n; --k >= 0;) {
+            double t = LU[p][k];
+            LU[p][k] = LU[j][k];
+            LU[j][k] = t;
+          }
+          int k = piv[p];
+          piv[p] = piv[j];
+          piv[j] = k;
+          pivsign = -pivsign;
+        }
+
+        // Compute multipliers.
+
+        if (j < m & LU[j][j] != 0.0)
+          for (int i = m; --i > j;)
+            LU[i][j] /= LU[j][j];
+      }
+    }
+
+    /* ------------------------
+       default Methods
+     * ------------------------ */
+
+    /**
+     * Solve A*X = B
+     * 
+     * @param b
+     *        A Matrix with as many rows as A and any number of columns.
+     * @param n 
+     * @return X so that L*U*X = B(piv,:) or null for wrong size or singular matrix
+     */
+
+    protected Matrix solve(Matrix b, int n) {
+      for (int j = 0; j < n; j++)
+        if (LU[j][j] == 0)
+          return null; // matrix is singular
+
+      // Copy right hand side with pivoting
+      int nx = b.n;
+      Matrix x = b.getMatrixSelected(piv, nx);
+      double[][] a = x.a;
+
+      // Solve L*Y = B(piv,:)
+      for (int k = 0; k < n; k++)
+        for (int i = k + 1; i < n; i++)
+          for (int j = 0; j < nx; j++)
+            a[i][j] -= a[k][j] * LU[i][k];
+
+      // Solve U*X = Y;
+      for (int k = n; --k >= 0;) {
+        for (int j = nx; --j >= 0;)
+          a[k][j] /= LU[k][k];
+        for (int i = k; --i >= 0;)
+          for (int j = nx; --j >= 0;)
+            a[i][j] -= a[k][j] * LU[i][k];
+      }
+      return x;
+    }
+  }
+
+}