diff --git a/src/atoms/affine/hstack.c b/src/atoms/affine/hstack.c
index a0e0a591..c5c5dc18 100644
--- a/src/atoms/affine/hstack.c
+++ b/src/atoms/affine/hstack.c
@@ -63,7 +63,7 @@ static void jacobian_init_impl(expr *node)
 
     /* precompute sparsity pattern of this node's jacobian */
     int row_offset = 0;
-    A->nnz = 0;
+    int cursor = 0;
 
     for (int i = 0; i < hnode->n_args; i++)
     {
@@ -71,25 +71,26 @@ static void jacobian_init_impl(expr *node)
         CSR_matrix *B = child->jacobian->to_csr(child->jacobian);
 
         /* copy columns */
-        memcpy(A->i + A->nnz, B->i, B->nnz * sizeof(int));
+        memcpy(A->i + cursor, B->i, B->nnz * sizeof(int));
 
         /* set row pointers */
         for (int r = 0; r < child->size; r++)
         {
-            A->p[row_offset + r] = A->nnz + B->p[r];
+            A->p[row_offset + r] = cursor + B->p[r];
         }
 
-        A->nnz += B->nnz;
+        cursor += B->nnz;
         row_offset += child->size;
     }
-    A->p[node->size] = A->nnz;
+    A->p[node->size] = cursor;
+    assert(cursor == A->nnz);
     node->jacobian = new_sparse_matrix(A);
 }
 
 static void eval_jacobian(expr *node)
 {
     hstack_expr *hnode = (hstack_expr *) node;
-    node->jacobian->nnz = 0;
+    int cursor = 0;
 
     for (int i = 0; i < hnode->n_args; i++)
     {
@@ -97,9 +98,9 @@ static void eval_jacobian(expr *node)
         child->eval_jacobian(child);
         /* to_csr needed for stacked_pd */
         CSR_matrix *child_csr = child->jacobian->to_csr(child->jacobian);
-        memcpy(node->jacobian->x + node->jacobian->nnz, child_csr->x,
+        memcpy(node->jacobian->x + cursor, child_csr->x,
                child_csr->nnz * sizeof(double));
-        node->jacobian->nnz += child_csr->nnz;
+        cursor += child_csr->nnz;
     }
 }
 
@@ -121,7 +122,6 @@ static void wsum_hess_init_impl(expr *node)
     hnode->CSR_work = new_CSR_matrix(node->n_vars, node->n_vars, nnz_ub);
 
     /* fill sparsity pattern */
-    H->nnz = 0;
     for (int i = 0; i < hnode->n_args; i++)
     {
         matrix *child_hess = hnode->args[i]->wsum_hess;
diff --git a/src/atoms/bivariate_restricted_dom/quad_over_lin.c b/src/atoms/bivariate_restricted_dom/quad_over_lin.c
index e6c5b897..b4460f91 100644
--- a/src/atoms/bivariate_restricted_dom/quad_over_lin.c
+++ b/src/atoms/bivariate_restricted_dom/quad_over_lin.c
@@ -94,24 +94,25 @@ static void jacobian_init_impl(expr *node)
         CSR_matrix *jac = new_CSR_matrix(1, node->n_vars, nonzero_cols + 1);
 
         /* precompute column indices */
-        jac->nnz = 0;
+        int cursor = 0;
         for (int j = 0; j < node->n_vars; j++)
         {
             if (col_nz[j])
             {
-                jac->i[jac->nnz] = j;
-                jac->nnz++;
+                jac->i[cursor] = j;
+                cursor++;
             }
         }
-        assert(nonzero_cols == jac->nnz);
+        assert(nonzero_cols == cursor);
 
         sp_free(col_nz);
 
         /* insert y variable index at correct position */
-        insert_idx(y->var_id, jac->i, jac->nnz);
-        jac->nnz += 1;
+        insert_idx(y->var_id, jac->i, cursor);
+        cursor += 1;
+        assert(cursor == jac->nnz);
         jac->p[0] = 0;
-        jac->p[1] = jac->nnz;
+        jac->p[1] = cursor;
 
         /* find position where y should be inserted */
         node->work->iwork = (int *) sp_malloc(sizeof(int));
diff --git a/src/problem.c b/src/problem.c
index 09cc14d5..f0bf7997 100644
--- a/src/problem.c
+++ b/src/problem.c
@@ -531,9 +531,7 @@ void problem_jacobian(problem *prob)
         nnz_offset += c->jacobian->nnz;
     }
 
-    /* update actual nnz (may be less than allocated) */
-    J->nnz = nnz_offset;
-
+    assert(nnz_offset == J->nnz);
     prob->jacobian_called = true;
     clock_gettime(CLOCK_MONOTONIC, &timer.end);
     prob->stats.time_eval_jacobian += GET_ELAPSED_SECONDS(timer);
diff --git a/src/utils/CSR_sum.c b/src/utils/CSR_sum.c
index 19f99a3d..1392b353 100644
--- a/src/utils/CSR_sum.c
+++ b/src/utils/CSR_sum.c
@@ -29,7 +29,7 @@ void sum_csr_alloc(const CSR_matrix *A, const CSR_matrix *B, CSR_matrix *C)
     /* A and B must be different from C */
     assert(A != C && B != C);
 
-    C->nnz = 0;
+    int cursor = 0;
 
     for (int row = 0; row < A->m; row++)
     {
@@ -37,48 +37,49 @@ void sum_csr_alloc(const CSR_matrix *A, const CSR_matrix *B, CSR_matrix *C)
         int a_end = A->p[row + 1];
         int b_ptr = B->p[row];
         int b_end = B->p[row + 1];
-        C->p[row] = C->nnz;
+        C->p[row] = cursor;
 
         /* Merge while both have elements (only column indices) */
         while (a_ptr < a_end && b_ptr < b_end)
         {
             if (A->i[a_ptr] < B->i[b_ptr])
             {
-                C->i[C->nnz] = A->i[a_ptr];
+                C->i[cursor] = A->i[a_ptr];
                 a_ptr++;
             }
             else if (B->i[b_ptr] < A->i[a_ptr])
             {
-                C->i[C->nnz] = B->i[b_ptr];
+                C->i[cursor] = B->i[b_ptr];
                 b_ptr++;
             }
             else
             {
-                C->i[C->nnz] = A->i[a_ptr];
+                C->i[cursor] = A->i[a_ptr];
                 a_ptr++;
                 b_ptr++;
             }
-            C->nnz++;
+            cursor++;
         }
 
         /* Copy remaining elements from A */
         if (a_ptr < a_end)
         {
             int a_remaining = a_end - a_ptr;
-            memcpy(C->i + C->nnz, A->i + a_ptr, a_remaining * sizeof(int));
-            C->nnz += a_remaining;
+            memcpy(C->i + cursor, A->i + a_ptr, a_remaining * sizeof(int));
+            cursor += a_remaining;
         }
 
         /* Copy remaining elements from B */
         if (b_ptr < b_end)
         {
             int b_remaining = b_end - b_ptr;
-            memcpy(C->i + C->nnz, B->i + b_ptr, b_remaining * sizeof(int));
-            C->nnz += b_remaining;
+            memcpy(C->i + cursor, B->i + b_ptr, b_remaining * sizeof(int));
+            cursor += b_remaining;
         }
     }
 
-    C->p[A->m] = C->nnz;
+    C->p[A->m] = cursor;
+    C->nnz = cursor;
 }
 
 void sum_csr_fill_values(const CSR_matrix *A, const CSR_matrix *B, CSR_matrix *C)
@@ -157,7 +158,7 @@ void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
 
     C->n = A->n;
     C->p[0] = 0;
-    C->nnz = 0;
+    int cursor = 0;
 
     int *cols = iwork;
     int *col_to_pos = iwork;
@@ -171,7 +172,7 @@ void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
         // Build sparsity pattern of the row resulting from summing
         // the block of rows from A
         // -----------------------------------------------------------------
-        C->p[block] = C->nnz;
+        C->p[block] = cursor;
         int count = 0;
         for (int row = start_row; row < end_row; row++)
         {
@@ -191,14 +192,14 @@ void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
             int col = cols[t];
             if (t == 0 || col != prev_col)
             {
-                C->i[C->nnz + unique_nnz] = col;
+                C->i[cursor + unique_nnz] = col;
                 prev_col = col;
                 unique_nnz++;
             }
         }
 
-        C->nnz += unique_nnz;
-        C->p[block + 1] = C->nnz;
+        cursor += unique_nnz;
+        C->p[block + 1] = cursor;
 
         // -----------------------------------------------------------------
         //         Build idx_map for all entries in this block
@@ -217,6 +218,8 @@ void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
             }
         }
     }
+
+    C->nnz = cursor;
 }
 
 /* iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz */
@@ -226,7 +229,7 @@ void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
     assert(C->m == row_spacing);
     C->n = A->n;
     C->p[0] = 0;
-    C->nnz = 0;
+    int cursor = 0;
 
     int *cols = iwork;
     int *col_to_pos = iwork;
@@ -237,7 +240,7 @@ void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
         // Build sparsity pattern of the row resulting from summing
         // evenly spaced rows from A
         // -----------------------------------------------------------------
-        C->p[C_row] = C->nnz;
+        C->p[C_row] = cursor;
         int count = 0;
         for (int row = C_row; row < A->m; row += row_spacing)
         {
@@ -257,14 +260,14 @@ void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
             int col = cols[t];
             if (t == 0 || col != prev_col)
             {
-                C->i[C->nnz + unique_nnz] = col;
+                C->i[cursor + unique_nnz] = col;
                 prev_col = col;
                 unique_nnz++;
             }
         }
 
-        C->nnz += unique_nnz;
-        C->p[C_row + 1] = C->nnz;
+        cursor += unique_nnz;
+        C->p[C_row + 1] = cursor;
 
         // -----------------------------------------------------------------
         //         Build idx_map for all entries in evenly spaced rows
@@ -283,6 +286,8 @@ void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
             }
         }
     }
+
+    C->nnz = cursor;
 }
 
 void accumulator(const double *vals, int nnz, const int *idx_map, double *out)