Function psor 

pub fn psor(
    a: &[f64],
    b: &[f64],
    c: &[f64],
    d: &[f64],
    floor: Option<&[f64]>,
    cap: Option<&[f64]>,
    omega: f64,
    tol: f64,
    max_iter: usize,
) -> PsorResult

Solve A x = d projected onto floor <= x <= cap by projected SOR, where A is tridiagonal in the same layout as thomas_algorithm: a is the sub-diagonal (a[i-1] multiplies x[i-1] in row i), b the diagonal and c the super-diagonal.

omega in (0, 2) is the relaxation factor (1 = projected Gauss-Seidel; ~1.2-1.6 typically accelerates diffusion operators). Convergence requires the usual SOR conditions (diagonally dominant or symmetric positive definite A), which theta-scheme matrices satisfy. The iteration stops when the largest update falls below tol.