1 | #
|
---|
2 | # Solves a*x = b for x, using LU decomposition.
|
---|
3 | #
|
---|
4 | module 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
|
---|
84 | end
|
---|