home | career | drupal | java | mac | mysql | perl | scala | uml | unix

# Java example source code file (LUDecomposition.java)

This example Java source code file (LUDecomposition.java) is included in the alvinalexander.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

## Java - Java tags/keywords

array2drowrealmatrix, arrayrealvector, decompositionsolver, default_too_small, ludecomposition, realmatrix, singularmatrixexception, solver

## The LUDecomposition.java Java example source code

```/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements.  See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License.  You may obtain a copy of the License at
*
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

package org.apache.commons.math3.linear;

import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.util.FastMath;

/**
* Calculates the LUP-decomposition of a square matrix.
* <p>The LUP-decomposition of a matrix A consists of three matrices L, U and
* P that satisfy: P×A = L×U. L is lower triangular (with unit
* diagonal terms), U is upper triangular and P is a permutation matrix. All
* matrices are m×m.</p>
* <p>As shown by the presence of the P matrix, this decomposition is
* implemented using partial pivoting.</p>
* <p>This class is based on the class with similar name from the
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA library.
* <ul>
*   <li>a {@link #getP() getP} method has been added,
*   <li>the {@code det} method has been renamed as {@link #getDeterminant()
*   getDeterminant},</li>
*   <li>the {@code getDoublePivot} method has been removed (but the int based
*   {@link #getPivot() getPivot} method has been kept),</li>
*   <li>the {@code solve} and {@code isNonSingular} methods have been replaced
*   by a {@link #getSolver() getSolver} method and the equivalent methods
*   provided by the returned {@link DecompositionSolver}.</li>
* </ul>
*
* @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld
* @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia
* @since 2.0 (changed to concrete class in 3.0)
*/
public class LUDecomposition {
/** Default bound to determine effective singularity in LU decomposition. */
private static final double DEFAULT_TOO_SMALL = 1e-11;
/** Entries of LU decomposition. */
private final double[][] lu;
/** Pivot permutation associated with LU decomposition. */
private final int[] pivot;
/** Parity of the permutation associated with the LU decomposition. */
private boolean even;
/** Singularity indicator. */
private boolean singular;
/** Cached value of L. */
private RealMatrix cachedL;
/** Cached value of U. */
private RealMatrix cachedU;
/** Cached value of P. */
private RealMatrix cachedP;

/**
* Calculates the LU-decomposition of the given matrix.
* This constructor uses 1e-11 as default value for the singularity
* threshold.
*
* @param matrix Matrix to decompose.
* @throws NonSquareMatrixException if matrix is not square.
*/
public LUDecomposition(RealMatrix matrix) {
this(matrix, DEFAULT_TOO_SMALL);
}

/**
* Calculates the LU-decomposition of the given matrix.
* @param matrix The matrix to decompose.
* @param singularityThreshold threshold (based on partial row norm)
* under which a matrix is considered singular
* @throws NonSquareMatrixException if matrix is not square
*/
public LUDecomposition(RealMatrix matrix, double singularityThreshold) {
if (!matrix.isSquare()) {
throw new NonSquareMatrixException(matrix.getRowDimension(),
matrix.getColumnDimension());
}

final int m = matrix.getColumnDimension();
lu = matrix.getData();
pivot = new int[m];
cachedL = null;
cachedU = null;
cachedP = null;

// Initialize permutation array and parity
for (int row = 0; row < m; row++) {
pivot[row] = row;
}
even     = true;
singular = false;

// Loop over columns
for (int col = 0; col < m; col++) {

// upper
for (int row = 0; row < col; row++) {
final double[] luRow = lu[row];
double sum = luRow[col];
for (int i = 0; i < row; i++) {
sum -= luRow[i] * lu[i][col];
}
luRow[col] = sum;
}

// lower
int max = col; // permutation row
double largest = Double.NEGATIVE_INFINITY;
for (int row = col; row < m; row++) {
final double[] luRow = lu[row];
double sum = luRow[col];
for (int i = 0; i < col; i++) {
sum -= luRow[i] * lu[i][col];
}
luRow[col] = sum;

// maintain best permutation choice
if (FastMath.abs(sum) > largest) {
largest = FastMath.abs(sum);
max = row;
}
}

// Singularity check
if (FastMath.abs(lu[max][col]) < singularityThreshold) {
singular = true;
return;
}

// Pivot if necessary
if (max != col) {
double tmp = 0;
final double[] luMax = lu[max];
final double[] luCol = lu[col];
for (int i = 0; i < m; i++) {
tmp = luMax[i];
luMax[i] = luCol[i];
luCol[i] = tmp;
}
int temp = pivot[max];
pivot[max] = pivot[col];
pivot[col] = temp;
even = !even;
}

// Divide the lower elements by the "winning" diagonal elt.
final double luDiag = lu[col][col];
for (int row = col + 1; row < m; row++) {
lu[row][col] /= luDiag;
}
}
}

/**
* Returns the matrix L of the decomposition.
* <p>L is a lower-triangular matrix
* @return the L matrix (or null if decomposed matrix is singular)
*/
public RealMatrix getL() {
if ((cachedL == null) && !singular) {
final int m = pivot.length;
cachedL = MatrixUtils.createRealMatrix(m, m);
for (int i = 0; i < m; ++i) {
final double[] luI = lu[i];
for (int j = 0; j < i; ++j) {
cachedL.setEntry(i, j, luI[j]);
}
cachedL.setEntry(i, i, 1.0);
}
}
return cachedL;
}

/**
* Returns the matrix U of the decomposition.
* <p>U is an upper-triangular matrix
* @return the U matrix (or null if decomposed matrix is singular)
*/
public RealMatrix getU() {
if ((cachedU == null) && !singular) {
final int m = pivot.length;
cachedU = MatrixUtils.createRealMatrix(m, m);
for (int i = 0; i < m; ++i) {
final double[] luI = lu[i];
for (int j = i; j < m; ++j) {
cachedU.setEntry(i, j, luI[j]);
}
}
}
return cachedU;
}

/**
* Returns the P rows permutation matrix.
* <p>P is a sparse matrix with exactly one element set to 1.0 in
* each row and each column, all other elements being set to 0.0.</p>
* <p>The positions of the 1 elements are given by the {@link #getPivot()
* pivot permutation vector}.</p>
* @return the P rows permutation matrix (or null if decomposed matrix is singular)
* @see #getPivot()
*/
public RealMatrix getP() {
if ((cachedP == null) && !singular) {
final int m = pivot.length;
cachedP = MatrixUtils.createRealMatrix(m, m);
for (int i = 0; i < m; ++i) {
cachedP.setEntry(i, pivot[i], 1.0);
}
}
return cachedP;
}

/**
* Returns the pivot permutation vector.
* @return the pivot permutation vector
* @see #getP()
*/
public int[] getPivot() {
return pivot.clone();
}

/**
* Return the determinant of the matrix
* @return determinant of the matrix
*/
public double getDeterminant() {
if (singular) {
return 0;
} else {
final int m = pivot.length;
double determinant = even ? 1 : -1;
for (int i = 0; i < m; i++) {
determinant *= lu[i][i];
}
return determinant;
}
}

/**
* Get a solver for finding the A × X = B solution in exact linear
* sense.
* @return a solver
*/
public DecompositionSolver getSolver() {
return new Solver(lu, pivot, singular);
}

/** Specialized solver. */
private static class Solver implements DecompositionSolver {

/** Entries of LU decomposition. */
private final double[][] lu;

/** Pivot permutation associated with LU decomposition. */
private final int[] pivot;

/** Singularity indicator. */
private final boolean singular;

/**
* Build a solver from decomposed matrix.
* @param lu entries of LU decomposition
* @param pivot pivot permutation associated with LU decomposition
* @param singular singularity indicator
*/
private Solver(final double[][] lu, final int[] pivot, final boolean singular) {
this.lu       = lu;
this.pivot    = pivot;
this.singular = singular;
}

/** {@inheritDoc} */
public boolean isNonSingular() {
return !singular;
}

/** {@inheritDoc} */
public RealVector solve(RealVector b) {
final int m = pivot.length;
if (b.getDimension() != m) {
throw new DimensionMismatchException(b.getDimension(), m);
}
if (singular) {
throw new SingularMatrixException();
}

final double[] bp = new double[m];

// Apply permutations to b
for (int row = 0; row < m; row++) {
bp[row] = b.getEntry(pivot[row]);
}

// Solve LY = b
for (int col = 0; col < m; col++) {
final double bpCol = bp[col];
for (int i = col + 1; i < m; i++) {
bp[i] -= bpCol * lu[i][col];
}
}

// Solve UX = Y
for (int col = m - 1; col >= 0; col--) {
bp[col] /= lu[col][col];
final double bpCol = bp[col];
for (int i = 0; i < col; i++) {
bp[i] -= bpCol * lu[i][col];
}
}

return new ArrayRealVector(bp, false);
}

/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b) {

final int m = pivot.length;
if (b.getRowDimension() != m) {
throw new DimensionMismatchException(b.getRowDimension(), m);
}
if (singular) {
throw new SingularMatrixException();
}

final int nColB = b.getColumnDimension();

// Apply permutations to b
final double[][] bp = new double[m][nColB];
for (int row = 0; row < m; row++) {
final double[] bpRow = bp[row];
final int pRow = pivot[row];
for (int col = 0; col < nColB; col++) {
bpRow[col] = b.getEntry(pRow, col);
}
}

// Solve LY = b
for (int col = 0; col < m; col++) {
final double[] bpCol = bp[col];
for (int i = col + 1; i < m; i++) {
final double[] bpI = bp[i];
final double luICol = lu[i][col];
for (int j = 0; j < nColB; j++) {
bpI[j] -= bpCol[j] * luICol;
}
}
}

// Solve UX = Y
for (int col = m - 1; col >= 0; col--) {
final double[] bpCol = bp[col];
final double luDiag = lu[col][col];
for (int j = 0; j < nColB; j++) {
bpCol[j] /= luDiag;
}
for (int i = 0; i < col; i++) {
final double[] bpI = bp[i];
final double luICol = lu[i][col];
for (int j = 0; j < nColB; j++) {
bpI[j] -= bpCol[j] * luICol;
}
}
}

return new Array2DRowRealMatrix(bp, false);
}

/**
* Get the inverse of the decomposed matrix.
*
* @return the inverse matrix.
* @throws SingularMatrixException if the decomposed matrix is singular.
*/
public RealMatrix getInverse() {
return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
}
}
}
```

## Other Java examples (source code examples)

Here is a short list of links related to this Java LUDecomposition.java source code file:

Copyright 1998-2019 Alvin Alexander, alvinalexander.com