import numpy as np def echelonize(a, eps = 1e-8): m = a.shape[0] n = a.shape[1] i = 0 j = 0 while i < m and j < n: # print("Iteration:", i, j) # print(a) # Find maximal element in j-th column, starting from i-th row indMax = i vMax = abs(a[i, j]) for k in range(i+1, m): v = abs(a[k, j]) if v > vMax: vMax = v indMax = k if vMax <= eps: for k in range(i, m): a[k, j] = 0. j += 1 continue # print("max. element in column", i, "is in row", indMax) if indMax != i: # Swap rows i, indMax a[[i, indMax]] = a[[indMax, i]] a[indMax] = -a[indMax] # To save a sign of determinant # print("swapping rows", i, indMax) # print(a) # Annihilate j-th column starting from i+1-th row r = (-1./a[i, j]) for k in range(i+1, m): a[k] += a[i]*(a[k, j]*r) a[k, j] = 0. # print("annihilating column", j) # print(a) i += 1 j += 1 return i def reversePass(a, r, eps = 1e-8): '''a is the matrix in the echelonized form, r is the rank of matrix a. Reduce the matrix a to diagonal form with units at the beginnings of rows and seroes in the columns above the first nonzero elements of rows''' m = a.shape[0] # Number of rows n = a.shape[1] # Number of columns for i in range(r-1, -1, -1): # Find the first nonzero element in i-th row j0 = n for j in range(n): if abs(a[i, j]) > eps: j0 = j break if j0 >= n: assert "Incorrect echelon form of matrix" # multiplyMatrixRow(a, i, 1./a[i][j0]) a[i] *= (1./a[i, j0]) for k in range(i): # addMatrixRows(a, k, i, -a[k][j0]) a[k] -= a[i]*a[k, j0] def inverseMatrix(a, eps = 1e-8): '''Return the inverse matrix to the matrix a of real numbers. The matrix a is represented as 2-dimensional numpy-array''' # To do... pass def solveLinearSystem(a, b, eps = 1e-8): '''Solve a linear system a*x = b, where the matrix a is represented by a 2-dimensional numpy-array of real numbers, and b is a linear numpy-array or any list/tuple''' # To do... pass def determinant(a, eps = 1e-8): m = a.shape[0] n = a.shape[1] if m != n: raise ValueError("Determinant of non-square matrix") b = a.copy() r = echelonize(b, eps) if r < m: return 0. d = 1. for i in range(m): d *= b[i, i] return d