LU Decomposition : Algorithms « Collections Data Structure « Java

Home
Java
1.2D Graphics GUI
2.3D
3.Advanced Graphics
4.Ant
5.Apache Common
6.Chart
7.Class
8.Collections Data Structure
9.Data Type
10.Database SQL JDBC
11.Design Pattern
12.Development Class
13.EJB3
14.Email
15.Event
16.File Input Output
17.Game
18.Generics
19.GWT
20.Hibernate
21.I18N
22.J2EE
23.J2ME
24.JavaFX
25.JDK 6
26.JDK 7
27.JNDI LDAP
28.JPA
29.JSP
30.JSTL
31.Language Basics
32.Network Protocol
33.PDF RTF
34.Reflection
35.Regular Expressions
36.Scripting
37.Security
38.Servlets
39.Spring
40.Swing Components
41.Swing JFC
42.SWT JFace Eclipse
43.Threads
44.Tiny Application
45.Velocity
46.Web Services SOA
47.XML
Java » Collections Data Structure » Algorithms 




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][j1.0;
        else {
          L[i][j0.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][j0.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 != ? 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 inputthrows 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((Doublev.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
java2s.com  | Contact Us | Privacy Policy
Copyright 2009 - 12 Demo Source and Support. All rights reserved.
All other trademarks are property of their respective owners.