Gauss–Jordan Elimination

gj-elim.rb

require 'bigdecimal'

def gj_elim(els)
  h = 0
  k = 0
  m = els.size
  n = els.first.size

  # form row-echelon matrix
  while h < m && k < n
    unless els[h][k].zero?
      [*(h + 1)...m].each do |i|
        f = BigDecimal(els[i][k]) / els[h][k]
        els[i][k..] = (k...n).map { |j| j == k ? 0 : els[i][j] - els[h][j] * f }
      end
      h += 1
    end

    k += 1
  end

  # back-substitution
  [*0...m].reverse.each do |i|
    t = els[i][n - 1]
    [*(i + 1)...m].each do |j|
      t -= els[i][j] * els[j][n - 1]
      els[i][j] = 0
    end
    els[i][n - 1] = t / els[i][i]
    els[i][i] = 1
  end

  els.map { |rows| rows.map(&:to_f) }
end