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;
}
}
}