LU decomposition with partial pivoting

Gaussian elimination, with its multipliers stored in L and the eliminated rows in U, factors P A = L U. The script reuses the factorization to solve multiple right-hand sides and compute the determinant, both far cheaper than refactoring from scratch each time.

import numpy as np

from mathematicskit.linalg import lu_decompose, lu_det, lu_solve

Factor once, solve for several right-hand sides

A = np.array([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]])
result = lu_decompose(A)
print("P @ A == L @ U:", np.allclose(result.P @ A, result.L @ result.U))

for b in (np.array([5.0, 12.0, 24.0]), np.array([1.0, 0.0, -1.0])):
    x = lu_solve(result, b)
    print(f"b={b} -> x={np.round(x, 6)}, residual={np.linalg.norm(A @ x - b):.2e}")

P @ A == L @ U: True
b=[ 5. 12. 24.] -> x=[ 1.5  3.  -1. ], residual=8.88e-16
b=[ 1.  0. -1.] -> x=[ 1.5 -2.5  0.5], residual=4.97e-16

Determinant via the pivoted triangular factors

det_lu = lu_det(A)
det_numpy = np.linalg.det(A)
print(f"det via LU: {det_lu:.10f}, numpy.linalg.det: {det_numpy:.10f}")

det via LU: 4.0000000000, numpy.linalg.det: 4.0000000000