diff --git a/src/sparse/cs_matrix_ops.rs b/src/sparse/cs_matrix_ops.rs
index 07119c19..34a68f9b 100644
--- a/src/sparse/cs_matrix_ops.rs
+++ b/src/sparse/cs_matrix_ops.rs
@@ -241,7 +241,11 @@ where
                 &mut res,
             );
 
-            for p in res.data.p[j]..nz {
+            // Keep the output sorted.
+            let range = res.data.p[j]..nz;
+            res.data.i[range.clone()].sort();
+
+            for p in range {
                 res.data.vals[p] = workspace[res.data.i[p]]
             }
         }
diff --git a/tests/sparse/cs_ops.rs b/tests/sparse/cs_ops.rs
index fa98fdb3..49dfc2bc 100644
--- a/tests/sparse/cs_ops.rs
+++ b/tests/sparse/cs_ops.rs
@@ -47,8 +47,8 @@ fn cs_mat_mul() {
 #[test]
 fn cs_mat_add() {
     let m1 = Matrix4x5::new(
-        4.0, 1.0, 4.0, 0.0, 9.0,
-        5.0, 6.0, 0.0, 8.0, 10.0,
+        4.0, 1.0, 4.0, 0.0, 0.0,
+        5.0, 6.0, 0.0, 8.0, 0.0,
         9.0, 10.0, 11.0, 12.0, 0.0,
         0.0, 0.0, 1.0, 0.0, 10.0
     );