Project: BoneJ
package org.doube.jama;
 
import ij.IJ; 
 
import java.text.NumberFormat; 
import java.text.DecimalFormat; 
import java.text.DecimalFormatSymbols; 
import java.util.Locale; 
import java.io.PrintWriter; 
import java.io.BufferedReader; 
import java.io.StreamTokenizer; 
 
/**
 * Jama = Java Matrix class. 
 * <P> 
 * The Java Matrix Class provides the fundamental operations of numerical linear 
 * algebra. Various constructors create Matrices from two dimensional arrays of 
 * double precision floating point numbers. Various "gets" and "sets" provide 
 * access to submatrices and matrix elements. Several methods implement basic 
 * matrix arithmetic, including matrix addition and multiplication, matrix 
 * norms, and element-by-element array operations. Methods for reading and 
 * printing matrices are also included. All the operations in this version of 
 * the Matrix Class involve real matrices. Complex matrices may be handled in a 
 * future version. 
 * <P> 
 * Five fundamental matrix decompositions, which consist of pairs or triples of 
 * matrices, permutation vectors, and the like, produce results in five 
 * decomposition classes. These decompositions are accessed by the Matrix class 
 * to compute solutions of simultaneous linear equations, determinants, inverses 
 * and other matrix functions. The five decompositions are: 
 * <P> 
 * <UL> 
 * <LI>Cholesky Decomposition of symmetric, positive definite matrices. 
 * <LI>LU Decomposition of rectangular matrices. 
 * <LI>QR Decomposition of rectangular matrices. 
 * <LI>Singular Value Decomposition of rectangular matrices. 
 * <LI>Eigenvalue Decomposition of both symmetric and nonsymmetric square 
 * matrices. 
 * </UL> 
 * <DL> 
 * <DT><B>Example of use:</B></DT> 
 * <P> 
 * <DD>Solve a linear system A x = b and compute the residual norm, ||b - A x||. 
 * <P> 
 *  
 * <PRE> 
 * double[][] vals = { { 1., 2., 3 }, { 4., 5., 6. }, { 7., 8., 10. } }; 
 * Matrix A = new Matrix(vals); 
 * Matrix b = Matrix.random(3, 1); 
 * Matrix x = A.solve(b); 
 * Matrix r = A.times(x).minus(b); 
 * double rnorm = r.normInf(); 
 * </PRE> 
 *  
 * </DD> 
 * </DL> 
 *  
 * @author The MathWorks, Inc. and the National Institute of Standards and 
 *         Technology. 
 * @version 5 August 1998 
 */
 
 
@SuppressWarnings("serial"
public class Matrix implements Cloneable, java.io.Serializable { 
 
 /*
  * ------------------------ Class variables ------------------------ 
  */
 
 
 /**
  * Array for internal storage of elements. 
  *  
  * @serial internal array storage. 
  */
 
 private double[][] A; 
 
 /**
  * Row and column dimensions. 
  *  
  * @serial row dimension. 
  * @serial column dimension. 
  */
 
 private int m, n; 
 
 /*
  * ------------------------ Constructors ------------------------ 
  */
 
 
 /**
  * Construct an m-by-n matrix of zeros. 
  *  
  * @param m 
  *            Number of rows. 
  * @param n 
  *            Number of colums. 
  */
 
 
 public Matrix(int m, int n) { 
  this.m = m; 
  this.n = n; 
  A = new double[m][n]; 
 
 
 /**
  * Construct an m-by-n constant matrix. 
  *  
  * @param m 
  *            Number of rows. 
  * @param n 
  *            Number of colums. 
  * @param s 
  *            Fill the matrix with this scalar value. 
  */
 
 
 public Matrix(int m, int n, double s) { 
  this.m = m; 
  this.n = n; 
  A = new double[m][n]; 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = s; 
   
  
 
 
 /**
  * Construct a matrix from a 2-D array. 
  *  
  * @param A 
  *            Two-dimensional array of doubles. 
  * @exception IllegalArgumentException 
  *                All rows must have the same length 
  * @see #constructWithCopy 
  */
 
 
 public Matrix(double[][] A) { 
  m = A.length; 
  n = A[0].length; 
  for (int i = 0; i < m; i++) { 
   if (A[i].length != n) { 
    throw new IllegalArgumentException( 
      "All rows must have the same length."); 
   
  
  this.A = A; 
 
 
 /**
  * Construct a matrix quickly without checking arguments. 
  *  
  * @param A 
  *            Two-dimensional array of doubles. 
  * @param m 
  *            Number of rows. 
  * @param n 
  *            Number of colums. 
  */
 
 
 public Matrix(double[][] A, int m, int n) { 
  this.A = A; 
  this.m = m; 
  this.n = n; 
 
 
 /**
  * Construct a matrix from a one-dimensional packed array 
  *  
  * @param vals 
  *            One-dimensional array of doubles, packed by columns (ala 
  *            Fortran). 
  * @param m 
  *            Number of rows. 
  * @exception IllegalArgumentException 
  *                Array length must be a multiple of m. 
  */
 
 
 public Matrix(double vals[], int m) { 
  this.m = m; 
  n = (m != 0 ? vals.length / m : 0); 
  if (m * n != vals.length) { 
   throw new IllegalArgumentException( 
     "Array length must be a multiple of m."); 
  
  A = new double[m][n]; 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = vals[i + j * m]; 
   
  
 
 
 /*
  * ------------------------ Public Methods ------------------------ 
  */
 
 
 /**
  * Construct a matrix from a copy of a 2-D array. 
  *  
  * @param A 
  *            Two-dimensional array of doubles. 
  * @exception IllegalArgumentException 
  *                All rows must have the same length 
  */
 
 
 public static Matrix constructWithCopy(double[][] A) { 
  int m = A.length; 
  int n = A[0].length; 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   if (A[i].length != n) { 
    throw new IllegalArgumentException( 
      "All rows must have the same length."); 
   
   for (int j = 0; j < n; j++) { 
    C[i][j] = A[i][j]; 
   
  
  return X; 
 
 
 /**
  * Make a deep copy of a matrix 
  */
 
 
 public Matrix copy() { 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = A[i][j]; 
   
  
  return X; 
 
 
 /**
  * Clone the Matrix object. 
  */
 
 
 public Object clone() { 
  return this.copy(); 
 
 
 /**
  * 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[][] C = new double[m][n]; 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = A[i][j]; 
   
  
  return C; 
 
 
 /**
  * Make a one-dimensional column packed copy of the internal array. 
  *  
  * @return Matrix elements packed in a one-dimensional array by columns. 
  */
 
 
 public double[] getColumnPackedCopy() { 
  double[] vals = new double[m * n]; 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    vals[i + j * m] = A[i][j]; 
   
  
  return vals; 
 
 
 /**
  * Make a one-dimensional row packed copy of the internal array. 
  *  
  * @return Matrix elements packed in a one-dimensional array by rows. 
  */
 
 
 public double[] getRowPackedCopy() { 
  double[] vals = new double[m * n]; 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    vals[i * n + j] = A[i][j]; 
   
  
  return vals; 
 
 
 /**
  * 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; 
 
 
 /**
  * Get a single element. 
  *  
  * @param i 
  *            Row index. 
  * @param j 
  *            Column index. 
  * @return A(i,j) 
  * @exception ArrayIndexOutOfBoundsException 
  */
 
 
 public double get(int i, int j) { 
  return A[i][j]; 
 
 
 /**
  * Get a submatrix. 
  *  
  * @param i0 
  *            Initial row index 
  * @param i1 
  *            Final row index 
  * @param j0 
  *            Initial column index 
  * @param j1 
  *            Final column index 
  * @return A(i0:i1,j0:j1) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public Matrix getMatrix(int i0, int i1, int j0, int j1) { 
  Matrix X = new Matrix(i1 - i0 + 1, j1 - j0 + 1); 
  double[][] B = X.getArray(); 
  try { 
   for (int i = i0; i <= i1; i++) { 
    for (int j = j0; j <= j1; j++) { 
     B[i - i0][j - j0] = A[i][j]; 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
  return X; 
 
 
 /**
  * Get a submatrix. 
  *  
  * @param r 
  *            Array of row indices. 
  * @param c 
  *            Array of column indices. 
  * @return A(r(:),c(:)) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public Matrix getMatrix(int[] r, int[] c) { 
  Matrix X = new Matrix(r.length, c.length); 
  double[][] B = X.getArray(); 
  try { 
   for (int i = 0; i < r.length; i++) { 
    for (int j = 0; j < c.length; j++) { 
     B[i][j] = A[r[i]][c[j]]; 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
  return X; 
 
 
 /**
  * Get a submatrix. 
  *  
  * @param i0 
  *            Initial row index 
  * @param i1 
  *            Final row index 
  * @param c 
  *            Array of column indices. 
  * @return A(i0:i1,c(:)) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public Matrix getMatrix(int i0, int i1, int[] c) { 
  Matrix X = new Matrix(i1 - i0 + 1, c.length); 
  double[][] B = X.getArray(); 
  try { 
   for (int i = i0; i <= i1; i++) { 
    for (int j = 0; j < c.length; j++) { 
     B[i - i0][j] = A[i][c[j]]; 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
  return X; 
 
 
 /**
  * Get a submatrix. 
  *  
  * @param r 
  *            Array of row indices. 
  * @param i0 
  *            Initial column index 
  * @param i1 
  *            Final column index 
  * @return A(r(:),j0:j1) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public Matrix getMatrix(int[] r, int j0, int j1) { 
  Matrix X = new Matrix(r.length, j1 - j0 + 1); 
  double[][] B = X.getArray(); 
  try { 
   for (int i = 0; i < r.length; i++) { 
    for (int j = j0; j <= j1; j++) { 
     B[i][j - j0] = A[r[i]][j]; 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
  return X; 
 
 
 /**
  * Set a single element. 
  *  
  * @param i 
  *            Row index. 
  * @param j 
  *            Column index. 
  * @param s 
  *            A(i,j). 
  * @exception ArrayIndexOutOfBoundsException 
  */
 
 
 public void set(int i, int j, double s) { 
  A[i][j] = s; 
 
 
 /**
  * Set a submatrix. 
  *  
  * @param i0 
  *            Initial row index 
  * @param i1 
  *            Final row index 
  * @param j0 
  *            Initial column index 
  * @param j1 
  *            Final column index 
  * @param X 
  *            A(i0:i1,j0:j1) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public void setMatrix(int i0, int i1, int j0, int j1, Matrix X) { 
  try { 
   for (int i = i0; i <= i1; i++) { 
    for (int j = j0; j <= j1; j++) { 
     A[i][j] = X.get(i - i0, j - j0); 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
 
 
 /**
  * Set a submatrix. 
  *  
  * @param r 
  *            Array of row indices. 
  * @param c 
  *            Array of column indices. 
  * @param X 
  *            A(r(:),c(:)) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public void setMatrix(int[] r, int[] c, Matrix X) { 
  try { 
   for (int i = 0; i < r.length; i++) { 
    for (int j = 0; j < c.length; j++) { 
     A[r[i]][c[j]] = X.get(i, j); 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
 
 
 /**
  * Set a submatrix. 
  *  
  * @param r 
  *            Array of row indices. 
  * @param j0 
  *            Initial column index 
  * @param j1 
  *            Final column index 
  * @param X 
  *            A(r(:),j0:j1) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public void setMatrix(int[] r, int j0, int j1, Matrix X) { 
  try { 
   for (int i = 0; i < r.length; i++) { 
    for (int j = j0; j <= j1; j++) { 
     A[r[i]][j] = X.get(i, j - j0); 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
 
 
 /**
  * Set a submatrix. 
  *  
  * @param i0 
  *            Initial row index 
  * @param i1 
  *            Final row index 
  * @param c 
  *            Array of column indices. 
  * @param X 
  *            A(i0:i1,c(:)) 
  * @exception ArrayIndexOutOfBoundsException 
  *                Submatrix indices 
  */
 
 
 public void setMatrix(int i0, int i1, int[] c, Matrix X) { 
  try { 
   for (int i = i0; i <= i1; i++) { 
    for (int j = 0; j < c.length; j++) { 
     A[i][c[j]] = X.get(i - i0, j); 
    
   
  catch (ArrayIndexOutOfBoundsException e) { 
   throw new ArrayIndexOutOfBoundsException("Submatrix indices"); 
  
 
 
 /**
  * Matrix transpose. 
  *  
  * @return A' 
  */
 
 
 public Matrix transpose() { 
  Matrix X = new Matrix(n, m); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[j][i] = A[i][j]; 
   
  
  return X; 
 
 
 /**
  * One norm 
  *  
  * @return maximum column sum. 
  */
 
 
 public double norm1() { 
  double f = 0
  for (int j = 0; j < n; j++) { 
   double s = 0
   for (int i = 0; i < m; i++) { 
    s += Math.abs(A[i][j]); 
   
   f = Math.max(f, s); 
  
  return f; 
 
 
 /**
  * Two norm 
  *  
  * @return maximum singular value. 
  */
 
 
 public double norm2() { 
  return (new SingularValueDecomposition(this).norm2()); 
 
 
 /**
  * Infinity norm 
  *  
  * @return maximum row sum. 
  */
 
 
 public double normInf() { 
  double f = 0
  for (int i = 0; i < m; i++) { 
   double s = 0
   for (int j = 0; j < n; j++) { 
    s += Math.abs(A[i][j]); 
   
   f = Math.max(f, s); 
  
  return f; 
 
 
 /**
  * Frobenius norm 
  *  
  * @return sqrt of sum of squares of all elements. 
  */
 
 
 public double normF() { 
  double f = 0
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    f = Maths.hypot(f, A[i][j]); 
   
  
  return f; 
 
 
 /**
  * Unary minus 
  *  
  * @return -A 
  */
 
 
 public Matrix uminus() { 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = -A[i][j]; 
   
  
  return X; 
 
 
 /**
  * C = A + B 
  *  
  * @param B 
  *            another matrix 
  * @return A + B 
  */
 
 
 public Matrix plus(Matrix B) { 
  checkMatrixDimensions(B); 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = A[i][j] + B.A[i][j]; 
   
  
  return X; 
 
 
 /**
  * A = A + B 
  *  
  * @param B 
  *            another matrix 
  * @return A + B 
  */
 
 
 public Matrix plusEquals(Matrix B) { 
  checkMatrixDimensions(B); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = A[i][j] + B.A[i][j]; 
   
  
  return this
 
 
 /**
  * C = A - B 
  *  
  * @param B 
  *            another matrix 
  * @return A - B 
  */
 
 
 public Matrix minus(Matrix B) { 
  checkMatrixDimensions(B); 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = A[i][j] - B.A[i][j]; 
   
  
  return X; 
 
 
 /**
  * A = A - B 
  *  
  * @param B 
  *            another matrix 
  * @return A - B 
  */
 
 
 public Matrix minusEquals(Matrix B) { 
  checkMatrixDimensions(B); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = A[i][j] - B.A[i][j]; 
   
  
  return this
 
 
 /**
  * Element-by-element multiplication, C = A.*B 
  *  
  * @param B 
  *            another matrix 
  * @return A.*B 
  */
 
 
 public Matrix arrayTimes(Matrix B) { 
  checkMatrixDimensions(B); 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = A[i][j] * B.A[i][j]; 
   
  
  return X; 
 
 
 /**
  * Element-by-element multiplication in place, A = A.*B 
  *  
  * @param B 
  *            another matrix 
  * @return A.*B 
  */
 
 
 public Matrix arrayTimesEquals(Matrix B) { 
  checkMatrixDimensions(B); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = A[i][j] * B.A[i][j]; 
   
  
  return this
 
 
 /**
  * Element-by-element right division, C = A./B 
  *  
  * @param B 
  *            another matrix 
  * @return A./B 
  */
 
 
 public Matrix arrayRightDivide(Matrix B) { 
  checkMatrixDimensions(B); 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = A[i][j] / B.A[i][j]; 
   
  
  return X; 
 
 
 /**
  * Element-by-element right division in place, A = A./B 
  *  
  * @param B 
  *            another matrix 
  * @return A./B 
  */
 
 
 public Matrix arrayRightDivideEquals(Matrix B) { 
  checkMatrixDimensions(B); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = A[i][j] / B.A[i][j]; 
   
  
  return this
 
 
 /**
  * Element-by-element left division, C = A.\B 
  *  
  * @param B 
  *            another matrix 
  * @return A.\B 
  */
 
 
 public Matrix arrayLeftDivide(Matrix B) { 
  checkMatrixDimensions(B); 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = B.A[i][j] / A[i][j]; 
   
  
  return X; 
 
 
 /**
  * Element-by-element left division in place, A = A.\B 
  *  
  * @param B 
  *            another matrix 
  * @return A.\B 
  */
 
 
 public Matrix arrayLeftDivideEquals(Matrix B) { 
  checkMatrixDimensions(B); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = B.A[i][j] / A[i][j]; 
   
  
  return this
 
 
 /**
  * Multiply a matrix by a scalar, C = s*A 
  *  
  * @param s 
  *            scalar 
  * @return s*A 
  */
 
 
 public Matrix times(double s) { 
  Matrix X = new Matrix(m, n); 
  double[][] C = X.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    C[i][j] = s * A[i][j]; 
   
  
  return X; 
 
 
 /**
  * Multiply a matrix by a scalar in place, A = s*A 
  *  
  * @param s 
  *            scalar 
  * @return replace A by s*A 
  */
 
 
 public Matrix timesEquals(double s) { 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    A[i][j] = s * A[i][j]; 
   
  
  return this
 
 
 /**
  * Linear algebraic matrix multiplication, A * B 
  *  
  * @param B 
  *            another matrix 
  * @return Matrix product, A * B 
  * @exception IllegalArgumentException 
  *                Matrix inner dimensions must agree. 
  */
 
 
 public Matrix times(Matrix B) { 
  if (B.m != n) { 
   throw new IllegalArgumentException( 
     "Matrix inner dimensions must agree."); 
  
  Matrix X = new Matrix(m, B.n); 
  double[][] C = X.getArray(); 
  double[] Bcolj = new double[n]; 
  for (int j = 0; j < B.n; j++) { 
   for (int k = 0; k < n; k++) { 
    Bcolj[k] = B.A[k][j]; 
   
   for (int i = 0; i < m; i++) { 
    double[] Arowi = A[i]; 
    double s = 0
    for (int k = 0; k < n; k++) { 
     s += Arowi[k] * Bcolj[k]; 
    
    C[i][j] = s; 
   
  
  return X; 
 
 
 /**
  * Calculate the dot product of two column or row vectors 
  *  
  * @param B 
  * @return the dot product 
  */
 
 public double dotProduct(Matrix B) { 
 
  if (m != B.m || n != B.n) { 
   throw new IllegalArgumentException("Matrix dimensions must agree"); 
  
  if (m != 1 && n != 1) { 
   throw new IllegalArgumentException( 
     "Matrix is not a row or column vector"); 
  
 
  double dot = 0
 
  // column vector 
  if (n == 1) { 
   for (int i = 0; i < m; i++) { 
    dot += A[i][0] * B.A[i][0]; 
   
  else if (m == 1) { 
   // row vector 
   for (int i = 0; i < n; i++) { 
    dot += A[0][i] * B.A[0][i]; 
   
  
  return dot; 
 
 
 /**
  * Calculate the cross product of 3-dimensional vectors 
  *  
  * @param B 
  * @return 
  */
 
 public Matrix crossProduct(Matrix B) { 
  if (m != B.m || n != B.n) { 
   throw new IllegalArgumentException("Matrix dimensions must agree"); 
  
  Matrix X = new Matrix(m, n); 
  if (m == 1 && n == 3) { 
   X.A[0][0] = A[0][1] * B.A[0][2] - A[0][2] * B.A[0][1]; 
   X.A[0][1] = A[0][2] * B.A[0][0] - A[0][0] * B.A[0][2]; 
   X.A[0][2] = A[0][0] * B.A[0][1] - A[0][1] * B.A[0][0]; 
  else if (m == 3 && n == 1) { 
   X.A[0][0] = A[1][0] * B.A[2][0] - A[2][0] * B.A[1][0]; 
   X.A[1][0] = A[2][0] * B.A[0][0] - A[0][0] * B.A[2][0]; 
   X.A[2][0] = A[0][0] * B.A[1][0] - A[1][0] * B.A[0][0]; 
  else { 
   throw new IllegalArgumentException( 
     "Matrix must be a 3-element row or column vector"); 
  
  return X; 
 
 
 /**
  * LU Decomposition 
  *  
  * @return LUDecomposition 
  * @see LUDecomposition 
  */
 
 
 public LUDecomposition lu() { 
  return new LUDecomposition(this); 
 
 
 /**
  * QR Decomposition 
  *  
  * @return QRDecomposition 
  * @see QRDecomposition 
  */
 
 
 public QRDecomposition qr() { 
  return new QRDecomposition(this); 
 
 
 /**
  * Cholesky Decomposition 
  *  
  * @return CholeskyDecomposition 
  * @see CholeskyDecomposition 
  */
 
 
 public CholeskyDecomposition chol() { 
  return new CholeskyDecomposition(this); 
 
 
 /**
  * Singular Value Decomposition 
  *  
  * @return SingularValueDecomposition 
  * @see SingularValueDecomposition 
  */
 
 
 public SingularValueDecomposition svd() { 
  return new SingularValueDecomposition(this); 
 
 
 /**
  * Eigenvalue Decomposition 
  *  
  * @return EigenvalueDecomposition 
  * @see EigenvalueDecomposition 
  */
 
 
 public EigenvalueDecomposition eig() { 
  return new EigenvalueDecomposition(this); 
 
 
 /**
  * Solve A*X = B 
  *  
  * @param B 
  *            right hand side 
  * @return solution if A is square, least squares solution otherwise 
  */
 
 
 public Matrix solve(Matrix B) { 
  return (m == n ? (new LUDecomposition(this)).solve(B) 
    : (new QRDecomposition(this)).solve(B)); 
 
 
 /**
  * Solve X*A = B, which is also A'*X' = B' 
  *  
  * @param B 
  *            right hand side 
  * @return solution if A is square, least squares solution otherwise. 
  */
 
 
 public Matrix solveTranspose(Matrix B) { 
  return transpose().solve(B.transpose()); 
 
 
 /**
  * Matrix inverse or pseudoinverse 
  *  
  * @return inverse(A) if A is square, pseudoinverse otherwise. 
  */
 
 
 public Matrix inverse() { 
  return solve(identity(m, m)); 
 
 
 /**
  * Matrix determinant 
  *  
  * @return determinant 
  */
 
 
 public double det() { 
  return new LUDecomposition(this).det(); 
 
 
 /**
  * Matrix rank 
  *  
  * @return effective numerical rank, obtained from SVD. 
  */
 
 
 public int rank() { 
  return new SingularValueDecomposition(this).rank(); 
 
 
 /**
  * Matrix condition (2 norm) 
  *  
  * @return ratio of largest to smallest singular value. 
  */
 
 
 public double cond() { 
  return new SingularValueDecomposition(this).cond(); 
 
 
 /**
  * Matrix trace. 
  *  
  * @return sum of the diagonal elements. 
  */
 
 
 public double trace() { 
  double t = 0
  for (int i = 0; i < Math.min(m, n); i++) { 
   t += A[i][i]; 
  
  return t; 
 
 
 /**
  * Generate matrix with random elements 
  *  
  * @param m 
  *            Number of rows. 
  * @param n 
  *            Number of colums. 
  * @return An m-by-n matrix with uniformly distributed random elements. 
  */
 
 
 public static Matrix random(int m, int n) { 
  Matrix A = new Matrix(m, n); 
  double[][] X = A.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    X[i][j] = Math.random(); 
   
  
  return A; 
 
 
 /**
  * Generate identity matrix 
  *  
  * @param m 
  *            Number of rows. 
  * @param n 
  *            Number of colums. 
  * @return An m-by-n matrix with ones on the diagonal and zeros elsewhere. 
  */
 
 
 public static Matrix identity(int m, int n) { 
  Matrix A = new Matrix(m, n); 
  double[][] X = A.getArray(); 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    X[i][j] = (i == j ? 1.0 : 0.0); 
   
  
  return A; 
 
 
 /**
  * Print the matrix to stdout. Line the elements up in columns with a 
  * Fortran-like 'Fw.d' style format. 
  *  
  * @param w 
  *            Column width. 
  * @param d 
  *            Number of digits after the decimal. 
  */
 
 
 public void print(int w, int d) { 
  print(new PrintWriter(System.out, true), w, d); 
 
 
 /**
  * Print the matrix to the output stream. Line the elements up in columns 
  * with a Fortran-like 'Fw.d' style format. 
  *  
  * @param output 
  *            Output stream. 
  * @param w 
  *            Column width. 
  * @param d 
  *            Number of digits after the decimal. 
  */
 
 
 public void print(PrintWriter output, int w, int d) { 
  DecimalFormat format = new DecimalFormat(); 
  format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US)); 
  format.setMinimumIntegerDigits(1); 
  format.setMaximumFractionDigits(d); 
  format.setMinimumFractionDigits(d); 
  format.setGroupingUsed(false); 
  print(output, format, w + 2); 
 
 
 /**
  * Print the matrix to stdout. Line the elements up in columns. Use the 
  * format object, and right justify within columns of width characters. Note 
  * that is the matrix is to be read back in, you probably will want to use a 
  * NumberFormat that is set to US Locale. 
  *  
  * @param format 
  *            A Formatting object for individual elements. 
  * @param width 
  *            Field width for each column. 
  * @see java.text.DecimalFormat#setDecimalFormatSymbols 
  */
 
 
 public void print(NumberFormat format, int width) { 
  print(new PrintWriter(System.out, true), format, width); 
 
 
 // DecimalFormat is a little disappointing coming from Fortran or C's 
 // printf. 
 // Since it doesn't pad on the left, the elements will come out different 
 // widths. Consequently, we'll pass the desired column width in as an 
 // argument and do the extra padding ourselves. 
 
 /**
  * Print the matrix to the output stream. Line the elements up in columns. 
  * Use the format object, and right justify within columns of width 
  * characters. Note that is the matrix is to be read back in, you probably 
  * will want to use a NumberFormat that is set to US Locale. 
  *  
  * @param output 
  *            the output stream. 
  * @param format 
  *            A formatting object to format the matrix elements 
  * @param width 
  *            Column width. 
  * @see java.text.DecimalFormat#setDecimalFormatSymbols 
  */
 
 
 public void print(PrintWriter output, NumberFormat format, int width) { 
  output.println(); // start on new line. 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    String s = format.format(A[i][j]); // format the number 
    int padding = Math.max(1, width - s.length()); // At _least_ 1 
    // space 
    for (int k = 0; k < padding; k++) 
     output.print(' '); 
    output.print(s); 
   
   output.println(); 
  
  output.println(); // end with blank line. 
 
 
 /**
  * Print Matrix to ImageJ log window 
  */
 
 public void printToIJLog() { 
  printToIJLog(""); 
  return
 
 
 /**
  * Print the Matrix to the ImageJ log 
  *  
  * @param title 
  *            Title of the Matrix 
  */
 
 public void printToIJLog(String title) { 
  if (!title.isEmpty()) 
   IJ.log(title); 
  int nCols = getColumnDimension(); 
  int nRows = getRowDimension(); 
  double[][] eVal = getArrayCopy(); 
  for (int r = 0; r < nRows; r++) { 
   String row = "||"
   for (int c = 0; c < nCols; c++) { 
    row = row + IJ.d2s(eVal[r][c], 3) + "|"
   
   row = row + "|"
   IJ.log(row); 
  
  IJ.log(""); 
  return
 
 
 /**
  * Read a matrix from a stream. The format is the same the print method, so 
  * printed matrices can be read back in (provided they were printed using US 
  * Locale). Elements are separated by whitespace, all the elements for each 
  * row appear on a single line, the last row is followed by a blank line. 
  *  
  * @param input 
  *            the input stream. 
  */
 
 
 @SuppressWarnings("unchecked"
 public static Matrix read(BufferedReader input) throws java.io.IOException { 
  StreamTokenizer tokenizer = new StreamTokenizer(input); 
 
  // Although StreamTokenizer will parse numbers, it doesn't recognize 
  // scientific notation (E or D); however, Double.valueOf does. 
  // The strategy here is to disable StreamTokenizer's number parsing. 
  // We'll only get whitespace delimited words, EOL's and EOF's. 
  // These words should all be numbers, for Double.valueOf to parse. 
 
  tokenizer.resetSyntax(); 
  tokenizer.wordChars(0255); 
  tokenizer.whitespaceChars(0, ' '); 
  tokenizer.eolIsSignificant(true); 
  java.util.Vector v = new java.util.Vector(); 
 
  // Ignore initial empty lines 
  while (tokenizer.nextToken() == StreamTokenizer.TT_EOL) 
   
  if (tokenizer.ttype == StreamTokenizer.TT_EOF) 
   throw new java.io.IOException("Unexpected EOF on matrix read."); 
  do { 
   v.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st 
   // row. 
  while (tokenizer.nextToken() == StreamTokenizer.TT_WORD); 
 
  int n = v.size(); // Now we've got the number of columns! 
  double row[] = new double[n]; 
  for (int j = 0; j < n; j++) 
   // extract the elements of the 1st row. 
   row[j] = ((Double) v.elementAt(j)).doubleValue(); 
  v.removeAllElements(); 
  v.addElement(row); // Start storing rows instead of columns. 
  while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) { 
   // While non-empty lines 
   v.addElement(row = new double[n]); 
   int j = 0
   do { 
    if (j >= n) 
     throw new java.io.IOException("Row " + v.size() 
       " is too long."); 
    row[j++] = Double.valueOf(tokenizer.sval).doubleValue(); 
   while (tokenizer.nextToken() == StreamTokenizer.TT_WORD); 
   if (j < n) 
    throw new java.io.IOException("Row " + v.size() 
      " is too short."); 
  
  int m = v.size(); // Now we've got the number of rows. 
  double[][] A = new double[m][]; 
  v.copyInto(A); // copy the rows out of the vector 
  return new Matrix(A); 
 
 
 /**
  * Create an m * n Matrix filled with 1 
  *  
  * @param m 
  *            number of rows 
  * @param n 
  *            number of columns 
  * @return m * n Matrix filled with 1 
  */
 
 public static Matrix ones(final int m, final int n) { 
  double[][] ones = new double[m][n]; 
  for (int i = 0; i < m; i++) { 
   for (int j = 0; j < n; j++) { 
    ones[i][j] = 1
   
  
  return new Matrix(ones); 
 
 
 /**
  * Create an n * n square identity matrix with 1 on the diagonal and 0 
  * elsewhere 
  *  
  * @param n 
  *            square matrix dimension 
  * @return n * n identity matrix 
  */
 
 public static Matrix eye(final int n) { 
  return eye(n, n); 
 
 
 /**
  * Create an m * n identity matrix 
  *  
  * @param m 
  * @param n 
  * @return 
  */
 
 public static Matrix eye(int m, int n) { 
  double[][] eye = new double[m][n]; 
  final int min = Math.min(m, n); 
  for (int i = 0; i < min; i++) { 
   eye[i][i] = 1
  
  return new Matrix(eye); 
 
 
 /**
  * Get the diagonal of the matrix as a column vector 
  *  
  * @return Column vector containing diagonal 
  */
 
 public Matrix diag() { 
  final int min = Math.min(m, n); 
  double[][] diag = new double[min][1]; 
  for (int i = 0; i < min; i++) { 
   diag[i][0] = get(i, i); 
  
  return new Matrix(diag); 
 
 
 /**
  * Check if a 3 x 3 Matrix is right handed. If the matrix is a rotation 
  * matrix, then right-handedness implies rotation only, while 
  * left-handedness implies a reflection will be performed. 
  *  
  * @return true if the Matrix is right handed, false if it is left handed 
  */
 
 public boolean isRightHanded() { 
  if (m != 3 || n != 3) { 
   throw new IllegalArgumentException(); 
  
 
  final double x0 = get(00); 
  final double x1 = get(10); 
  final double x2 = get(20); 
  final double y0 = get(01); 
  final double y1 = get(11); 
  final double y2 = get(21); 
  final double z0 = get(02); 
  final double z1 = get(12); 
  final double z2 = get(22); 
 
  final double c0 = x1 * y2 - x2 * y1; 
  final double c1 = x2 * y0 - x0 * y2; 
  final double c2 = x0 * y1 - x1 * y0; 
 
  final double dot = c0 * z0 + c1 * z1 + c2 * z2; 
 
  if (dot > 0
   return true
  else 
   return false
 
 
 /**
  * Check if a rotation matrix will flip the direction of the z component of 
  * the original 
  *  
  * @return true if the rotation matrix will cause z-flipping 
  */
 
 public boolean isZFlipped() { 
  final double x2 = get(20); 
  final double y2 = get(21); 
  final double z2 = get(22); 
 
  final double dot = x2 + y2 + z2; 
 
  if (dot < 0
   return true
  else 
   return false
 
 
 /** Check if size(A) == size(B) **/ 
 
 private void checkMatrixDimensions(Matrix B) { 
  if (B.m != m || B.n != n) { 
   throw new IllegalArgumentException("Matrix dimensions must agree."); 
  
 
}