ggggggg

import numpy as np

def pivotedOuterProductLU(A):
    """
    Perform an LU factorization with *column* pivoting, using an
    outer-product (rank-1) update approach.
    Returns L, U, and P, where P is the column permutation matrix,
    so that A*P = L*U.
    """
    A = A.astype(float)        # make a floating copy
    n = A.shape[0]
    U = A.copy()
    L = np.eye(n, dtype=float)
    P = np.eye(n, dtype=float)
    
    for k in range(n-1):
        
        # 1) Pivot search in row k
        pivot_col_relative = np.argmax(np.abs(U[k, k:]))  # index in [0..(n-k-1)]
        pivot_col = pivot_col_relative + k
        
        # 2) Swap columns if needed
        if pivot_col != k:
            # Swap columns k and pivot_col in U
            U[:, [k, pivot_col]] = U[:, [pivot_col, k]]
            
            # Swap columns k and pivot_col in P
            P[:, [k, pivot_col]] = P[:, [pivot_col, k]]
            
            # Swap columns in L (below row k only matters, but simpler to swap all):
            if k > 0:
                L[:, [k, pivot_col]] = L[:, [pivot_col, k]]
        
        # 3) Use U[k,k] as pivot
        pivot_val = U[k, k]
        
        # Guard against a zero (or extremely small) pivot if needed
        if abs(pivot_val) < 1e-15:
            raise ValueError("Encountered zero pivot; factorization fails.")
        
        # 4) Compute multipliers in L
        for i in range(k+1, n):
            L[i, k] = U[i, k] / pivot_val
        
        # 5) Rank-1 update of trailing submatrix
        for i in range(k+1, n):
            for j in range(k+1, n):
                U[i, j] -= L[i, k] * U[k, j]
    
    return L, U, P

# --------------- TEST on the given 4x4 matrix ---------------
A_test = np.array([
    [ 4,  3,  2,  1],
    [ 8,  8,  5,  2],
    [16, 12, 10,  5],
    [32, 24, 20, 11]
], dtype=float)

L, U, P = pivotedOuterProductLU(A_test)

print("Matrix A:")
print(A_test, "\n")

print("L:")
print(L, "\n")

print("U:")
print(U, "\n")

print("P (column permutation matrix):")
print(P, "\n")

# Check that A * P ≈ L * U
lhs = A_test @ P
rhs = L @ U

print("Check A*P vs L*U (they should be very close):")
print("A*P =\n", lhs, "\n")
print("L*U =\n", rhs, "\n")

# If we compare lhs and rhs, they should match (up to floating-point rounding).
diff = np.linalg.norm(lhs - rhs)
print(f"Difference norm: {diff:e}")