Function spgemm 

pub fn spgemm(
    a: &CsrMatrix<f32>,
    b: &CsrMatrix<f32>,
) -> Result<CsrMatrix<f32>, SparseError>

Sparse matrix-matrix multiply: C = A * B (both CSR).

Uses Gustavson’s algorithm: for each row of A, scatter-gather into a dense workspace, then compress into CSR.

Errors

Returns error if A.cols() != B.rows().