Function compute_residual 

pub fn compute_residual(
    matrix: &dyn Matrix,
    x: &[Precision],
    b: &[Precision],
    residual: &mut [Precision],
) -> Result<()>

Compute residual vector: r = A*x - b