Gaussian elimination akin to the verion on rosettacode.org

gaussian_elimination.py

# This is a modified version of the python implementation found on
# https://rosettacode.org/wiki/Gaussian_elimination#Python
#
# NOTE: I observe the rosettacode code has a bug for other dimensions of B than in the example (some index error)

import copy
from fractions import Fraction
        
 
def gauss_rosettacode(a, b):
    a = copy.deepcopy(a)
    b = copy.deepcopy(b)
    n = len(a)
    p = len(b[0])
    det = 1
    for i in range(n - 1):
        k = i
        for j in range(i + 1, n):
            if abs(a[j][i]) > abs(a[k][i]):
                k = j
        if k != i:
            try:
                a[i], a[k] = a[k], a[i]
            except IndexError:
                pass
            try:
                b[i], b[k] = b[k], b[i]
            except IndexError:
                pass
            det = -det
 
        for j in range(i + 1, n):
            t = a[j][i]/a[i][i]
            for k in range(i + 1, n):
                a[j][k] -= t*a[i][k]
            for k in range(p):
                b[j][k] -= t*b[i][k]
 
    for i in range(n - 1, -1, -1):
        for j in range(i + 1, n):
            t = a[i][j]
            for k in range(p):
                b[i][k] -= t*b[j][k]
        t = 1/a[i][i]
        det *= a[i][i]
        for j in range(p):
            b[i][j] *= t
            
    return b, det
    

def gaussian_elimination(a, b):
    pivots_a = []
    
    n = len(a)
    m = len(b)
    p = len(b[0])
    
    for i in range(n - 1):
        
        k = i
        for j in range(i + 1, n):
            if abs(a[j][i]) > abs(a[k][i]):
                k = j
        if k != i:
            pivots_a = [-1] + pivots_a
            a[i], a[k] = a[k], a[i]
            b[i], b[k] = b[k], b[i]
 
        for j in range(i + 1, n):
            t = a[j][i] / a[i][i]
            for k in range(i + 1, n):
                a[j][k] -= a[i][k] * t
                
        for j in range(i + 1, m):
            for k in range(m):
                b[j][k] -= b[i][k] * t
 
    for i in range(n - 1, -1, -1):
        aii = a[i][i]
        for j in range(i + 1, n):
            aij = a[i][j]
            b[i] = [bij - aij * b[j][k] for k, bij in enumerate(b[i])]
        b[i] = [bij / aii for bij in b[i]]

        pivots_a.append(aii)

    from functools import reduce
    det_a = reduce(lambda x, y: x * y, pivots_a) 
    
    return b, det_a


def matrix_map(fun, mat):
    def map_fun(row):
        return list(map(fun, row))
    return list(map(map_fun, mat))
    
    
def matrix_multiply(scalar, mat):
    def multiply_by_scalar(elem):
        return scalar * elem
    return matrix_map(multiply_by_scalar, mat)
 
 
def print_matrix(x):
    print()
    for row in x:
        print(row)


if __name__ == '__main__':
    """
    The 'gaussian_elimination' function takes two matrices A and B, with A a square matrix 
    and it return the determinant of A and a matrix x such that 
    A * x = B
    Note: Can be used to invert A by setting B to the identity.
    """
    
    A = [[2, 9, 4], [7, 5, 3], [6, 1, 8]]
    B = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    
    x, deta_a = gauss_rosettacode(A, B)
    print_matrix(x)
    print(deta_a)

    A = matrix_map(Fraction, A)
    B = matrix_map(Fraction, B)
    
    x, _deta_a = gauss_rosettacode(A, B)
    print_matrix(x)
 
    x, deta_a = gaussian_elimination(A, B)
    print_matrix(x)
    print(deta_a)
    
    print_matrix(matrix_multiply(deta_a, x))