summaryrefslogtreecommitdiff
path: root/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
diff options
context:
space:
mode:
authorshigek <shigek@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2003-07-18 15:26:37 +0000
committershigek <shigek@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2003-07-18 15:26:37 +0000
commit1bb3071bde8344427a5fe3f5998dc491d2630d2b (patch)
tree8310d06d0fd319977dbc2b560d230c48f48f03da /ext/bigdecimal/lib/bigdecimal/ludcmp.rb
parente242cf9defb5442ef223535abe93399352cf139e (diff)
As discussed in ruby-dev ML:
lib directory moved. util.rb created instead of bigdecimal-rational.rb git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@4093 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
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