Mark Bishop

Analytical Chemistry and

Signal Acquisition/Analysis



/**
 * Copyright 2009 Mark Bishop
 *  This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version. 
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details: http://www.gnu.org/licenses
 
    The author makes no warranty for the accuracy, completeness, safety, or usefulness
    of any information provided and does not represent that its use would not infringe
    privately owned right.
 */
 
import java.util.ArrayList;
import java.util.List;
 
public class CpuLU {
	/**
	 * 
	 * @param A
	 *            A[m,m]
	 * @return A Inverse such that A(A Inverse) = I[m]
	 */
	public static float[][] Inverse(float[][] A) {
		int n = A.length;
		float[][] AInv = new float[n][n];
		float[][] A1 = Utilities.CloneMatrix(A);
 
		List<Object> LUPQ = GaussUpperTriangularFullPivot(A1);
		A1 = (float[][]) LUPQ.get(0);
		int[] P = (int[]) LUPQ.get(1);
		int[] Q = (int[]) LUPQ.get(2);
 
		// Ax = b, solve for each column x[j] AInverse[j] and b = I[j] (one
		// column at a
		// time)
		float[][] identity = CpuMatrixOps.Eye(n);
		float[] bSub = new float[n];
		float[] xSub = new float[n];
		for (int i = 0; i < n; i++) {
			bSub = identity[i];
			xSub = CompactLU_SolveForX(A1, bSub, P, Q);
			for (int j = 0; j < n; j++) {
				AInv[j][i] = xSub[j];
			}
		}
		return AInv;
	}
 
	/**
	 * Full pivot upper triangular Gaussian elimination of leftmost square of A
	 * 
	 * @param A
	 *            (Will be overwritten!)
	 * @return A List<Object> LUPQ where:
	 * 
	 *         The first object is a full-pivot LU in compact encoded form
	 *         stored as a float[][]
	 * 
	 *         The second object is a row permutation vector stored as an int[]
	 * 
	 *         The third object is the column permutation vector stored as an
	 *         int[]
	 */
	public static List<Object> GaussUpperTriangularFullPivot(float[][] A) {
 
		int m = A.length;
		int n = A[0].length;
		int i = 0, j = 0;
 
		// initialize permutation vectors
		int[] P = new int[m];
		int[] Q = new int[n];
		for (i = 0; i < m; i++) {
			P[i] = i;
			Q[i] = i;
		}
 
		for (int k = 0; k < m; k++) {
			int exchangeRow = k;
			int exchangeColumn = k;
 
			// Find the element with the highest absolute value in the pivot
			// submatrix and record its row and column
			float max = Math.abs(A[k][k]);
			for (i = k; i < m; i++) {
				for (j = k; j < m; j++) {
					if (Math.abs(A[i][j]) > max) {
						max = Math.abs(A[i][j]);
						exchangeRow = i;
						exchangeColumn = j;
					}
				}
			}
 
			// exchange rows
			if (exchangeRow != k) {
				for (j = 0; j < n; j++) {
					// for A
					float exchangeStorage = A[k][j];
					A[k][j] = A[exchangeRow][j];
					A[exchangeRow][j] = exchangeStorage;
				}
				// for P
				int exchangeIndex = P[k];
				P[k] = P[exchangeRow];
				P[exchangeRow] = exchangeIndex;
			}
 
			// Exchange columns to bring element with highest absolute value to
			// the
			// principal diagonal
			if (exchangeColumn != k) {
				for (i = 0; i < m; i++) {
					// for A
					float exchangeStorage = A[i][k];
					A[i][k] = A[i][exchangeColumn];
					A[i][exchangeColumn] = exchangeStorage;
				}
				// for P
				int exchangeIndex = Q[k];
				Q[k] = Q[exchangeColumn];
				Q[exchangeColumn] = exchangeIndex;
			}
 
			if (A[k][k] == 0) {
				return null; // singular matrix
			}
 
			for (i = k + 1; i < m; i++) {
				float rowFactor = A[i][k] / A[k][k];
				for (j = k; j < n; j++) {
					A[i][j] = A[i][j] - rowFactor * A[k][j];
				}
				A[i][k] = rowFactor;
			}
		}
		List<Object> APQ = new ArrayList<Object>();
		APQ.add(A);
		APQ.add(P);
		APQ.add(Q);
		return APQ;
	}
 
	/**
	 * Given PAQ = LU and b, Solve for x (Ax = b)
	 * 
	 * @param LU
	 *            Compact form LU
	 * @param b
	 * @param P
	 *            row permutation (can be null if no row exchanges)
	 * @param Q
	 *            column permutation (can be null if no column exchanges)
	 * @return x
	 */
	public static float[] CompactLU_SolveForX(float[][] LU, float[] b, int[] P,
			int[] Q) {
		float[] x = b.clone();
		int n = LU.length;
		if (P != null) {
			float[] tempb = new float[n];
			for (int i = 0; i < n; i++) {
				tempb[i] = x[P[i]];
			}
			x = tempb; // x now holds Pb
		}
		float temp = 0;
 
		// forward calculate for Ly=Pb
		for (int i = 1; i < n; i++) {
			for (int j = 0; j <= i - 1; j++) {
				temp = temp - LU[i][j] * x[j];
			}
			x[i] = temp + x[i];
			temp = 0;
		}
		// x now holds y
		// back calculate for Uz=y
		x = CpuMatrixOps.BackCalculateX(LU, x);
		// x now holds z = x before permutation is applied
		if (Q == null) {
			return x;
		} else {
			float[] tempb = new float[n];
			for (int i = 0; i < n; i++) {
				tempb[Q[i]] = x[i];
			}
			x = tempb;
			return x;
		}
	}
 
	/**
	 * Parse a compact form LU into L and U
	 * 
	 * @param LU
	 *            An LU factorization in compact form
	 * @return result = L and U as an array of float[][] of length 2 where
	 *         result[0] = L and result[1] = U
	 */
	public static float[][][] ParseLU(float[][] LU) {
		int n = LU.length;
		float[][][] result = new float[2][][];
		float[][] L = CpuMatrixOps.Eye(n);
		float[][] U = new float[n][n];
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				if (i <= j) {
					U[i][j] = LU[i][j];
				} else {
					L[i][j] = LU[i][j];
				}
			}
		}
		result[0] = L;
		result[1] = U;
		return result;
	}
}
 
 
ID may be exposed