Instantly share code, notes, and snippets.

@sreiter /LUSolver.txt Secret
Created Jan 24, 2018

Embed
What would you like to do?
LU solver with pivotization for VRL Studio
package eu.mihosoft.vrl.user;
// LU Solver with pivotization
@ComponentInfo(name="LUSolver", category="Custom")
public class LUSolver implements java.io.Serializable {
private static final long serialVersionUID=1L;
private transient double[][] L;
private transient double[][] U;
private transient double[] x;
private transient int[] perm;
private transient int n;
private void init(int size)
{
n = size;
L = new double[n][n];
U = new double[n][n];
perm = new int[n];
}
@MethodInfo(noGUI=true)
void computeLU (double[][] A)
{
if(A.length != n)
init(A.length)
for(int i = 0; i < n; ++i)
perm[i] = i;
for (int j = 0; j < n; j++) {
L[j][j] = 1;
// check for necessary permutation.
// For better qualities we first search for an entry bigger than 'diagThreshold'.
// If no such entry is found, we search for an entry != 0.
double diagThreshold = 1.0e-6
for(int outerTry = 0; outerTry < 2; ++outerTry){
for(int itry = j; itry < n; ++itry) {
for (int i = 0; i <= j; i++) {
double s1 = 0;
for (int k = 0; k < i; k++)
s1 += U[k][j] * L[i][k];
U[i][j] = A[perm[i]][j] - s1;
}
if(Math.abs(U[j][j]) > diagThreshold)
break; // we found a good row
else if(itry + 1 < n){
// continue search with permutated matrix. Swap entries (itry+1) and j in perm.
int tmp = perm[itry + 1]
perm[itry + 1] = perm[j]
perm[j] = tmp
}
}
// after all rows below 'j' have been checked, we reduce the threshold to 0 and check again
diagThreshold = 0;
}
for (int i = j; i < n; i++) {
double s2 = 0;
for (int k = 0; k < j; k++)
s2 += U[k][j] * L[i][k];
L[i][j] = (A[perm[i]][j] - s2) / U[j][j];
}
}
}
@MethodInfo(noGUI=true)
void applyLU (double[] xOut, double[] b)
{
// solve forward
for (int i = 0; i < n; i++) {
xOut[i] = b[perm[i]];
for (int k = 0; k < i; k++) {
xOut[i] -= L[i][k]*xOut[k]
}
}
// solve backward
for (int i = n-1; i >= 0; i--) {
for (int k = i+1; k < n; k++) {
xOut[i] -= U[i][k]*xOut[k]
}
xOut[i] /= U[i][i];
}
}
@MethodInfo(noGUI=true)
void solve (double[] xOut, double[][] A, double[] b)
{
computeLU(A)
applyLU (xOut, b)
}
@MethodInfo(hide=false, valueName="x", valueOptions="serialization=false")
double[] solve (
@ParamInfo(name="A", options="serialization=false") double[][] A,
@ParamInfo(name="b", options="serialization=false") double[] b)
{
computeLU(A)
double[] xOut = new double[b.length]
applyLU (xOut, b)
return xOut
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment