lu_decomposition.py

# lower–upper (LU) decomposition - https://en.wikipedia.org/wiki/LU_decomposition
import numpy

def LUDecompose (table):
    # Table that contains our data
    # Table has to be a square array so we need to check first
    rows,columns=numpy.shape(table)
    L=numpy.zeros((rows,columns))
    U=numpy.zeros((rows,columns))
    if rows!=columns:
        return []
    for i in range (columns):
        for j in range(i-1):
            sum=0
            for k in range (j-1):
                sum+=L[i][k]*U[k][j]
            L[i][j]=(table[i][j]-sum)/U[j][j]
        L[i][i]=1
        for j in range(i-1,columns):
            sum1=0
            for k in range(i-1):
                sum1+=L[i][k]*U[k][j]
            U[i][j]=table[i][j]-sum1
    return L,U

if __name__ == "__main__":
    matrix =numpy.array([[2,-2,1],
                         [0,1,2],
                         [5,3,1]])
    L,U = LUDecompose(matrix)
    print(L)
    print(U)