Method: LUSolve.ludecomp

Defined in:
lib/bigdecimal/ludcmp.rb

.ludecomp(a, n, zero = 0, one = 1) ⇒ Object

Performs LU decomposition of the n by n matrix a.



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
 
# File 'lib/bigdecimal/ludcmp.rb', line 13

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