Do not use nnz as cursor (just write it once)

src/atoms/affine/hstack.c

@@ -63,43 +63,44 @@ 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++)
  {
  expr *child = hnode->args[i];
  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++)
  {
  expr *child = hnode->args[i];
  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;

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));

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);

src/utils/CSR_sum.c

@@ -29,56 +29,57 @@ 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++)
  {
  int a_ptr = A->p[row];
  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)