summaryrefslogtreecommitdiff log msg author committer range
path: root/ext/bigdecimal/lib/ludcmp.rb
blob: c36f0dea5b6b4e6b32609a8c233bfa54d62144da (plain)
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 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