pub fn compute_residue(
    a: &CsrMatrix<f64>,
    x: &DVector<f64>,
    b: &DVector<f64>,
    residue: &mut DVector<f64>
)