Module: LUSolve

Included in:
Newton
Defined in:
lib/bigdecimal/ludcmp.rb

Overview

Solves a*x = b for x, using LU decomposition.

Class Method Summary collapse

Class Method Details

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

Performs LU decomposition of the n by n matrix a.



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

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

.lusolve(a, b, ps, zero = 0.0) ⇒ Object

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.



66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
 
# File 'lib/bigdecimal/ludcmp.rb', line 66

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