Function sor_solve 

pub fn sor_solve(
    a: &CsrMatrix,
    b: &[f64],
    omega: f64,
    tol: f64,
    max_iter: usize,
) -> (Vec<f64>, u32)

Successive Over-Relaxation (SOR) solver. Use omega = 1 for Gauss-Seidel.

Returns (x, n_iter).