- * 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. - *
- * 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; - } - } - -} +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. + *
+ * 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. + *
+ * 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; + } + } + +}