import numpy as np def echelonize(a, eps=1e-8): """Reduce a matrix to row echelon form. Return the rank of matrix""" m, n = a.shape # number of rows, columns i = 0; j = 0 while i < m and j < n: # Find the maximal element in j-th column, # starting from i-th row k = i; maxVal = abs(a[i, j]) for l in range(i+1, m): if abs(a[l, j]) > maxVal: maxVal = abs(a[l, j]) k = l if maxVal <= eps: # Zero column for l in range(i, m): a[l, j] = 0. j += 1 continue if k != i: a[[i, k]] = a[[k, i]] # Swap rows i, k. To save a determinant, a[k] *= (-1.) # change a sign of one row. r = a[i, j] for l in range(i+1, m): a[l] -= a[i]*(a[l, j]/r) a[l, j] = 0. i += 1 j += 1 return i # Return the rank of matrix = number of nonzero rows def determinant(a): if a.shape[0] != a.shape[1]: raise ValueError("determinant of non-square matrix") b = a.copy() r = echelonize(b) # r is the rank of the matrix if r < b.shape[0]: # print("Singular matrix") return 0. # Determinant of singular matrix = 0 d = 1. for i in range(b.shape[0]): d *= b[i, i] return d def reversePass(a, r, eps = 1e-8): '''a is an echelonized matrix, r is its rank. Annihilate elements in columns above the first nonzero elements of rows, and normalize rows so that first nonzero elements = 1. Example: [2 3 0 1] [1 0 -7.5 0] [0 1 5 2] --> [0 1 5 0] [0 0 0 1] [0 0 0 1] ''' m, n = a.shape assert r <= m for i in range(r-1, -1, -1): j = 0 while j < n and abs(a[i, j]) <= eps: j += 1 if j >= n: continue assert abs(a[i, j]) > eps a[i] *= 1./a[i, j] for k in range(i): a[k] -= a[i]*a[k, j] def invMatrix(a, eps = 1e-8): """Return an inverse matrix of a matrix a""" pass # To do... def solveLinearSystem(a, b, eps = 1e-8): """Solve a linear system a x = b where a is non-singular square matrix and b is array of free terms. Return an array x""" pass # To do...