kotovalexarian-likes-github/ruby--ruby
1
0
Fork 0
mirror of https://github.com/ruby/ruby.git synced 2022-11-09 12:17:21 -05:00
Code Releases Activity

* lib/matrix: alias {row|column}_size to {row|column}_count and use the latter.

Browse source 
[Bug #7369] [ruby-core:49409]

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@38300 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
This commit is contained in:
marcandre 2012-12-10 16:53:57 +00:00
parent 4c02cff191
commit 8aac5f48fc
4 changed files with 124 additions and 116 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
ChangeLog
View file

@ -1,3 +1,9 @@
Tue Dec 11 01:53:37 2012 Marc-Andre Lafortune <ruby-core@marc-andre.ca>
* lib/matrix: alias {row|column}_size to {row|column}_count and use
the latter.
[Bug #7369] [ruby-core:49409]
Tue Dec 11 00:26:58 2012 Shugo Maeda <shugo@ruby-lang.org> Tue Dec 11 00:26:58 2012 Shugo Maeda <shugo@ruby-lang.org>
* fix the behavior when a module is included into a refinement. * fix the behavior when a module is included into a refinement.

176
lib/matrix.rb
View file

@ -36,7 +36,7 @@ end
# * Matrix.[](*rows) # * Matrix.[](*rows)
# * Matrix.rows(rows, copy = true) # * Matrix.rows(rows, copy = true)
# * Matrix.columns(columns) # * Matrix.columns(columns)
# * Matrix.build(row_size, column_size, &block) # * Matrix.build(row_count, column_count, &block)
# * Matrix.diagonal(*values) # * Matrix.diagonal(*values)
# * Matrix.scalar(n, value) # * Matrix.scalar(n, value)
# * Matrix.identity(n) # * Matrix.identity(n)
@ -48,8 +48,8 @@ end
# #
# To access Matrix elements/columns/rows/submatrices/properties: # To access Matrix elements/columns/rows/submatrices/properties:
# * #[](i, j) # * #[](i, j)
# * #row_size # * #row_count (row_size)
# * #column_size # * #column_count (column_size)
# * #row(i) # * #row(i)
# * #column(j) # * #column(j)
# * #collect # * #collect
@ -172,7 +172,7 @@ class Matrix
end end
# #
# Creates a matrix of size +row_size+ x +column_size+. # Creates a matrix of size +row_count+ x +column_count+.
# It fills the values by calling the given block, # It fills the values by calling the given block,
# passing the current row and column. # passing the current row and column.
# Returns an enumerator if no block is given. # Returns an enumerator if no block is given.
@ -182,17 +182,17 @@ class Matrix
# m = Matrix.build(3) { rand } # m = Matrix.build(3) { rand }
# => a 3x3 matrix with random elements # => a 3x3 matrix with random elements
# #
def Matrix.build(row_size, column_size = row_size) def Matrix.build(row_count, column_count = row_count)
row_size = CoercionHelper.coerce_to_int(row_size) row_count = CoercionHelper.coerce_to_int(row_count)
column_size = CoercionHelper.coerce_to_int(column_size) column_count = CoercionHelper.coerce_to_int(column_count)
raise ArgumentError if row_size < 0 || column_size < 0 raise ArgumentError if row_count < 0 || column_count < 0
return to_enum :build, row_size, column_size unless block_given? return to_enum :build, row_count, column_count unless block_given?
rows = Array.new(row_size) do |i| rows = Array.new(row_count) do |i|
Array.new(column_size) do |j| Array.new(column_count) do |j|
yield i, j yield i, j
end end
end end
new rows, column_size new rows, column_count
end end
# #
@ -243,9 +243,9 @@ class Matrix
# => 0 0 # => 0 0
# 0 0 # 0 0
# #
def Matrix.zero(row_size, column_size = row_size) def Matrix.zero(row_count, column_count = row_count)
rows = Array.new(row_size){Array.new(column_size, 0)} rows = Array.new(row_count){Array.new(column_count, 0)}
new rows, column_size new rows, column_count
end end
# #
@ -273,8 +273,8 @@ class Matrix
end end
# #
# Creates a empty matrix of +row_size+ x +column_size+. # Creates a empty matrix of +row_count+ x +column_count+.
# At least one of +row_size+ or +column_size+ must be 0. # At least one of +row_count+ or +column_count+ must be 0.
# #
# m = Matrix.empty(2, 0) # m = Matrix.empty(2, 0)
# m == Matrix[ [], [] ] # m == Matrix[ [], [] ]
@ -285,26 +285,26 @@ class Matrix
# m * n # m * n
# => Matrix[[0, 0, 0], [0, 0, 0]] # => Matrix[[0, 0, 0], [0, 0, 0]]
# #
def Matrix.empty(row_size = 0, column_size = 0) def Matrix.empty(row_count = 0, column_count = 0)
Matrix.Raise ArgumentError, "One size must be 0" if column_size != 0 && row_size != 0 Matrix.Raise ArgumentError, "One size must be 0" if column_count != 0 && row_count != 0
Matrix.Raise ArgumentError, "Negative size" if column_size < 0 || row_size < 0 Matrix.Raise ArgumentError, "Negative size" if column_count < 0 || row_count < 0
new([[]]*row_size, column_size) new([[]]*row_count, column_count)
end end
# #
# Matrix.new is private; use Matrix.rows, columns, [], etc... to create. # Matrix.new is private; use Matrix.rows, columns, [], etc... to create.
# #
def initialize(rows, column_size = rows[0].size) def initialize(rows, column_count = rows[0].size)
# No checking is done at this point. rows must be an Array of Arrays. # No checking is done at this point. rows must be an Array of Arrays.
# column_size must be the size of the first row, if there is one, # column_count must be the size of the first row, if there is one,
# otherwise it *must* be specified and can be any integer >= 0 # otherwise it *must* be specified and can be any integer >= 0
@rows = rows @rows = rows
@column_size = column_size @column_count = column_count
end end
def new_matrix(rows, column_size = rows[0].size) # :nodoc: def new_matrix(rows, column_count = rows[0].size) # :nodoc:
self.class.send(:new, rows, column_size) # bypass privacy of Matrix.new self.class.send(:new, rows, column_count) # bypass privacy of Matrix.new
end end
private :new_matrix private :new_matrix
@ -327,14 +327,16 @@ class Matrix
# #
# Returns the number of rows. # Returns the number of rows.
# #
def row_size def row_count
@rows.size @rows.size
end end
alias_method :row_size, :row_count
# #
# Returns the number of columns. # Returns the number of columns.
# #
attr_reader :column_size attr_reader :column_count
alias_method :column_size, :column_count
# #
# Returns row vector number +i+ of the matrix as a Vector (starting at 0 like # Returns row vector number +i+ of the matrix as a Vector (starting at 0 like
@ -356,14 +358,14 @@ class Matrix
# #
def column(j) # :yield: e def column(j) # :yield: e
if block_given? if block_given?
return self if j >= column_size || j < -column_size return self if j >= column_count || j < -column_count
row_size.times do |i| row_count.times do |i|
yield @rows[i][j] yield @rows[i][j]
end end
self self
else else
return nil if j >= column_size || j < -column_size return nil if j >= column_count || j < -column_count
col = Array.new(row_size) {|i| col = Array.new(row_count) {|i|
@rows[i][j] @rows[i][j]
} }
Vector.elements(col, false) Vector.elements(col, false)
@ -380,7 +382,7 @@ class Matrix
def collect(&block) # :yield: e def collect(&block) # :yield: e
return to_enum(:collect) unless block_given? return to_enum(:collect) unless block_given?
rows = @rows.collect{|row| row.collect(&block)} rows = @rows.collect{|row| row.collect(&block)}
new_matrix rows, column_size new_matrix rows, column_count
end end
alias map collect alias map collect
@ -402,7 +404,7 @@ class Matrix
# #
def each(which = :all) # :yield: e def each(which = :all) # :yield: e
return to_enum :each, which unless block_given? return to_enum :each, which unless block_given?
last = column_size - 1 last = column_count - 1
case which case which
when :all when :all
block = Proc.new block = Proc.new
@ -415,7 +417,7 @@ class Matrix
end end
when :off_diagonal when :off_diagonal
@rows.each_with_index do |row, row_index| @rows.each_with_index do |row, row_index|
column_size.times do |col_index| column_count.times do |col_index|
yield row[col_index] unless row_index == col_index yield row[col_index] unless row_index == col_index
end end
end end
@ -427,7 +429,7 @@ class Matrix
end end
when :strict_lower when :strict_lower
@rows.each_with_index do |row, row_index| @rows.each_with_index do |row, row_index|
[row_index, column_size].min.times do |col_index| [row_index, column_count].min.times do |col_index|
yield row[col_index] yield row[col_index]
end end
end end
@ -463,7 +465,7 @@ class Matrix
# #
def each_with_index(which = :all) # :yield: e, row, column def each_with_index(which = :all) # :yield: e, row, column
return to_enum :each_with_index, which unless block_given? return to_enum :each_with_index, which unless block_given?
last = column_size - 1 last = column_count - 1
case which case which
when :all when :all
@rows.each_with_index do |row, row_index| @rows.each_with_index do |row, row_index|
@ -477,7 +479,7 @@ class Matrix
end end
when :off_diagonal when :off_diagonal
@rows.each_with_index do |row, row_index| @rows.each_with_index do |row, row_index|
column_size.times do |col_index| column_count.times do |col_index|
yield row[col_index], row_index, col_index unless row_index == col_index yield row[col_index], row_index, col_index unless row_index == col_index
end end
end end
@ -489,7 +491,7 @@ class Matrix
end end
when :strict_lower when :strict_lower
@rows.each_with_index do |row, row_index| @rows.each_with_index do |row, row_index|
[row_index, column_size].min.times do |col_index| [row_index, column_count].min.times do |col_index|
yield row[col_index], row_index, col_index yield row[col_index], row_index, col_index
end end
end end
@ -552,39 +554,39 @@ class Matrix
# #
# Like Array#[], negative indices count backward from the end of the # Like Array#[], negative indices count backward from the end of the
# row or column (-1 is the last element). Returns nil if the starting # row or column (-1 is the last element). Returns nil if the starting
# row or column is greater than row_size or column_size respectively. # row or column is greater than row_count or column_count respectively.
# #
def minor(*param) def minor(*param)
case param.size case param.size
when 2 when 2
row_range, col_range = param row_range, col_range = param
from_row = row_range.first from_row = row_range.first
from_row += row_size if from_row < 0 from_row += row_count if from_row < 0
to_row = row_range.end to_row = row_range.end
to_row += row_size if to_row < 0 to_row += row_count if to_row < 0
to_row += 1 unless row_range.exclude_end? to_row += 1 unless row_range.exclude_end?
size_row = to_row - from_row size_row = to_row - from_row
from_col = col_range.first from_col = col_range.first
from_col += column_size if from_col < 0 from_col += column_count if from_col < 0
to_col = col_range.end to_col = col_range.end
to_col += column_size if to_col < 0 to_col += column_count if to_col < 0
to_col += 1 unless col_range.exclude_end? to_col += 1 unless col_range.exclude_end?
size_col = to_col - from_col size_col = to_col - from_col
when 4 when 4
from_row, size_row, from_col, size_col = param from_row, size_row, from_col, size_col = param
return nil if size_row < 0 || size_col < 0 return nil if size_row < 0 || size_col < 0
from_row += row_size if from_row < 0 from_row += row_count if from_row < 0
from_col += column_size if from_col < 0 from_col += column_count if from_col < 0
else else
Matrix.Raise ArgumentError, param.inspect Matrix.Raise ArgumentError, param.inspect
end end
return nil if from_row > row_size || from_col > column_size || from_row < 0 || from_col < 0 return nil if from_row > row_count || from_col > column_count || from_row < 0 || from_col < 0
rows = @rows[from_row, size_row].collect{|row| rows = @rows[from_row, size_row].collect{|row|
row[from_col, size_col] row[from_col, size_col]
} }
new_matrix rows, [column_size - from_col, size_col].min new_matrix rows, [column_count - from_col, size_col].min
end end
#-- #--
@ -605,7 +607,7 @@ class Matrix
# or the number of columns is 0. # or the number of columns is 0.
# #
def empty? def empty?
column_size == 0 || row_size == 0 column_count == 0 || row_count == 0
end end
# #
@ -651,9 +653,9 @@ class Matrix
def orthogonal? def orthogonal?
Matrix.Raise ErrDimensionMismatch unless square? Matrix.Raise ErrDimensionMismatch unless square?
rows.each_with_index do |row, i| rows.each_with_index do |row, i|
column_size.times do |j| column_count.times do |j|
s = 0 s = 0
row_size.times do |k| row_count.times do |k|
s += row[k] * rows[k][j] s += row[k] * rows[k][j]
end end
return false unless s == (i == j ? 1 : 0) return false unless s == (i == j ? 1 : 0)
@ -668,7 +670,7 @@ class Matrix
# #
def permutation? def permutation?
Matrix.Raise ErrDimensionMismatch unless square? Matrix.Raise ErrDimensionMismatch unless square?
cols = Array.new(column_size) cols = Array.new(column_count)
rows.each_with_index do |row, i| rows.each_with_index do |row, i|
found = false found = false
row.each_with_index do |e, j| row.each_with_index do |e, j|
@ -709,7 +711,7 @@ class Matrix
# Returns +true+ is this is a square matrix. # Returns +true+ is this is a square matrix.
# #
def square? def square?
column_size == row_size column_count == row_count
end end
# #
@ -731,9 +733,9 @@ class Matrix
def unitary? def unitary?
Matrix.Raise ErrDimensionMismatch unless square? Matrix.Raise ErrDimensionMismatch unless square?
rows.each_with_index do |row, i| rows.each_with_index do |row, i|
column_size.times do |j| column_count.times do |j|
s = 0 s = 0
row_size.times do |k| row_count.times do |k|
s += row[k].conj * rows[k][j] s += row[k].conj * rows[k][j]
end end
return false unless s == (i == j ? 1 : 0) return false unless s == (i == j ? 1 : 0)
@ -765,13 +767,13 @@ class Matrix
# #
def ==(other) def ==(other)
return false unless Matrix === other && return false unless Matrix === other &&
column_size == other.column_size # necessary for empty matrices column_count == other.column_count # necessary for empty matrices
rows == other.rows rows == other.rows
end end
def eql?(other) def eql?(other)
return false unless Matrix === other && return false unless Matrix === other &&
column_size == other.column_size # necessary for empty matrices column_count == other.column_count # necessary for empty matrices
rows.eql? other.rows rows.eql? other.rows
end end
@ -781,7 +783,7 @@ class Matrix
# There should be no good reason to do this since Matrices are immutable. # There should be no good reason to do this since Matrices are immutable.
# #
def clone def clone
new_matrix @rows.map(&:dup), column_size new_matrix @rows.map(&:dup), column_count
end end
# #
@ -807,22 +809,22 @@ class Matrix
rows = @rows.collect {|row| rows = @rows.collect {|row|
row.collect {|e| e * m } row.collect {|e| e * m }
} }
return new_matrix rows, column_size return new_matrix rows, column_count
when Vector when Vector
m = self.class.column_vector(m) m = self.class.column_vector(m)
r = self * m r = self * m
return r.column(0) return r.column(0)
when Matrix when Matrix
Matrix.Raise ErrDimensionMismatch if column_size != m.row_size Matrix.Raise ErrDimensionMismatch if column_count != m.row_count
rows = Array.new(row_size) {|i| rows = Array.new(row_count) {|i|
Array.new(m.column_size) {|j| Array.new(m.column_count) {|j|
(0 ... column_size).inject(0) do |vij, k| (0 ... column_count).inject(0) do |vij, k|
vij + self[i, k] * m[k, j] vij + self[i, k] * m[k, j]
end end
} }
} }
return new_matrix rows, m.column_size return new_matrix rows, m.column_count
else else
return apply_through_coercion(m, __method__) return apply_through_coercion(m, __method__)
end end
@ -845,14 +847,14 @@ class Matrix
return apply_through_coercion(m, __method__) return apply_through_coercion(m, __method__)
end end
Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count and column_count == m.column_count
rows = Array.new(row_size) {|i| rows = Array.new(row_count) {|i|
Array.new(column_size) {|j| Array.new(column_count) {|j|
self[i, j] + m[i, j] self[i, j] + m[i, j]
} }
} }
new_matrix rows, column_size new_matrix rows, column_count
end end
# #
@ -872,14 +874,14 @@ class Matrix
return apply_through_coercion(m, __method__) return apply_through_coercion(m, __method__)
end end
Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count and column_count == m.column_count
rows = Array.new(row_size) {|i| rows = Array.new(row_count) {|i|
Array.new(column_size) {|j| Array.new(column_count) {|j|
self[i, j] - m[i, j] self[i, j] - m[i, j]
} }
} }
new_matrix rows, column_size new_matrix rows, column_count
end end
# #
@ -894,7 +896,7 @@ class Matrix
rows = @rows.collect {|row| rows = @rows.collect {|row|
row.collect {|e| e / other } row.collect {|e| e / other }
} }
return new_matrix rows, column_size return new_matrix rows, column_count
when Matrix when Matrix
return self * other.inverse return self * other.inverse
else else
@ -910,12 +912,12 @@ class Matrix
# #
def inverse def inverse
Matrix.Raise ErrDimensionMismatch unless square? Matrix.Raise ErrDimensionMismatch unless square?
self.class.I(row_size).send(:inverse_from, self) self.class.I(row_count).send(:inverse_from, self)
end end
alias inv inverse alias inv inverse
def inverse_from(src) # :nodoc: def inverse_from(src) # :nodoc:
last = row_size - 1 last = row_count - 1
a = src.to_a a = src.to_a
0.upto(last) do |k| 0.upto(last) do |k|
@ -974,7 +976,7 @@ class Matrix
x = self x = self
if other <= 0 if other <= 0
x = self.inverse x = self.inverse
return self.class.identity(self.column_size) if other == 0 return self.class.identity(self.column_count) if other == 0
other = -other other = -other
end end
z = nil z = nil
@ -1008,7 +1010,7 @@ class Matrix
def determinant def determinant
Matrix.Raise ErrDimensionMismatch unless square? Matrix.Raise ErrDimensionMismatch unless square?
m = @rows m = @rows
case row_size case row_count
# Up to 4x4, give result using Laplacian expansion by minors. # Up to 4x4, give result using Laplacian expansion by minors.
# This will typically be faster, as well as giving good results # This will typically be faster, as well as giving good results
# in case of Floats # in case of Floats
@ -1057,7 +1059,7 @@ class Matrix
# intermediate results with better precision. # intermediate results with better precision.
# #
def determinant_bareiss def determinant_bareiss
size = row_size size = row_count
last = size - 1 last = size - 1
a = to_a a = to_a
no_pivot = Proc.new{ return 0 } no_pivot = Proc.new{ return 0 }
@ -1106,8 +1108,8 @@ class Matrix
# We currently use Bareiss' multistep integer-preserving gaussian elimination # We currently use Bareiss' multistep integer-preserving gaussian elimination
# (see comments on determinant) # (see comments on determinant)
a = to_a a = to_a
last_column = column_size - 1 last_column = column_count - 1
last_row = row_size - 1 last_row = row_count - 1
pivot_row = 0 pivot_row = 0
previous_pivot = 1 previous_pivot = 1
0.upto(last_column) do |k| 0.upto(last_column) do |k|
@ -1152,7 +1154,7 @@ class Matrix
# #
def trace def trace
Matrix.Raise ErrDimensionMismatch unless square? Matrix.Raise ErrDimensionMismatch unless square?
(0...column_size).inject(0) do |tr, i| (0...column_count).inject(0) do |tr, i|
tr + @rows[i][i] tr + @rows[i][i]
end end
end end
@ -1169,8 +1171,8 @@ class Matrix
# 2 4 6 # 2 4 6
# #
def transpose def transpose
return self.class.empty(column_size, 0) if row_size.zero? return self.class.empty(column_count, 0) if row_count.zero?
new_matrix @rows.transpose, row_size new_matrix @rows.transpose, row_count
end end
alias t transpose alias t transpose
@ -1286,7 +1288,7 @@ class Matrix
# Returns an array of the row vectors of the matrix. See Vector. # Returns an array of the row vectors of the matrix. See Vector.
# #
def row_vectors def row_vectors
Array.new(row_size) {|i| Array.new(row_count) {|i|
row(i) row(i)
} }
end end
@ -1295,7 +1297,7 @@ class Matrix
# Returns an array of the column vectors of the matrix. See Vector. # Returns an array of the column vectors of the matrix. See Vector.
# #
def column_vectors def column_vectors
Array.new(column_size) {|i| Array.new(column_count) {|i|
column(i) column(i)
} }
end end
@ -1331,7 +1333,7 @@ class Matrix
# #
def to_s def to_s
if empty? if empty?
"#{self.class}.empty(#{row_size}, #{column_size})" "#{self.class}.empty(#{row_count}, #{column_count})"
else else
"#{self.class}[" + @rows.collect{|row| "#{self.class}[" + @rows.collect{|row|
"[" + row.collect{|e| e.to_s}.join(", ") + "]" "[" + row.collect{|e| e.to_s}.join(", ") + "]"
@ -1344,7 +1346,7 @@ class Matrix
# #
def inspect def inspect
if empty? if empty?
"#{self.class}.empty(#{row_size}, #{column_size})" "#{self.class}.empty(#{row_count}, #{column_count})"
else else
"#{self.class}#{@rows.inspect}" "#{self.class}#{@rows.inspect}"
end end

2
lib/matrix/eigenvalue_decomposition.rb
View file

@ -20,7 +20,7 @@ class Matrix
# @v: Array for internal storage of eigenvectors. # @v: Array for internal storage of eigenvectors.
# @h: Array for internal storage of nonsymmetric Hessenberg form. # @h: Array for internal storage of nonsymmetric Hessenberg form.
raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix) raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
@size = a.row_size @size = a.row_count
@d = Array.new(@size, 0) @d = Array.new(@size, 0)
@e = Array.new(@size, 0) @e = Array.new(@size, 0)

56
lib/matrix/lup_decomposition.rb
View file

@ -19,7 +19,7 @@ class Matrix
include Matrix::ConversionHelper include Matrix::ConversionHelper
def l def l
Matrix.build(@row_size, @col_size) do |i, j| Matrix.build(@row_count, @column_count) do |i, j|
if (i > j) if (i > j)
@lu[i][j] @lu[i][j]
elsif (i == j) elsif (i == j)
@ -33,7 +33,7 @@ class Matrix
# Returns the upper triangular factor +U+ # Returns the upper triangular factor +U+
def u def u
Matrix.build(@col_size, @col_size) do |i, j| Matrix.build(@column_count, @column_count) do |i, j|
if (i <= j) if (i <= j)
@lu[i][j] @lu[i][j]
else else
@ -45,9 +45,9 @@ class Matrix
# Returns the permutation matrix +P+ # Returns the permutation matrix +P+
def p def p
rows = Array.new(@row_size){Array.new(@col_size, 0)} rows = Array.new(@row_count){Array.new(@column_count, 0)}
@pivots.each_with_index{|p, i| rows[i][p] = 1} @pivots.each_with_index{|p, i| rows[i][p] = 1}
Matrix.send :new, rows, @col_size Matrix.send :new, rows, @column_count
end end
# Returns +L+, +U+, +P+ in an array # Returns +L+, +U+, +P+ in an array
@ -64,7 +64,7 @@ class Matrix
# Returns +true+ if +U+, and hence +A+, is singular. # Returns +true+ if +U+, and hence +A+, is singular.
def singular? () def singular? ()
@col_size.times do |j| @column_count.times do |j|
if (@lu[j][j] == 0) if (@lu[j][j] == 0)
return true return true
end end
@ -76,11 +76,11 @@ class Matrix
# from the factorization. # from the factorization.
def det def det
if (@row_size != @col_size) if (@row_count != @column_count)
Matrix.Raise Matrix::ErrDimensionMismatch unless square? Matrix.Raise Matrix::ErrDimensionMismatch unless square?
end end
d = @pivot_sign d = @pivot_sign
@col_size.times do |j| @column_count.times do |j|
d *= @lu[j][j] d *= @lu[j][j]
end end
d d
@ -96,24 +96,24 @@ class Matrix
Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular." Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."
end end
if b.is_a? Matrix if b.is_a? Matrix
if (b.row_size != @row_size) if (b.row_count != @row_count)
Matrix.Raise Matrix::ErrDimensionMismatch Matrix.Raise Matrix::ErrDimensionMismatch
end end
# Copy right hand side with pivoting # Copy right hand side with pivoting
nx = b.column_size nx = b.column_count
m = @pivots.map{|row| b.row(row).to_a} m = @pivots.map{|row| b.row(row).to_a}
# Solve L*Y = P*b # Solve L*Y = P*b
@col_size.times do |k| @column_count.times do |k|
(k+1).upto(@col_size-1) do |i| (k+1).upto(@column_count-1) do |i|
nx.times do |j| nx.times do |j|
m[i][j] -= m[k][j]*@lu[i][k] m[i][j] -= m[k][j]*@lu[i][k]
end end
end end
end end
# Solve U*m = Y # Solve U*m = Y
(@col_size-1).downto(0) do |k| (@column_count-1).downto(0) do |k|
nx.times do |j| nx.times do |j|
m[k][j] = m[k][j].quo(@lu[k][k]) m[k][j] = m[k][j].quo(@lu[k][k])
end end
@ -126,7 +126,7 @@ class Matrix
Matrix.send :new, m, nx Matrix.send :new, m, nx
else # same algorithm, specialized for simpler case of a vector else # same algorithm, specialized for simpler case of a vector
b = convert_to_array(b) b = convert_to_array(b)
if (b.size != @row_size) if (b.size != @row_count)
Matrix.Raise Matrix::ErrDimensionMismatch Matrix.Raise Matrix::ErrDimensionMismatch
end end
@ -134,13 +134,13 @@ class Matrix
m = b.values_at(*@pivots) m = b.values_at(*@pivots)
# Solve L*Y = P*b # Solve L*Y = P*b
@col_size.times do |k| @column_count.times do |k|
(k+1).upto(@col_size-1) do |i| (k+1).upto(@column_count-1) do |i|
m[i] -= m[k]*@lu[i][k] m[i] -= m[k]*@lu[i][k]
end end
end end
# Solve U*m = Y # Solve U*m = Y
(@col_size-1).downto(0) do |k| (@column_count-1).downto(0) do |k|
m[k] = m[k].quo(@lu[k][k]) m[k] = m[k].quo(@lu[k][k])
k.times do |i| k.times do |i|
m[i] -= m[k]*@lu[i][k] m[i] -= m[k]*@lu[i][k]
@ -154,28 +154,28 @@ class Matrix
raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix) raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
# Use a "left-looking", dot-product, Crout/Doolittle algorithm. # Use a "left-looking", dot-product, Crout/Doolittle algorithm.
@lu = a.to_a @lu = a.to_a
@row_size = a.row_size @row_count = a.row_count
@col_size = a.column_size @column_count = a.column_count
@pivots = Array.new(@row_size) @pivots = Array.new(@row_count)
@row_size.times do |i| @row_count.times do |i|
@pivots[i] = i @pivots[i] = i
end end
@pivot_sign = 1 @pivot_sign = 1
lu_col_j = Array.new(@row_size) lu_col_j = Array.new(@row_count)
# Outer loop. # Outer loop.
@col_size.times do |j| @column_count.times do |j|
# Make a copy of the j-th column to localize references. # Make a copy of the j-th column to localize references.
@row_size.times do |i| @row_count.times do |i|
lu_col_j[i] = @lu[i][j] lu_col_j[i] = @lu[i][j]
end end
# Apply previous transformations. # Apply previous transformations.
@row_size.times do |i| @row_count.times do |i|
lu_row_i = @lu[i] lu_row_i = @lu[i]
# Most of the time is spent in the following dot product. # Most of the time is spent in the following dot product.
@ -192,13 +192,13 @@ class Matrix
# Find pivot and exchange if necessary. # Find pivot and exchange if necessary.
p = j p = j
(j+1).upto(@row_size-1) do |i| (j+1).upto(@row_count-1) do |i|
if (lu_col_j[i].abs > lu_col_j[p].abs) if (lu_col_j[i].abs > lu_col_j[p].abs)
p = i p = i
end end
end end
if (p != j) if (p != j)
@col_size.times do |k| @column_count.times do |k|
t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t
end end
k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k
@ -207,8 +207,8 @@ class Matrix
# Compute multipliers. # Compute multipliers.
if (j < @row_size && @lu[j][j] != 0) if (j < @row_count && @lu[j][j] != 0)
(j+1).upto(@row_size-1) do |i| (j+1).upto(@row_count-1) do |i|
@lu[i][j] = @lu[i][j].quo(@lu[j][j]) @lu[i][j] = @lu[i][j].quo(@lu[j][j])
end end
end end