Gaussian (Gauss-Jordan) Elimination in Ruby

I had to do some Gaussian elimination for an assignment. So here's Gaussian elimination in Ruby:

```# Performs an in-place Gaussian elimination on an NxN matrix 'matrix' (2D array
# of Numeric objects) and an N-element vector 'vector.' (array of N Numerics).
def gaussianElimination(matrix, vector)
0.upto(matrix.length - 2) do |pivotIdx|
# Find the best pivot. This is the one who has the largest absolute value
# relative to his row (scaled partial pivoting). This step can be omitted
# to improve speed at the cost of increased error.
maxRelVal = 0
maxIdx = pivotIdx
(pivotIdx).upto(matrix.length - 1) do |row|
relVal = matrix[row][pivotIdx] / matrix[row].map{ |x| x.abs }.max
if relVal >= maxRelVal
maxRelVal = relVal
maxIdx = row
end
end

# Swap the best pivot row into place.
matrix[pivotIdx], matrix[maxIdx] = matrix[maxIdx], matrix[pivotIdx]
vector[pivotIdx], vector[maxIdx] = vector[maxIdx], vector[pivotIdx]

pivot = matrix[pivotIdx][pivotIdx]
# Loop over each row below the pivot row.
(pivotIdx+1).upto(matrix.length - 1) do |row|
# Find factor so that [this row] = [this row] - factor*[pivot row]
# leaves 0 in the pivot column.
factor = matrix[row][pivotIdx]/pivot
# We know it will be zero.
matrix[row][pivotIdx] = 0.0
# Compute [this row] = [this row] - factor*[pivot row] for the other cols.
(pivotIdx+1).upto(matrix[row].length - 1) do |col|
matrix[row][col] -= factor*matrix[pivotIdx][col]
end
vector[row] -= factor*vector[pivotIdx]
end
end

return [matrix,vector]
end

# Assumes 'matrix' is in row echelon form.
def backSubstitution(matrix, vector)
(matrix.length - 1).downto( 0 ) do |row|
tail = vector[row]
(row+1).upto(matrix.length - 1) do |col|
tail -= matrix[row][col] * vector[col]
matrix[row][col] = 0.0
end
vector[row] = tail / matrix[row][row]
matrix[row][row] = 1.0
end
end

# Example usage:
require 'pp'

# A system of equations: matrix * X = vector
matrix =
[
[1.0, 1.0, 1.0, 1.0],
[0.0, 1.0, 2.0, 3.0],
[1.0, 2.0, 4.0, 8.0],
[0.0, 1.0, 4.0, 12.0],
]
vector = [1.0, 0.0, 2.0, 0.0]

# Create a backup for verification.
vector_backup= vector.dup

# Gaussian elemination to put the system in row echelon form.
gaussianElimination(matrix, vector)
# Back-substitution to solve the system.
backSubstitution(matrix, vector)

# Print the result.
pp matrix
pp vector

# Verify the result.
pass = true

0.upto(matrix_backup.length - 1) do  |eqn|
sum = 0
0.upto(matrix_backup[eqn].length - 1) do |term|
sum += matrix_backup[eqn][term] * vector[term]
end
if (sum - vector_backup[eqn]).abs > 0.0000000001
pass = false
break
end
end

if pass
puts "Verification PASSED."
else
puts "Verification FAILED."
end
```