Class Matrix::LUPDecomposition
In: lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb
Parent: Object

For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit lower triangular matrix L, an n-by-n upper triangular matrix U, and a m-by-m permutation matrix P so that L*U = P*A. If m < n, then L is m-by-m and U is m-by-n.

The LUP decomposition with pivoting always exists, even if the matrix is singular, so the constructor will never fail. The primary use of the LU decomposition is in the solution of square systems of simultaneous linear equations. This will fail if singular? returns true.

Methods

det   determinant   l   new   p   singular?   solve   to_a   to_ary   u  

Included Modules

Matrix::ConversionHelper

Attributes

pivots  [R]  Returns the pivoting indices

Public Class methods

[Source]

     # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 153
153:     def initialize a
154:       raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
155:       # Use a "left-looking", dot-product, Crout/Doolittle algorithm.
156:       @lu = a.to_a
157:       @row_size = a.row_size
158:       @col_size = a.column_size
159:       @pivots = Array.new(@row_size)
160:       @row_size.times do |i|
161:          @pivots[i] = i
162:       end
163:       @pivot_sign = 1
164:       lu_col_j = Array.new(@row_size)
165: 
166:       # Outer loop.
167: 
168:       @col_size.times do |j|
169: 
170:         # Make a copy of the j-th column to localize references.
171: 
172:         @row_size.times do |i|
173:           lu_col_j[i] = @lu[i][j]
174:         end
175: 
176:         # Apply previous transformations.
177: 
178:         @row_size.times do |i|
179:           lu_row_i = @lu[i]
180: 
181:           # Most of the time is spent in the following dot product.
182: 
183:           kmax = [i, j].min
184:           s = 0
185:           kmax.times do |k|
186:             s += lu_row_i[k]*lu_col_j[k]
187:           end
188: 
189:           lu_row_i[j] = lu_col_j[i] -= s
190:         end
191: 
192:         # Find pivot and exchange if necessary.
193: 
194:         p = j
195:         (j+1).upto(@row_size-1) do |i|
196:           if (lu_col_j[i].abs > lu_col_j[p].abs)
197:             p = i
198:           end
199:         end
200:         if (p != j)
201:           @col_size.times do |k|
202:             t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t
203:           end
204:           k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k
205:           @pivot_sign = -@pivot_sign
206:         end
207: 
208:         # Compute multipliers.
209: 
210:         if (j < @row_size && @lu[j][j] != 0)
211:           (j+1).upto(@row_size-1) do |i|
212:             @lu[i][j] = @lu[i][j].quo(@lu[j][j])
213:           end
214:         end
215:       end
216:     end

Public Instance methods

Returns the determinant of A, calculated efficiently from the factorization.

[Source]

    # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 78
78:     def det
79:       if (@row_size != @col_size)
80:         Matrix.Raise Matrix::ErrDimensionMismatch unless square?
81:       end
82:       d = @pivot_sign
83:       @col_size.times do |j|
84:         d *= @lu[j][j]
85:       end
86:       d
87:     end
determinant()

Alias for det

[Source]

    # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 21
21:     def l
22:       Matrix.build(@row_size, @col_size) do |i, j|
23:         if (i > j)
24:           @lu[i][j]
25:         elsif (i == j)
26:           1
27:         else
28:           0
29:         end
30:       end
31:     end

Returns the permutation matrix P

[Source]

    # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 47
47:     def p
48:       rows = Array.new(@row_size){Array.new(@col_size, 0)}
49:       @pivots.each_with_index{|p, i| rows[i][p] = 1}
50:       Matrix.send :new, rows, @col_size
51:     end

Returns true if U, and hence A, is singular.

[Source]

    # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 66
66:     def singular? ()
67:       @col_size.times do |j|
68:         if (@lu[j][j] == 0)
69:           return true
70:         end
71:       end
72:       false
73:     end

Returns m so that A*m = b, or equivalently so that L*U*m = P*b b can be a Matrix or a Vector

[Source]

     # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 94
 94:     def solve b
 95:       if (singular?)
 96:         Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."
 97:       end
 98:       if b.is_a? Matrix
 99:         if (b.row_size != @row_size)
100:           Matrix.Raise Matrix::ErrDimensionMismatch
101:         end
102: 
103:         # Copy right hand side with pivoting
104:         nx = b.column_size
105:         m = @pivots.map{|row| b.row(row).to_a}
106: 
107:         # Solve L*Y = P*b
108:         @col_size.times do |k|
109:           (k+1).upto(@col_size-1) do |i|
110:             nx.times do |j|
111:               m[i][j] -= m[k][j]*@lu[i][k]
112:             end
113:           end
114:         end
115:         # Solve U*m = Y
116:         (@col_size-1).downto(0) do |k|
117:           nx.times do |j|
118:             m[k][j] = m[k][j].quo(@lu[k][k])
119:           end
120:           k.times do |i|
121:             nx.times do |j|
122:               m[i][j] -= m[k][j]*@lu[i][k]
123:             end
124:           end
125:         end
126:         Matrix.send :new, m, nx
127:       else # same algorithm, specialized for simpler case of a vector
128:         b = convert_to_array(b)
129:         if (b.size != @row_size)
130:           Matrix.Raise Matrix::ErrDimensionMismatch
131:         end
132: 
133:         # Copy right hand side with pivoting
134:         m = b.values_at(*@pivots)
135: 
136:         # Solve L*Y = P*b
137:         @col_size.times do |k|
138:           (k+1).upto(@col_size-1) do |i|
139:             m[i] -= m[k]*@lu[i][k]
140:           end
141:         end
142:         # Solve U*m = Y
143:         (@col_size-1).downto(0) do |k|
144:           m[k] = m[k].quo(@lu[k][k])
145:           k.times do |i|
146:             m[i] -= m[k]*@lu[i][k]
147:           end
148:         end
149:         Vector.elements(m, false)
150:       end
151:     end
to_a()

Alias for to_ary

Returns L, U, P in an array

[Source]

    # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 55
55:     def to_ary
56:       [l, u, p]
57:     end

Returns the upper triangular factor U

[Source]

    # File lib/backports/1.9.2/stdlib/matrix/lup_decomposition.rb, line 35
35:     def u
36:       Matrix.build(@col_size, @col_size) do |i, j|
37:         if (i <= j)
38:           @lu[i][j]
39:         else
40:           0
41:         end
42:       end
43:     end

[Validate]