Module: LUSolve

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

Overview

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

Instance Method Summary collapse

Instance Method Details

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

Performs LU decomposition of the n by n matrix a.



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

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.



62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
# File 'lib/bigdecimal/ludcmp.rb', line 62

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