# frozen_string_literal: false

[0] Fix | Delete

class Matrix

[1] Fix | Delete

# Adapted from JAMA: http://math.nist.gov/javanumerics/jama/

[2] Fix | Delete

[3] Fix | Delete

#

[4] Fix | Delete

# For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n

[5] Fix | Delete

# unit lower triangular matrix L, an n-by-n upper triangular matrix U,

[6] Fix | Delete

# and a m-by-m permutation matrix P so that L*U = P*A.

[7] Fix | Delete

# If m < n, then L is m-by-m and U is m-by-n.

[8] Fix | Delete

#

[9] Fix | Delete

# The LUP decomposition with pivoting always exists, even if the matrix is

[10] Fix | Delete

# singular, so the constructor will never fail. The primary use of the

[11] Fix | Delete

# LU decomposition is in the solution of square systems of simultaneous

[12] Fix | Delete

# linear equations. This will fail if singular? returns true.

[13] Fix | Delete

#

[14] Fix | Delete

[15] Fix | Delete

class LUPDecomposition

[16] Fix | Delete

# Returns the lower triangular factor +L+

[17] Fix | Delete

[18] Fix | Delete

include Matrix::ConversionHelper

[19] Fix | Delete

[20] Fix | Delete

def l

[21] Fix | Delete

Matrix.build(@row_count, [@column_count, @row_count].min) do |i, j|

[22] Fix | Delete

if (i > j)

[23] Fix | Delete

@lu[i][j]

[24] Fix | Delete

elsif (i == j)

[25] Fix | Delete

1

[26] Fix | Delete

else

[27] Fix | Delete

0

[28] Fix | Delete

end

[29] Fix | Delete

end

[30] Fix | Delete

end

[31] Fix | Delete

[32] Fix | Delete

# Returns the upper triangular factor +U+

[33] Fix | Delete

[34] Fix | Delete

def u

[35] Fix | Delete

Matrix.build([@column_count, @row_count].min, @column_count) do |i, j|

[36] Fix | Delete

if (i <= j)

[37] Fix | Delete

@lu[i][j]

[38] Fix | Delete

else

[39] Fix | Delete

0

[40] Fix | Delete

end

[41] Fix | Delete

end

[42] Fix | Delete

end

[43] Fix | Delete

[44] Fix | Delete

# Returns the permutation matrix +P+

[45] Fix | Delete

[46] Fix | Delete

def p

[47] Fix | Delete

rows = Array.new(@row_count){Array.new(@row_count, 0)}

[48] Fix | Delete

@pivots.each_with_index{|p, i| rows[i][p] = 1}

[49] Fix | Delete

Matrix.send :new, rows, @row_count

[50] Fix | Delete

end

[51] Fix | Delete

[52] Fix | Delete

# Returns +L+, +U+, +P+ in an array

[53] Fix | Delete

[54] Fix | Delete

def to_ary

[55] Fix | Delete

[l, u, p]

[56] Fix | Delete

end

[57] Fix | Delete

alias_method :to_a, :to_ary

[58] Fix | Delete

[59] Fix | Delete

# Returns the pivoting indices

[60] Fix | Delete

[61] Fix | Delete

attr_reader :pivots

[62] Fix | Delete

[63] Fix | Delete

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

[64] Fix | Delete

[65] Fix | Delete

def singular?

[66] Fix | Delete

@column_count.times do |j|

[67] Fix | Delete

if (@lu[j][j] == 0)

[68] Fix | Delete

return true

[69] Fix | Delete

end

[70] Fix | Delete

end

[71] Fix | Delete

false

[72] Fix | Delete

end

[73] Fix | Delete

[74] Fix | Delete

# Returns the determinant of +A+, calculated efficiently

[75] Fix | Delete

# from the factorization.

[76] Fix | Delete

[77] Fix | Delete

def det

[78] Fix | Delete

if (@row_count != @column_count)

[79] Fix | Delete

Matrix.Raise Matrix::ErrDimensionMismatch

[80] Fix | Delete

end

[81] Fix | Delete

d = @pivot_sign

[82] Fix | Delete

@column_count.times do |j|

[83] Fix | Delete

d *= @lu[j][j]

[84] Fix | Delete

end

[85] Fix | Delete

d

[86] Fix | Delete

end

[87] Fix | Delete

alias_method :determinant, :det

[88] Fix | Delete

[89] Fix | Delete

# Returns +m+ so that <tt>A*m = b</tt>,

[90] Fix | Delete

# or equivalently so that <tt>L*U*m = P*b</tt>

[91] Fix | Delete

# +b+ can be a Matrix or a Vector

[92] Fix | Delete

[93] Fix | Delete

def solve b

[94] Fix | Delete

if (singular?)

[95] Fix | Delete

Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."

[96] Fix | Delete

end

[97] Fix | Delete

if b.is_a? Matrix

[98] Fix | Delete

if (b.row_count != @row_count)

[99] Fix | Delete

Matrix.Raise Matrix::ErrDimensionMismatch

[100] Fix | Delete

end

[101] Fix | Delete

[102] Fix | Delete

# Copy right hand side with pivoting

[103] Fix | Delete

nx = b.column_count

[104] Fix | Delete

m = @pivots.map{|row| b.row(row).to_a}

[105] Fix | Delete

[106] Fix | Delete

# Solve L*Y = P*b

[107] Fix | Delete

@column_count.times do |k|

[108] Fix | Delete

(k+1).upto(@column_count-1) do |i|

[109] Fix | Delete

nx.times do |j|

[110] Fix | Delete

m[i][j] -= m[k][j]*@lu[i][k]

[111] Fix | Delete

end

[112] Fix | Delete

end

[113] Fix | Delete

end

[114] Fix | Delete

# Solve U*m = Y

[115] Fix | Delete

(@column_count-1).downto(0) do |k|

[116] Fix | Delete

nx.times do |j|

[117] Fix | Delete

m[k][j] = m[k][j].quo(@lu[k][k])

[118] Fix | Delete

end

[119] Fix | Delete

k.times do |i|

[120] Fix | Delete

nx.times do |j|

[121] Fix | Delete

m[i][j] -= m[k][j]*@lu[i][k]

[122] Fix | Delete

end

[123] Fix | Delete

end

[124] Fix | Delete

end

[125] Fix | Delete

Matrix.send :new, m, nx

[126] Fix | Delete

else # same algorithm, specialized for simpler case of a vector

[127] Fix | Delete

b = convert_to_array(b)

[128] Fix | Delete

if (b.size != @row_count)

[129] Fix | Delete

Matrix.Raise Matrix::ErrDimensionMismatch

[130] Fix | Delete

end

[131] Fix | Delete

[132] Fix | Delete

# Copy right hand side with pivoting

[133] Fix | Delete

m = b.values_at(*@pivots)

[134] Fix | Delete

[135] Fix | Delete

# Solve L*Y = P*b

[136] Fix | Delete

@column_count.times do |k|

[137] Fix | Delete

(k+1).upto(@column_count-1) do |i|

[138] Fix | Delete

m[i] -= m[k]*@lu[i][k]

[139] Fix | Delete

end

[140] Fix | Delete

end

[141] Fix | Delete

# Solve U*m = Y

[142] Fix | Delete

(@column_count-1).downto(0) do |k|

[143] Fix | Delete

m[k] = m[k].quo(@lu[k][k])

[144] Fix | Delete

k.times do |i|

[145] Fix | Delete

m[i] -= m[k]*@lu[i][k]

[146] Fix | Delete

end

[147] Fix | Delete

end

[148] Fix | Delete

Vector.elements(m, false)

[149] Fix | Delete

end

[150] Fix | Delete

end

[151] Fix | Delete

[152] Fix | Delete

def initialize a

[153] Fix | Delete

raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)

[154] Fix | Delete

# Use a "left-looking", dot-product, Crout/Doolittle algorithm.

[155] Fix | Delete

@lu = a.to_a

[156] Fix | Delete

@row_count = a.row_count

[157] Fix | Delete

@column_count = a.column_count

[158] Fix | Delete

@pivots = Array.new(@row_count)

[159] Fix | Delete

@row_count.times do |i|

[160] Fix | Delete

@pivots[i] = i

[161] Fix | Delete

end

[162] Fix | Delete

@pivot_sign = 1

[163] Fix | Delete

lu_col_j = Array.new(@row_count)

[164] Fix | Delete

[165] Fix | Delete

# Outer loop.

[166] Fix | Delete

[167] Fix | Delete

@column_count.times do |j|

[168] Fix | Delete

[169] Fix | Delete

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

[170] Fix | Delete

[171] Fix | Delete

@row_count.times do |i|

[172] Fix | Delete

lu_col_j[i] = @lu[i][j]

[173] Fix | Delete

end

[174] Fix | Delete

[175] Fix | Delete

# Apply previous transformations.

[176] Fix | Delete

[177] Fix | Delete

@row_count.times do |i|

[178] Fix | Delete

lu_row_i = @lu[i]

[179] Fix | Delete

[180] Fix | Delete

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

[181] Fix | Delete

[182] Fix | Delete

kmax = [i, j].min

[183] Fix | Delete

s = 0

[184] Fix | Delete

kmax.times do |k|

[185] Fix | Delete

s += lu_row_i[k]*lu_col_j[k]

[186] Fix | Delete

end

[187] Fix | Delete

[188] Fix | Delete

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

[189] Fix | Delete

end

[190] Fix | Delete

[191] Fix | Delete

# Find pivot and exchange if necessary.

[192] Fix | Delete

[193] Fix | Delete

p = j

[194] Fix | Delete

(j+1).upto(@row_count-1) do |i|

[195] Fix | Delete

if (lu_col_j[i].abs > lu_col_j[p].abs)

[196] Fix | Delete

p = i

[197] Fix | Delete

end

[198] Fix | Delete

end

[199] Fix | Delete

if (p != j)

[200] Fix | Delete

@column_count.times do |k|

[201] Fix | Delete

t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t

[202] Fix | Delete

end

[203] Fix | Delete

k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k

[204] Fix | Delete

@pivot_sign = -@pivot_sign

[205] Fix | Delete

end

[206] Fix | Delete

[207] Fix | Delete

# Compute multipliers.

[208] Fix | Delete

[209] Fix | Delete

if (j < @row_count && @lu[j][j] != 0)

[210] Fix | Delete

(j+1).upto(@row_count-1) do |i|

[211] Fix | Delete

@lu[i][j] = @lu[i][j].quo(@lu[j][j])

[212] Fix | Delete

end

[213] Fix | Delete

end

[214] Fix | Delete

end

[215] Fix | Delete

end

[216] Fix | Delete

end

[217] Fix | Delete

end

[218] Fix | Delete

[219] Fix | Delete