/tools/Ruby/lib/ruby/1.8/bigdecimal/ludcmp.rb

 1#
 2# Solves a*x = b for x, using LU decomposition.
 3#
 4module LUSolve
 5  # Performs LU decomposition of the n by n matrix a.
 6  def ludecomp(a,n,zero=0,one=1)
 7    prec = BigDecimal.limit(nil)
 8    ps     = []
 9    scales = []
10    for i in 0...n do  # pick up largest(abs. val.) element in each row.
11      ps <<= i
12      nrmrow  = zero
13      ixn = i*n
14      for j in 0...n do
15         biggst = a[ixn+j].abs
16         nrmrow = biggst if biggst>nrmrow
17      end
18      if nrmrow>zero then
19         scales <<= one.div(nrmrow,prec)
20      else 
21         raise "Singular matrix"
22      end
23    end
24    n1          = n - 1
25    for k in 0...n1 do # Gaussian elimination with partial pivoting.
26      biggst  = zero;
27      for i in k...n do
28         size = a[ps[i]*n+k].abs*scales[ps[i]]
29         if size>biggst then
30            biggst = size
31            pividx  = i
32         end
33      end
34      raise "Singular matrix" if biggst<=zero
35      if pividx!=k then
36        j = ps[k]
37        ps[k] = ps[pividx]
38        ps[pividx] = j
39      end
40      pivot   = a[ps[k]*n+k]
41      for i in (k+1)...n do
42        psin = ps[i]*n
43        a[psin+k] = mult = a[psin+k].div(pivot,prec)
44        if mult!=zero then
45           pskn = ps[k]*n
46           for j in (k+1)...n do
47             a[psin+j] -= mult.mult(a[pskn+j],prec)
48           end
49        end
50      end
51    end
52    raise "Singular matrix" if a[ps[n1]*n+n1] == zero
53    ps
54  end
55
56  # Solves a*x = b for x, using LU decomposition.
57  #
58  # a is a matrix, b is a constant vector, x is the solution vector.
59  #
60  # ps is the pivot, a vector which indicates the permutation of rows performed
61  # during LU decomposition.
62  def lusolve(a,b,ps,zero=0.0)
63    prec = BigDecimal.limit(nil)
64    n = ps.size
65    x = []
66    for i in 0...n do
67      dot = zero
68      psin = ps[i]*n
69      for j in 0...i do
70        dot = a[psin+j].mult(x[j],prec) + dot
71      end
72      x <<= b[ps[i]] - dot
73    end
74    (n-1).downto(0) do |i|
75       dot = zero
76       psin = ps[i]*n
77       for j in (i+1)...n do
78         dot = a[psin+j].mult(x[j],prec) + dot
79       end
80       x[i]  = (x[i]-dot).div(a[psin+i],prec)
81    end
82    x
83  end
84end