pub fn lu_solve(l: &Matrix, u: &Matrix, piv: &[usize], b: &[f64]) -> Vec<f64>
Solve LUx = Pb using forward/back substitution.