1 files changed, 84 insertions, 0 deletions
diff --git a/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/ludcmp.rb b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
new file mode 100644
index 0000000000..8f4888725e
--- /dev/null
+++ b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
@@ -0,0 +1,84 @@
+#
+# Solves a*x = b for x, using LU decomposition.
+#
+module LUSolve
+  # Performs LU decomposition of the n by n matrix a.
+  def ludecomp(a,n,zero=0,one=1)
+    prec = BigDecimal.limit(nil)
+    ps     = []
+    scales = []
+    for i in 0...n do  # pick up largest(abs. val.) element in each row.
+      ps <<= i
+      nrmrow  = zero
+      ixn = i*n
+      for j in 0...n do
+         biggst = a[ixn+j].abs
+         nrmrow = biggst if biggst>nrmrow
+      end
+      if nrmrow>zero then
+         scales <<= one.div(nrmrow,prec)
+      else 
+         raise "Singular matrix"
+      end
+    end
+    n1          = n - 1
+    for k in 0...n1 do # Gaussian elimination with partial pivoting.
+      biggst  = zero;
+      for i in k...n do
+         size = a[ps[i]*n+k].abs*scales[ps[i]]
+         if size>biggst then
+            biggst = size
+            pividx  = i
+         end
+      end
+      raise "Singular matrix" if biggst<=zero
+      if pividx!=k then
+        j = ps[k]
+        ps[k] = ps[pividx]
+        ps[pividx] = j
+      end
+      pivot   = a[ps[k]*n+k]
+      for i in (k+1)...n do
+        psin = ps[i]*n
+        a[psin+k] = mult = a[psin+k].div(pivot,prec)
+        if mult!=zero then
+           pskn = ps[k]*n
+           for j in (k+1)...n do
+             a[psin+j] -= mult.mult(a[pskn+j],prec)
+           end
+        end
+      end
+    end
+    raise "Singular matrix" if a[ps[n1]*n+n1] == zero
+    ps
+  end
+
+  # Solves a*x = b for x, using LU decomposition.
+  #
+  # 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.
+  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