a is a matrix, b is a constant vector, x is the solution vector.
ps is the pivot, a vector which indicates the permutation of rows performed
during LU decomposition.
# File ext/bigdecimal/lib/bigdecimal/ludcmp.rb, line 67
def lusolve(a,b,ps,zero=0.0)
prec = BigDecimal.limit(nil)
n = ps.size
x = []
for i in 0...n do
dot = zero
psin = ps[i]*n
for j in 0...i do
dot = a[psin+j].mult(x[j],prec) + dot
end
x <<= b[ps[i]] - dot
end
(n-1).downto(0) do |i|
dot = zero
psin = ps[i]*n
for j in (i+1)...n do
dot = a[psin+j].mult(x[j],prec) + dot
end
x[i] = (x[i]-dot).div(a[psin+i],prec)
end
x
end
end