Support for the Ruby 2.4 series has ended. See here for reference.

### In Files

• matrix/lup_decomposition.rb

Quicksearch

# Matrix::LUPDecomposition

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.

### Attributes

pivots[R]

Returns the pivoting indices

### Public Class Methods

new(a) click to toggle source
```
# File matrix/lup_decomposition.rb, line 154
def initialize a
raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
# Use a "left-looking", dot-product, Crout/Doolittle algorithm.
@lu = a.to_a
@row_count = a.row_count
@column_count = a.column_count
@pivots = Array.new(@row_count)
@row_count.times do |i|
@pivots[i] = i
end
@pivot_sign = 1
lu_col_j = Array.new(@row_count)

# Outer loop.

@column_count.times do |j|

# Make a copy of the j-th column to localize references.

@row_count.times do |i|
lu_col_j[i] = @lu[i][j]
end

# Apply previous transformations.

@row_count.times do |i|
lu_row_i = @lu[i]

# Most of the time is spent in the following dot product.

kmax = [i, j].min
s = 0
kmax.times do |k|
s += lu_row_i[k]*lu_col_j[k]
end

lu_row_i[j] = lu_col_j[i] -= s
end

# Find pivot and exchange if necessary.

p = j
(j+1).upto(@row_count-1) do |i|
if (lu_col_j[i].abs > lu_col_j[p].abs)
p = i
end
end
if (p != j)
@column_count.times do |k|
t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t
end
k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k
@pivot_sign = -@pivot_sign
end

# Compute multipliers.

if (j < @row_count && @lu[j][j] != 0)
(j+1).upto(@row_count-1) do |i|
@lu[i][j] = @lu[i][j].quo(@lu[j][j])
end
end
end
end
```

### Public Instance Methods

det() click to toggle source

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

```
# File matrix/lup_decomposition.rb, line 79
def det
if (@row_count != @column_count)
Matrix.Raise Matrix::ErrDimensionMismatch
end
d = @pivot_sign
@column_count.times do |j|
d *= @lu[j][j]
end
d
end
```
Also aliased as: determinant
determinant() click to toggle source
Alias for: det
l() click to toggle source
```
# File matrix/lup_decomposition.rb, line 22
def l
Matrix.build(@row_count, [@column_count, @row_count].min) do |i, j|
if (i > j)
@lu[i][j]
elsif (i == j)
1
else
0
end
end
end
```
p() click to toggle source

Returns the permutation matrix `P`

```
# File matrix/lup_decomposition.rb, line 48
def p
rows = Array.new(@row_count){Array.new(@row_count, 0)}
@pivots.each_with_index{|p, i| rows[i][p] = 1}
Matrix.send :new, rows, @row_count
end
```
singular?() click to toggle source

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

```
# File matrix/lup_decomposition.rb, line 67
def singular?
@column_count.times do |j|
if (@lu[j][j] == 0)
return true
end
end
false
end
```
solve(b) click to toggle source

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

```
# File matrix/lup_decomposition.rb, line 95
def solve b
if (singular?)
Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."
end
if b.is_a? Matrix
if (b.row_count != @row_count)
Matrix.Raise Matrix::ErrDimensionMismatch
end

# Copy right hand side with pivoting
nx = b.column_count
m = @pivots.map{|row| b.row(row).to_a}

# Solve L*Y = P*b
@column_count.times do |k|
(k+1).upto(@column_count-1) do |i|
nx.times do |j|
m[i][j] -= m[k][j]*@lu[i][k]
end
end
end
# Solve U*m = Y
(@column_count-1).downto(0) do |k|
nx.times do |j|
m[k][j] = m[k][j].quo(@lu[k][k])
end
k.times do |i|
nx.times do |j|
m[i][j] -= m[k][j]*@lu[i][k]
end
end
end
Matrix.send :new, m, nx
else # same algorithm, specialized for simpler case of a vector
b = convert_to_array(b)
if (b.size != @row_count)
Matrix.Raise Matrix::ErrDimensionMismatch
end

# Copy right hand side with pivoting
m = b.values_at(*@pivots)

# Solve L*Y = P*b
@column_count.times do |k|
(k+1).upto(@column_count-1) do |i|
m[i] -= m[k]*@lu[i][k]
end
end
# Solve U*m = Y
(@column_count-1).downto(0) do |k|
m[k] = m[k].quo(@lu[k][k])
k.times do |i|
m[i] -= m[k]*@lu[i][k]
end
end
Vector.elements(m, false)
end
end
```
to_a() click to toggle source
Alias for: to_ary
to_ary() click to toggle source

Returns `L`, `U`, `P` in an array

```
# File matrix/lup_decomposition.rb, line 56
def to_ary
[l, u, p]
end
```
Also aliased as: to_a
u() click to toggle source

Returns the upper triangular factor `U`

```
# File matrix/lup_decomposition.rb, line 36
def u
Matrix.build([@column_count, @row_count].min, @column_count) do |i, j|
if (i <= j)
@lu[i][j]
else
0
end
end
end
```