/*
* Copyright 2011-2013, by Vladimir Kostyukov, Mike Anderson and Contributors.
*
* This file is adapted from the la4j project (http://la4j.org)
*
* Licensed 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.
*
* Contributor(s): -
*
*/
package mikera.matrixx.decompose.impl.lu;
import mikera.vectorz.Vector;
import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.decompose.ILUPResult;
import mikera.matrixx.impl.PermutationMatrix;
/**
* Performs a standard LUP decomposition
*
* @author Mike
*
*/
public class SimpleLUP {
public static ILUPResult decompose(AMatrix matrix) {
return decomposeLUPInternal(Matrix.create(matrix));
}
public static ILUPResult decomposeLUP(Matrix matrix) {
return decomposeLUPInternal(matrix.clone());
}
/**
* Performs a LUP decomposition on the given matrix.
*
* Warning: Destructively modifies the source matrix
* @param m
* @return
*/
private static ILUPResult decomposeLUPInternal(Matrix m) {
if (!m.isSquare()) {
throw new IllegalArgumentException("Wrong matrix size: " + "not square");
}
int n = m.rowCount();
PermutationMatrix p = PermutationMatrix.createIdentity(n);
for (int j = 0; j < n; j++) {
Vector jcolumn = m.getColumn(j).toVector();
for (int i = 0; i < n; i++) {
int kmax = Math.min(i, j);
double s = 0.0;
for (int k = 0; k < kmax; k++) {
s += m.get(i, k) * jcolumn.unsafeGet(k);
}
jcolumn.set(i, jcolumn.unsafeGet(i) - s);
m.set(i, j, jcolumn.unsafeGet(i));
}
int biggest = j;
for (int i = j + 1; i < n; i++) {
if (Math.abs(jcolumn.unsafeGet(i)) > Math.abs(jcolumn.unsafeGet(biggest)))
biggest = i;
}
if (biggest != j) {
m.swapRows(biggest, j);
p.swapRows(biggest, j);
}
if ((j < n) && (m.get(j, j) != 0.0)) {
for (int i = j + 1; i < n; i++) {
m.set(i, j, m.get(i, j) / m.get(j, j));
}
}
}
Matrix l = Matrix.create(n, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < i; j++) {
l.unsafeSet(i, j, m.get(i, j));
}
l.unsafeSet(i, i, 1.0);
}
// clear low elements to ensure upper triangle only is populated
Matrix u = m;
for (int i = 0; i < n; i++) {
for (int j = 0; j < i; j++) {
u.unsafeSet(i, j, 0.0);
}
}
return new LUPResult (l, u, p);
}
}