summaryrefslogtreecommitdiff
path: root/trunk/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
diff options
context:
space:
mode:
Diffstat (limited to 'trunk/ext/bigdecimal/lib/bigdecimal/ludcmp.rb')
-rw-r--r--trunk/ext/bigdecimal/lib/bigdecimal/ludcmp.rb84
1 files changed, 0 insertions, 84 deletions
diff --git a/trunk/ext/bigdecimal/lib/bigdecimal/ludcmp.rb b/trunk/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
deleted file mode 100644
index 8f4888725e..0000000000
--- a/trunk/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
+++ /dev/null
@@ -1,84 +0,0 @@
-#
-# 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