Skip to content

Instantly share code, notes, and snippets.

@jamiees2
Created June 2, 2015 23:51
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jamiees2/ba3490416212ebbbbde8 to your computer and use it in GitHub Desktop.
Save jamiees2/ba3490416212ebbbbde8 to your computer and use it in GitHub Desktop.
A fast LU decomposition algorithm, along with computing the determinant
def matrixMul(A, B):
TB = list(zip(*B))
return [[sum(ea*eb for ea,eb in zip(a,b)) for b in TB] for a in A]
def pivotize(m):
"""Creates the pivoting matrix for m."""
n = len(m)
ID = [[float(i == j) for i in range(n)] for j in range(n)]
r = 0
for j in range(n):
row = max(range(j, n), key=lambda i: abs(m[i][j]))
if j != row:
ID[j], ID[row] = ID[row], ID[j]
r += 1
return ID, r
def lu(A):
"""Decomposes a nxn matrix A by PA=LU and returns L, U and P."""
n = len(A)
L = [[0.0] * n for i in range(n)]
U = [[0.0] * n for i in range(n)]
P, r = pivotize(A)
A2 = matrixMul(P, A)
for j in range(n):
L[j][j] = 1.0
for i in range(j+1):
s1 = sum(U[k][j] * L[i][k] for k in range(i))
U[i][j] = A2[i][j] - s1
for i in range(j, n):
s2 = sum(U[k][j] * L[i][k] for k in range(j))
L[i][j] = (A2[i][j] - s2) / U[j][j]
return (L, U, P, r)
def trace(m):
n = len(m)
r = 1
for i in range(n):
if len(m[i]) <= i:
break
r *= m[i][i]
return r
def det(m):
l, u, p, r = lu(m)
return (-1)**r * trace(l) * trace(u)
mat = [[1, 3, 5], [2, 4, 7], [1, 1, 0]]
print(det(mat))
@jamiees2
Copy link
Author

jamiees2 commented Jun 2, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment