AAM: Algorithms for Applied Mathematics

def palu(self, debug=False):
    m,n = self.dims()
    if m!=n: raise ValueError("Cannot PA=LU decompose nonsquare matrix.")
    # Initialize grids
    grid = [[entry for entry in row] for row in self._data] # Skim the data
    pmat = [[1 if r==c else 0 for c in range(n)] for r in range(m)]
    debug_out = []
    # Start at upper left
    r,c = 0,0
    while (r < m) and (c < n):
        # Check if rows need swapped
        debug_out.append( str(grid) )
        cur_max = abs(grid[r][c])
        piv = r
        for r_prime in range(r+1, m):
            if abs(grid[r_prime][c]) > cur_max:
                cur_max = abs(grid[r_prime][c])
                piv = r_prime
        # if cur_max == 0, shift right.
        if cur_max == 0:
            c += 1
            continue
            # otherwise check if rows need swapped
        elif piv != r:
            debug_out.append(f"Swap R({r})<->R({r_prime})")
            # interchange elements of lists
            grid[piv],grid[r] = grid[r],grid[piv]
            pmat[piv],pmat[r] = pmat[r],pmat[piv]
        # No scaling
        # Elimination
        for r_prime in range(r+1, m):
            if grid[r_prime][c] == 0: continue
            alpha = grid[r_prime][c] / grid[r][c]
            # Since we want to maintain the lower-triangular part for L, we only alter
            # columns c+1 through n of grid[r_prime]
            grid[r_prime][c+1:] = [-alpha*grid[r][c_prime] + grid[r_prime][c_prime] for c_prime in range(c+1,n)]
            # Then store the alpha value in grid[r_prime][c] to pull out for L later
            grid[r_prime][c] = alpha
            debug_out.append(f"Elim ({-alpha})*R({r}) + R({r_prime})")
        r += 1
        c += 1
    P = AlgoMatrix(pmat)
    L = AlgoMatrix([[1 if r==c else (grid[r][c] if r > c else 0) for c in range(n)] for r in range(m)])
    U = AlgoMatrix([[grid[r][c] if r <= c else 0 for c in range(n)] for r in range(m)])
    if debug: print("\n".join(debug_out))
    return P,L,U