summaryrefslogtreecommitdiff
path: root/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
diff options
context:
space:
mode:
Diffstat (limited to 'ext/bigdecimal/lib/bigdecimal/ludcmp.rb')
-rw-r--r--ext/bigdecimal/lib/bigdecimal/ludcmp.rb75
1 files changed, 75 insertions, 0 deletions
diff --git a/ext/bigdecimal/lib/bigdecimal/ludcmp.rb b/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
new file mode 100644
index 00000000000..c36f0dea5b6
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
@@ -0,0 +1,75 @@
+#
+# ludcmp.rb
+#
+module LUSolve
+ def ludecomp(a,n,zero=0.0,one=1.0)
+ 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/nrmrow
+ 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]/pivot
+ if mult!=zero then
+ pskn = ps[k]*n
+ for j in (k+1)...n do
+ a[psin+j] -= mult*a[pskn+j]
+ end
+ end
+ end
+ end
+ raise "Singular matrix" if a[ps[n1]*n+n1] == zero
+ ps
+ end
+
+ def lusolve(a,b,ps,zero=0.0)
+ 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]*x[j] + 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]*x[j] + dot
+ end
+ x[i] = (x[i]-dot)/a[psin+i]
+ end
+ x
+ end
+end