def initialize a
raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
@lu = a.to_a
@row_size = a.row_size
@col_size = a.column_size
@pivots = Array.new(@row_size)
@row_size.times do |i|
@pivots[i] = i
end
@pivot_sign = 1
lu_col_j = Array.new(@row_size)
@col_size.times do |j|
@row_size.times do |i|
lu_col_j[i] = @lu[i][j]
end
@row_size.times do |i|
lu_row_i = @lu[i]
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
p = j
(j+1).upto(@row_size-1) do |i|
if (lu_col_j[i].abs > lu_col_j[p].abs)
p = i
end
end
if (p != j)
@col_size.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
if (j < @row_size && @lu[j][j] != 0)
(j+1).upto(@row_size-1) do |i|
@lu[i][j] = @lu[i][j].quo(@lu[j][j])
end
end
end
end