content: package aima.core.util.math;
  
  import java.io.BufferedReader;
  import java.io.PrintWriter;
  import java.io.StreamTokenizer;
  import java.text.DecimalFormat;
  import java.text.DecimalFormatSymbols;
  import java.text.NumberFormat;
  import java.util.List;
  import java.util.Locale;
  
  /**
   * 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
   */
  public class Matrix implements Cloneable, java.io.Serializable {
  	private static final long serialVersionUID = 1L;
  
  	/*
  	 * ------------------------ Class variables ------------------------
  	 */
  
  	/**
  	 * Array for internal storage of elements.
  	 * 
  	 * @serial internal array storage.
  	 */
  	private final double[][] A;
  
  	/**
  	 * Row and column dimensions.
  	 * 
  	 * @serial row dimension.
  	 * @serial column dimension.
  	 */
  	private final int m, n;
  
  	/*
  	 * ------------------------ Constructors ------------------------
  	 */
  
  	/** Construct a diagonal Matrix from the given List of doubles */
  	public static Matrix createDiagonalMatrix(List<Double> values) {
  		Matrix m = new Matrix(values.size(), values.size(), 0);
  		for (int i = 0; i < values.size(); i++) {
  			m.set(i, i, values.get(i));
  		}
  		return m;
  	}
  
  	/**
  	 * 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.
  	 */
  
  	@Override
  	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 j0
  	 *            Initial column index
  	 * @param j1
  	 *            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;
  	}
  
  	/**
  	 * 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) {
  		// assumed m == n
  		return new LUDecomposition(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.
  	}
  
  	/**
  	 * 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(0, 255);
  		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);
  	}
  
  	@Override
  	public String toString() {
  		StringBuffer buf = new StringBuffer();
  		for (int i = 0; i < getRowDimension(); i++) {
  
  			for (int j = 0; j < getColumnDimension(); j++) {
  				buf.append(get(i, j));
  				buf.append(" ");
  			}
  			buf.append("\n");
  		}
  
  		return buf.toString();
  	}
  
  	/*
  	 * ------------------------ Private Methods ------------------------
  	 */
  
  	/** 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.");
  		}
  	}
  }
  