LU Decomposition : Algorithms « Collections Data Structure « Java





LU Decomposition

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

/**
 * 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.
 */
public class LUDecomposition implements java.io.Serializable {
  private static final long serialVersionUID = 1L;

  /*
   * ------------------------ Class variables ------------------------
   */

  /**
   * Array for internal storage of decomposition.
   * 
   * @serial internal array storage.
   */
  private final double[][] LU;

  /**
   * Row and column dimensions, and pivot sign.
   * 
   * @serial column dimension.
   * @serial row dimension.
   * @serial pivot sign.
   */
  private final int m, n;

  private int pivsign;

  /**
   * Internal storage of pivot vector.
   * 
   * @serial pivot vector.
   */
  private final int[] piv;

  /*
   * ------------------------ Constructor ------------------------
   */

  /**
   * LU Decomposition, a structure to access L, U and piv.
   * 
   * @param A
   *            Rectangular matrix
   */
  public LUDecomposition(Matrix A) {

    // Use a "left-looking", dot-product, Crout/Doolittle algorithm.

    LU = A.getArrayCopy();
    m = A.getRowDimension();
    n = A.getColumnDimension();
    piv = new int[m];
    for (int i = 0; i < m; i++) {
      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 = 0; i < m; i++) {
        LUcolj[i] = LU[i][j];
      }

      // Apply previous transformations.

      for (int i = 0; i < m; i++) {
        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 = 0; k < kmax; k++) {
          s += LUrowi[k] * LUcolj[k];
        }

        LUrowi[j] = LUcolj[i] -= s;
      }

      // Find pivot and exchange if necessary.

      int p = j;
      for (int i = j + 1; i < m; i++) {
        if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
          p = i;
        }
      }
      if (p != j) {
        for (int k = 0; k < n; k++) {
          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 = j + 1; i < m; i++) {
          LU[i][j] /= LU[j][j];
        }
      }
    }
  }

  /*
   * ------------------------ Temporary, experimental code.
   * ------------------------\
   * 
   * \ LU Decomposition, computed by Gaussian elimination. <P> This
   * constructor computes L and U with the "daxpy"-based elimination algorithm
   * used in LINPACK and MATLAB. In Java, we suspect the dot-product, Crout
   * algorithm will be faster. We have temporarily included this constructor
   * until timing experiments confirm this suspicion. <P> @param A Rectangular
   * matrix @param linpackflag Use Gaussian elimination. Actual value ignored.
   * 
   * @return Structure to access L, U and piv. \
   * 
   * public LUDecomposition (Matrix A, int linpackflag) { // Initialize. LU =
   * A.getArrayCopy(); m = A.getRowDimension(); n = A.getColumnDimension();
   * piv = new int[m]; for (int i = 0; i < m; i++) { piv[i] = i; } pivsign =
   * 1; // Main loop. for (int k = 0; k < n; k++) { // Find pivot. int p = k;
   * for (int i = k+1; i < m; i++) { if (Math.abs(LU[i][k]) >
   * Math.abs(LU[p][k])) { p = i; } } // Exchange if necessary. if (p != k) {
   * for (int j = 0; j < n; j++) { double t = LU[p][j]; LU[p][j] = LU[k][j];
   * LU[k][j] = t; } int t = piv[p]; piv[p] = piv[k]; piv[k] = t; pivsign =
   * -pivsign; } // Compute multipliers and eliminate k-th column. if
   * (LU[k][k] != 0.0) { for (int i = k+1; i < m; i++) { LU[i][k] /= LU[k][k];
   * for (int j = k+1; j < n; j++) { LU[i][j] -= LU[i][k]LU[k][j]; } } } } } \
   * ------------------------ End of temporary code. ------------------------
   */

  /*
   * ------------------------ Public Methods ------------------------
   */

  /**
   * Is the matrix nonsingular?
   * 
   * @return true if U, and hence A, is nonsingular.
   */
  public boolean isNonsingular() {
    for (int j = 0; j < n; j++) {
      if (LU[j][j] == 0)
        return false;
    }
    return true;
  }

  /**
   * Return lower triangular factor
   * 
   * @return L
   */
  public Matrix getL() {
    Matrix X = new Matrix(m, n);
    double[][] L = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if (i > j) {
          L[i][j] = LU[i][j];
        } else if (i == j) {
          L[i][j] = 1.0;
        } else {
          L[i][j] = 0.0;
        }
      }
    }
    return X;
  }

  /**
   * Return upper triangular factor
   * 
   * @return U
   */
  public Matrix getU() {
    Matrix X = new Matrix(n, n);
    double[][] U = X.getArray();
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        if (i <= j) {
          U[i][j] = LU[i][j];
        } else {
          U[i][j] = 0.0;
        }
      }
    }
    return X;
  }

  /**
   * Return pivot permutation vector
   * 
   * @return piv
   */
  public int[] getPivot() {
    int[] p = new int[m];
    for (int i = 0; i < m; i++) {
      p[i] = piv[i];
    }
    return p;
  }

  /**
   * Return pivot permutation vector as a one-dimensional double array
   * 
   * @return (double) piv
   */
  public double[] getDoublePivot() {
    double[] vals = new double[m];
    for (int i = 0; i < m; i++) {
      vals[i] = piv[i];
    }
    return vals;
  }

  /**
   * Determinant
   * 
   * @return det(A)
   * @exception IllegalArgumentException
   *                Matrix must be square
   */
  public double det() {
    if (m != n) {
      throw new IllegalArgumentException("Matrix must be square.");
    }
    double d = pivsign;
    for (int j = 0; j < n; j++) {
      d *= LU[j][j];
    }
    return d;
  }

  /**
   * Solve A*X = B
   * 
   * @param B
   *            A Matrix with as many rows as A and any number of columns.
   * @return X so that L*U*X = B(piv,:)
   * @exception IllegalArgumentException
   *                Matrix row dimensions must agree.
   * @exception RuntimeException
   *                Matrix is singular.
   */
  public Matrix solve(Matrix B) {
    if (B.getRowDimension() != m) {
      throw new IllegalArgumentException(
          "Matrix row dimensions must agree.");
    }
    if (!this.isNonsingular()) {
      throw new RuntimeException("Matrix is singular.");
    }

    // Copy right hand side with pivoting
    int nx = B.getColumnDimension();
    Matrix Xmat = B.getMatrix(piv, 0, nx - 1);
    double[][] X = Xmat.getArray();

    // 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++) {
          X[i][j] -= X[k][j] * LU[i][k];
        }
      }
    }
    // Solve U*X = Y;
    for (int k = n - 1; k >= 0; k--) {
      for (int j = 0; j < nx; j++) {
        X[k][j] /= LU[k][k];
      }
      for (int i = 0; i < k; i++) {
        for (int j = 0; j < nx; j++) {
          X[i][j] -= X[k][j] * LU[i][k];
        }
      }
    }
    return Xmat;
  }
}


/**
 * 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
 */
 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.
   */

  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.");
    }
  }
}

   
    
    
    
    
    
    
    
  










Related examples in the same category

1.AnagramsAnagrams
2.Hanoi puzzleHanoi puzzle
3.FibonacciFibonacci
4.Sieve Sieve
5.Find connections using a depth-first searchFind connections using a depth-first search
6.Find connections using hill climbing.
7.Find optimal solution using least-cost
8.Find the lost keysFind the lost keys
9.Compute the area of a triangle using Heron's FormulaCompute the area of a triangle using Heron's Formula
10.Compute prime numbers
11.Print a table of fahrenheit and celsius temperatures 1
12.Print a table of fahrenheit and celsius temperatures 2
13.Print a table of Fahrenheit and Celsius temperatures 3Print a table of Fahrenheit and Celsius temperatures 3
14.Soundex - the Soundex Algorithm, as described by KnuthSoundex - the Soundex Algorithm, as described by Knuth
15.A programmable Finite State Machine implementation.
16.An extendable Graph datastructure.
17.Utilities for flop (floating-point operation) counting.
18.Reverse Polish Notation
19.Permutator test
20.implements the LZF lossless data compression algorithm
21.Linear Interpolation
22.Utility class for generating the k-subsets of the numbers 0 to n
23.VersionVersion