Skip to content

Instantly share code, notes, and snippets.

@bryan-lunt

bryan-lunt/lu_fact.py

Created September 17, 2021 19:36
Show Gist options
  • Save bryan-lunt/428033e51ea0491bc39018142303b3de to your computer and use it in GitHub Desktop.
Save bryan-lunt/428033e51ea0491bc39018142303b3de to your computer and use it in GitHub Desktop.
Download ZIP
Simple LU factorization examples in Python
Raw
lu_fact.py
import numpy as np
def LU_recursive(A):
"""LU factorization without pivoting. Will fail if any diagonal is 0.
Recursive implementation.
"""
L = np.zeros_like(A,dtype=np.float_)
U = np.eye(*A.shape,dtype=np.float_)
a = A[0,0]
L[0,0] = 1.0
U[0,0] = a
if A.shape[0] > 1:
L[1:,0] = A[1:,0]
U[0,1:] = A[0,1:]
L[1:,0] *= np.power(a,-1.0)
L[1:,1:] = A[1:,1:] - np.outer(L[1:,0],A[0,1:])
L[1:,1:], U[1:,1:] = LU_recursive(L[1:,1:])
return L,U
def LU_loop(A):
"""LU factorization without pivoting. Will fail if any diagonal is 0.
Loop implementation.
"""
L = np.copy(A)
U = np.eye(*A.shape,dtype=np.float_)
i = 0
M = A.shape[0]
while( True ):
a = L[i,i]
L[i,i] = 1.0
U[i,i] = a
if(i == M-1):
break
else:
U[i,i+1:] = L[i,i+1:]
L[i,i+1:] = 0.0
L[i+1:,i] *= np.power(a,-1.0)
U[i+1:,i] = 0.0
L[i+1:,i+1:] -= np.outer(L[i+1:,i],U[i,i+1:])
i=i+1
return L,U
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment