module LUSolve

def lusolve(a,b,ps,zero=0.0)

during LU decomposition.
ps is the pivot, a vector which indicates the permutation of rows performed

a is a matrix, b is a constant vector, x is the solution vector.

Solves a*x = b for x, using LU decomposition. 
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