class Matrix

The Matrix class represents a mathematical matrix. It provides methods for creating matrices, operating on them arithmetically and algebraically, and determining their mathematical properties such as trace, rank, inverse, determinant, or eigensystem.

Constants

SELECTORS

Attributes

rows[R]

instance creations

Public Class Methods

I(n)
Alias for: identity
[](*rows) click to toggle source

Creates a matrix where each argument is a row.

Matrix[ [25, 93], [-1, 66] ]
   =>  25 93
       -1 66
# File lib/matrix.rb, line 51
def Matrix.[](*rows)
  rows(rows, false)
end
build(row_count, column_count = row_count) { |i, j| ... } click to toggle source

Creates a matrix of size row_count x column_count. It fills the values by calling the given block, passing the current row and column. Returns an enumerator if no block is given.

m = Matrix.build(2, 4) {|row, col| col - row }
  => Matrix[[0, 1, 2, 3], [-1, 0, 1, 2]]
m = Matrix.build(3) { rand }
  => a 3x3 matrix with random elements
# File lib/matrix.rb, line 96
def Matrix.build(row_count, column_count = row_count)
  row_count = CoercionHelper.coerce_to_int(row_count)
  column_count = CoercionHelper.coerce_to_int(column_count)
  raise ArgumentError if row_count < 0 || column_count < 0
  return to_enum :build, row_count, column_count unless block_given?
  rows = Array.new(row_count) do |i|
    Array.new(column_count) do |j|
      yield i, j
    end
  end
  new rows, column_count
end
column_vector(column) click to toggle source

Creates a single-column matrix where the values of that column are as given in column.

Matrix.column_vector([4,5,6])
  => 4
     5
     6
# File lib/matrix.rb, line 182
def Matrix.column_vector(column)
  column = convert_to_array(column)
  new [column].transpose, 1
end
columns(columns) click to toggle source

Creates a matrix using columns as an array of column vectors.

Matrix.columns([[25, 93], [-1, 66]])
   =>  25 -1
       93 66
# File lib/matrix.rb, line 81
def Matrix.columns(columns)
  rows(columns, false).transpose
end
combine(*matrices) { |map{|m| m}| ... } click to toggle source

Create a matrix by combining matrices entrywise, using the given block

x = Matrix[[6, 6], [4, 4]]
y = Matrix[[1, 2], [3, 4]]
Matrix.combine(x, y) {|a, b| a - b} # => Matrix[[5, 4], [1, 0]]
# File lib/matrix.rb, line 259
def Matrix.combine(*matrices)
  return to_enum(__method__, *matrices) unless block_given?

  return Matrix.empty if matrices.empty?
  matrices.map!(&CoercionHelper.method(:coerce_to_matrix))
  x = matrices.first
  matrices.each do |m|
    Matrix.Raise ErrDimensionMismatch unless x.row_count == m.row_count && x.column_count == m.column_count
  end

  rows = Array.new(x.row_count) do |i|
    Array.new(x.column_count) do |j|
      yield matrices.map{|m| m[i,j]}
    end
  end
  new rows, x.column_count
end
diagonal(*values) click to toggle source

Creates a matrix where the diagonal elements are composed of values.

Matrix.diagonal(9, 5, -3)
  =>  9  0  0
      0  5  0
      0  0 -3
# File lib/matrix.rb, line 116
def Matrix.diagonal(*values)
  size = values.size
  return Matrix.empty if size == 0
  rows = Array.new(size) {|j|
    row = Array.new(size, 0)
    row[j] = values[j]
    row
  }
  new rows
end
empty(row_count = 0, column_count = 0) click to toggle source

Creates a empty matrix of row_count x column_count. At least one of row_count or column_count must be 0.

m = Matrix.empty(2, 0)
m == Matrix[ [], [] ]
  => true
n = Matrix.empty(0, 3)
n == Matrix.columns([ [], [], [] ])
  => true
m * n
  => Matrix[[0, 0, 0], [0, 0, 0]]
# File lib/matrix.rb, line 200
def Matrix.empty(row_count = 0, column_count = 0)
  raise ArgumentError, "One size must be 0" if column_count != 0 && row_count != 0
  raise ArgumentError, "Negative size" if column_count < 0 || row_count < 0

  new([[]]*row_count, column_count)
end
hstack(x, *matrices) click to toggle source

Create a matrix by stacking matrices horizontally

x = Matrix[[1, 2], [3, 4]]
y = Matrix[[5, 6], [7, 8]]
Matrix.hstack(x, y) # => Matrix[[1, 2, 5, 6], [3, 4, 7, 8]]
# File lib/matrix.rb, line 235
def Matrix.hstack(x, *matrices)
  x = CoercionHelper.coerce_to_matrix(x)
  result = x.send(:rows).map(&:dup)
  total_column_count = x.column_count
  matrices.each do |m|
    m = CoercionHelper.coerce_to_matrix(m)
    if m.row_count != x.row_count
      raise ErrDimensionMismatch, "The given matrices must have #{x.row_count} rows, but one has #{m.row_count}"
    end
    result.each_with_index do |row, i|
      row.concat m.send(:rows)[i]
    end
    total_column_count += m.column_count
  end
  new result, total_column_count
end
identity(n) click to toggle source

Creates an n by n identity matrix.

Matrix.identity(2)
  => 1 0
     0 1
# File lib/matrix.rb, line 144
def Matrix.identity(n)
  scalar(n, 1)
end
Also aliased as: unit, I
new(rows, column_count = rows[0].size) click to toggle source

::new is private; use ::rows, columns, [], etc… to create.

# File lib/matrix.rb, line 284
def initialize(rows, column_count = rows[0].size)
  # No checking is done at this point. rows must be an Array of Arrays.
  # 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
  @rows = rows
  @column_count = column_count
end
row_vector(row) click to toggle source

Creates a single-row matrix where the values of that row are as given in row.

Matrix.row_vector([4,5,6])
  => 4 5 6
# File lib/matrix.rb, line 169
def Matrix.row_vector(row)
  row = convert_to_array(row)
  new [row]
end
rows(rows, copy = true) click to toggle source

Creates a matrix where rows is an array of arrays, each of which is a row of the matrix. If the optional argument copy is false, use the given arrays as the internal structure of the matrix without copying.

Matrix.rows([[25, 93], [-1, 66]])
   =>  25 93
       -1 66
# File lib/matrix.rb, line 63
def Matrix.rows(rows, copy = true)
  rows = convert_to_array(rows, copy)
  rows.map! do |row|
    convert_to_array(row, copy)
  end
  size = (rows[0] || []).size
  rows.each do |row|
    raise ErrDimensionMismatch, "row size differs (#{row.size} should be #{size})" unless row.size == size
  end
  new rows, size
end
scalar(n, value) click to toggle source

Creates an n by n diagonal matrix where each diagonal element is value.

Matrix.scalar(2, 5)
  => 5 0
     0 5
# File lib/matrix.rb, line 134
def Matrix.scalar(n, value)
  diagonal(*Array.new(n, value))
end
unit(n)
Alias for: identity
vstack(x, *matrices) click to toggle source

Create a matrix by stacking matrices vertically

x = Matrix[[1, 2], [3, 4]]
y = Matrix[[5, 6], [7, 8]]
Matrix.vstack(x, y) # => Matrix[[1, 2], [3, 4], [5, 6], [7, 8]]
# File lib/matrix.rb, line 214
def Matrix.vstack(x, *matrices)
  x = CoercionHelper.coerce_to_matrix(x)
  result = x.send(:rows).map(&:dup)
  matrices.each do |m|
    m = CoercionHelper.coerce_to_matrix(m)
    if m.column_count != x.column_count
      raise ErrDimensionMismatch, "The given matrices must have #{x.column_count} columns, but one has #{m.column_count}"
    end
    result.concat(m.send(:rows))
  end
  new result, x.column_count
end
zero(row_count, column_count = row_count) click to toggle source

Creates a zero matrix.

Matrix.zero(2)
  => 0 0
     0 0
# File lib/matrix.rb, line 158
def Matrix.zero(row_count, column_count = row_count)
  rows = Array.new(row_count){Array.new(column_count, 0)}
  new rows, column_count
end

Public Instance Methods

*(m) click to toggle source

Matrix multiplication.

Matrix[[2,4], [6,8]] * Matrix.identity(2)
  => 2 4
     6 8
# File lib/matrix.rb, line 1019
def *(m) # m is matrix or vector or number
  case(m)
  when Numeric
    rows = @rows.collect {|row|
      row.collect {|e| e * m }
    }
    return new_matrix rows, column_count
  when Vector
    m = self.class.column_vector(m)
    r = self * m
    return r.column(0)
  when Matrix
    Matrix.Raise ErrDimensionMismatch if column_count != m.row_count

    rows = Array.new(row_count) {|i|
      Array.new(m.column_count) {|j|
        (0 ... column_count).inject(0) do |vij, k|
          vij + self[i, k] * m[k, j]
        end
      }
    }
    return new_matrix rows, m.column_count
  else
    return apply_through_coercion(m, __method__)
  end
end
**(other) click to toggle source

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
# File lib/matrix.rb, line 1196
def **(other)
  case other
  when Integer
    x = self
    if other <= 0
      x = self.inverse
      return self.class.identity(self.column_count) if other == 0
      other = -other
    end
    z = nil
    loop do
      z = z ? z * x : x if other[0] == 1
      return z if (other >>= 1).zero?
      x *= x
    end
  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
end
+(m) click to toggle source

Matrix addition.

Matrix.scalar(2,5) + Matrix[[1,0], [-4,7]]
  =>  6  0
     -4 12
# File lib/matrix.rb, line 1052
def +(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "+", self.class, m.class
  when Vector
    m = self.class.column_vector(m)
  when Matrix
  else
    return apply_through_coercion(m, __method__)
  end

  Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count && column_count == m.column_count

  rows = Array.new(row_count) {|i|
    Array.new(column_count) {|j|
      self[i, j] + m[i, j]
    }
  }
  new_matrix rows, column_count
end
+@() click to toggle source
# File lib/matrix.rb, line 1219
def +@
  self
end
-(m) click to toggle source

Matrix subtraction.

Matrix[[1,5], [4,2]] - Matrix[[9,3], [-4,1]]
  => -8  2
      8  1
# File lib/matrix.rb, line 1079
def -(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "-", self.class, m.class
  when Vector
    m = self.class.column_vector(m)
  when Matrix
  else
    return apply_through_coercion(m, __method__)
  end

  Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count && column_count == m.column_count

  rows = Array.new(row_count) {|i|
    Array.new(column_count) {|j|
      self[i, j] - m[i, j]
    }
  }
  new_matrix rows, column_count
end
-@() click to toggle source
# File lib/matrix.rb, line 1223
def -@
  collect {|e| -e }
end
/(other) click to toggle source

Matrix division (multiplication by the inverse).

Matrix[[7,6], [3,9]] / Matrix[[2,9], [3,1]]
  => -7  1
     -3 -6
# File lib/matrix.rb, line 1106
def /(other)
  case other
  when Numeric
    rows = @rows.collect {|row|
      row.collect {|e| e / other }
    }
    return new_matrix rows, column_count
  when Matrix
    return self * other.inverse
  else
    return apply_through_coercion(other, __method__)
  end
end
==(other) click to toggle source

Returns true if and only if the two matrices contain equal elements.

# File lib/matrix.rb, line 982
def ==(other)
  return false unless Matrix === other &&
                      column_count == other.column_count # necessary for empty matrices
  rows == other.rows
end
[](i, j) click to toggle source

Returns element (i,j) of the matrix. That is: row i, column j.

# File lib/matrix.rb, line 299
def [](i, j)
  @rows.fetch(i){return nil}[j]
end
Also aliased as: element, component
matrix[range, range] = matrix/element click to toggle source
matrix[range, integer] = vector/column_matrix/element
matrix[integer, range] = vector/row_matrix/element
matrix[integer, integer] = element

Set element or elements of matrix.

# File lib/matrix.rb, line 313
def []=(i, j, v)
  raise FrozenError, "can't modify frozen Matrix" if frozen?
  rows = check_range(i, :row) or row = check_int(i, :row)
  columns = check_range(j, :column) or column = check_int(j, :column)
  if rows && columns
    set_row_and_col_range(rows, columns, v)
  elsif rows
    set_row_range(rows, column, v)
  elsif columns
    set_col_range(row, columns, v)
  else
    set_value(row, column, v)
  end
end
Also aliased as: set_element, set_component
adjugate() click to toggle source

Returns the adjugate of the matrix.

Matrix[ [7,6],[3,9] ].adjugate
  => 9 -6
     -3 7
# File lib/matrix.rb, line 755
def adjugate
  Matrix.Raise ErrDimensionMismatch unless square?
  Matrix.build(row_count, column_count) do |row, column|
    cofactor(column, row)
  end
end
antisymmetric?() click to toggle source

Returns true if this is an antisymmetric matrix. Raises an error if matrix is not square.

# File lib/matrix.rb, line 934
def antisymmetric?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:upper) do |e, row, col|
    return false unless e == -rows[col][row]
  end
  true
end
check_int(val, direction) click to toggle source
# File lib/matrix.rb, line 338
        def check_int(val, direction)
  count = direction == :row ? row_count : column_count
  CoercionHelper.check_int(val, count, direction)
end
coerce(other) click to toggle source

The coerce method provides support for Ruby type coercion. This coercion mechanism is used by Ruby to handle mixed-type numeric operations: it is intended to find a compatible common type between the two operands of the operator. See also Numeric#coerce.

# File lib/matrix.rb, line 1532
def coerce(other)
  case other
  when Numeric
    return Scalar.new(other), self
  else
    raise TypeError, "#{self.class} can't be coerced into #{other.class}"
  end
end
cofactor(row, column) click to toggle source

Returns the (row, column) cofactor which is obtained by multiplying the first minor by (-1)**(row + column).

Matrix.diagonal(9, 5, -3, 4).cofactor(1, 1)
  => -108
# File lib/matrix.rb, line 740
def cofactor(row, column)
  raise RuntimeError, "cofactor of empty matrix is not defined" if empty?
  Matrix.Raise ErrDimensionMismatch unless square?

  det_of_minor = first_minor(row, column).determinant
  det_of_minor * (-1) ** (row + column)
end
collect(which = :all) { |e| ... } click to toggle source

Returns a matrix that is the result of iteration of the given block over all elements of the matrix. 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] ].collect { |e| e**2 }

    => 1  4
       9 16
# File lib/matrix.rb, line 470
def collect(which = :all, &block) # :yield: e
  return to_enum(:collect, which) unless block_given?
  dup.collect!(which, &block)
end
collect!(which = :all) { |e| ... } click to toggle source

Invokes the given block for each element of matrix, replacing the element with the value returned by the block. 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

# File lib/matrix.rb, line 488
def collect!(which = :all)
  return to_enum(:collect!, which) unless block_given?
  raise FrozenError, "can't modify frozen Matrix" if frozen?
  each_with_index(which){ |e, row_index, col_index| @rows[row_index][col_index] = yield e }
end
column(j) { |e| ... } click to toggle source

Returns column vector number j of the matrix as a Vector (starting at 0 like an array). When a block is given, the elements of that vector are iterated.

# File lib/matrix.rb, line 439
def column(j) # :yield: e
  if block_given?
    return self if j >= column_count || j < -column_count
    row_count.times do |i|
      yield @rows[i][j]
    end
    self
  else
    return nil if j >= column_count || j < -column_count
    col = Array.new(row_count) {|i|
      @rows[i][j]
    }
    Vector.elements(col, false)
  end
end
column_vectors() click to toggle source

Returns an array of the column vectors of the matrix. See Vector.

# File lib/matrix.rb, line 1553
def column_vectors
  Array.new(column_count) {|i|
    column(i)
  }
end
combine(*matrices, &block) click to toggle source
# File lib/matrix.rb, line 277
def combine(*matrices, &block)
  Matrix.combine(self, *matrices, &block)
end
component(i, j)
Alias for: []
conjugate() click to toggle source

Returns the conjugate of the matrix.

Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
  => 1+2i   i  0
        1   2  3
Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].conjugate
  => 1-2i  -i  0
        1   2  3
# File lib/matrix.rb, line 1478
def conjugate
  collect(&:conjugate)
end
determinant() click to toggle source

Returns the determinant of the matrix.

Beware that using Float values can yield erroneous results because of their lack of precision. Consider using exact types like Rational or BigDecimal instead.

Matrix[[7,6], [3,9]].determinant
  => 45
# File lib/matrix.rb, line 1241
def determinant
  Matrix.Raise ErrDimensionMismatch unless square?
  m = @rows
  case row_count
    # Up to 4x4, give result using Laplacian expansion by minors.
    # This will typically be faster, as well as giving good results
    # in case of Floats
  when 0
    +1
  when 1
    + m[0][0]
  when 2
    + m[0][0] * m[1][1] - m[0][1] * m[1][0]
  when 3
    m0, m1, m2 = m
    + m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \
    - m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \
    + m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0]
  when 4
    m0, m1, m2, m3 = m
    + m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \
    - m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \
    + m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \
    - m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \
    + m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \
    - m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \
    + m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \
    - m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \
    + m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \
    - m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \
    + m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \
    - m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0]
  else
    # For bigger matrices, use an efficient and general algorithm.
    # Currently, we use the Gauss-Bareiss algorithm
    determinant_bareiss
  end
end
determinant_bareiss() click to toggle source
# File lib/matrix.rb, line 1292
        def determinant_bareiss
  size = row_count
  last = size - 1
  a = to_a
  no_pivot = Proc.new{ return 0 }
  sign = +1
  pivot = 1
  size.times do |k|
    previous_pivot = pivot
    if (pivot = a[k][k]) == 0
      switch = (k+1 ... size).find(no_pivot) {|row|
        a[row][k] != 0
      }
      a[switch], a[k] = a[k], a[switch]
      pivot = a[k][k]
      sign = -sign
    end
    (k+1).upto(last) do |i|
      ai = a[i]
      (k+1).upto(last) do |j|
        ai[j] =  (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot
      end
    end
  end
  sign * pivot
end
determinant_e() click to toggle source

deprecated; use #determinant

# File lib/matrix.rb, line 1322
def determinant_e
  warn "Matrix#determinant_e is deprecated; use #determinant", uplevel: 1
  determinant
end
diagonal?() click to toggle source

Returns true if this is a diagonal matrix. Raises an error if matrix is not square.

# File lib/matrix.rb, line 801
def diagonal?
  Matrix.Raise ErrDimensionMismatch unless square?
  each(:off_diagonal).all?(&:zero?)
end
each(which = :all) { |e| ... } click to toggle source

Yields all elements of the matrix, starting with those of the first row, or returns an Enumerator if 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
    

    Matrix[ [1,2], [3,4] ].each(:strict_lower).to_a # => [3]

# File lib/matrix.rb, line 517
def each(which = :all) # :yield: e
  return to_enum :each, which unless block_given?
  last = column_count - 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_count.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_count].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
    raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end
each_with_index(which = :all) { |e, row, column| ... } click to toggle source

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
  # => Prints:
  #    1 at 0, 0
  #    2 at 0, 1
  #    3 at 1, 0
  #    4 at 1, 1
# File lib/matrix.rb, line 578
def each_with_index(which = :all) # :yield: e, row, column
  return to_enum :each_with_index, which unless block_given?
  last = column_count - 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_count.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_count].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
    raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end
eigensystem() click to toggle source

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
# File lib/matrix.rb, line 1445
def eigensystem
  EigenvalueDecomposition.new(self)
end
element(i, j)
Alias for: []
elements_to_f() click to toggle source

Deprecated.

Use map(&:to_f)

# File lib/matrix.rb, line 1576
def elements_to_f
  warn "Matrix#elements_to_f is deprecated, use map(&:to_f)", uplevel: 1
  map(&:to_f)
end
elements_to_i() click to toggle source

Deprecated.

Use map(&:to_i)

# File lib/matrix.rb, line 1584
def elements_to_i
  warn "Matrix#elements_to_i is deprecated, use map(&:to_i)", uplevel: 1
  map(&:to_i)
end
elements_to_r() click to toggle source

Deprecated.

Use map(&:to_r)

# File lib/matrix.rb, line 1592
def elements_to_r
  warn "Matrix#elements_to_r is deprecated, use map(&:to_r)", uplevel: 1
  map(&:to_r)
end
empty?() click to toggle source

Returns true if this is an empty matrix, i.e. if the number of rows or the number of columns is 0.

# File lib/matrix.rb, line 810
def empty?
  column_count == 0 || row_count == 0
end
eql?(other) click to toggle source
# File lib/matrix.rb, line 988
def eql?(other)
  return false unless Matrix === other &&
                      column_count == other.column_count # necessary for empty matrices
  rows.eql? other.rows
end
first_minor(row, column) click to toggle source

Returns the submatrix obtained by deleting the specified row and column.

Matrix.diagonal(9, 5, -3, 4).first_minor(1, 2)
  => 9 0 0
     0 0 0
     0 0 4
# File lib/matrix.rb, line 713
def first_minor(row, column)
  raise RuntimeError, "first_minor of empty matrix is not defined" if empty?

  unless 0 <= row && row < row_count
    raise ArgumentError, "invalid row (#{row.inspect} for 0..#{row_count - 1})"
  end

  unless 0 <= column && column < column_count
    raise ArgumentError, "invalid column (#{column.inspect} for 0..#{column_count - 1})"
  end

  arrays = to_a
  arrays.delete_at(row)
  arrays.each do |array|
    array.delete_at(column)
  end

  new_matrix arrays, column_count - 1
end
freeze() click to toggle source
Calls superclass method Object#freeze
# File lib/matrix.rb, line 496
def freeze
  @rows.freeze
  super
end
hadamard_product(m) click to toggle source

Hadamard product

Matrix[[1,2], [3,4]].hadamard_product(Matrix[[1,2], [3,2]])
  => 1  4
     9  8
# File lib/matrix.rb, line 1126
def hadamard_product(m)
  combine(m){|a, b| a * b}
end
hash() click to toggle source

Returns a hash-code for the matrix.

# File lib/matrix.rb, line 1005
def hash
  @rows.hash
end
hermitian?() click to toggle source

Returns true if this is an hermitian matrix. Raises an error if matrix is not square.

# File lib/matrix.rb, line 818
def hermitian?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:upper).all? do |e, row, col|
    e == rows[col][row].conj
  end
end
hstack(*matrices) click to toggle source

Returns a new matrix resulting by stacking horizontally the receiver with the given matrices

x = Matrix[[1, 2], [3, 4]]
y = Matrix[[5, 6], [7, 8]]
x.hstack(y) # => Matrix[[1, 2, 5, 6], [3, 4, 7, 8]]
# File lib/matrix.rb, line 1336
def hstack(*matrices)
  self.class.hstack(self, *matrices)
end
imaginary() click to toggle source

Returns the imaginary part of the matrix.

Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
  => 1+2i  i  0
        1  2  3
Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].imaginary
  =>   2i  i  0
        0  0  0
# File lib/matrix.rb, line 1492
def imaginary
  collect(&:imaginary)
end
index(value, selector = :all) → [row, column] click to toggle source
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]
# File lib/matrix.rb, line 641
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
initialize_copy(m) click to toggle source
Calls superclass method
# File lib/matrix.rb, line 997
        def initialize_copy(m)
  super
  @rows = @rows.map(&:dup) unless frozen?
end
inspect() click to toggle source

Overrides Object#inspect

# File lib/matrix.rb, line 1617
def inspect
  if empty?
    "#{self.class}.empty(#{row_count}, #{column_count})"
  else
    "#{self.class}#{@rows.inspect}"
  end
end
inverse() click to toggle source

Returns the inverse of the matrix.

Matrix[[-1, -1], [0, -1]].inverse
  => -1  1
      0 -1
# File lib/matrix.rb, line 1137
def inverse
  Matrix.Raise ErrDimensionMismatch unless square?
  self.class.I(row_count).send(:inverse_from, self)
end
laplace_expansion(row: nil, column: nil) click to toggle source

Returns the Laplace expansion along given row or column.

Matrix[[7,6], [3,9]].laplace_expansion(column: 1)
 => 45

Matrix[[Vector[1, 0], Vector[0, 1]], [2, 3]].laplace_expansion(row: 0)
 => Vector[3, -2]
# File lib/matrix.rb, line 772
def laplace_expansion(row: nil, column: nil)
  num = row || column

  if !num || (row && column)
    raise ArgumentError, "exactly one the row or column arguments must be specified"
  end

  Matrix.Raise ErrDimensionMismatch unless square?
  raise RuntimeError, "laplace_expansion of empty matrix is not defined" if empty?

  unless 0 <= num && num < row_count
    raise ArgumentError, "invalid num (#{num.inspect} for 0..#{row_count - 1})"
  end

  send(row ? :row : :column, num).map.with_index { |e, k|
    e * cofactor(*(row ? [num, k] : [k,num]))
  }.inject(:+)
end
lower_triangular?() click to toggle source

Returns true if this is a lower triangular matrix.

# File lib/matrix.rb, line 828
def lower_triangular?
  each(:strict_upper).all?(&:zero?)
end
lup() click to toggle source

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 == p * a      # => true
a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)]
# File lib/matrix.rb, line 1460
def lup
  LUPDecomposition.new(self)
end
minor(*param) click to toggle source

Returns a section of the matrix. The parameters are either:

  • start_row, nrows, start_col, ncols; OR

  • row_range, col_range

Matrix.diagonal(9, 5, -3).minor(0..1, 0..2)
  => 9 0 0
     0 5 0

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 is greater than #row_count or column_count respectively.

# File lib/matrix.rb, line 672
def minor(*param)
  case param.size
  when 2
    row_range, col_range = param
    from_row = row_range.first
    from_row += row_count if from_row < 0
    to_row = row_range.end
    to_row += row_count if to_row < 0
    to_row += 1 unless row_range.exclude_end?
    size_row = to_row - from_row

    from_col = col_range.first
    from_col += column_count if from_col < 0
    to_col = col_range.end
    to_col += column_count if to_col < 0
    to_col += 1 unless col_range.exclude_end?
    size_col = to_col - from_col
  when 4
    from_row, size_row, from_col, size_col = param
    return nil if size_row < 0 || size_col < 0
    from_row += row_count if from_row < 0
    from_col += column_count if from_col < 0
  else
    raise ArgumentError, param.inspect
  end

  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|
    row[from_col, size_col]
  }
  new_matrix rows, [column_count - from_col, size_col].min
end
normal?() click to toggle source

Returns true if this is a normal matrix. Raises an error if matrix is not square.

# File lib/matrix.rb, line 836
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
orthogonal?() click to toggle source

Returns true if this is an orthogonal matrix Raises an error if matrix is not square.

# File lib/matrix.rb, line 854
def orthogonal?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row, i|
    column_count.times do |j|
      s = 0
      row_count.times do |k|
        s += row[k] * rows[k][j]
      end
      return false unless s == (i == j ? 1 : 0)
    end
  end
  true
end
permutation?() click to toggle source

Returns true if this is a permutation matrix Raises an error if matrix is not square.

# File lib/matrix.rb, line 872
def permutation?
  Matrix.Raise ErrDimensionMismatch unless square?
  cols = Array.new(column_count)
  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
rank() click to toggle source

Returns the rank of the matrix. Beware that using Float values can yield erroneous results because of their lack of precision. Consider using exact types like Rational or BigDecimal instead.

Matrix[[7,6], [3,9]].rank
  => 2
# File lib/matrix.rb, line 1349
def rank
  # We currently use Bareiss' multistep integer-preserving gaussian elimination
  # (see comments on determinant)
  a = to_a
  last_column = column_count - 1
  last_row = row_count - 1
  pivot_row = 0
  previous_pivot = 1
  0.upto(last_column) do |k|
    switch_row = (pivot_row .. last_row).find {|row|
      a[row][k] != 0
    }
    if switch_row
      a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
      pivot = a[pivot_row][k]
      (pivot_row+1).upto(last_row) do |i|
         ai = a[i]
         (k+1).upto(last_column) do |j|
           ai[j] =  (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
         end
       end
      pivot_row += 1
      previous_pivot = pivot
    end
  end
  pivot_row
end
rank_e() click to toggle source

deprecated; use #rank

# File lib/matrix.rb, line 1380
def rank_e
  warn "Matrix#rank_e is deprecated; use #rank", uplevel: 1
  rank
end
real() click to toggle source

Returns the real part of the matrix.

Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
  => 1+2i  i  0
        1  2  3
Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].real
  =>    1  0  0
        1  2  3
# File lib/matrix.rb, line 1506
def real
  collect(&:real)
end
real?() click to toggle source

Returns true if all entries of the matrix are real.

# File lib/matrix.rb, line 893
def real?
  all?(&:real?)
end
rect() click to toggle source

Returns an array containing matrices corresponding to the real and imaginary parts of the matrix

m.rect == [m.real, m.imag] # ==> true for all matrices m

# File lib/matrix.rb, line 1516
def rect
  [real, imag]
end
regular?() click to toggle source

Returns true if this is a regular (i.e. non-singular) matrix.

# File lib/matrix.rb, line 900
def regular?
  not singular?
end
round(ndigits=0) click to toggle source

Returns a matrix with entries rounded to the given precision (see Float#round)

# File lib/matrix.rb, line 1388
def round(ndigits=0)
  map{|e| e.round(ndigits)}
end
row(i) { |e| ... } click to toggle source

Returns row vector number i of the matrix as a Vector (starting at 0 like an array). When a block is given, the elements of that vector are iterated.

# File lib/matrix.rb, line 425
def row(i, &block) # :yield: e
  if block_given?
    @rows.fetch(i){return self}.each(&block)
    self
  else
    Vector.elements(@rows.fetch(i){return nil})
  end
end
row_count() click to toggle source

Returns the number of rows.

# File lib/matrix.rb, line 410
def row_count
  @rows.size
end
row_vectors() click to toggle source

Returns an array of the row vectors of the matrix. See Vector.

# File lib/matrix.rb, line 1544
def row_vectors
  Array.new(row_count) {|i|
    row(i)
  }
end
set_col_range(row, col_range, value) click to toggle source
# File lib/matrix.rb, line 394
        def set_col_range(row, col_range, value)
  value = if value.is_a?(Vector)
    value.to_a
  elsif value.is_a?(Matrix)
    Matrix.Raise ErrDimensionMismatch unless value.row_count == 1
    value.row(0).to_a
  else
    Array.new(col_range.size, value)
  end
  Matrix.Raise ErrDimensionMismatch unless col_range.size == value.size
  @rows[row][col_range] = value
end
set_column_vector(row_range, col, value) click to toggle source
# File lib/matrix.rb, line 387
        def set_column_vector(row_range, col, value)
  value.each_with_index do |e, index|
    r = row_range.begin + index
    @rows[r][col] = e
  end
end
set_row_and_col_range(row_range, col_range, value) click to toggle source
# File lib/matrix.rb, line 349
        def set_row_and_col_range(row_range, col_range, value)
  if value.is_a?(Matrix)
    if row_range.size != value.row_count || col_range.size != value.column_count
      raise ErrDimensionMismatch, [
        'Expected a Matrix of dimensions',
        "#{row_range.size}x#{col_range.size}",
        'got',
        "#{value.row_count}x#{value.column_count}",
      ].join(' ')
    end
    source = value.instance_variable_get :@rows
    row_range.each_with_index do |row, i|
      @rows[row][col_range] = source[i]
    end
  elsif value.is_a?(Vector)
    raise ErrDimensionMismatch, 'Expected a Matrix or a value, got a Vector'
  else
    value_to_set = Array.new(col_range.size, value)
    row_range.each do |i|
      @rows[i][col_range] = value_to_set
    end
  end
end
set_row_range(row_range, col, value) click to toggle source
# File lib/matrix.rb, line 373
        def set_row_range(row_range, col, value)
  if value.is_a?(Vector)
    Matrix.Raise ErrDimensionMismatch unless row_range.size == value.size
    set_column_vector(row_range, col, value)
  elsif value.is_a?(Matrix)
    Matrix.Raise ErrDimensionMismatch unless value.column_count == 1
    value = value.column(0)
    Matrix.Raise ErrDimensionMismatch unless row_range.size == value.size
    set_column_vector(row_range, col, value)
  else
    @rows[row_range].each{|e| e[col] = value }
  end
end
set_value(row, col, value) click to toggle source
# File lib/matrix.rb, line 343
        def set_value(row, col, value)
  raise ErrDimensionMismatch, "Expected a a value, got a #{value.class}" if value.respond_to?(:to_matrix)

  @rows[row][col] = value
end
singular?() click to toggle source

Returns true if this is a singular matrix.

# File lib/matrix.rb, line 907
def singular?
  determinant == 0
end
square?() click to toggle source

Returns true if this is a square matrix.

# File lib/matrix.rb, line 914
def square?
  column_count == row_count
end
symmetric?() click to toggle source

Returns true if this is a symmetric matrix. Raises an error if matrix is not square.

# File lib/matrix.rb, line 922
def symmetric?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:strict_upper) do |e, row, col|
    return false if e != rows[col][row]
  end
  true
end
to_a() click to toggle source

Returns an array of arrays that describe the rows of the matrix.

# File lib/matrix.rb, line 1569
def to_a
  @rows.collect(&:dup)
end
to_matrix() click to toggle source

Explicit conversion to a Matrix. Returns self

# File lib/matrix.rb, line 1562
def to_matrix
  self
end
to_s() click to toggle source

Overrides Object#to_s

# File lib/matrix.rb, line 1604
def to_s
  if empty?
    "#{self.class}.empty(#{row_count}, #{column_count})"
  else
    "#{self.class}[" + @rows.collect{|row|
      "[" + row.collect{|e| e.to_s}.join(", ") + "]"
    }.join(", ")+"]"
  end
end
trace() click to toggle source

Returns the trace (sum of diagonal elements) of the matrix.

Matrix[[7,6], [3,9]].trace
  => 16
# File lib/matrix.rb, line 1397
def trace
  Matrix.Raise ErrDimensionMismatch unless square?
  (0...column_count).inject(0) do |tr, i|
    tr + @rows[i][i]
  end
end
transpose() click to toggle source

Returns the transpose of the matrix.

Matrix[[1,2], [3,4], [5,6]]
  => 1 2
     3 4
     5 6
Matrix[[1,2], [3,4], [5,6]].transpose
  => 1 3 5
     2 4 6
# File lib/matrix.rb, line 1415
def transpose
  return self.class.empty(column_count, 0) if row_count.zero?
  new_matrix @rows.transpose, row_count
end
unitary?() click to toggle source

Returns true if this is a unitary matrix Raises an error if matrix is not square.

# File lib/matrix.rb, line 947
def unitary?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row, i|
    column_count.times do |j|
      s = 0
      row_count.times do |k|
        s += row[k].conj * rows[k][j]
      end
      return false unless s == (i == j ? 1 : 0)
    end
  end
  true
end
upper_triangular?() click to toggle source

Returns true if this is an upper triangular matrix.

# File lib/matrix.rb, line 964
def upper_triangular?
  each(:strict_lower).all?(&:zero?)
end
vstack(*matrices) click to toggle source

Returns a new matrix resulting by stacking vertically the receiver with the given matrices

x = Matrix[[1, 2], [3, 4]]
y = Matrix[[5, 6], [7, 8]]
x.vstack(y) # => Matrix[[1, 2], [3, 4], [5, 6], [7, 8]]
# File lib/matrix.rb, line 1429
def vstack(*matrices)
  self.class.vstack(self, *matrices)
end
zero?() click to toggle source

Returns true if this is a matrix with only zero elements

# File lib/matrix.rb, line 971
def zero?
  all?(&:zero?)
end

Private Instance Methods

check_range(val, direction) click to toggle source

Returns range or nil

# File lib/matrix.rb, line 332
          def check_range(val, direction)
    return unless val.is_a?(Range)
    count = direction == :row ? row_count : column_count
    CoercionHelper.check_range(val, count, direction)
  end

  private def check_int(val, direction)
    count = direction == :row ? row_count : column_count
    CoercionHelper.check_int(val, count, direction)
  end

  private def set_value(row, col, value)
    raise ErrDimensionMismatch, "Expected a a value, got a #{value.class}" if value.respond_to?(:to_matrix)

    @rows[row][col] = value
  end

  private def set_row_and_col_range(row_range, col_range, value)
    if value.is_a?(Matrix)
      if row_range.size != value.row_count || col_range.size != value.column_count
        raise ErrDimensionMismatch, [
          'Expected a Matrix of dimensions',
          "#{row_range.size}x#{col_range.size}",
          'got',
          "#{value.row_count}x#{value.column_count}",
        ].join(' ')
      end
      source = value.instance_variable_get :@rows
      row_range.each_with_index do |row, i|
        @rows[row][col_range] = source[i]
      end
    elsif value.is_a?(Vector)
      raise ErrDimensionMismatch, 'Expected a Matrix or a value, got a Vector'
    else
      value_to_set = Array.new(col_range.size, value)
      row_range.each do |i|
        @rows[i][col_range] = value_to_set
      end
    end
  end

  private def set_row_range(row_range, col, value)
    if value.is_a?(Vector)
      Matrix.Raise ErrDimensionMismatch unless row_range.size == value.size
      set_column_vector(row_range, col, value)
    elsif value.is_a?(Matrix)
      Matrix.Raise ErrDimensionMismatch unless value.column_count == 1
      value = value.column(0)
      Matrix.Raise ErrDimensionMismatch unless row_range.size == value.size
      set_column_vector(row_range, col, value)
    else
      @rows[row_range].each{|e| e[col] = value }
    end
  end

  private def set_column_vector(row_range, col, value)
    value.each_with_index do |e, index|
      r = row_range.begin + index
      @rows[r][col] = e
    end
  end

  private def set_col_range(row, col_range, value)
    value = if value.is_a?(Vector)
      value.to_a
    elsif value.is_a?(Matrix)
      Matrix.Raise ErrDimensionMismatch unless value.row_count == 1
      value.row(0).to_a
    else
      Array.new(col_range.size, value)
    end
    Matrix.Raise ErrDimensionMismatch unless col_range.size == value.size
    @rows[row][col_range] = value
  end

  #
  # Returns the number of rows.
  #
  def row_count
    @rows.size
  end

  alias_method :row_size, :row_count
  #
  # Returns the number of columns.
  #
  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
  # an array).  When a block is given, the elements of that vector are iterated.
  #
  def row(i, &block) # :yield: e
    if block_given?
      @rows.fetch(i){return self}.each(&block)
      self
    else
      Vector.elements(@rows.fetch(i){return nil})
    end
  end

  #
  # Returns column vector number +j+ of the matrix as a Vector (starting at 0
  # like an array).  When a block is given, the elements of that vector are
  # iterated.
  #
  def column(j) # :yield: e
    if block_given?
      return self if j >= column_count || j < -column_count
      row_count.times do |i|
        yield @rows[i][j]
      end
      self
    else
      return nil if j >= column_count || j < -column_count
      col = Array.new(row_count) {|i|
        @rows[i][j]
      }
      Vector.elements(col, false)
    end
  end

  #
  # Returns a matrix that is the result of iteration of the given block over all
  # elements of the matrix.
  # 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] ].collect { |e| e**2 }
  #     => 1  4
  #        9 16
  #
  def collect(which = :all, &block) # :yield: e
    return to_enum(:collect, which) unless block_given?
    dup.collect!(which, &block)
  end
  alias_method :map, :collect

  #
  # Invokes the given block for each element of matrix, replacing the element with the value
  # returned by the block.
  # 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
  #
  def collect!(which = :all)
    return to_enum(:collect!, which) unless block_given?
    raise FrozenError, "can't modify frozen Matrix" if frozen?
    each_with_index(which){ |e, row_index, col_index| @rows[row_index][col_index] = yield e }
  end

  alias map! collect!

  def freeze
    @rows.freeze
    super
  end

  #
  # Yields all elements of the matrix, starting with those of the first row,
  # or returns an Enumerator if 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
  #   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_count - 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_count.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_count].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
      raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
    end
    self
  end

  #
  # 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
  #     # => 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_count - 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_count.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_count].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
      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
  # *  row_range, col_range
  #
  #   Matrix.diagonal(9, 5, -3).minor(0..1, 0..2)
  #     => 9 0 0
  #        0 5 0
  #
  # 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 is greater than row_count or column_count respectively.
  #
  def minor(*param)
    case param.size
    when 2
      row_range, col_range = param
      from_row = row_range.first
      from_row += row_count if from_row < 0
      to_row = row_range.end
      to_row += row_count if to_row < 0
      to_row += 1 unless row_range.exclude_end?
      size_row = to_row - from_row

      from_col = col_range.first
      from_col += column_count if from_col < 0
      to_col = col_range.end
      to_col += column_count if to_col < 0
      to_col += 1 unless col_range.exclude_end?
      size_col = to_col - from_col
    when 4
      from_row, size_row, from_col, size_col = param
      return nil if size_row < 0 || size_col < 0
      from_row += row_count if from_row < 0
      from_col += column_count if from_col < 0
    else
      raise ArgumentError, param.inspect
    end

    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|
      row[from_col, size_col]
    }
    new_matrix rows, [column_count - from_col, size_col].min
  end

  #
  # Returns the submatrix obtained by deleting the specified row and column.
  #
  #   Matrix.diagonal(9, 5, -3, 4).first_minor(1, 2)
  #     => 9 0 0
  #        0 0 0
  #        0 0 4
  #
  def first_minor(row, column)
    raise RuntimeError, "first_minor of empty matrix is not defined" if empty?

    unless 0 <= row && row < row_count
      raise ArgumentError, "invalid row (#{row.inspect} for 0..#{row_count - 1})"
    end

    unless 0 <= column && column < column_count
      raise ArgumentError, "invalid column (#{column.inspect} for 0..#{column_count - 1})"
    end

    arrays = to_a
    arrays.delete_at(row)
    arrays.each do |array|
      array.delete_at(column)
    end

    new_matrix arrays, column_count - 1
  end

  #
  # Returns the (row, column) cofactor which is obtained by multiplying
  # the first minor by (-1)**(row + column).
  #
  #   Matrix.diagonal(9, 5, -3, 4).cofactor(1, 1)
  #     => -108
  #
  def cofactor(row, column)
    raise RuntimeError, "cofactor of empty matrix is not defined" if empty?
    Matrix.Raise ErrDimensionMismatch unless square?

    det_of_minor = first_minor(row, column).determinant
    det_of_minor * (-1) ** (row + column)
  end

  #
  # Returns the adjugate of the matrix.
  #
  #   Matrix[ [7,6],[3,9] ].adjugate
  #     => 9 -6
  #        -3 7
  #
  def adjugate
    Matrix.Raise ErrDimensionMismatch unless square?
    Matrix.build(row_count, column_count) do |row, column|
      cofactor(column, row)
    end
  end

  #
  # Returns the Laplace expansion along given row or column.
  #
  #    Matrix[[7,6], [3,9]].laplace_expansion(column: 1)
  #     => 45
  #
  #    Matrix[[Vector[1, 0], Vector[0, 1]], [2, 3]].laplace_expansion(row: 0)
  #     => Vector[3, -2]
  #
  #
  def laplace_expansion(row: nil, column: nil)
    num = row || column

    if !num || (row && column)
      raise ArgumentError, "exactly one the row or column arguments must be specified"
    end

    Matrix.Raise ErrDimensionMismatch unless square?
    raise RuntimeError, "laplace_expansion of empty matrix is not defined" if empty?

    unless 0 <= num && num < row_count
      raise ArgumentError, "invalid num (#{num.inspect} for 0..#{row_count - 1})"
    end

    send(row ? :row : :column, num).map.with_index { |e, k|
      e * cofactor(*(row ? [num, k] : [k,num]))
    }.inject(:+)
  end
  alias_method :cofactor_expansion, :laplace_expansion


  #--
  # TESTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
  #++

  #
  # Returns +true+ if 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.
  #
  def empty?
    column_count == 0 || row_count == 0
  end

  #
  # Returns +true+ if this is an hermitian matrix.
  # Raises an error if matrix is not square.
  #
  def hermitian?
    Matrix.Raise ErrDimensionMismatch unless square?
    each_with_index(:upper).all? do |e, row, col|
      e == rows[col][row].conj
    end
  end

  #
  # Returns +true+ if this is a lower triangular matrix.
  #
  def lower_triangular?
    each(:strict_upper).all?(&:zero?)
  end

  #
  # Returns +true+ if 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+ if 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_count.times do |j|
        s = 0
        row_count.times do |k|
          s += row[k] * rows[k][j]
        end
        return false unless s == (i == j ? 1 : 0)
      end
    end
    true
  end

  #
  # Returns +true+ if this is a permutation matrix
  # Raises an error if matrix is not square.
  #
  def permutation?
    Matrix.Raise ErrDimensionMismatch unless square?
    cols = Array.new(column_count)
    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?
    all?(&:real?)
  end

  #
  # Returns +true+ if this is a regular (i.e. non-singular) matrix.
  #
  def regular?
    not singular?
  end

  #
  # Returns +true+ if this is a singular matrix.
  #
  def singular?
    determinant == 0
  end

  #
  # Returns +true+ if this is a square matrix.
  #
  def square?
    column_count == row_count
  end

  #
  # Returns +true+ if 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) do |e, row, col|
      return false if e != rows[col][row]
    end
    true
  end

  #
  # Returns +true+ if this is an antisymmetric matrix.
  # Raises an error if matrix is not square.
  #
  def antisymmetric?
    Matrix.Raise ErrDimensionMismatch unless square?
    each_with_index(:upper) do |e, row, col|
      return false unless e == -rows[col][row]
    end
    true
  end
  alias_method :skew_symmetric?, :antisymmetric?

  #
  # Returns +true+ if 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_count.times do |j|
        s = 0
        row_count.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+ if this is an upper triangular matrix.
  #
  def upper_triangular?
    each(:strict_lower).all?(&:zero?)
  end

  #
  # Returns +true+ if this is a matrix with only zero elements
  #
  def zero?
    all?(&:zero?)
  end

  #--
  # OBJECT METHODS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
  #++

  #
  # Returns +true+ if and only if the two matrices contain equal elements.
  #
  def ==(other)
    return false unless Matrix === other &&
                        column_count == other.column_count # necessary for empty matrices
    rows == other.rows
  end

  def eql?(other)
    return false unless Matrix === other &&
                        column_count == other.column_count # necessary for empty matrices
    rows.eql? other.rows
  end

  #
  # Called for dup & clone.
  #
  private def initialize_copy(m)
    super
    @rows = @rows.map(&:dup) unless frozen?
  end

  #
  # Returns a hash-code for the matrix.
  #
  def hash
    @rows.hash
  end

  #--
  # ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
  #++

  #
  # Matrix multiplication.
  #   Matrix[[2,4], [6,8]] * Matrix.identity(2)
  #     => 2 4
  #        6 8
  #
  def *(m) # m is matrix or vector or number
    case(m)
    when Numeric
      rows = @rows.collect {|row|
        row.collect {|e| e * m }
      }
      return new_matrix rows, column_count
    when Vector
      m = self.class.column_vector(m)
      r = self * m
      return r.column(0)
    when Matrix
      Matrix.Raise ErrDimensionMismatch if column_count != m.row_count

      rows = Array.new(row_count) {|i|
        Array.new(m.column_count) {|j|
          (0 ... column_count).inject(0) do |vij, k|
            vij + self[i, k] * m[k, j]
          end
        }
      }
      return new_matrix rows, m.column_count
    else
      return apply_through_coercion(m, __method__)
    end
  end

  #
  # Matrix addition.
  #   Matrix.scalar(2,5) + Matrix[[1,0], [-4,7]]
  #     =>  6  0
  #        -4 12
  #
  def +(m)
    case m
    when Numeric
      Matrix.Raise ErrOperationNotDefined, "+", self.class, m.class
    when Vector
      m = self.class.column_vector(m)
    when Matrix
    else
      return apply_through_coercion(m, __method__)
    end

    Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count && column_count == m.column_count

    rows = Array.new(row_count) {|i|
      Array.new(column_count) {|j|
        self[i, j] + m[i, j]
      }
    }
    new_matrix rows, column_count
  end

  #
  # Matrix subtraction.
  #   Matrix[[1,5], [4,2]] - Matrix[[9,3], [-4,1]]
  #     => -8  2
  #         8  1
  #
  def -(m)
    case m
    when Numeric
      Matrix.Raise ErrOperationNotDefined, "-", self.class, m.class
    when Vector
      m = self.class.column_vector(m)
    when Matrix
    else
      return apply_through_coercion(m, __method__)
    end

    Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count && column_count == m.column_count

    rows = Array.new(row_count) {|i|
      Array.new(column_count) {|j|
        self[i, j] - m[i, j]
      }
    }
    new_matrix rows, column_count
  end

  #
  # Matrix division (multiplication by the inverse).
  #   Matrix[[7,6], [3,9]] / Matrix[[2,9], [3,1]]
  #     => -7  1
  #        -3 -6
  #
  def /(other)
    case other
    when Numeric
      rows = @rows.collect {|row|
        row.collect {|e| e / other }
      }
      return new_matrix rows, column_count
    when Matrix
      return self * other.inverse
    else
      return apply_through_coercion(other, __method__)
    end
  end

  #
  # Hadamard product
  #    Matrix[[1,2], [3,4]].hadamard_product(Matrix[[1,2], [3,2]])
  #      => 1  4
  #         9  8
  #
  def hadamard_product(m)
    combine(m){|a, b| a * b}
  end
  alias_method :entrywise_product, :hadamard_product

  #
  # Returns the inverse of the matrix.
  #   Matrix[[-1, -1], [0, -1]].inverse
  #     => -1  1
  #         0 -1
  #
  def inverse
    Matrix.Raise ErrDimensionMismatch unless square?
    self.class.I(row_count).send(:inverse_from, self)
  end
  alias_method :inv, :inverse

  private def inverse_from(src) # :nodoc:
    last = row_count - 1
    a = src.to_a

    0.upto(last) do |k|
      i = k
      akk = a[k][k].abs
      (k+1).upto(last) do |j|
        v = a[j][k].abs
        if v > akk
          i = j
          akk = v
        end
      end
      Matrix.Raise ErrNotRegular if akk == 0
      if i != k
        a[i], a[k] = a[k], a[i]
        @rows[i], @rows[k] = @rows[k], @rows[i]
      end
      akk = a[k][k]

      0.upto(last) do |ii|
        next if ii == k
        q = a[ii][k].quo(akk)
        a[ii][k] = 0

        (k + 1).upto(last) do |j|
          a[ii][j] -= a[k][j] * q
        end
        0.upto(last) do |j|
          @rows[ii][j] -= @rows[k][j] * q
        end
      end

      (k+1).upto(last) do |j|
        a[k][j] = a[k][j].quo(akk)
      end
      0.upto(last) do |j|
        @rows[k][j] = @rows[k][j].quo(akk)
      end
    end
    self
  end

  #
  # 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
  #
  def **(other)
    case other
    when Integer
      x = self
      if other <= 0
        x = self.inverse
        return self.class.identity(self.column_count) if other == 0
        other = -other
      end
      z = nil
      loop do
        z = z ? z * x : x if other[0] == 1
        return z if (other >>= 1).zero?
        x *= x
      end
    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
  end

  def +@
    self
  end

  def -@
    collect {|e| -e }
  end

  #--
  # MATRIX FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
  #++

  #
  # Returns the determinant of the matrix.
  #
  # Beware that using Float values can yield erroneous results
  # because of their lack of precision.
  # Consider using exact types like Rational or BigDecimal instead.
  #
  #   Matrix[[7,6], [3,9]].determinant
  #     => 45
  #
  def determinant
    Matrix.Raise ErrDimensionMismatch unless square?
    m = @rows
    case row_count
      # Up to 4x4, give result using Laplacian expansion by minors.
      # This will typically be faster, as well as giving good results
      # in case of Floats
    when 0
      +1
    when 1
      + m[0][0]
    when 2
      + m[0][0] * m[1][1] - m[0][1] * m[1][0]
    when 3
      m0, m1, m2 = m
      + m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \
      - m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \
      + m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0]
    when 4
      m0, m1, m2, m3 = m
      + m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \
      - m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \
      + m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \
      - m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \
      + m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \
      - m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \
      + m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \
      - m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \
      + m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \
      - m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \
      + m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \
      - m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0]
    else
      # For bigger matrices, use an efficient and general algorithm.
      # Currently, we use the Gauss-Bareiss algorithm
      determinant_bareiss
    end
  end
  alias_method :det, :determinant

  #
  # Private. Use Matrix#determinant
  #
  # Returns the determinant of the matrix, using
  # Bareiss' multistep integer-preserving gaussian elimination.
  # It has the same computational cost order O(n^3) as standard Gaussian elimination.
  # Intermediate results are fraction free and of lower complexity.
  # A matrix of Integers will have thus intermediate results that are also Integers,
  # with smaller bignums (if any), while a matrix of Float will usually have
  # intermediate results with better precision.
  #
  private def determinant_bareiss
    size = row_count
    last = size - 1
    a = to_a
    no_pivot = Proc.new{ return 0 }
    sign = +1
    pivot = 1
    size.times do |k|
      previous_pivot = pivot
      if (pivot = a[k][k]) == 0
        switch = (k+1 ... size).find(no_pivot) {|row|
          a[row][k] != 0
        }
        a[switch], a[k] = a[k], a[switch]
        pivot = a[k][k]
        sign = -sign
      end
      (k+1).upto(last) do |i|
        ai = a[i]
        (k+1).upto(last) do |j|
          ai[j] =  (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot
        end
      end
    end
    sign * pivot
  end

  #
  # deprecated; use Matrix#determinant
  #
  def determinant_e
    warn "Matrix#determinant_e is deprecated; use #determinant", uplevel: 1
    determinant
  end
  alias_method :det_e, :determinant_e

  #
  # Returns a new matrix resulting by stacking horizontally
  # the receiver with the given matrices
  #
  #   x = Matrix[[1, 2], [3, 4]]
  #   y = Matrix[[5, 6], [7, 8]]
  #   x.hstack(y) # => Matrix[[1, 2, 5, 6], [3, 4, 7, 8]]
  #
  def hstack(*matrices)
    self.class.hstack(self, *matrices)
  end

  #
  # Returns the rank of the matrix.
  # Beware that using Float values can yield erroneous results
  # because of their lack of precision.
  # Consider using exact types like Rational or BigDecimal instead.
  #
  #   Matrix[[7,6], [3,9]].rank
  #     => 2
  #
  def rank
    # We currently use Bareiss' multistep integer-preserving gaussian elimination
    # (see comments on determinant)
    a = to_a
    last_column = column_count - 1
    last_row = row_count - 1
    pivot_row = 0
    previous_pivot = 1
    0.upto(last_column) do |k|
      switch_row = (pivot_row .. last_row).find {|row|
        a[row][k] != 0
      }
      if switch_row
        a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
        pivot = a[pivot_row][k]
        (pivot_row+1).upto(last_row) do |i|
           ai = a[i]
           (k+1).upto(last_column) do |j|
             ai[j] =  (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
           end
         end
        pivot_row += 1
        previous_pivot = pivot
      end
    end
    pivot_row
  end

  #
  # deprecated; use Matrix#rank
  #
  def rank_e
    warn "Matrix#rank_e is deprecated; use #rank", uplevel: 1
    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.
  #   Matrix[[7,6], [3,9]].trace
  #     => 16
  #
  def trace
    Matrix.Raise ErrDimensionMismatch unless square?
    (0...column_count).inject(0) do |tr, i|
      tr + @rows[i][i]
    end
  end
  alias_method :tr, :trace

  #
  # Returns the transpose of the matrix.
  #   Matrix[[1,2], [3,4], [5,6]]
  #     => 1 2
  #        3 4
  #        5 6
  #   Matrix[[1,2], [3,4], [5,6]].transpose
  #     => 1 3 5
  #        2 4 6
  #
  def transpose
    return self.class.empty(column_count, 0) if row_count.zero?
    new_matrix @rows.transpose, row_count
  end
  alias_method :t, :transpose

  #
  # Returns a new matrix resulting by stacking vertically
  # the receiver with the given matrices
  #
  #   x = Matrix[[1, 2], [3, 4]]
  #   y = Matrix[[5, 6], [7, 8]]
  #   x.vstack(y) # => Matrix[[1, 2], [3, 4], [5, 6], [7, 8]]
  #
  def vstack(*matrices)
    self.class.vstack(self, *matrices)
  end

  #--
  # 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_method :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 == p * a      # => true
  #   a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)]
  #
  def lup
    LUPDecomposition.new(self)
  end
  alias_method :lup_decomposition, :lup

  #--
  # COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  #++

  #
  # Returns the conjugate of the matrix.
  #   Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
  #     => 1+2i   i  0
  #           1   2  3
  #   Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].conjugate
  #     => 1-2i  -i  0
  #           1   2  3
  #
  def conjugate
    collect(&:conjugate)
  end
  alias_method :conj, :conjugate

  #
  # Returns the imaginary part of the matrix.
  #   Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
  #     => 1+2i  i  0
  #           1  2  3
  #   Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].imaginary
  #     =>   2i  i  0
  #           0  0  0
  #
  def imaginary
    collect(&:imaginary)
  end
  alias_method :imag, :imaginary

  #
  # Returns the real part of the matrix.
  #   Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]]
  #     => 1+2i  i  0
  #           1  2  3
  #   Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].real
  #     =>    1  0  0
  #           1  2  3
  #
  def real
    collect(&:real)
  end

  #
  # Returns an array containing matrices corresponding to the real and imaginary
  # parts of the matrix
  #
  # m.rect == [m.real, m.imag]  # ==> true for all matrices m
  #
  def rect
    [real, imag]
  end
  alias_method :rectangular, :rect

  #--
  # CONVERTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
  #++

  #
  # The coerce method provides support for Ruby type coercion.
  # This coercion mechanism is used by Ruby to handle mixed-type
  # numeric operations: it is intended to find a compatible common
  # type between the two operands of the operator.
  # See also Numeric#coerce.
  #
  def coerce(other)
    case other
    when Numeric
      return Scalar.new(other), self
    else
      raise TypeError, "#{self.class} can't be coerced into #{other.class}"
    end
  end

  #
  # Returns an array of the row vectors of the matrix.  See Vector.
  #
  def row_vectors
    Array.new(row_count) {|i|
      row(i)
    }
  end

  #
  # Returns an array of the column vectors of the matrix.  See Vector.
  #
  def column_vectors
    Array.new(column_count) {|i|
      column(i)
    }
  end

  #
  # Explicit conversion to a Matrix. Returns self
  #
  def to_matrix
    self
  end

  #
  # Returns an array of arrays that describe the rows of the matrix.
  #
  def to_a
    @rows.collect(&:dup)
  end

  # Deprecated.
  #
  # Use map(&:to_f)
  def elements_to_f
    warn "Matrix#elements_to_f is deprecated, use map(&:to_f)", uplevel: 1
    map(&:to_f)
  end

  # Deprecated.
  #
  # Use map(&:to_i)
  def elements_to_i
    warn "Matrix#elements_to_i is deprecated, use map(&:to_i)", uplevel: 1
    map(&:to_i)
  end

  # Deprecated.
  #
  # Use map(&:to_r)
  def elements_to_r
    warn "Matrix#elements_to_r is deprecated, use map(&:to_r)", uplevel: 1
    map(&:to_r)
  end

  #--
  # PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
  #++

  #
  # Overrides Object#to_s
  #
  def to_s
    if empty?
      "#{self.class}.empty(#{row_count}, #{column_count})"
    else
      "#{self.class}[" + @rows.collect{|row|
        "[" + row.collect{|e| e.to_s}.join(", ") + "]"
      }.join(", ")+"]"
    end
  end

  #
  # Overrides Object#inspect
  #
  def inspect
    if empty?
      "#{self.class}.empty(#{row_count}, #{column_count})"
    else
      "#{self.class}#{@rows.inspect}"
    end
  end

  # Private helper modules

  module ConversionHelper # :nodoc:
    #
    # Converts the obj to an Array. If copy is set to true
    # a copy of obj will be made if necessary.
    #
    private def convert_to_array(obj, copy = false) # :nodoc:
      case obj
      when Array
        copy ? obj.dup : obj
      when Vector
        obj.to_a
      else
        begin
          converted = obj.to_ary
        rescue Exception => e
          raise TypeError, "can't convert #{obj.class} into an Array (#{e.message})"
        end
        raise TypeError, "#{obj.class}#to_ary should return an Array" unless converted.is_a? Array
        converted
      end
    end
  end

  extend ConversionHelper

  module CoercionHelper # :nodoc:
    #
    # Applies the operator +oper+ with argument +obj+
    # through coercion of +obj+
    #
    private def apply_through_coercion(obj, oper)
      coercion = obj.coerce(self)
      raise TypeError unless coercion.is_a?(Array) && coercion.length == 2
      coercion[0].public_send(oper, coercion[1])
    rescue
      raise TypeError, "#{obj.inspect} can't be coerced into #{self.class}"
    end

    #
    # Helper method to coerce a value into a specific class.
    # Raises a TypeError if the coercion fails or the returned value
    # is not of the right class.
    # (from Rubinius)
    #
    def self.coerce_to(obj, cls, meth) # :nodoc:
      return obj if obj.kind_of?(cls)
      raise TypeError, "Expected a #{cls} but got a #{obj.class}" unless obj.respond_to? meth
      begin
        ret = obj.__send__(meth)
      rescue Exception => e
        raise TypeError, "Coercion error: #{obj.inspect}.#{meth} => #{cls} failed:\n" \
                         "(#{e.message})"
      end
      raise TypeError, "Coercion error: obj.#{meth} did NOT return a #{cls} (was #{ret.class})" unless ret.kind_of? cls
      ret
    end

    def self.coerce_to_int(obj)
      coerce_to(obj, Integer, :to_int)
    end

    def self.coerce_to_matrix(obj)
      coerce_to(obj, Matrix, :to_matrix)
    end

    # Returns `nil` for non Ranges
    # Checks range validity, return canonical range with 0 <= begin <= end < count
    def self.check_range(val, count, kind)
      canonical = (val.begin + (val.begin < 0 ? count : 0))..
                  (val.end ? val.end + (val.end < 0 ? count : 0) - (val.exclude_end? ? 1 : 0)
                           : count - 1)
      unless 0 <= canonical.begin && canonical.begin <= canonical.end && canonical.end < count
        raise IndexError, "given range #{val} is outside of #{kind} dimensions: 0...#{count}"
      end
      canonical
    end

    def self.check_int(val, count, kind)
      val = CoercionHelper.coerce_to_int(val)
      if val >= count || val < -count
        raise IndexError, "given #{kind} #{val} is outside of #{-count}...#{count}"
      end
      val
    end
  end

  include CoercionHelper

  # Private CLASS

  class Scalar < Numeric # :nodoc:
    include ExceptionForMatrix
    include CoercionHelper

    def initialize(value)
      @value = value
    end

    # ARITHMETIC
    def +(other)
      case other
      when Numeric
        Scalar.new(@value + other)
      when Vector, Matrix
        Scalar.Raise ErrOperationNotDefined, "+", @value.class, other.class
      else
        apply_through_coercion(other, __method__)
      end
    end

    def -(other)
      case other
      when Numeric
        Scalar.new(@value - other)
      when Vector, Matrix
        Scalar.Raise ErrOperationNotDefined, "-", @value.class, other.class
      else
        apply_through_coercion(other, __method__)
      end
    end

    def *(other)
      case other
      when Numeric
        Scalar.new(@value * other)
      when Vector, Matrix
        other.collect{|e| @value * e}
      else
        apply_through_coercion(other, __method__)
      end
    end

    def /(other)
      case other
      when Numeric
        Scalar.new(@value / other)
      when Vector
        Scalar.Raise ErrOperationNotDefined, "/", @value.class, other.class
      when Matrix
        self * other.inverse
      else
        apply_through_coercion(other, __method__)
      end
    end

    def **(other)
      case other
      when Numeric
        Scalar.new(@value ** other)
      when Vector
        Scalar.Raise ErrOperationNotDefined, "**", @value.class, other.class
      when Matrix
        #other.powered_by(self)
        Scalar.Raise ErrOperationNotImplemented, "**", @value.class, other.class
      else
        apply_through_coercion(other, __method__)
      end
    end
  end

end
set_component(i, j, v)
Alias for: []=
set_element(i, j, v)
Alias for: []=