LUDecomposition.py

"""
newA
[[1 2 1]
 [1 1 1]
 [2 4 2]]

L, U = LU([(3-1,2-1,0),(3-1,2-1,3)], newA)
L= [[ 1.  0.  0.]
 [-0.  1. -0.]
 [ 0.  3.  1.]]
U= [[ 1.  2.  1.]
 [ 1.  1.  1.]
 [-1.  1. -1.]]
LU= [[1. 2. 1.]
 [1. 1. 1.]
 [2. 4. 2.]]
A= [[1 2 1]
 [1 1 1]
 [2 4 2]]
"""

import numpy as np
from numpy.linalg import inv 

def LU(steps, A):
    """
    steps should be a list of tuples, which is in the format like (row2, row1, mul).
    (row2, row1, mul) means that row2 - mul * row1,
    where row2 and row1 are the row id of matrix A, starting from 0.
    """
    assert len(A.shape) == 2 and A.shape[0] == A.shape[1]  # assume A is a square matrix
    
    # calculate intermediate matrices for debugging
    mats = []
    for s in steps:
        assert len(s) == 3  # ensure the correct format
        mat = np.identity(A.shape[0])
        mat[s[0], s[1]] = -s[2]
        mats.append(mat)
    
    # calculate L and U
    invL = np.identity(A.shape[0])
    for m in mats:
        invL = invL.dot(m)
    U = invL.dot(A)
    L = inv(invL)
    
    print("L=", L)
    print("U=", U)
    print("LU=", L.dot(U))
    print("A=", A)
    
    return L, U