Merged FieldLUDecomposition and FieldLUDecompositionImpl (see MATH-662).
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1174537 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
parent
81821bc466
commit
4709a37a4f
|
@ -17,49 +17,192 @@
|
||||||
|
|
||||||
package org.apache.commons.math.linear;
|
package org.apache.commons.math.linear;
|
||||||
|
|
||||||
|
import java.lang.reflect.Array;
|
||||||
|
|
||||||
|
import org.apache.commons.math.Field;
|
||||||
import org.apache.commons.math.FieldElement;
|
import org.apache.commons.math.FieldElement;
|
||||||
|
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* An interface to classes that implement an algorithm to calculate the
|
* Calculates the LUP-decomposition of a square matrix.
|
||||||
* LU-decomposition of a real matrix.
|
* <p>The LUP-decomposition of a matrix A consists of three matrices
|
||||||
* <p>The LU-decomposition of matrix A is a set of three matrices: P, L and U
|
* L, U and P that satisfy: PA = LU, L is lower triangular, and U is
|
||||||
* such that P×A = L×U. P is a rows permutation matrix that is used
|
* upper triangular and P is a permutation matrix. All matrices are
|
||||||
* to rearrange the rows of A before so that it can be decomposed. L is a lower
|
* m×m.</p>
|
||||||
* triangular matrix with unit diagonal terms and U is an upper triangular matrix.</p>
|
* <p>Since {@link FieldElement field elements} do not provide an ordering
|
||||||
* <p>This interface is based on the class with similar name from the
|
* operator, the permutation matrix is computed here only in order to avoid
|
||||||
|
* a zero pivot element, no attempt is done to get the largest pivot
|
||||||
|
* element.</p>
|
||||||
|
* <p>This class is based on the class with similar name from the
|
||||||
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
|
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
|
||||||
* <ul>
|
* <ul>
|
||||||
* <li>a {@link #getP() getP} method has been added,</li>
|
* <li>a {@link #getP() getP} method has been added,</li>
|
||||||
* <li>the <code>det</code> method has been renamed as {@link #getDeterminant()
|
* <li>the {@code det} method has been renamed as {@link #getDeterminant()
|
||||||
* getDeterminant},</li>
|
* getDeterminant},</li>
|
||||||
* <li>the <code>getDoublePivot</code> method has been removed (but the int based
|
* <li>the {@code getDoublePivot} method has been removed (but the int based
|
||||||
* {@link #getPivot() getPivot} method has been kept),</li>
|
* {@link #getPivot() getPivot} method has been kept),</li>
|
||||||
* <li>the <code>solve</code> and <code>isNonSingular</code> methods have been replaced
|
* <li>the {@code solve} and {@code isNonSingular} methods have been replaced
|
||||||
* by a {@link #getSolver() getSolver} method and the equivalent methods provided by
|
* by a {@link #getSolver() getSolver} method and the equivalent methods
|
||||||
* the returned {@link DecompositionSolver}.</li>
|
* provided by the returned {@link DecompositionSolver}.</li>
|
||||||
* </ul>
|
* </ul>
|
||||||
*
|
*
|
||||||
* @param <T> the type of the field elements
|
* @param <T> the type of the field elements
|
||||||
* @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
|
* @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
|
||||||
* @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
|
* @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
|
||||||
* @version $Id$
|
* @version $Id$
|
||||||
* @since 2.0
|
* @since 2.0 (changed to concrete class in 3.0)
|
||||||
*/
|
*/
|
||||||
public interface FieldLUDecomposition<T extends FieldElement<T>> {
|
public class FieldLUDecomposition<T extends FieldElement<T>> {
|
||||||
|
|
||||||
|
/** Field to which the elements belong. */
|
||||||
|
private final Field<T> field;
|
||||||
|
|
||||||
|
/** Entries of LU decomposition. */
|
||||||
|
private T[][] lu;
|
||||||
|
|
||||||
|
/** Pivot permutation associated with LU decomposition. */
|
||||||
|
private int[] pivot;
|
||||||
|
|
||||||
|
/** Parity of the permutation associated with the LU decomposition. */
|
||||||
|
private boolean even;
|
||||||
|
|
||||||
|
/** Singularity indicator. */
|
||||||
|
private boolean singular;
|
||||||
|
|
||||||
|
/** Cached value of L. */
|
||||||
|
private FieldMatrix<T> cachedL;
|
||||||
|
|
||||||
|
/** Cached value of U. */
|
||||||
|
private FieldMatrix<T> cachedU;
|
||||||
|
|
||||||
|
/** Cached value of P. */
|
||||||
|
private FieldMatrix<T> cachedP;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Calculates the LU-decomposition of the given matrix.
|
||||||
|
* @param matrix The matrix to decompose.
|
||||||
|
* @throws NonSquareMatrixException if matrix is not square
|
||||||
|
*/
|
||||||
|
public FieldLUDecomposition(FieldMatrix<T> matrix) {
|
||||||
|
if (!matrix.isSquare()) {
|
||||||
|
throw new NonSquareMatrixException(matrix.getRowDimension(),
|
||||||
|
matrix.getColumnDimension());
|
||||||
|
}
|
||||||
|
|
||||||
|
final int m = matrix.getColumnDimension();
|
||||||
|
field = matrix.getField();
|
||||||
|
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++) {
|
||||||
|
|
||||||
|
T sum = field.getZero();
|
||||||
|
|
||||||
|
// upper
|
||||||
|
for (int row = 0; row < col; row++) {
|
||||||
|
final T[] luRow = lu[row];
|
||||||
|
sum = luRow[col];
|
||||||
|
for (int i = 0; i < row; i++) {
|
||||||
|
sum = sum.subtract(luRow[i].multiply(lu[i][col]));
|
||||||
|
}
|
||||||
|
luRow[col] = sum;
|
||||||
|
}
|
||||||
|
|
||||||
|
// lower
|
||||||
|
int nonZero = col; // permutation row
|
||||||
|
for (int row = col; row < m; row++) {
|
||||||
|
final T[] luRow = lu[row];
|
||||||
|
sum = luRow[col];
|
||||||
|
for (int i = 0; i < col; i++) {
|
||||||
|
sum = sum.subtract(luRow[i].multiply(lu[i][col]));
|
||||||
|
}
|
||||||
|
luRow[col] = sum;
|
||||||
|
|
||||||
|
if (lu[nonZero][col].equals(field.getZero())) {
|
||||||
|
// try to select a better permutation choice
|
||||||
|
++nonZero;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Singularity check
|
||||||
|
if (nonZero >= m) {
|
||||||
|
singular = true;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Pivot if necessary
|
||||||
|
if (nonZero != col) {
|
||||||
|
T tmp = field.getZero();
|
||||||
|
for (int i = 0; i < m; i++) {
|
||||||
|
tmp = lu[nonZero][i];
|
||||||
|
lu[nonZero][i] = lu[col][i];
|
||||||
|
lu[col][i] = tmp;
|
||||||
|
}
|
||||||
|
int temp = pivot[nonZero];
|
||||||
|
pivot[nonZero] = pivot[col];
|
||||||
|
pivot[col] = temp;
|
||||||
|
even = !even;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Divide the lower elements by the "winning" diagonal elt.
|
||||||
|
final T luDiag = lu[col][col];
|
||||||
|
for (int row = col + 1; row < m; row++) {
|
||||||
|
final T[] luRow = lu[row];
|
||||||
|
luRow[col] = luRow[col].divide(luDiag);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns the matrix L of the decomposition.
|
* Returns the matrix L of the decomposition.
|
||||||
* <p>L is an lower-triangular matrix</p>
|
* <p>L is a lower-triangular matrix</p>
|
||||||
* @return the L matrix (or null if decomposed matrix is singular)
|
* @return the L matrix (or null if decomposed matrix is singular)
|
||||||
*/
|
*/
|
||||||
FieldMatrix<T> getL();
|
public FieldMatrix<T> getL() {
|
||||||
|
if ((cachedL == null) && !singular) {
|
||||||
|
final int m = pivot.length;
|
||||||
|
cachedL = new Array2DRowFieldMatrix<T>(field, m, m);
|
||||||
|
for (int i = 0; i < m; ++i) {
|
||||||
|
final T[] luI = lu[i];
|
||||||
|
for (int j = 0; j < i; ++j) {
|
||||||
|
cachedL.setEntry(i, j, luI[j]);
|
||||||
|
}
|
||||||
|
cachedL.setEntry(i, i, field.getOne());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return cachedL;
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns the matrix U of the decomposition.
|
* Returns the matrix U of the decomposition.
|
||||||
* <p>U is an upper-triangular matrix</p>
|
* <p>U is an upper-triangular matrix</p>
|
||||||
* @return the U matrix (or null if decomposed matrix is singular)
|
* @return the U matrix (or null if decomposed matrix is singular)
|
||||||
*/
|
*/
|
||||||
FieldMatrix<T> getU();
|
public FieldMatrix<T> getU() {
|
||||||
|
if ((cachedU == null) && !singular) {
|
||||||
|
final int m = pivot.length;
|
||||||
|
cachedU = new Array2DRowFieldMatrix<T>(field, m, m);
|
||||||
|
for (int i = 0; i < m; ++i) {
|
||||||
|
final T[] luI = lu[i];
|
||||||
|
for (int j = i; j < m; ++j) {
|
||||||
|
cachedU.setEntry(i, j, luI[j]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return cachedU;
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns the P rows permutation matrix.
|
* Returns the P rows permutation matrix.
|
||||||
|
@ -70,25 +213,240 @@ public interface FieldLUDecomposition<T extends FieldElement<T>> {
|
||||||
* @return the P rows permutation matrix (or null if decomposed matrix is singular)
|
* @return the P rows permutation matrix (or null if decomposed matrix is singular)
|
||||||
* @see #getPivot()
|
* @see #getPivot()
|
||||||
*/
|
*/
|
||||||
FieldMatrix<T> getP();
|
public FieldMatrix<T> getP() {
|
||||||
|
if ((cachedP == null) && !singular) {
|
||||||
|
final int m = pivot.length;
|
||||||
|
cachedP = new Array2DRowFieldMatrix<T>(field, m, m);
|
||||||
|
for (int i = 0; i < m; ++i) {
|
||||||
|
cachedP.setEntry(i, pivot[i], field.getOne());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return cachedP;
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns the pivot permutation vector.
|
* Returns the pivot permutation vector.
|
||||||
* @return the pivot permutation vector
|
* @return the pivot permutation vector
|
||||||
* @see #getP()
|
* @see #getP()
|
||||||
*/
|
*/
|
||||||
int[] getPivot();
|
public int[] getPivot() {
|
||||||
|
return pivot.clone();
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Return the determinant of the matrix
|
* Return the determinant of the matrix.
|
||||||
* @return determinant of the matrix
|
* @return determinant of the matrix
|
||||||
*/
|
*/
|
||||||
T getDeterminant();
|
public T getDeterminant() {
|
||||||
|
if (singular) {
|
||||||
|
return field.getZero();
|
||||||
|
} else {
|
||||||
|
final int m = pivot.length;
|
||||||
|
T determinant = even ? field.getOne() : field.getZero().subtract(field.getOne());
|
||||||
|
for (int i = 0; i < m; i++) {
|
||||||
|
determinant = determinant.multiply(lu[i][i]);
|
||||||
|
}
|
||||||
|
return determinant;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Get a solver for finding the A × X = B solution in exact linear sense.
|
* Get a solver for finding the A × X = B solution in exact linear sense.
|
||||||
* @return a solver
|
* @return a solver
|
||||||
*/
|
*/
|
||||||
FieldDecompositionSolver<T> getSolver();
|
public FieldDecompositionSolver<T> getSolver() {
|
||||||
|
return new Solver<T>(field, lu, pivot, singular);
|
||||||
|
}
|
||||||
|
|
||||||
|
/** Specialized solver. */
|
||||||
|
private static class Solver<T extends FieldElement<T>> implements FieldDecompositionSolver<T> {
|
||||||
|
|
||||||
|
/** Field to which the elements belong. */
|
||||||
|
private final Field<T> field;
|
||||||
|
|
||||||
|
/** Entries of LU decomposition. */
|
||||||
|
private final T[][] lu;
|
||||||
|
|
||||||
|
/** Pivot permutation associated with LU decomposition. */
|
||||||
|
private final int[] pivot;
|
||||||
|
|
||||||
|
/** Singularity indicator. */
|
||||||
|
private final boolean singular;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Build a solver from decomposed matrix.
|
||||||
|
* @param field field to which the matrix elements belong
|
||||||
|
* @param lu entries of LU decomposition
|
||||||
|
* @param pivot pivot permutation associated with LU decomposition
|
||||||
|
* @param singular singularity indicator
|
||||||
|
*/
|
||||||
|
private Solver(final Field<T> field, final T[][] lu,
|
||||||
|
final int[] pivot, final boolean singular) {
|
||||||
|
this.field = field;
|
||||||
|
this.lu = lu;
|
||||||
|
this.pivot = pivot;
|
||||||
|
this.singular = singular;
|
||||||
|
}
|
||||||
|
|
||||||
|
/** {@inheritDoc} */
|
||||||
|
public boolean isNonSingular() {
|
||||||
|
return !singular;
|
||||||
|
}
|
||||||
|
|
||||||
|
/** {@inheritDoc} */
|
||||||
|
public FieldVector<T> solve(FieldVector<T> b) {
|
||||||
|
try {
|
||||||
|
return solve((ArrayFieldVector<T>) b);
|
||||||
|
} catch (ClassCastException cce) {
|
||||||
|
|
||||||
|
final int m = pivot.length;
|
||||||
|
if (b.getDimension() != m) {
|
||||||
|
throw new DimensionMismatchException(b.getDimension(), m);
|
||||||
|
}
|
||||||
|
if (singular) {
|
||||||
|
throw new SingularMatrixException();
|
||||||
|
}
|
||||||
|
|
||||||
|
@SuppressWarnings("unchecked") // field is of type T
|
||||||
|
final T[] bp = (T[]) Array.newInstance(field.getZero().getClass(), 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 T bpCol = bp[col];
|
||||||
|
for (int i = col + 1; i < m; i++) {
|
||||||
|
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Solve UX = Y
|
||||||
|
for (int col = m - 1; col >= 0; col--) {
|
||||||
|
bp[col] = bp[col].divide(lu[col][col]);
|
||||||
|
final T bpCol = bp[col];
|
||||||
|
for (int i = 0; i < col; i++) {
|
||||||
|
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return new ArrayFieldVector<T>(field, bp, false);
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/** Solve the linear equation A × X = B.
|
||||||
|
* <p>The A matrix is implicit here. It is </p>
|
||||||
|
* @param b right-hand side of the equation A × X = B
|
||||||
|
* @return a vector X such that A × X = B
|
||||||
|
* @throws DimensionMismatchException if the matrices dimensions do not match.
|
||||||
|
* @throws SingularMatrixException if the decomposed matrix is singular.
|
||||||
|
*/
|
||||||
|
public ArrayFieldVector<T> solve(ArrayFieldVector<T> b) {
|
||||||
|
final int m = pivot.length;
|
||||||
|
if (b.data.length != m) {
|
||||||
|
throw new DimensionMismatchException(b.data.length, m);
|
||||||
|
}
|
||||||
|
if (singular) {
|
||||||
|
throw new SingularMatrixException();
|
||||||
|
}
|
||||||
|
|
||||||
|
@SuppressWarnings("unchecked")
|
||||||
|
// field is of type T
|
||||||
|
final T[] bp = (T[]) Array.newInstance(field.getZero().getClass(),
|
||||||
|
m);
|
||||||
|
|
||||||
|
// Apply permutations to b
|
||||||
|
for (int row = 0; row < m; row++) {
|
||||||
|
bp[row] = b.data[pivot[row]];
|
||||||
|
}
|
||||||
|
|
||||||
|
// Solve LY = b
|
||||||
|
for (int col = 0; col < m; col++) {
|
||||||
|
final T bpCol = bp[col];
|
||||||
|
for (int i = col + 1; i < m; i++) {
|
||||||
|
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Solve UX = Y
|
||||||
|
for (int col = m - 1; col >= 0; col--) {
|
||||||
|
bp[col] = bp[col].divide(lu[col][col]);
|
||||||
|
final T bpCol = bp[col];
|
||||||
|
for (int i = 0; i < col; i++) {
|
||||||
|
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return new ArrayFieldVector<T>(bp, false);
|
||||||
|
}
|
||||||
|
|
||||||
|
/** {@inheritDoc} */
|
||||||
|
public FieldMatrix<T> solve(FieldMatrix<T> 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
|
||||||
|
@SuppressWarnings("unchecked") // field is of type T
|
||||||
|
final T[][] bp = (T[][]) Array.newInstance(field.getZero().getClass(), new int[] { m, nColB });
|
||||||
|
for (int row = 0; row < m; row++) {
|
||||||
|
final T[] 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 T[] bpCol = bp[col];
|
||||||
|
for (int i = col + 1; i < m; i++) {
|
||||||
|
final T[] bpI = bp[i];
|
||||||
|
final T luICol = lu[i][col];
|
||||||
|
for (int j = 0; j < nColB; j++) {
|
||||||
|
bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Solve UX = Y
|
||||||
|
for (int col = m - 1; col >= 0; col--) {
|
||||||
|
final T[] bpCol = bp[col];
|
||||||
|
final T luDiag = lu[col][col];
|
||||||
|
for (int j = 0; j < nColB; j++) {
|
||||||
|
bpCol[j] = bpCol[j].divide(luDiag);
|
||||||
|
}
|
||||||
|
for (int i = 0; i < col; i++) {
|
||||||
|
final T[] bpI = bp[i];
|
||||||
|
final T luICol = lu[i][col];
|
||||||
|
for (int j = 0; j < nColB; j++) {
|
||||||
|
bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return new Array2DRowFieldMatrix<T>(field, bp, false);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/** {@inheritDoc} */
|
||||||
|
public FieldMatrix<T> getInverse() {
|
||||||
|
final int m = pivot.length;
|
||||||
|
final T one = field.getOne();
|
||||||
|
FieldMatrix<T> identity = new Array2DRowFieldMatrix<T>(field, m, m);
|
||||||
|
for (int i = 0; i < m; ++i) {
|
||||||
|
identity.setEntry(i, i, one);
|
||||||
|
}
|
||||||
|
return solve(identity);
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
|
@ -1,411 +0,0 @@
|
||||||
/*
|
|
||||||
* 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
|
|
||||||
*
|
|
||||||
* http://www.apache.org/licenses/LICENSE-2.0
|
|
||||||
*
|
|
||||||
* 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.math.linear;
|
|
||||||
|
|
||||||
import java.lang.reflect.Array;
|
|
||||||
|
|
||||||
import org.apache.commons.math.Field;
|
|
||||||
import org.apache.commons.math.FieldElement;
|
|
||||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
|
||||||
|
|
||||||
/**
|
|
||||||
* 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: PA = LU, L is lower triangular, and U is
|
|
||||||
* upper triangular and P is a permutation matrix. All matrices are
|
|
||||||
* m×m.</p>
|
|
||||||
* <p>Since {@link FieldElement field elements} do not provide an ordering
|
|
||||||
* operator, the permutation matrix is computed here only in order to avoid
|
|
||||||
* a zero pivot element, no attempt is done to get the largest pivot element.</p>
|
|
||||||
*
|
|
||||||
* @param <T> the type of the field elements
|
|
||||||
* @version $Id$
|
|
||||||
* @since 2.0
|
|
||||||
*/
|
|
||||||
public class FieldLUDecompositionImpl<T extends FieldElement<T>> implements FieldLUDecomposition<T> {
|
|
||||||
|
|
||||||
/** Field to which the elements belong. */
|
|
||||||
private final Field<T> field;
|
|
||||||
|
|
||||||
/** Entries of LU decomposition. */
|
|
||||||
private T lu[][];
|
|
||||||
|
|
||||||
/** Pivot permutation associated with LU decomposition */
|
|
||||||
private int[] pivot;
|
|
||||||
|
|
||||||
/** Parity of the permutation associated with the LU decomposition */
|
|
||||||
private boolean even;
|
|
||||||
|
|
||||||
/** Singularity indicator. */
|
|
||||||
private boolean singular;
|
|
||||||
|
|
||||||
/** Cached value of L. */
|
|
||||||
private FieldMatrix<T> cachedL;
|
|
||||||
|
|
||||||
/** Cached value of U. */
|
|
||||||
private FieldMatrix<T> cachedU;
|
|
||||||
|
|
||||||
/** Cached value of P. */
|
|
||||||
private FieldMatrix<T> cachedP;
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Calculates the LU-decomposition of the given matrix.
|
|
||||||
* @param matrix The matrix to decompose.
|
|
||||||
* @throws NonSquareMatrixException if matrix is not square
|
|
||||||
*/
|
|
||||||
public FieldLUDecompositionImpl(FieldMatrix<T> matrix) {
|
|
||||||
if (!matrix.isSquare()) {
|
|
||||||
throw new NonSquareMatrixException(matrix.getRowDimension(),
|
|
||||||
matrix.getColumnDimension());
|
|
||||||
}
|
|
||||||
|
|
||||||
final int m = matrix.getColumnDimension();
|
|
||||||
field = matrix.getField();
|
|
||||||
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++) {
|
|
||||||
|
|
||||||
T sum = field.getZero();
|
|
||||||
|
|
||||||
// upper
|
|
||||||
for (int row = 0; row < col; row++) {
|
|
||||||
final T[] luRow = lu[row];
|
|
||||||
sum = luRow[col];
|
|
||||||
for (int i = 0; i < row; i++) {
|
|
||||||
sum = sum.subtract(luRow[i].multiply(lu[i][col]));
|
|
||||||
}
|
|
||||||
luRow[col] = sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
// lower
|
|
||||||
int nonZero = col; // permutation row
|
|
||||||
for (int row = col; row < m; row++) {
|
|
||||||
final T[] luRow = lu[row];
|
|
||||||
sum = luRow[col];
|
|
||||||
for (int i = 0; i < col; i++) {
|
|
||||||
sum = sum.subtract(luRow[i].multiply(lu[i][col]));
|
|
||||||
}
|
|
||||||
luRow[col] = sum;
|
|
||||||
|
|
||||||
if (lu[nonZero][col].equals(field.getZero())) {
|
|
||||||
// try to select a better permutation choice
|
|
||||||
++nonZero;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Singularity check
|
|
||||||
if (nonZero >= m) {
|
|
||||||
singular = true;
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Pivot if necessary
|
|
||||||
if (nonZero != col) {
|
|
||||||
T tmp = field.getZero();
|
|
||||||
for (int i = 0; i < m; i++) {
|
|
||||||
tmp = lu[nonZero][i];
|
|
||||||
lu[nonZero][i] = lu[col][i];
|
|
||||||
lu[col][i] = tmp;
|
|
||||||
}
|
|
||||||
int temp = pivot[nonZero];
|
|
||||||
pivot[nonZero] = pivot[col];
|
|
||||||
pivot[col] = temp;
|
|
||||||
even = !even;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Divide the lower elements by the "winning" diagonal elt.
|
|
||||||
final T luDiag = lu[col][col];
|
|
||||||
for (int row = col + 1; row < m; row++) {
|
|
||||||
final T[] luRow = lu[row];
|
|
||||||
luRow[col] = luRow[col].divide(luDiag);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public FieldMatrix<T> getL() {
|
|
||||||
if ((cachedL == null) && !singular) {
|
|
||||||
final int m = pivot.length;
|
|
||||||
cachedL = new Array2DRowFieldMatrix<T>(field, m, m);
|
|
||||||
for (int i = 0; i < m; ++i) {
|
|
||||||
final T[] luI = lu[i];
|
|
||||||
for (int j = 0; j < i; ++j) {
|
|
||||||
cachedL.setEntry(i, j, luI[j]);
|
|
||||||
}
|
|
||||||
cachedL.setEntry(i, i, field.getOne());
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return cachedL;
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public FieldMatrix<T> getU() {
|
|
||||||
if ((cachedU == null) && !singular) {
|
|
||||||
final int m = pivot.length;
|
|
||||||
cachedU = new Array2DRowFieldMatrix<T>(field, m, m);
|
|
||||||
for (int i = 0; i < m; ++i) {
|
|
||||||
final T[] luI = lu[i];
|
|
||||||
for (int j = i; j < m; ++j) {
|
|
||||||
cachedU.setEntry(i, j, luI[j]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return cachedU;
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public FieldMatrix<T> getP() {
|
|
||||||
if ((cachedP == null) && !singular) {
|
|
||||||
final int m = pivot.length;
|
|
||||||
cachedP = new Array2DRowFieldMatrix<T>(field, m, m);
|
|
||||||
for (int i = 0; i < m; ++i) {
|
|
||||||
cachedP.setEntry(i, pivot[i], field.getOne());
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return cachedP;
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public int[] getPivot() {
|
|
||||||
return pivot.clone();
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public T getDeterminant() {
|
|
||||||
if (singular) {
|
|
||||||
return field.getZero();
|
|
||||||
} else {
|
|
||||||
final int m = pivot.length;
|
|
||||||
T determinant = even ? field.getOne() : field.getZero().subtract(field.getOne());
|
|
||||||
for (int i = 0; i < m; i++) {
|
|
||||||
determinant = determinant.multiply(lu[i][i]);
|
|
||||||
}
|
|
||||||
return determinant;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public FieldDecompositionSolver<T> getSolver() {
|
|
||||||
return new Solver<T>(field, lu, pivot, singular);
|
|
||||||
}
|
|
||||||
|
|
||||||
/** Specialized solver. */
|
|
||||||
private static class Solver<T extends FieldElement<T>> implements FieldDecompositionSolver<T> {
|
|
||||||
|
|
||||||
/** Field to which the elements belong. */
|
|
||||||
private final Field<T> field;
|
|
||||||
|
|
||||||
/** Entries of LU decomposition. */
|
|
||||||
private final T lu[][];
|
|
||||||
|
|
||||||
/** Pivot permutation associated with LU decomposition. */
|
|
||||||
private final int[] pivot;
|
|
||||||
|
|
||||||
/** Singularity indicator. */
|
|
||||||
private final boolean singular;
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Build a solver from decomposed matrix.
|
|
||||||
* @param field field to which the matrix elements belong
|
|
||||||
* @param lu entries of LU decomposition
|
|
||||||
* @param pivot pivot permutation associated with LU decomposition
|
|
||||||
* @param singular singularity indicator
|
|
||||||
*/
|
|
||||||
private Solver(final Field<T> field, final T[][] lu,
|
|
||||||
final int[] pivot, final boolean singular) {
|
|
||||||
this.field = field;
|
|
||||||
this.lu = lu;
|
|
||||||
this.pivot = pivot;
|
|
||||||
this.singular = singular;
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public boolean isNonSingular() {
|
|
||||||
return !singular;
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public FieldVector<T> solve(FieldVector<T> b) {
|
|
||||||
try {
|
|
||||||
return solve((ArrayFieldVector<T>) b);
|
|
||||||
} catch (ClassCastException cce) {
|
|
||||||
|
|
||||||
final int m = pivot.length;
|
|
||||||
if (b.getDimension() != m) {
|
|
||||||
throw new DimensionMismatchException(b.getDimension(), m);
|
|
||||||
}
|
|
||||||
if (singular) {
|
|
||||||
throw new SingularMatrixException();
|
|
||||||
}
|
|
||||||
|
|
||||||
@SuppressWarnings("unchecked") // field is of type T
|
|
||||||
final T[] bp = (T[]) Array.newInstance(field.getZero().getClass(), 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 T bpCol = bp[col];
|
|
||||||
for (int i = col + 1; i < m; i++) {
|
|
||||||
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Solve UX = Y
|
|
||||||
for (int col = m - 1; col >= 0; col--) {
|
|
||||||
bp[col] = bp[col].divide(lu[col][col]);
|
|
||||||
final T bpCol = bp[col];
|
|
||||||
for (int i = 0; i < col; i++) {
|
|
||||||
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return new ArrayFieldVector<T>(field, bp, false);
|
|
||||||
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/** Solve the linear equation A × X = B.
|
|
||||||
* <p>The A matrix is implicit here. It is </p>
|
|
||||||
* @param b right-hand side of the equation A × X = B
|
|
||||||
* @return a vector X such that A × X = B
|
|
||||||
* @throws DimensionMismatchException if the matrices dimensions do not match.
|
|
||||||
* @throws SingularMatrixException if the decomposed matrix is singular.
|
|
||||||
*/
|
|
||||||
public ArrayFieldVector<T> solve(ArrayFieldVector<T> b) {
|
|
||||||
final int m = pivot.length;
|
|
||||||
if (b.data.length != m) {
|
|
||||||
throw new DimensionMismatchException(b.data.length, m);
|
|
||||||
}
|
|
||||||
if (singular) {
|
|
||||||
throw new SingularMatrixException();
|
|
||||||
}
|
|
||||||
|
|
||||||
@SuppressWarnings("unchecked")
|
|
||||||
// field is of type T
|
|
||||||
final T[] bp = (T[]) Array.newInstance(field.getZero().getClass(),
|
|
||||||
m);
|
|
||||||
|
|
||||||
// Apply permutations to b
|
|
||||||
for (int row = 0; row < m; row++) {
|
|
||||||
bp[row] = b.data[pivot[row]];
|
|
||||||
}
|
|
||||||
|
|
||||||
// Solve LY = b
|
|
||||||
for (int col = 0; col < m; col++) {
|
|
||||||
final T bpCol = bp[col];
|
|
||||||
for (int i = col + 1; i < m; i++) {
|
|
||||||
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Solve UX = Y
|
|
||||||
for (int col = m - 1; col >= 0; col--) {
|
|
||||||
bp[col] = bp[col].divide(lu[col][col]);
|
|
||||||
final T bpCol = bp[col];
|
|
||||||
for (int i = 0; i < col; i++) {
|
|
||||||
bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return new ArrayFieldVector<T>(bp, false);
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public FieldMatrix<T> solve(FieldMatrix<T> 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
|
|
||||||
@SuppressWarnings("unchecked") // field is of type T
|
|
||||||
final T[][] bp = (T[][]) Array.newInstance(field.getZero().getClass(), new int[] { m, nColB });
|
|
||||||
for (int row = 0; row < m; row++) {
|
|
||||||
final T[] 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 T[] bpCol = bp[col];
|
|
||||||
for (int i = col + 1; i < m; i++) {
|
|
||||||
final T[] bpI = bp[i];
|
|
||||||
final T luICol = lu[i][col];
|
|
||||||
for (int j = 0; j < nColB; j++) {
|
|
||||||
bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Solve UX = Y
|
|
||||||
for (int col = m - 1; col >= 0; col--) {
|
|
||||||
final T[] bpCol = bp[col];
|
|
||||||
final T luDiag = lu[col][col];
|
|
||||||
for (int j = 0; j < nColB; j++) {
|
|
||||||
bpCol[j] = bpCol[j].divide(luDiag);
|
|
||||||
}
|
|
||||||
for (int i = 0; i < col; i++) {
|
|
||||||
final T[] bpI = bp[i];
|
|
||||||
final T luICol = lu[i][col];
|
|
||||||
for (int j = 0; j < nColB; j++) {
|
|
||||||
bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return new Array2DRowFieldMatrix<T>(field, bp, false);
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/** {@inheritDoc} */
|
|
||||||
public FieldMatrix<T> getInverse() {
|
|
||||||
final int m = pivot.length;
|
|
||||||
final T one = field.getOne();
|
|
||||||
FieldMatrix<T> identity = new Array2DRowFieldMatrix<T>(field, m, m);
|
|
||||||
for (int i = 0; i < m; ++i) {
|
|
||||||
identity.setEntry(i, i, one);
|
|
||||||
}
|
|
||||||
return solve(identity);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
|
@ -26,7 +26,7 @@ import org.apache.commons.math.linear.Array2DRowFieldMatrix;
|
||||||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||||
import org.apache.commons.math.linear.ArrayFieldVector;
|
import org.apache.commons.math.linear.ArrayFieldVector;
|
||||||
import org.apache.commons.math.linear.FieldDecompositionSolver;
|
import org.apache.commons.math.linear.FieldDecompositionSolver;
|
||||||
import org.apache.commons.math.linear.FieldLUDecompositionImpl;
|
import org.apache.commons.math.linear.FieldLUDecomposition;
|
||||||
import org.apache.commons.math.linear.FieldMatrix;
|
import org.apache.commons.math.linear.FieldMatrix;
|
||||||
import org.apache.commons.math.linear.MatrixUtils;
|
import org.apache.commons.math.linear.MatrixUtils;
|
||||||
import org.apache.commons.math.linear.QRDecomposition;
|
import org.apache.commons.math.linear.QRDecomposition;
|
||||||
|
@ -154,7 +154,7 @@ public class AdamsNordsieckTransformer {
|
||||||
// compute exact coefficients
|
// compute exact coefficients
|
||||||
FieldMatrix<BigFraction> bigP = buildP(nSteps);
|
FieldMatrix<BigFraction> bigP = buildP(nSteps);
|
||||||
FieldDecompositionSolver<BigFraction> pSolver =
|
FieldDecompositionSolver<BigFraction> pSolver =
|
||||||
new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver();
|
new FieldLUDecomposition<BigFraction>(bigP).getSolver();
|
||||||
|
|
||||||
BigFraction[] u = new BigFraction[nSteps];
|
BigFraction[] u = new BigFraction[nSteps];
|
||||||
Arrays.fill(u, BigFraction.ONE);
|
Arrays.fill(u, BigFraction.ONE);
|
||||||
|
|
|
@ -434,8 +434,8 @@ public final class BlockFieldMatrixTest {
|
||||||
@Test
|
@Test
|
||||||
public void testTranspose() {
|
public void testTranspose() {
|
||||||
FieldMatrix<Fraction> m = new BlockFieldMatrix<Fraction>(testData);
|
FieldMatrix<Fraction> m = new BlockFieldMatrix<Fraction>(testData);
|
||||||
FieldMatrix<Fraction> mIT = new FieldLUDecompositionImpl<Fraction>(m).getSolver().getInverse().transpose();
|
FieldMatrix<Fraction> mIT = new FieldLUDecomposition<Fraction>(m).getSolver().getInverse().transpose();
|
||||||
FieldMatrix<Fraction> mTI = new FieldLUDecompositionImpl<Fraction>(m.transpose()).getSolver().getInverse();
|
FieldMatrix<Fraction> mTI = new FieldLUDecomposition<Fraction>(m.transpose()).getSolver().getInverse();
|
||||||
TestUtils.assertEquals(mIT, mTI);
|
TestUtils.assertEquals(mIT, mTI);
|
||||||
m = new BlockFieldMatrix<Fraction>(testData2);
|
m = new BlockFieldMatrix<Fraction>(testData2);
|
||||||
FieldMatrix<Fraction> mt = new BlockFieldMatrix<Fraction>(testData2T);
|
FieldMatrix<Fraction> mt = new BlockFieldMatrix<Fraction>(testData2T);
|
||||||
|
@ -532,7 +532,7 @@ public final class BlockFieldMatrixTest {
|
||||||
Assert.assertEquals(2, p.getRowDimension());
|
Assert.assertEquals(2, p.getRowDimension());
|
||||||
Assert.assertEquals(2, p.getColumnDimension());
|
Assert.assertEquals(2, p.getColumnDimension());
|
||||||
// Invert p
|
// Invert p
|
||||||
FieldMatrix<Fraction> pInverse = new FieldLUDecompositionImpl<Fraction>(p).getSolver().getInverse();
|
FieldMatrix<Fraction> pInverse = new FieldLUDecomposition<Fraction>(p).getSolver().getInverse();
|
||||||
Assert.assertEquals(2, pInverse.getRowDimension());
|
Assert.assertEquals(2, pInverse.getRowDimension());
|
||||||
Assert.assertEquals(2, pInverse.getColumnDimension());
|
Assert.assertEquals(2, pInverse.getColumnDimension());
|
||||||
|
|
||||||
|
@ -547,7 +547,7 @@ public final class BlockFieldMatrixTest {
|
||||||
new Fraction(1), new Fraction(-2), new Fraction(1)
|
new Fraction(1), new Fraction(-2), new Fraction(1)
|
||||||
};
|
};
|
||||||
Fraction[] solution;
|
Fraction[] solution;
|
||||||
solution = new FieldLUDecompositionImpl<Fraction>(coefficients)
|
solution = new FieldLUDecomposition<Fraction>(coefficients)
|
||||||
.getSolver()
|
.getSolver()
|
||||||
.solve(new ArrayFieldVector<Fraction>(constants, false)).toArray();
|
.solve(new ArrayFieldVector<Fraction>(constants, false)).toArray();
|
||||||
Assert.assertEquals(new Fraction(2).multiply(solution[0]).
|
Assert.assertEquals(new Fraction(2).multiply(solution[0]).
|
||||||
|
|
|
@ -24,7 +24,7 @@ import org.apache.commons.math.TestUtils;
|
||||||
import org.apache.commons.math.fraction.Fraction;
|
import org.apache.commons.math.fraction.Fraction;
|
||||||
import org.apache.commons.math.fraction.FractionField;
|
import org.apache.commons.math.fraction.FractionField;
|
||||||
|
|
||||||
public class FieldLUDecompositionImplTest {
|
public class FieldLUDecompositionTest {
|
||||||
private Fraction[][] testData = {
|
private Fraction[][] testData = {
|
||||||
{ new Fraction(1), new Fraction(2), new Fraction(3)},
|
{ new Fraction(1), new Fraction(2), new Fraction(3)},
|
||||||
{ new Fraction(2), new Fraction(5), new Fraction(3)},
|
{ new Fraction(2), new Fraction(5), new Fraction(3)},
|
||||||
|
@ -58,7 +58,7 @@ public class FieldLUDecompositionImplTest {
|
||||||
public void testDimensions() {
|
public void testDimensions() {
|
||||||
FieldMatrix<Fraction> matrix =
|
FieldMatrix<Fraction> matrix =
|
||||||
new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
||||||
FieldLUDecomposition<Fraction> LU = new FieldLUDecompositionImpl<Fraction>(matrix);
|
FieldLUDecomposition<Fraction> LU = new FieldLUDecomposition<Fraction>(matrix);
|
||||||
Assert.assertEquals(testData.length, LU.getL().getRowDimension());
|
Assert.assertEquals(testData.length, LU.getL().getRowDimension());
|
||||||
Assert.assertEquals(testData.length, LU.getL().getColumnDimension());
|
Assert.assertEquals(testData.length, LU.getL().getColumnDimension());
|
||||||
Assert.assertEquals(testData.length, LU.getU().getRowDimension());
|
Assert.assertEquals(testData.length, LU.getU().getRowDimension());
|
||||||
|
@ -73,7 +73,7 @@ public class FieldLUDecompositionImplTest {
|
||||||
public void testNonSquare() {
|
public void testNonSquare() {
|
||||||
try {
|
try {
|
||||||
// we don't use FractionField.getInstance() for testing purposes
|
// we don't use FractionField.getInstance() for testing purposes
|
||||||
new FieldLUDecompositionImpl<Fraction>(new Array2DRowFieldMatrix<Fraction>(new Fraction[][] {
|
new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(new Fraction[][] {
|
||||||
{ Fraction.ZERO, Fraction.ZERO },
|
{ Fraction.ZERO, Fraction.ZERO },
|
||||||
{ Fraction.ZERO, Fraction.ZERO },
|
{ Fraction.ZERO, Fraction.ZERO },
|
||||||
{ Fraction.ZERO, Fraction.ZERO }
|
{ Fraction.ZERO, Fraction.ZERO }
|
||||||
|
@ -88,14 +88,14 @@ public class FieldLUDecompositionImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testPAEqualLU() {
|
public void testPAEqualLU() {
|
||||||
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
||||||
FieldLUDecomposition<Fraction> lu = new FieldLUDecompositionImpl<Fraction>(matrix);
|
FieldLUDecomposition<Fraction> lu = new FieldLUDecomposition<Fraction>(matrix);
|
||||||
FieldMatrix<Fraction> l = lu.getL();
|
FieldMatrix<Fraction> l = lu.getL();
|
||||||
FieldMatrix<Fraction> u = lu.getU();
|
FieldMatrix<Fraction> u = lu.getU();
|
||||||
FieldMatrix<Fraction> p = lu.getP();
|
FieldMatrix<Fraction> p = lu.getP();
|
||||||
TestUtils.assertEquals(p.multiply(matrix), l.multiply(u));
|
TestUtils.assertEquals(p.multiply(matrix), l.multiply(u));
|
||||||
|
|
||||||
matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testDataMinus);
|
matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testDataMinus);
|
||||||
lu = new FieldLUDecompositionImpl<Fraction>(matrix);
|
lu = new FieldLUDecomposition<Fraction>(matrix);
|
||||||
l = lu.getL();
|
l = lu.getL();
|
||||||
u = lu.getU();
|
u = lu.getU();
|
||||||
p = lu.getP();
|
p = lu.getP();
|
||||||
|
@ -105,21 +105,21 @@ public class FieldLUDecompositionImplTest {
|
||||||
for (int i = 0; i < matrix.getRowDimension(); ++i) {
|
for (int i = 0; i < matrix.getRowDimension(); ++i) {
|
||||||
matrix.setEntry(i, i, Fraction.ONE);
|
matrix.setEntry(i, i, Fraction.ONE);
|
||||||
}
|
}
|
||||||
lu = new FieldLUDecompositionImpl<Fraction>(matrix);
|
lu = new FieldLUDecomposition<Fraction>(matrix);
|
||||||
l = lu.getL();
|
l = lu.getL();
|
||||||
u = lu.getU();
|
u = lu.getU();
|
||||||
p = lu.getP();
|
p = lu.getP();
|
||||||
TestUtils.assertEquals(p.multiply(matrix), l.multiply(u));
|
TestUtils.assertEquals(p.multiply(matrix), l.multiply(u));
|
||||||
|
|
||||||
matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), singular);
|
matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), singular);
|
||||||
lu = new FieldLUDecompositionImpl<Fraction>(matrix);
|
lu = new FieldLUDecomposition<Fraction>(matrix);
|
||||||
Assert.assertFalse(lu.getSolver().isNonSingular());
|
Assert.assertFalse(lu.getSolver().isNonSingular());
|
||||||
Assert.assertNull(lu.getL());
|
Assert.assertNull(lu.getL());
|
||||||
Assert.assertNull(lu.getU());
|
Assert.assertNull(lu.getU());
|
||||||
Assert.assertNull(lu.getP());
|
Assert.assertNull(lu.getP());
|
||||||
|
|
||||||
matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), bigSingular);
|
matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), bigSingular);
|
||||||
lu = new FieldLUDecompositionImpl<Fraction>(matrix);
|
lu = new FieldLUDecomposition<Fraction>(matrix);
|
||||||
Assert.assertFalse(lu.getSolver().isNonSingular());
|
Assert.assertFalse(lu.getSolver().isNonSingular());
|
||||||
Assert.assertNull(lu.getL());
|
Assert.assertNull(lu.getL());
|
||||||
Assert.assertNull(lu.getU());
|
Assert.assertNull(lu.getU());
|
||||||
|
@ -131,7 +131,7 @@ public class FieldLUDecompositionImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testLLowerTriangular() {
|
public void testLLowerTriangular() {
|
||||||
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
||||||
FieldMatrix<Fraction> l = new FieldLUDecompositionImpl<Fraction>(matrix).getL();
|
FieldMatrix<Fraction> l = new FieldLUDecomposition<Fraction>(matrix).getL();
|
||||||
for (int i = 0; i < l.getRowDimension(); i++) {
|
for (int i = 0; i < l.getRowDimension(); i++) {
|
||||||
Assert.assertEquals(Fraction.ONE, l.getEntry(i, i));
|
Assert.assertEquals(Fraction.ONE, l.getEntry(i, i));
|
||||||
for (int j = i + 1; j < l.getColumnDimension(); j++) {
|
for (int j = i + 1; j < l.getColumnDimension(); j++) {
|
||||||
|
@ -144,7 +144,7 @@ public class FieldLUDecompositionImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testUUpperTriangular() {
|
public void testUUpperTriangular() {
|
||||||
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
||||||
FieldMatrix<Fraction> u = new FieldLUDecompositionImpl<Fraction>(matrix).getU();
|
FieldMatrix<Fraction> u = new FieldLUDecomposition<Fraction>(matrix).getU();
|
||||||
for (int i = 0; i < u.getRowDimension(); i++) {
|
for (int i = 0; i < u.getRowDimension(); i++) {
|
||||||
for (int j = 0; j < i; j++) {
|
for (int j = 0; j < i; j++) {
|
||||||
Assert.assertEquals(Fraction.ZERO, u.getEntry(i, j));
|
Assert.assertEquals(Fraction.ZERO, u.getEntry(i, j));
|
||||||
|
@ -156,7 +156,7 @@ public class FieldLUDecompositionImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testPPermutation() {
|
public void testPPermutation() {
|
||||||
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
|
||||||
FieldMatrix<Fraction> p = new FieldLUDecompositionImpl<Fraction>(matrix).getP();
|
FieldMatrix<Fraction> p = new FieldLUDecomposition<Fraction>(matrix).getP();
|
||||||
|
|
||||||
FieldMatrix<Fraction> ppT = p.multiply(p.transpose());
|
FieldMatrix<Fraction> ppT = p.multiply(p.transpose());
|
||||||
FieldMatrix<Fraction> id =
|
FieldMatrix<Fraction> id =
|
||||||
|
@ -212,11 +212,11 @@ public class FieldLUDecompositionImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testSingular() {
|
public void testSingular() {
|
||||||
FieldLUDecomposition<Fraction> lu =
|
FieldLUDecomposition<Fraction> lu =
|
||||||
new FieldLUDecompositionImpl<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData));
|
new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData));
|
||||||
Assert.assertTrue(lu.getSolver().isNonSingular());
|
Assert.assertTrue(lu.getSolver().isNonSingular());
|
||||||
lu = new FieldLUDecompositionImpl<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), singular));
|
lu = new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), singular));
|
||||||
Assert.assertFalse(lu.getSolver().isNonSingular());
|
Assert.assertFalse(lu.getSolver().isNonSingular());
|
||||||
lu = new FieldLUDecompositionImpl<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), bigSingular));
|
lu = new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), bigSingular));
|
||||||
Assert.assertFalse(lu.getSolver().isNonSingular());
|
Assert.assertFalse(lu.getSolver().isNonSingular());
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -224,7 +224,7 @@ public class FieldLUDecompositionImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testMatricesValues1() {
|
public void testMatricesValues1() {
|
||||||
FieldLUDecomposition<Fraction> lu =
|
FieldLUDecomposition<Fraction> lu =
|
||||||
new FieldLUDecompositionImpl<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData));
|
new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData));
|
||||||
FieldMatrix<Fraction> lRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
|
FieldMatrix<Fraction> lRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
|
||||||
{ new Fraction(1), new Fraction(0), new Fraction(0) },
|
{ new Fraction(1), new Fraction(0), new Fraction(0) },
|
||||||
{ new Fraction(2), new Fraction(1), new Fraction(0) },
|
{ new Fraction(2), new Fraction(1), new Fraction(0) },
|
||||||
|
@ -265,7 +265,7 @@ public class FieldLUDecompositionImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testMatricesValues2() {
|
public void testMatricesValues2() {
|
||||||
FieldLUDecomposition<Fraction> lu =
|
FieldLUDecomposition<Fraction> lu =
|
||||||
new FieldLUDecompositionImpl<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), luData));
|
new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), luData));
|
||||||
FieldMatrix<Fraction> lRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
|
FieldMatrix<Fraction> lRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
|
||||||
{ new Fraction(1), new Fraction(0), new Fraction(0) },
|
{ new Fraction(1), new Fraction(0), new Fraction(0) },
|
||||||
{ new Fraction(3), new Fraction(1), new Fraction(0) },
|
{ new Fraction(3), new Fraction(1), new Fraction(0) },
|
|
@ -65,13 +65,13 @@ public class FieldLUSolverTest {
|
||||||
@Test
|
@Test
|
||||||
public void testSingular() {
|
public void testSingular() {
|
||||||
FieldDecompositionSolver<Fraction> solver;
|
FieldDecompositionSolver<Fraction> solver;
|
||||||
solver = new FieldLUDecompositionImpl<Fraction>(createFractionMatrix(testData))
|
solver = new FieldLUDecomposition<Fraction>(createFractionMatrix(testData))
|
||||||
.getSolver();
|
.getSolver();
|
||||||
Assert.assertTrue(solver.isNonSingular());
|
Assert.assertTrue(solver.isNonSingular());
|
||||||
solver = new FieldLUDecompositionImpl<Fraction>(createFractionMatrix(singular))
|
solver = new FieldLUDecomposition<Fraction>(createFractionMatrix(singular))
|
||||||
.getSolver();
|
.getSolver();
|
||||||
Assert.assertFalse(solver.isNonSingular());
|
Assert.assertFalse(solver.isNonSingular());
|
||||||
solver = new FieldLUDecompositionImpl<Fraction>(createFractionMatrix(bigSingular))
|
solver = new FieldLUDecomposition<Fraction>(createFractionMatrix(bigSingular))
|
||||||
.getSolver();
|
.getSolver();
|
||||||
Assert.assertFalse(solver.isNonSingular());
|
Assert.assertFalse(solver.isNonSingular());
|
||||||
}
|
}
|
||||||
|
@ -80,7 +80,7 @@ public class FieldLUSolverTest {
|
||||||
@Test
|
@Test
|
||||||
public void testSolveDimensionErrors() {
|
public void testSolveDimensionErrors() {
|
||||||
FieldDecompositionSolver<Fraction> solver;
|
FieldDecompositionSolver<Fraction> solver;
|
||||||
solver = new FieldLUDecompositionImpl<Fraction>(createFractionMatrix(testData))
|
solver = new FieldLUDecomposition<Fraction>(createFractionMatrix(testData))
|
||||||
.getSolver();
|
.getSolver();
|
||||||
FieldMatrix<Fraction> b = createFractionMatrix(new int[2][2]);
|
FieldMatrix<Fraction> b = createFractionMatrix(new int[2][2]);
|
||||||
try {
|
try {
|
||||||
|
@ -101,7 +101,7 @@ public class FieldLUSolverTest {
|
||||||
@Test
|
@Test
|
||||||
public void testSolveSingularityErrors() {
|
public void testSolveSingularityErrors() {
|
||||||
FieldDecompositionSolver solver;
|
FieldDecompositionSolver solver;
|
||||||
solver = new FieldLUDecompositionImpl(createFractionMatrix(singular))
|
solver = new FieldLUDecomposition(createFractionMatrix(singular))
|
||||||
.getSolver();
|
.getSolver();
|
||||||
FieldMatrix b = createFractionMatrix(new int[2][2]);
|
FieldMatrix b = createFractionMatrix(new int[2][2]);
|
||||||
try {
|
try {
|
||||||
|
@ -122,7 +122,7 @@ public class FieldLUSolverTest {
|
||||||
@Test
|
@Test
|
||||||
public void testSolve() {
|
public void testSolve() {
|
||||||
FieldDecompositionSolver solver;
|
FieldDecompositionSolver solver;
|
||||||
solver = new FieldLUDecompositionImpl<Fraction>(createFractionMatrix(testData))
|
solver = new FieldLUDecomposition<Fraction>(createFractionMatrix(testData))
|
||||||
.getSolver();
|
.getSolver();
|
||||||
FieldMatrix<Fraction> b = createFractionMatrix(new int[][] {
|
FieldMatrix<Fraction> b = createFractionMatrix(new int[][] {
|
||||||
{ 1, 0 }, { 2, -5 }, { 3, 1 }
|
{ 1, 0 }, { 2, -5 }, { 3, 1 }
|
||||||
|
@ -172,6 +172,6 @@ public class FieldLUSolverTest {
|
||||||
}
|
}
|
||||||
|
|
||||||
private double getDeterminant(final FieldMatrix<Fraction> m) {
|
private double getDeterminant(final FieldMatrix<Fraction> m) {
|
||||||
return new FieldLUDecompositionImpl<Fraction>(m).getDeterminant().doubleValue();
|
return new FieldLUDecomposition<Fraction>(m).getDeterminant().doubleValue();
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
|
@ -297,8 +297,8 @@ public final class FieldMatrixImplTest {
|
||||||
@Test
|
@Test
|
||||||
public void testTranspose() {
|
public void testTranspose() {
|
||||||
FieldMatrix<Fraction> m = new Array2DRowFieldMatrix<Fraction>(testData);
|
FieldMatrix<Fraction> m = new Array2DRowFieldMatrix<Fraction>(testData);
|
||||||
FieldMatrix<Fraction> mIT = new FieldLUDecompositionImpl<Fraction>(m).getSolver().getInverse().transpose();
|
FieldMatrix<Fraction> mIT = new FieldLUDecomposition<Fraction>(m).getSolver().getInverse().transpose();
|
||||||
FieldMatrix<Fraction> mTI = new FieldLUDecompositionImpl<Fraction>(m.transpose()).getSolver().getInverse();
|
FieldMatrix<Fraction> mTI = new FieldLUDecomposition<Fraction>(m.transpose()).getSolver().getInverse();
|
||||||
TestUtils.assertEquals(mIT, mTI);
|
TestUtils.assertEquals(mIT, mTI);
|
||||||
m = new Array2DRowFieldMatrix<Fraction>(testData2);
|
m = new Array2DRowFieldMatrix<Fraction>(testData2);
|
||||||
FieldMatrix<Fraction> mt = new Array2DRowFieldMatrix<Fraction>(testData2T);
|
FieldMatrix<Fraction> mt = new Array2DRowFieldMatrix<Fraction>(testData2T);
|
||||||
|
@ -395,7 +395,7 @@ public final class FieldMatrixImplTest {
|
||||||
Assert.assertEquals(2, p.getRowDimension());
|
Assert.assertEquals(2, p.getRowDimension());
|
||||||
Assert.assertEquals(2, p.getColumnDimension());
|
Assert.assertEquals(2, p.getColumnDimension());
|
||||||
// Invert p
|
// Invert p
|
||||||
FieldMatrix<Fraction> pInverse = new FieldLUDecompositionImpl<Fraction>(p).getSolver().getInverse();
|
FieldMatrix<Fraction> pInverse = new FieldLUDecomposition<Fraction>(p).getSolver().getInverse();
|
||||||
Assert.assertEquals(2, pInverse.getRowDimension());
|
Assert.assertEquals(2, pInverse.getRowDimension());
|
||||||
Assert.assertEquals(2, pInverse.getColumnDimension());
|
Assert.assertEquals(2, pInverse.getColumnDimension());
|
||||||
|
|
||||||
|
@ -410,7 +410,7 @@ public final class FieldMatrixImplTest {
|
||||||
new Fraction(1), new Fraction(-2), new Fraction(1)
|
new Fraction(1), new Fraction(-2), new Fraction(1)
|
||||||
};
|
};
|
||||||
Fraction[] solution;
|
Fraction[] solution;
|
||||||
solution = new FieldLUDecompositionImpl<Fraction>(coefficients)
|
solution = new FieldLUDecomposition<Fraction>(coefficients)
|
||||||
.getSolver()
|
.getSolver()
|
||||||
.solve(new ArrayFieldVector<Fraction>(constants, false)).toArray();
|
.solve(new ArrayFieldVector<Fraction>(constants, false)).toArray();
|
||||||
Assert.assertEquals(new Fraction(2).multiply(solution[0]).
|
Assert.assertEquals(new Fraction(2).multiply(solution[0]).
|
||||||
|
|
|
@ -291,8 +291,8 @@ public class SparseFieldMatrixTest {
|
||||||
@Test
|
@Test
|
||||||
public void testTranspose() {
|
public void testTranspose() {
|
||||||
FieldMatrix<Fraction> m = createSparseMatrix(testData);
|
FieldMatrix<Fraction> m = createSparseMatrix(testData);
|
||||||
FieldMatrix<Fraction> mIT = new FieldLUDecompositionImpl<Fraction>(m).getSolver().getInverse().transpose();
|
FieldMatrix<Fraction> mIT = new FieldLUDecomposition<Fraction>(m).getSolver().getInverse().transpose();
|
||||||
FieldMatrix<Fraction> mTI = new FieldLUDecompositionImpl<Fraction>(m.transpose()).getSolver().getInverse();
|
FieldMatrix<Fraction> mTI = new FieldLUDecomposition<Fraction>(m.transpose()).getSolver().getInverse();
|
||||||
assertClose("inverse-transpose", mIT, mTI, normTolerance);
|
assertClose("inverse-transpose", mIT, mTI, normTolerance);
|
||||||
m = createSparseMatrix(testData2);
|
m = createSparseMatrix(testData2);
|
||||||
FieldMatrix<Fraction> mt = createSparseMatrix(testData2T);
|
FieldMatrix<Fraction> mt = createSparseMatrix(testData2T);
|
||||||
|
@ -387,7 +387,7 @@ public class SparseFieldMatrixTest {
|
||||||
Assert.assertEquals(2, p.getRowDimension());
|
Assert.assertEquals(2, p.getRowDimension());
|
||||||
Assert.assertEquals(2, p.getColumnDimension());
|
Assert.assertEquals(2, p.getColumnDimension());
|
||||||
// Invert p
|
// Invert p
|
||||||
FieldMatrix<Fraction> pInverse = new FieldLUDecompositionImpl<Fraction>(p).getSolver().getInverse();
|
FieldMatrix<Fraction> pInverse = new FieldLUDecomposition<Fraction>(p).getSolver().getInverse();
|
||||||
Assert.assertEquals(2, pInverse.getRowDimension());
|
Assert.assertEquals(2, pInverse.getRowDimension());
|
||||||
Assert.assertEquals(2, pInverse.getColumnDimension());
|
Assert.assertEquals(2, pInverse.getColumnDimension());
|
||||||
|
|
||||||
|
@ -397,7 +397,7 @@ public class SparseFieldMatrixTest {
|
||||||
FieldMatrix<Fraction> coefficients = createSparseMatrix(coefficientsData);
|
FieldMatrix<Fraction> coefficients = createSparseMatrix(coefficientsData);
|
||||||
Fraction[] constants = { new Fraction(1), new Fraction(-2), new Fraction(1) };
|
Fraction[] constants = { new Fraction(1), new Fraction(-2), new Fraction(1) };
|
||||||
Fraction[] solution;
|
Fraction[] solution;
|
||||||
solution = new FieldLUDecompositionImpl<Fraction>(coefficients)
|
solution = new FieldLUDecomposition<Fraction>(coefficients)
|
||||||
.getSolver()
|
.getSolver()
|
||||||
.solve(new ArrayFieldVector<Fraction>(constants, false)).toArray();
|
.solve(new ArrayFieldVector<Fraction>(constants, false)).toArray();
|
||||||
Assert.assertEquals((new Fraction(2).multiply((solution[0])).add(new Fraction(3).multiply(solution[1])).subtract(new Fraction(2).multiply(solution[2]))).doubleValue(),
|
Assert.assertEquals((new Fraction(2).multiply((solution[0])).add(new Fraction(3).multiply(solution[1])).subtract(new Fraction(2).multiply(solution[2]))).doubleValue(),
|
||||||
|
|
Loading…
Reference in New Issue