Browse files

Imported matrix.rb from MRI r33295.

Since many of the failing Matrix specs are guarded with ruby_bug, we either
must fail them or fix them even though we are targetting 1.9.2. The fixes
would not seem to cause problems for Matrix users, but if they do, we can
revisit whether to import from 1.9.2 and fail the specs or update the specs
to use ruby_version_is rather than ruby_bug. There is a fine line between
a version difference and a bug. I tend to favor only using ruby_bug when
the value computed is clearly in error or there is a hang or segfault.
  • Loading branch information...
1 parent a636518 commit 6f7146be5d15f01bed99f1753b99254748efc539 @brixen brixen committed Sep 19, 2011
Showing with 1,487 additions and 54 deletions.
  1. +383 −54 lib/19/matrix.rb
  2. +886 −0 lib/19/matrix/eigenvalue_decomposition.rb
  3. +218 −0 lib/19/matrix/lup_decomposition.rb
View
437 lib/19/matrix.rb
@@ -56,14 +56,25 @@ module ExceptionForMatrix # :nodoc:
# * <tt> #map </tt>
# * <tt> #each </tt>
# * <tt> #each_with_index </tt>
+# * <tt> #find_index </tt>
# * <tt> #minor(*param) </tt>
#
# Properties of a matrix:
+# * <tt> #diagonal? </tt>
# * <tt> #empty? </tt>
+# * <tt> #hermitian? </tt>
+# * <tt> #lower_triangular? </tt>
+# * <tt> #normal? </tt>
+# * <tt> #orthogonal? </tt>
+# * <tt> #permutation? </tt>
# * <tt> #real? </tt>
# * <tt> #regular? </tt>
# * <tt> #singular? </tt>
# * <tt> #square? </tt>
+# * <tt> #symmetric? </tt>
+# * <tt> #unitary? </tt>
+# * <tt> #upper_triangular? </tt>
+# * <tt> #zero? </tt>
#
# Matrix arithmetic:
# * <tt> *(m) </tt>
@@ -78,11 +89,18 @@ module ExceptionForMatrix # :nodoc:
# * <tt> #determinant </tt>
# * <tt> #det </tt>
# * <tt> #rank </tt>
+# * <tt> #round </tt>
# * <tt> #trace </tt>
# * <tt> #tr </tt>
# * <tt> #transpose </tt>
# * <tt> #t </tt>
#
+# Matrix decompositions:
+# * <tt> #eigen </tt>
+# * <tt> #eigensystem </tt>
+# * <tt> #lup </tt>
+# * <tt> #lup_decomposition </tt>
+#
# Complex arithmetic:
# * <tt> conj </tt>
# * <tt> conjugate </tt>
@@ -105,6 +123,8 @@ module ExceptionForMatrix # :nodoc:
class Matrix
include Enumerable
include ExceptionForMatrix
+ autoload :EigenvalueDecomposition, "matrix/eigenvalue_decomposition"
+ autoload :LUPDecomposition, "matrix/lup_decomposition"
# instance creations
private_class_method :new
@@ -118,7 +138,7 @@ class Matrix
# -1 66
#
def Matrix.[](*rows)
- Matrix.rows(rows, false)
+ rows(rows, false)
end
#
@@ -148,7 +168,7 @@ def Matrix.rows(rows, copy = true)
# 93 66
#
def Matrix.columns(columns)
- Matrix.rows(columns, false).transpose
+ rows(columns, false).transpose
end
#
@@ -200,7 +220,7 @@ def Matrix.diagonal(*values)
# 0 5
#
def Matrix.scalar(n, value)
- Matrix.diagonal(*Array.new(n, value))
+ diagonal(*Array.new(n, value))
end
#
@@ -210,21 +230,22 @@ def Matrix.scalar(n, value)
# 0 1
#
def Matrix.identity(n)
- Matrix.scalar(n, 1)
+ scalar(n, 1)
end
class << Matrix
alias unit identity
alias I identity
end
#
- # Creates an +n+ by +n+ zero matrix.
+ # Creates a zero matrix.
# Matrix.zero(2)
# => 0 0
# 0 0
#
- def Matrix.zero(n)
- Matrix.scalar(n, 0)
+ def Matrix.zero(row_size, column_size = row_size)
+ rows = Array.new(row_size){Array.new(column_size, 0)}
+ new rows, column_size
end
#
@@ -283,7 +304,7 @@ def initialize(rows, column_size = rows[0].size)
end
def new_matrix(rows, column_size = rows[0].size) # :nodoc:
- Matrix.send(:new, rows, column_size) # bypass privacy of Matrix.new
+ self.class.send(:new, rows, column_size) # bypass privacy of Matrix.new
end
private :new_matrix
@@ -365,40 +386,161 @@ def collect(&block) # :yield: e
#
# Yields all elements of the matrix, starting with those of the first row,
- # or returns an Enumerator is no block given
+ # or returns an Enumerator is no block given.
+ # Elements can be restricted by passing an argument:
+ # * :all (default): yields all elements
+ # * :diagonal: yields only elements on the diagonal
+ # * :off_diagonal: yields all elements except on the diagonal
+ # * :lower: yields only elements on or below the diagonal
+ # * :strict_lower: yields only elements below the diagonal
+ # * :strict_upper: yields only elements above the diagonal
+ # * :upper: yields only elements on or above the diagonal
+ #
# Matrix[ [1,2], [3,4] ].each { |e| puts e }
# # => prints the numbers 1 to 4
- #
- def each(&block) # :yield: e
- return to_enum(:each) unless block_given?
- @rows.each do |row|
- row.each(&block)
+ # Matrix[ [1,2], [3,4] ].each(:strict_lower).to_a # => [3]
+ #
+ def each(which = :all) # :yield: e
+ return to_enum :each, which unless block_given?
+ last = column_size - 1
+ case which
+ when :all
+ block = Proc.new
+ @rows.each do |row|
+ row.each(&block)
+ end
+ when :diagonal
+ @rows.each_with_index do |row, row_index|
+ yield row.fetch(row_index){return self}
+ end
+ when :off_diagonal
+ @rows.each_with_index do |row, row_index|
+ column_size.times do |col_index|
+ yield row[col_index] unless row_index == col_index
+ end
+ end
+ when :lower
+ @rows.each_with_index do |row, row_index|
+ 0.upto([row_index, last].min) do |col_index|
+ yield row[col_index]
+ end
+ end
+ when :strict_lower
+ @rows.each_with_index do |row, row_index|
+ [row_index, column_size].min.times do |col_index|
+ yield row[col_index]
+ end
+ end
+ when :strict_upper
+ @rows.each_with_index do |row, row_index|
+ (row_index+1).upto(last) do |col_index|
+ yield row[col_index]
+ end
+ end
+ when :upper
+ @rows.each_with_index do |row, row_index|
+ row_index.upto(last) do |col_index|
+ yield row[col_index]
+ end
+ end
+ else
+ Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
end
self
end
#
- # Yields all elements of the matrix, starting with those of the first row,
- # along with the row index and column index,
- # or returns an Enumerator is no block given
+ # Same as #each, but the row index and column index in addition to the element
+ #
# Matrix[ [1,2], [3,4] ].each_with_index do |e, row, col|
# puts "#{e} at #{row}, #{col}"
# end
- # # => 1 at 0, 0
- # # => 2 at 0, 1
- # # => 3 at 1, 0
- # # => 4 at 1, 1
- #
- def each_with_index(&block) # :yield: e, row, column
- return to_enum(:each_with_index) unless block_given?
- @rows.each_with_index do |row, row_index|
- row.each_with_index do |e, col_index|
- yield e, row_index, col_index
+ # # => Prints:
+ # # 1 at 0, 0
+ # # 2 at 0, 1
+ # # 3 at 1, 0
+ # # 4 at 1, 1
+ #
+ def each_with_index(which = :all) # :yield: e, row, column
+ return to_enum :each_with_index, which unless block_given?
+ last = column_size - 1
+ case which
+ when :all
+ @rows.each_with_index do |row, row_index|
+ row.each_with_index do |e, col_index|
+ yield e, row_index, col_index
+ end
+ end
+ when :diagonal
+ @rows.each_with_index do |row, row_index|
+ yield row.fetch(row_index){return self}, row_index, row_index
+ end
+ when :off_diagonal
+ @rows.each_with_index do |row, row_index|
+ column_size.times do |col_index|
+ yield row[col_index], row_index, col_index unless row_index == col_index
+ end
+ end
+ when :lower
+ @rows.each_with_index do |row, row_index|
+ 0.upto([row_index, last].min) do |col_index|
+ yield row[col_index], row_index, col_index
+ end
+ end
+ when :strict_lower
+ @rows.each_with_index do |row, row_index|
+ [row_index, column_size].min.times do |col_index|
+ yield row[col_index], row_index, col_index
+ end
+ end
+ when :strict_upper
+ @rows.each_with_index do |row, row_index|
+ (row_index+1).upto(last) do |col_index|
+ yield row[col_index], row_index, col_index
+ end
+ end
+ when :upper
+ @rows.each_with_index do |row, row_index|
+ row_index.upto(last) do |col_index|
+ yield row[col_index], row_index, col_index
+ end
end
+ else
+ Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
end
self
end
+ SELECTORS = {all: true, diagonal: true, off_diagonal: true, lower: true, strict_lower: true, strict_upper: true, upper: true}.freeze
+ #
+ # :call-seq:
+ # index(value, selector = :all) -> [row, column]
+ # index(selector = :all){ block } -> [row, column]
+ # index(selector = :all) -> an_enumerator
+ #
+ # The index method is specialized to return the index as [row, column]
+ # It also accepts an optional +selector+ argument, see #each for details.
+ #
+ # Matrix[ [1,2], [3,4] ].index(&:even?) # => [0, 1]
+ # Matrix[ [1,1], [1,1] ].index(1, :strict_lower) # => [1, 0]
+ #
+ def index(*args)
+ raise ArgumentError, "wrong number of arguments(#{args.size} for 0-2)" if args.size > 2
+ which = (args.size == 2 || SELECTORS.include?(args.last)) ? args.pop : :all
+ return to_enum :find_index, which, *args unless block_given? || args.size == 1
+ if args.size == 1
+ value = args.first
+ each_with_index(which) do |e, row_index, col_index|
+ return row_index, col_index if e == value
+ end
+ else
+ each_with_index(which) do |e, row_index, col_index|
+ return row_index, col_index if yield e
+ end
+ end
+ nil
+ end
+ alias_method :find_index, :index
#
# Returns a section of the matrix. The parameters are either:
# * start_row, nrows, start_col, ncols; OR
@@ -450,6 +592,15 @@ def minor(*param)
#++
#
+ # Returns +true+ is this is a diagonal matrix.
+ # Raises an error if matrix is not square.
+ #
+ def diagonal?
+ Matrix.Raise ErrDimensionMismatch unless square?
+ each(:off_diagonal).all?(&:zero?)
+ end
+
+ #
# Returns +true+ if this is an empty matrix, i.e. if the number of rows
# or the number of columns is 0.
#
@@ -458,6 +609,82 @@ def empty?
end
#
+ # Returns +true+ is this is an hermitian matrix.
+ # Raises an error if matrix is not square.
+ #
+ def hermitian?
+ Matrix.Raise ErrDimensionMismatch unless square?
+ each_with_index(:strict_upper).all? do |e, row, col|
+ e == rows[col][row].conj
+ end
+ end
+
+ #
+ # Returns +true+ is this is a lower triangular matrix.
+ #
+ def lower_triangular?
+ each(:strict_upper).all?(&:zero?)
+ end
+
+ #
+ # Returns +true+ is this is a normal matrix.
+ # Raises an error if matrix is not square.
+ #
+ def normal?
+ Matrix.Raise ErrDimensionMismatch unless square?
+ rows.each_with_index do |row_i, i|
+ rows.each_with_index do |row_j, j|
+ s = 0
+ rows.each_with_index do |row_k, k|
+ s += row_i[k] * row_j[k].conj - row_k[i].conj * row_k[j]
+ end
+ return false unless s == 0
+ end
+ end
+ true
+ end
+
+ #
+ # Returns +true+ is this is an orthogonal matrix
+ # Raises an error if matrix is not square.
+ #
+ def orthogonal?
+ Matrix.Raise ErrDimensionMismatch unless square?
+ rows.each_with_index do |row, i|
+ column_size.times do |j|
+ s = 0
+ row_size.times do |k|
+ s += row[k] * rows[k][j]
+ end
+ return false unless s == (i == j ? 1 : 0)
+ end
+ end
+ true
+ end
+
+ #
+ # Returns +true+ is this is a permutation matrix
+ # Raises an error if matrix is not square.
+ #
+ def permutation?
+ Matrix.Raise ErrDimensionMismatch unless square?
+ cols = Array.new(column_size)
+ rows.each_with_index do |row, i|
+ found = false
+ row.each_with_index do |e, j|
+ if e == 1
+ return false if found || cols[j]
+ found = cols[j] = true
+ elsif e != 0
+ return false
+ end
+ end
+ return false unless found
+ end
+ true
+ end
+
+ #
# Returns +true+ if all entries of the matrix are real.
#
def real?
@@ -485,6 +712,49 @@ def square?
column_size == row_size
end
+ #
+ # Returns +true+ is this is a symmetric matrix.
+ # Raises an error if matrix is not square.
+ #
+ def symmetric?
+ Matrix.Raise ErrDimensionMismatch unless square?
+ each_with_index(:strict_upper).all? do |e, row, col|
+ e == rows[col][row]
+ end
+ end
+
+ #
+ # Returns +true+ is this is a unitary matrix
+ # Raises an error if matrix is not square.
+ #
+ def unitary?
+ Matrix.Raise ErrDimensionMismatch unless square?
+ rows.each_with_index do |row, i|
+ column_size.times do |j|
+ s = 0
+ row_size.times do |k|
+ s += row[k].conj * rows[k][j]
+ end
+ return false unless s == (i == j ? 1 : 0)
+ end
+ end
+ true
+ end
+
+ #
+ # Returns +true+ is this is an upper triangular matrix.
+ #
+ def upper_triangular?
+ each(:strict_lower).all?(&:zero?)
+ end
+
+ #
+ # Returns +true+ is this is a matrix with only zero elements
+ #
+ def zero?
+ all?(&:zero?)
+ end
+
#--
# OBJECT METHODS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
#++
@@ -538,7 +808,7 @@ def *(m) # m is matrix or vector or number
}
return new_matrix rows, column_size
when Vector
- m = Matrix.column_vector(m)
+ m = self.class.column_vector(m)
r = self * m
return r.column(0)
when Matrix
@@ -568,7 +838,7 @@ def +(m)
when Numeric
Matrix.Raise ErrOperationNotDefined, "+", self.class, m.class
when Vector
- m = Matrix.column_vector(m)
+ m = self.class.column_vector(m)
when Matrix
else
return apply_through_coercion(m, __method__)
@@ -595,7 +865,7 @@ def -(m)
when Numeric
Matrix.Raise ErrOperationNotDefined, "-", self.class, m.class
when Vector
- m = Matrix.column_vector(m)
+ m = self.class.column_vector(m)
when Matrix
else
return apply_through_coercion(m, __method__)
@@ -639,7 +909,7 @@ def /(other)
#
def inverse
Matrix.Raise ErrDimensionMismatch unless square?
- Matrix.I(row_size).send(:inverse_from, self)
+ self.class.I(row_size).send(:inverse_from, self)
end
alias inv inverse
@@ -689,8 +959,10 @@ def inverse_from(src) # :nodoc:
private :inverse_from
#
- # Matrix exponentiation. Currently implemented for integer powers only.
+ # Matrix exponentiation.
# Equivalent to multiplying the matrix by itself N times.
+ # Non integer exponents will be handled by diagonalizing the matrix.
+ #
# Matrix[[7,6], [3,9]] ** 2
# => 67 96
# 48 99
@@ -701,7 +973,7 @@ def ** (other)
x = self
if other <= 0
x = self.inverse
- return Matrix.identity(self.column_size) if other == 0
+ return self.class.identity(self.column_size) if other == 0
other = -other
end
z = nil
@@ -710,8 +982,9 @@ def ** (other)
return z if (other >>= 1).zero?
x *= x
end
- when Float, Rational
- Matrix.Raise ErrOperationNotImplemented, "**", self.class, other.class
+ when Numeric
+ v, d, v_inv = eigensystem
+ v * self.class.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv
else
Matrix.Raise ErrOperationNotDefined, "**", self.class, other.class
end
@@ -834,7 +1107,6 @@ def rank
a = to_a
last_column = column_size - 1
last_row = row_size - 1
- rank = 0
pivot_row = 0
previous_pivot = 1
0.upto(last_column) do |k|
@@ -865,6 +1137,12 @@ def rank_e
rank
end
+ # Returns a matrix with entries rounded to the given precision
+ # (see Float#round)
+ #
+ def round(ndigits=0)
+ map{|e| e.round(ndigits)}
+ end
#
# Returns the trace (sum of diagonal elements) of the matrix.
@@ -890,12 +1168,44 @@ def trace
# 2 4 6
#
def transpose
- return Matrix.empty(column_size, 0) if row_size.zero?
+ return self.class.empty(column_size, 0) if row_size.zero?
new_matrix @rows.transpose, row_size
end
alias t transpose
#--
+ # DECOMPOSITIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
+ #++
+
+ #
+ # Returns the Eigensystem of the matrix; see +EigenvalueDecomposition+.
+ # m = Matrix[[1, 2], [3, 4]]
+ # v, d, v_inv = m.eigensystem
+ # d.diagonal? # => true
+ # v.inv == v_inv # => true
+ # (v * d * v_inv).round(5) == m # => true
+ #
+ def eigensystem
+ EigenvalueDecomposition.new(self)
+ end
+ alias eigen eigensystem
+
+ #
+ # Returns the LUP decomposition of the matrix; see +LUPDecomposition+.
+ # a = Matrix[[1, 2], [3, 4]]
+ # l, u, p = a.lup
+ # l.lower_triangular? # => true
+ # u.upper_triangular? # => true
+ # p.permutation? # => true
+ # l * u == a * p # => true
+ # a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)]
+ #
+ def lup
+ LUPDecomposition.new(self)
+ end
+ alias lup_decomposition lup
+
+ #--
# COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#++
@@ -1020,9 +1330,9 @@ def elements_to_r
#
def to_s
if empty?
- "Matrix.empty(#{row_size}, #{column_size})"
+ "#{self.class}.empty(#{row_size}, #{column_size})"
else
- "Matrix[" + @rows.collect{|row|
+ "#{self.class}[" + @rows.collect{|row|
"[" + row.collect{|e| e.to_s}.join(", ") + "]"
}.join(", ")+"]"
end
@@ -1033,9 +1343,9 @@ def to_s
#
def inspect
if empty?
- "Matrix.empty(#{row_size}, #{column_size})"
+ "#{self.class}.empty(#{row_size}, #{column_size})"
else
- "Matrix#{@rows.inspect}"
+ "#{self.class}#{@rows.inspect}"
end
end
@@ -1207,8 +1517,11 @@ def ** (other)
# Vector functions:
# * <tt> #inner_product(v) </tt>
# * <tt> #collect </tt>
+# * <tt> #magnitude </tt>
# * <tt> #map </tt>
# * <tt> #map2(v) </tt>
+# * <tt> #norm </tt>
+# * <tt> #normalize </tt>
# * <tt> #r </tt>
# * <tt> #size </tt>
#
@@ -1237,7 +1550,7 @@ class Vector
# Vector[7, 4, ...]
#
def Vector.[](*array)
- new convert_to_array(array, copy = false)
+ new convert_to_array(array, false)
end
#
@@ -1341,7 +1654,7 @@ def eql?(other)
# Return a copy of the vector.
#
def clone
- Vector.elements(@elements)
+ self.class.elements(@elements)
end
#
@@ -1362,7 +1675,7 @@ def *(x)
case x
when Numeric
els = @elements.collect{|e| e * x}
- Vector.elements(els, false)
+ self.class.elements(els, false)
when Matrix
Matrix.column_vector(self) * x
when Vector
@@ -1382,7 +1695,7 @@ def +(v)
els = collect2(v) {|v1, v2|
v1 + v2
}
- Vector.elements(els, false)
+ self.class.elements(els, false)
when Matrix
Matrix.column_vector(self) + v
else
@@ -1400,7 +1713,7 @@ def -(v)
els = collect2(v) {|v1, v2|
v1 - v2
}
- Vector.elements(els, false)
+ self.class.elements(els, false)
when Matrix
Matrix.column_vector(self) - v
else
@@ -1415,7 +1728,7 @@ def /(x)
case x
when Numeric
els = @elements.collect{|e| e / x}
- Vector.elements(els, false)
+ self.class.elements(els, false)
when Matrix, Vector
Vector.Raise ErrOperationNotDefined, "/", self.class, x.class
else
@@ -1447,25 +1760,41 @@ def inner_product(v)
def collect(&block) # :yield: e
return to_enum(:collect) unless block_given?
els = @elements.collect(&block)
- Vector.elements(els, false)
+ self.class.elements(els, false)
end
alias map collect
#
+ # Returns the modulus (Pythagorean distance) of the vector.
+ # Vector[5,8,2].r => 9.643650761
+ #
+ def magnitude
+ Math.sqrt(@elements.inject(0) {|v, e| v + e*e})
+ end
+ alias r magnitude
+ alias norm magnitude
+
+ #
# Like Vector#collect2, but returns a Vector instead of an Array.
#
def map2(v, &block) # :yield: e1, e2
return to_enum(:map2, v) unless block_given?
els = collect2(v, &block)
- Vector.elements(els, false)
+ self.class.elements(els, false)
end
+ class ZeroVectorError < StandardError
+ end
#
- # Returns the modulus (Pythagorean distance) of the vector.
- # Vector[5,8,2].r => 9.643650761
+ # Returns a new vector with the same direction but with norm 1.
+ # v = Vector[5,8,2].normalize
+ # # => Vector[0.5184758473652127, 0.8295613557843402, 0.20739033894608505]
+ # v.norm => 1.0
#
- def r
- Math.sqrt(@elements.inject(0) {|v, e| v + e*e})
+ def normalize
+ n = magnitude
+ raise ZeroVectorError, "Zero vectors can not be normalized" if n == 0
+ self / n
end
#--
@@ -1532,6 +1861,6 @@ def to_s
# Overrides Object#inspect
#
def inspect
- str = "Vector"+@elements.inspect
+ "Vector" + @elements.inspect
end
end
View
886 lib/19/matrix/eigenvalue_decomposition.rb
@@ -0,0 +1,886 @@
+class Matrix
+ # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
+
+ # Eigenvalues and eigenvectors of a real matrix.
+ #
+ # Computes the eigenvalues and eigenvectors of a matrix A.
+ #
+ # If A is diagonalizable, this provides matrices V and D
+ # such that A = V*D*V.inv, where D is the diagonal matrix with entries
+ # equal to the eigenvalues and V is formed by the eigenvectors.
+ #
+ # If A is symmetric, then V is orthogonal and thus A = V*D*V.t
+
+ class EigenvalueDecomposition
+
+ # Constructs the eigenvalue decomposition for a square matrix +A+
+ #
+ def initialize(a)
+ # @d, @e: Arrays for internal storage of eigenvalues.
+ # @v: Array for internal storage of eigenvectors.
+ # @h: Array for internal storage of nonsymmetric Hessenberg form.
+ raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
+ @size = a.row_size
+ @d = Array.new(@size, 0)
+ @e = Array.new(@size, 0)
+
+ if (@symmetric = a.symmetric?)
+ @v = a.to_a
+ tridiagonalize
+ diagonalize
+ else
+ @v = Array.new(@size) { Array.new(@size, 0) }
+ @h = a.to_a
+ @ort = Array.new(@size, 0)
+ reduce_to_hessenberg
+ hessenberg_to_real_schur
+ end
+ end
+
+ # Returns the eigenvector matrix +V+
+ #
+ def eigenvector_matrix
+ Matrix.send :new, build_eigenvectors.transpose
+ end
+ alias v eigenvector_matrix
+
+ # Returns the inverse of the eigenvector matrix +V+
+ #
+ def eigenvector_matrix_inv
+ r = Matrix.send :new, build_eigenvectors
+ r = r.transpose.inverse unless @symmetric
+ r
+ end
+ alias v_inv eigenvector_matrix_inv
+
+ # Returns the eigenvalues in an array
+ #
+ def eigenvalues
+ values = @d.dup
+ @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0}
+ values
+ end
+
+ # Returns an array of the eigenvectors
+ #
+ def eigenvectors
+ build_eigenvectors.map{|ev| Vector.send :new, ev}
+ end
+
+ # Returns the block diagonal eigenvalue matrix +D+
+ #
+ def eigenvalue_matrix
+ Matrix.diagonal(*eigenvalues)
+ end
+ alias d eigenvalue_matrix
+
+ # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv]
+ #
+ def to_ary
+ [v, d, v_inv]
+ end
+ alias_method :to_a, :to_ary
+
+ private
+ def build_eigenvectors
+ # JAMA stores complex eigenvectors in a strange way
+ # See http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html
+ @e.each_with_index.map do |imag, i|
+ if imag == 0
+ Array.new(@size){|j| @v[j][i]}
+ elsif imag > 0
+ Array.new(@size){|j| Complex(@v[j][i], @v[j][i+1])}
+ else
+ Array.new(@size){|j| Complex(@v[j][i], -@v[j][i-1])}
+ end
+ end
+ end
+ # Complex scalar division.
+
+ def cdiv(xr, xi, yr, yi)
+ if (yr.abs > yi.abs)
+ r = yi/yr
+ d = yr + r*yi
+ [(xr + r*xi)/d, (xi - r*xr)/d]
+ else
+ r = yr/yi
+ d = yi + r*yr
+ [(r*xr + xi)/d, (r*xi - xr)/d]
+ end
+ end
+
+
+ # Symmetric Householder reduction to tridiagonal form.
+
+ def tridiagonalize
+
+ # This is derived from the Algol procedures tred2 by
+ # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ # Fortran subroutine in EISPACK.
+
+ @size.times do |j|
+ @d[j] = @v[@size-1][j]
+ end
+
+ # Householder reduction to tridiagonal form.
+
+ (@size-1).downto(0+1) do |i|
+
+ # Scale to avoid under/overflow.
+
+ scale = 0.0
+ h = 0.0
+ i.times do |k|
+ scale = scale + @d[k].abs
+ end
+ if (scale == 0.0)
+ @e[i] = @d[i-1]
+ i.times do |j|
+ @d[j] = @v[i-1][j]
+ @v[i][j] = 0.0
+ @v[j][i] = 0.0
+ end
+ else
+
+ # Generate Householder vector.
+
+ i.times do |k|
+ @d[k] /= scale
+ h += @d[k] * @d[k]
+ end
+ f = @d[i-1]
+ g = Math.sqrt(h)
+ if (f > 0)
+ g = -g
+ end
+ @e[i] = scale * g
+ h -= f * g
+ @d[i-1] = f - g
+ i.times do |j|
+ @e[j] = 0.0
+ end
+
+ # Apply similarity transformation to remaining columns.
+
+ i.times do |j|
+ f = @d[j]
+ @v[j][i] = f
+ g = @e[j] + @v[j][j] * f
+ (j+1).upto(i-1) do |k|
+ g += @v[k][j] * @d[k]
+ @e[k] += @v[k][j] * f
+ end
+ @e[j] = g
+ end
+ f = 0.0
+ i.times do |j|
+ @e[j] /= h
+ f += @e[j] * @d[j]
+ end
+ hh = f / (h + h)
+ i.times do |j|
+ @e[j] -= hh * @d[j]
+ end
+ i.times do |j|
+ f = @d[j]
+ g = @e[j]
+ j.upto(i-1) do |k|
+ @v[k][j] -= (f * @e[k] + g * @d[k])
+ end
+ @d[j] = @v[i-1][j]
+ @v[i][j] = 0.0
+ end
+ end
+ @d[i] = h
+ end
+
+ # Accumulate transformations.
+
+ 0.upto(@size-1-1) do |i|
+ @v[@size-1][i] = @v[i][i]
+ @v[i][i] = 1.0
+ h = @d[i+1]
+ if (h != 0.0)
+ 0.upto(i) do |k|
+ @d[k] = @v[k][i+1] / h
+ end
+ 0.upto(i) do |j|
+ g = 0.0
+ 0.upto(i) do |k|
+ g += @v[k][i+1] * @v[k][j]
+ end
+ 0.upto(i) do |k|
+ @v[k][j] -= g * @d[k]
+ end
+ end
+ end
+ 0.upto(i) do |k|
+ @v[k][i+1] = 0.0
+ end
+ end
+ @size.times do |j|
+ @d[j] = @v[@size-1][j]
+ @v[@size-1][j] = 0.0
+ end
+ @v[@size-1][@size-1] = 1.0
+ @e[0] = 0.0
+ end
+
+
+ # Symmetric tridiagonal QL algorithm.
+
+ def diagonalize
+ # This is derived from the Algol procedures tql2, by
+ # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ # Fortran subroutine in EISPACK.
+
+ 1.upto(@size-1) do |i|
+ @e[i-1] = @e[i]
+ end
+ @e[@size-1] = 0.0
+
+ f = 0.0
+ tst1 = 0.0
+ eps = Float::EPSILON
+ @size.times do |l|
+
+ # Find small subdiagonal element
+
+ tst1 = [tst1, @d[l].abs + @e[l].abs].max
+ m = l
+ while (m < @size) do
+ if (@e[m].abs <= eps*tst1)
+ break
+ end
+ m+=1
+ end
+
+ # If m == l, @d[l] is an eigenvalue,
+ # otherwise, iterate.
+
+ if (m > l)
+ iter = 0
+ begin
+ iter = iter + 1 # (Could check iteration count here.)
+
+ # Compute implicit shift
+
+ g = @d[l]
+ p = (@d[l+1] - g) / (2.0 * @e[l])
+ r = Math.hypot(p, 1.0)
+ if (p < 0)
+ r = -r
+ end
+ @d[l] = @e[l] / (p + r)
+ @d[l+1] = @e[l] * (p + r)
+ dl1 = @d[l+1]
+ h = g - @d[l]
+ (l+2).upto(@size-1) do |i|
+ @d[i] -= h
+ end
+ f += h
+
+ # Implicit QL transformation.
+
+ p = @d[m]
+ c = 1.0
+ c2 = c
+ c3 = c
+ el1 = @e[l+1]
+ s = 0.0
+ s2 = 0.0
+ (m-1).downto(l) do |i|
+ c3 = c2
+ c2 = c
+ s2 = s
+ g = c * @e[i]
+ h = c * p
+ r = Math.hypot(p, @e[i])
+ @e[i+1] = s * r
+ s = @e[i] / r
+ c = p / r
+ p = c * @d[i] - s * g
+ @d[i+1] = h + s * (c * g + s * @d[i])
+
+ # Accumulate transformation.
+
+ @size.times do |k|
+ h = @v[k][i+1]
+ @v[k][i+1] = s * @v[k][i] + c * h
+ @v[k][i] = c * @v[k][i] - s * h
+ end
+ end
+ p = -s * s2 * c3 * el1 * @e[l] / dl1
+ @e[l] = s * p
+ @d[l] = c * p
+
+ # Check for convergence.
+
+ end while (@e[l].abs > eps*tst1)
+ end
+ @d[l] = @d[l] + f
+ @e[l] = 0.0
+ end
+
+ # Sort eigenvalues and corresponding vectors.
+
+ 0.upto(@size-2) do |i|
+ k = i
+ p = @d[i]
+ (i+1).upto(@size-1) do |j|
+ if (@d[j] < p)
+ k = j
+ p = @d[j]
+ end
+ end
+ if (k != i)
+ @d[k] = @d[i]
+ @d[i] = p
+ @size.times do |j|
+ p = @v[j][i]
+ @v[j][i] = @v[j][k]
+ @v[j][k] = p
+ end
+ end
+ end
+ end
+
+ # Nonsymmetric reduction to Hessenberg form.
+
+ def reduce_to_hessenberg
+ # This is derived from the Algol procedures orthes and ortran,
+ # by Martin and Wilkinson, Handbook for Auto. Comp.,
+ # Vol.ii-Linear Algebra, and the corresponding
+ # Fortran subroutines in EISPACK.
+
+ low = 0
+ high = @size-1
+
+ (low+1).upto(high-1) do |m|
+
+ # Scale column.
+
+ scale = 0.0
+ m.upto(high) do |i|
+ scale = scale + @h[i][m-1].abs
+ end
+ if (scale != 0.0)
+
+ # Compute Householder transformation.
+
+ h = 0.0
+ high.downto(m) do |i|
+ @ort[i] = @h[i][m-1]/scale
+ h += @ort[i] * @ort[i]
+ end
+ g = Math.sqrt(h)
+ if (@ort[m] > 0)
+ g = -g
+ end
+ h -= @ort[m] * g
+ @ort[m] = @ort[m] - g
+
+ # Apply Householder similarity transformation
+ # @h = (I-u*u'/h)*@h*(I-u*u')/h)
+
+ m.upto(@size-1) do |j|
+ f = 0.0
+ high.downto(m) do |i|
+ f += @ort[i]*@h[i][j]
+ end
+ f = f/h
+ m.upto(high) do |i|
+ @h[i][j] -= f*@ort[i]
+ end
+ end
+
+ 0.upto(high) do |i|
+ f = 0.0
+ high.downto(m) do |j|
+ f += @ort[j]*@h[i][j]
+ end
+ f = f/h
+ m.upto(high) do |j|
+ @h[i][j] -= f*@ort[j]
+ end
+ end
+ @ort[m] = scale*@ort[m]
+ @h[m][m-1] = scale*g
+ end
+ end
+
+ # Accumulate transformations (Algol's ortran).
+
+ @size.times do |i|
+ @size.times do |j|
+ @v[i][j] = (i == j ? 1.0 : 0.0)
+ end
+ end
+
+ (high-1).downto(low+1) do |m|
+ if (@h[m][m-1] != 0.0)
+ (m+1).upto(high) do |i|
+ @ort[i] = @h[i][m-1]
+ end
+ m.upto(high) do |j|
+ g = 0.0
+ m.upto(high) do |i|
+ g += @ort[i] * @v[i][j]
+ end
+ # Double division avoids possible underflow
+ g = (g / @ort[m]) / @h[m][m-1]
+ m.upto(high) do |i|
+ @v[i][j] += g * @ort[i]
+ end
+ end
+ end
+ end
+ end
+
+
+
+ # Nonsymmetric reduction from Hessenberg to real Schur form.
+
+ def hessenberg_to_real_schur
+
+ # This is derived from the Algol procedure hqr2,
+ # by Martin and Wilkinson, Handbook for Auto. Comp.,
+ # Vol.ii-Linear Algebra, and the corresponding
+ # Fortran subroutine in EISPACK.
+
+ # Initialize
+
+ nn = @size
+ n = nn-1
+ low = 0
+ high = nn-1
+ eps = Float::EPSILON
+ exshift = 0.0
+ p=q=r=s=z=0
+
+ # Store roots isolated by balanc and compute matrix norm
+
+ norm = 0.0
+ nn.times do |i|
+ if (i < low || i > high)
+ @d[i] = @h[i][i]
+ @e[i] = 0.0
+ end
+ ([i-1, 0].max).upto(nn-1) do |j|
+ norm = norm + @h[i][j].abs
+ end
+ end
+
+ # Outer loop over eigenvalue index
+
+ iter = 0
+ while (n >= low) do
+
+ # Look for single small sub-diagonal element
+
+ l = n
+ while (l > low) do
+ s = @h[l-1][l-1].abs + @h[l][l].abs
+ if (s == 0.0)
+ s = norm
+ end
+ if (@h[l][l-1].abs < eps * s)
+ break
+ end
+ l-=1
+ end
+
+ # Check for convergence
+ # One root found
+
+ if (l == n)
+ @h[n][n] = @h[n][n] + exshift
+ @d[n] = @h[n][n]
+ @e[n] = 0.0
+ n-=1
+ iter = 0
+
+ # Two roots found
+
+ elsif (l == n-1)
+ w = @h[n][n-1] * @h[n-1][n]
+ p = (@h[n-1][n-1] - @h[n][n]) / 2.0
+ q = p * p + w
+ z = Math.sqrt(q.abs)
+ @h[n][n] = @h[n][n] + exshift
+ @h[n-1][n-1] = @h[n-1][n-1] + exshift
+ x = @h[n][n]
+
+ # Real pair
+
+ if (q >= 0)
+ if (p >= 0)
+ z = p + z
+ else
+ z = p - z
+ end
+ @d[n-1] = x + z
+ @d[n] = @d[n-1]
+ if (z != 0.0)
+ @d[n] = x - w / z
+ end
+ @e[n-1] = 0.0
+ @e[n] = 0.0
+ x = @h[n][n-1]
+ s = x.abs + z.abs
+ p = x / s
+ q = z / s
+ r = Math.sqrt(p * p+q * q)
+ p /= r
+ q /= r
+
+ # Row modification
+
+ (n-1).upto(nn-1) do |j|
+ z = @h[n-1][j]
+ @h[n-1][j] = q * z + p * @h[n][j]
+ @h[n][j] = q * @h[n][j] - p * z
+ end
+
+ # Column modification
+
+ 0.upto(n) do |i|
+ z = @h[i][n-1]
+ @h[i][n-1] = q * z + p * @h[i][n]
+ @h[i][n] = q * @h[i][n] - p * z
+ end
+
+ # Accumulate transformations
+
+ low.upto(high) do |i|
+ z = @v[i][n-1]
+ @v[i][n-1] = q * z + p * @v[i][n]
+ @v[i][n] = q * @v[i][n] - p * z
+ end
+
+ # Complex pair
+
+ else
+ @d[n-1] = x + p
+ @d[n] = x + p
+ @e[n-1] = z
+ @e[n] = -z
+ end
+ n -= 2
+ iter = 0
+
+ # No convergence yet
+
+ else
+
+ # Form shift
+
+ x = @h[n][n]
+ y = 0.0
+ w = 0.0
+ if (l < n)
+ y = @h[n-1][n-1]
+ w = @h[n][n-1] * @h[n-1][n]
+ end
+
+ # Wilkinson's original ad hoc shift
+
+ if (iter == 10)
+ exshift += x
+ low.upto(n) do |i|
+ @h[i][i] -= x
+ end
+ s = @h[n][n-1].abs + @h[n-1][n-2].abs
+ x = y = 0.75 * s
+ w = -0.4375 * s * s
+ end
+
+ # MATLAB's new ad hoc shift
+
+ if (iter == 30)
+ s = (y - x) / 2.0
+ s *= s + w
+ if (s > 0)
+ s = Math.sqrt(s)
+ if (y < x)
+ s = -s
+ end
+ s = x - w / ((y - x) / 2.0 + s)
+ low.upto(n) do |i|
+ @h[i][i] -= s
+ end
+ exshift += s
+ x = y = w = 0.964
+ end
+ end
+
+ iter = iter + 1 # (Could check iteration count here.)
+
+ # Look for two consecutive small sub-diagonal elements
+
+ m = n-2
+ while (m >= l) do
+ z = @h[m][m]
+ r = x - z
+ s = y - z
+ p = (r * s - w) / @h[m+1][m] + @h[m][m+1]
+ q = @h[m+1][m+1] - z - r - s
+ r = @h[m+2][m+1]
+ s = p.abs + q.abs + r.abs
+ p /= s
+ q /= s
+ r /= s
+ if (m == l)
+ break
+ end
+ if (@h[m][m-1].abs * (q.abs + r.abs) <
+ eps * (p.abs * (@h[m-1][m-1].abs + z.abs +
+ @h[m+1][m+1].abs)))
+ break
+ end
+ m-=1
+ end
+
+ (m+2).upto(n) do |i|
+ @h[i][i-2] = 0.0
+ if (i > m+2)
+ @h[i][i-3] = 0.0
+ end
+ end
+
+ # Double QR step involving rows l:n and columns m:n
+
+ m.upto(n-1) do |k|
+ notlast = (k != n-1)
+ if (k != m)
+ p = @h[k][k-1]
+ q = @h[k+1][k-1]
+ r = (notlast ? @h[k+2][k-1] : 0.0)
+ x = p.abs + q.abs + r.abs
+ if (x != 0.0)
+ p /= x
+ q /= x
+ r /= x
+ end
+ end
+ if (x == 0.0)
+ break
+ end
+ s = Math.sqrt(p * p + q * q + r * r)
+ if (p < 0)
+ s = -s
+ end
+ if (s != 0)
+ if (k != m)
+ @h[k][k-1] = -s * x
+ elsif (l != m)
+ @h[k][k-1] = -@h[k][k-1]
+ end
+ p += s
+ x = p / s
+ y = q / s
+ z = r / s
+ q /= p
+ r /= p
+
+ # Row modification
+
+ k.upto(nn-1) do |j|
+ p = @h[k][j] + q * @h[k+1][j]
+ if (notlast)
+ p += r * @h[k+2][j]
+ @h[k+2][j] = @h[k+2][j] - p * z
+ end
+ @h[k][j] = @h[k][j] - p * x
+ @h[k+1][j] = @h[k+1][j] - p * y
+ end
+
+ # Column modification
+
+ 0.upto([n, k+3].min) do |i|
+ p = x * @h[i][k] + y * @h[i][k+1]
+ if (notlast)
+ p += z * @h[i][k+2]
+ @h[i][k+2] = @h[i][k+2] - p * r
+ end
+ @h[i][k] = @h[i][k] - p
+ @h[i][k+1] = @h[i][k+1] - p * q
+ end
+
+ # Accumulate transformations
+
+ low.upto(high) do |i|
+ p = x * @v[i][k] + y * @v[i][k+1]
+ if (notlast)
+ p += z * @v[i][k+2]
+ @v[i][k+2] = @v[i][k+2] - p * r
+ end
+ @v[i][k] = @v[i][k] - p
+ @v[i][k+1] = @v[i][k+1] - p * q
+ end
+ end # (s != 0)
+ end # k loop
+ end # check convergence
+ end # while (n >= low)
+
+ # Backsubstitute to find vectors of upper triangular form
+
+ if (norm == 0.0)
+ return
+ end
+
+ (nn-1).downto(0) do |n|
+ p = @d[n]
+ q = @e[n]
+
+ # Real vector
+
+ if (q == 0)
+ l = n
+ @h[n][n] = 1.0
+ (n-1).downto(0) do |i|
+ w = @h[i][i] - p
+ r = 0.0
+ l.upto(n) do |j|
+ r += @h[i][j] * @h[j][n]
+ end
+ if (@e[i] < 0.0)
+ z = w
+ s = r
+ else
+ l = i
+ if (@e[i] == 0.0)
+ if (w != 0.0)
+ @h[i][n] = -r / w
+ else
+ @h[i][n] = -r / (eps * norm)
+ end
+
+ # Solve real equations
+
+ else
+ x = @h[i][i+1]
+ y = @h[i+1][i]
+ q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i]
+ t = (x * s - z * r) / q
+ @h[i][n] = t
+ if (x.abs > z.abs)
+ @h[i+1][n] = (-r - w * t) / x
+ else
+ @h[i+1][n] = (-s - y * t) / z
+ end
+ end
+
+ # Overflow control
+
+ t = @h[i][n].abs
+ if ((eps * t) * t > 1)
+ i.upto(n) do |j|
+ @h[j][n] = @h[j][n] / t
+ end
+ end
+ end
+ end
+
+ # Complex vector
+
+ elsif (q < 0)
+ l = n-1
+
+ # Last vector component imaginary so matrix is triangular
+
+ if (@h[n][n-1].abs > @h[n-1][n].abs)
+ @h[n-1][n-1] = q / @h[n][n-1]
+ @h[n-1][n] = -(@h[n][n] - p) / @h[n][n-1]
+ else
+ cdivr, cdivi = cdiv(0.0, -@h[n-1][n], @h[n-1][n-1]-p, q)
+ @h[n-1][n-1] = cdivr
+ @h[n-1][n] = cdivi
+ end
+ @h[n][n-1] = 0.0
+ @h[n][n] = 1.0
+ (n-2).downto(0) do |i|
+ ra = 0.0
+ sa = 0.0
+ l.upto(n) do |j|
+ ra = ra + @h[i][j] * @h[j][n-1]
+ sa = sa + @h[i][j] * @h[j][n]
+ end
+ w = @h[i][i] - p
+
+ if (@e[i] < 0.0)
+ z = w
+ r = ra
+ s = sa
+ else
+ l = i
+ if (@e[i] == 0)
+ cdivr, cdivi = cdiv(-ra, -sa, w, q)
+ @h[i][n-1] = cdivr
+ @h[i][n] = cdivi
+ else
+
+ # Solve complex equations
+
+ x = @h[i][i+1]
+ y = @h[i+1][i]
+ vr = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] - q * q
+ vi = (@d[i] - p) * 2.0 * q
+ if (vr == 0.0 && vi == 0.0)
+ vr = eps * norm * (w.abs + q.abs +
+ x.abs + y.abs + z.abs)
+ end
+ cdivr, cdivi = cdiv(x*r-z*ra+q*sa, x*s-z*sa-q*ra, vr, vi)
+ @h[i][n-1] = cdivr
+ @h[i][n] = cdivi
+ if (x.abs > (z.abs + q.abs))
+ @h[i+1][n-1] = (-ra - w * @h[i][n-1] + q * @h[i][n]) / x
+ @h[i+1][n] = (-sa - w * @h[i][n] - q * @h[i][n-1]) / x
+ else
+ cdivr, cdivi = cdiv(-r-y*@h[i][n-1], -s-y*@h[i][n], z, q)
+ @h[i+1][n-1] = cdivr
+ @h[i+1][n] = cdivi
+ end
+ end
+
+ # Overflow control
+
+ t = [@h[i][n-1].abs, @h[i][n].abs].max
+ if ((eps * t) * t > 1)
+ i.upto(n) do |j|
+ @h[j][n-1] = @h[j][n-1] / t
+ @h[j][n] = @h[j][n] / t
+ end
+ end
+ end
+ end
+ end
+ end
+
+ # Vectors of isolated roots
+
+ nn.times do |i|
+ if (i < low || i > high)
+ i.upto(nn-1) do |j|
+ @v[i][j] = @h[i][j]
+ end
+ end
+ end
+
+ # Back transformation to get eigenvectors of original matrix
+
+ (nn-1).downto(low) do |j|
+ low.upto(high) do |i|
+ z = 0.0
+ low.upto([j, high].min) do |k|
+ z += @v[i][k] * @h[k][j]
+ end
+ @v[i][j] = z
+ end
+ end
+ end
+
+ end
+end
View
218 lib/19/matrix/lup_decomposition.rb
@@ -0,0 +1,218 @@
+class Matrix
+ # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
+
+ #
+ # 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.
+ #
+
+ class LUPDecomposition
+ # Returns the lower triangular factor +L+
+
+ include Matrix::ConversionHelper
+
+ def l
+ Matrix.build(@row_size, @col_size) do |i, j|
+ if (i > j)
+ @lu[i][j]
+ elsif (i == j)
+ 1
+ else
+ 0
+ end
+ end
+ end
+
+ # Returns the upper triangular factor +U+
+
+ def u
+ Matrix.build(@col_size, @col_size) do |i, j|
+ if (i <= j)
+ @lu[i][j]
+ else
+ 0
+ end
+ end
+ end
+
+ # Returns the permutation matrix +P+
+
+ def p
+ rows = Array.new(@row_size){Array.new(@col_size, 0)}
+ @pivots.each_with_index{|p, i| rows[i][p] = 1}
+ Matrix.send :new, rows, @col_size
+ end
+
+ # Returns +L+, +U+, +P+ in an array
+
+ def to_ary
+ [l, u, p]
+ end
+ alias_method :to_a, :to_ary
+
+ # Returns the pivoting indices
+
+ attr_reader :pivots
+
+ # Returns +true+ if +U+, and hence +A+, is singular.
+
+ def singular? ()
+ @col_size.times do |j|
+ if (@lu[j][j] == 0)
+ return true
+ end
+ end
+ false
+ end
+
+ # Returns the determinant of +A+, calculated efficiently
+ # from the factorization.
+
+ def det
+ if (@row_size != @col_size)
+ Matrix.Raise Matrix::ErrDimensionMismatch unless square?
+ end
+ d = @pivot_sign
+ @col_size.times do |j|
+ d *= @lu[j][j]
+ end
+ d
+ end
+ alias_method :determinant, :det
+
+ # Returns +m+ so that <tt>A*m = b</tt>,
+ # or equivalently so that <tt>L*U*m = P*b</tt>
+ # +b+ can be a Matrix or a Vector
+
+ def solve b
+ if (singular?)
+ Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."
+ end
+ if b.is_a? Matrix
+ if (b.row_size != @row_size)
+ Matrix.Raise Matrix::ErrDimensionMismatch
+ end
+
+ # Copy right hand side with pivoting
+ nx = b.column_size
+ m = @pivots.map{|row| b.row(row).to_a}
+
+ # Solve L*Y = P*b
+ @col_size.times do |k|
+ (k+1).upto(@col_size-1) do |i|
+ nx.times do |j|
+ m[i][j] -= m[k][j]*@lu[i][k]
+ end
+ end
+ end
+ # Solve U*m = Y
+ (@col_size-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_size)
+ Matrix.Raise Matrix::ErrDimensionMismatch
+ end
+
+ # Copy right hand side with pivoting
+ m = b.values_at(*@pivots)
+
+ # Solve L*Y = P*b
+ @col_size.times do |k|
+ (k+1).upto(@col_size-1) do |i|
+ m[i] -= m[k]*@lu[i][k]
+ end
+ end
+ # Solve U*m = Y
+ (@col_size-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
+
+ 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_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)
+
+ # Outer loop.
+
+ @col_size.times do |j|
+
+ # Make a copy of the j-th column to localize references.
+
+ @row_size.times do |i|
+ lu_col_j[i] = @lu[i][j]
+ end
+
+ # Apply previous transformations.
+
+ @row_size.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_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
+
+ # Compute multipliers.
+
+ 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
+ end
+end

0 comments on commit 6f7146b

Please sign in to comment.