COO support (#47)

* basic coo support

* run formatter only on macos

* ran formatter

* added comment

* formatter

* updated version number

Author: dance858

.github/workflows/formatting.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-latest]
+        os: [macos-latest]
 
     steps:
       - uses: actions/checkout@v5

CMakeLists.txt

@@ -4,7 +4,7 @@ set(CMAKE_C_STANDARD 99)
 
 set(DIFF_ENGINE_VERSION_MAJOR 0)
 set(DIFF_ENGINE_VERSION_MINOR 1)
-set(DIFF_ENGINE_VERSION_PATCH 2)
+set(DIFF_ENGINE_VERSION_PATCH 3)
 set(DIFF_ENGINE_VERSION "${DIFF_ENGINE_VERSION_MAJOR}.${DIFF_ENGINE_VERSION_MINOR}.${DIFF_ENGINE_VERSION_PATCH}")
 add_compile_definitions(DIFF_ENGINE_VERSION="${DIFF_ENGINE_VERSION}")
 

include/problem.h

@@ -19,6 +19,7 @@
 #define PROBLEM_H
 
 #include "expr.h"
+#include "utils/COO_Matrix.h"
 #include "utils/CSR_Matrix.h"
 #include "utils/Timer.h"
 #include <stdbool.h>
@@ -45,21 +46,22 @@ typedef struct problem
     int n_vars;
     int total_constraint_size;
 
-    /* Allocated by new_problem */
+    /* allocated by new_problem */
     double *constraint_values;
     double *gradient_values;
 
-    /* Allocated by problem_init_derivatives */
+    /* allocated by problem_init_derivatives */
     CSR_Matrix *jacobian;
     CSR_Matrix *lagrange_hessian;
-    int *hess_idx_map; /* Maps all wsum_hess nnz to lagrange_hessian (obj +
-                          constraints) */
+    int *hess_idx_map; /* maps all wsum_hess nnz to lagrange_hessian */
+    COO_Matrix *jacobian_coo;
+    COO_Matrix *lagrange_hessian_coo; /* lower triangular part stored in COO */
 
     /* for the affine shortcut we keep track of the first time the jacobian and
      * hessian are called */
     bool jacobian_called;
 
-    /* Statistics for performance measurement */
+    /* statistics for performance measurement */
     Diff_engine_stats stats;
     bool verbose;
 } problem;
@@ -70,6 +72,8 @@ problem *new_problem(expr *objective, expr **constraints, int n_constraints,
 void problem_init_jacobian(problem *prob);
 void problem_init_hessian(problem *prob);
 void problem_init_derivatives(problem *prob);
+void problem_init_jacobian_coo(problem *prob);
+void problem_init_hessian_coo_lower_triangular(problem *prob);
 void free_problem(problem *prob);
 
 double problem_objective_forward(problem *prob, const double *u);

include/utils/COO_Matrix.h

@@ -0,0 +1,43 @@
+#ifndef COO_MATRIX_H
+#define COO_MATRIX_H
+
+#include "CSR_Matrix.h"
+
+/* COO (Coordinate) Sparse Matrix Format
+ *
+ * For an m x n matrix with nnz nonzeros:
+ * - rows: array of size nnz containing row indices
+ * - cols: array of size nnz containing column indices
+ * - x: array of size nnz containing values
+ * - value_map: array of size nnz mapping CSR entries to COO entries (for
+ * lower-triangular COO)
+ * - m: number of rows
+ * - n: number of columns
+ * - nnz: number of nonzero entries
+ */
+typedef struct COO_Matrix
+{
+    int *rows;
+    int *cols;
+    double *x;
+    int *value_map;
+    int m;
+    int n;
+    int nnz;
+} COO_Matrix;
+
+/* Construct a COO matrix from a CSR matrix */
+COO_Matrix *new_coo_matrix(const CSR_Matrix *A);
+
+/* Construct a COO matrix containing only the lower-triangular
+ * entries (col <= row) of a symmetric CSR matrix. Populates
+ * value_map so that refresh_lower_triangular_coo can update
+ * values without recomputing structure. */
+COO_Matrix *new_coo_matrix_lower_triangular(const CSR_Matrix *A);
+
+/* Refresh COO values from a new CSR value array using value_map */
+void refresh_lower_triangular_coo(COO_Matrix *coo, const double *vals);
+
+void free_coo_matrix(COO_Matrix *matrix);
+
+#endif /* COO_MATRIX_H */

src/problem.c

@@ -237,6 +237,27 @@ void problem_init_hessian(problem *prob)
     prob->stats.time_init_derivatives += GET_ELAPSED_SECONDS(timer);
 }
 
+void problem_init_jacobian_coo(problem *prob)
+{
+    problem_init_jacobian(prob);
+    Timer timer;
+    clock_gettime(CLOCK_MONOTONIC, &timer.start);
+    prob->jacobian_coo = new_coo_matrix(prob->jacobian);
+    clock_gettime(CLOCK_MONOTONIC, &timer.end);
+    prob->stats.time_init_derivatives += GET_ELAPSED_SECONDS(timer);
+}
+
+void problem_init_hessian_coo_lower_triangular(problem *prob)
+{
+    problem_init_hessian(prob);
+    Timer timer;
+    clock_gettime(CLOCK_MONOTONIC, &timer.start);
+    prob->lagrange_hessian_coo =
+        new_coo_matrix_lower_triangular(prob->lagrange_hessian);
+    clock_gettime(CLOCK_MONOTONIC, &timer.end);
+    prob->stats.time_init_derivatives += GET_ELAPSED_SECONDS(timer);
+}
+
 void problem_init_derivatives(problem *prob)
 {
     problem_init_jacobian(prob);
@@ -286,6 +307,7 @@ void free_problem(problem *prob)
     free(prob->gradient_values);
     free_csr_matrix(prob->jacobian);
     free_csr_matrix(prob->lagrange_hessian);
+    free_coo_matrix(prob->jacobian_coo);
     free(prob->hess_idx_map);
 
     /* Release expression references (decrements refcount) */

src/utils/COO_Matrix.c

@@ -0,0 +1,109 @@
+/*
+ * Copyright 2026 Daniel Cederberg and William Zhang
+ *
+ * This file is part of the DNLP-differentiation-engine project.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "utils/COO_Matrix.h"
+#include <stdlib.h>
+#include <string.h>
+
+COO_Matrix *new_coo_matrix(const CSR_Matrix *A)
+{
+    COO_Matrix *coo = (COO_Matrix *) malloc(sizeof(COO_Matrix));
+    coo->m = A->m;
+    coo->n = A->n;
+    coo->nnz = A->nnz;
+    coo->rows = (int *) malloc(A->nnz * sizeof(int));
+    coo->cols = (int *) malloc(A->nnz * sizeof(int));
+    coo->x = (double *) malloc(A->nnz * sizeof(double));
+    coo->value_map = NULL;
+
+    for (int r = 0; r < A->m; r++)
+    {
+        for (int k = A->p[r]; k < A->p[r + 1]; k++)
+        {
+            coo->rows[k] = r;
+        }
+    }
+
+    memcpy(coo->cols, A->i, A->nnz * sizeof(int));
+    memcpy(coo->x, A->x, A->nnz * sizeof(double));
+
+    return coo;
+}
+
+COO_Matrix *new_coo_matrix_lower_triangular(const CSR_Matrix *A)
+{
+    /* Pass 1: count lower-triangular entries (col <= row) */
+    int count = 0;
+    for (int r = 0; r < A->m; r++)
+    {
+        for (int k = A->p[r]; k < A->p[r + 1]; k++)
+        {
+            if (A->i[k] <= r)
+            {
+                count++;
+            }
+        }
+    }
+
+    COO_Matrix *coo = (COO_Matrix *) malloc(sizeof(COO_Matrix));
+    coo->m = A->m;
+    coo->n = A->n;
+    coo->nnz = count;
+    coo->rows = (int *) malloc(count * sizeof(int));
+    coo->cols = (int *) malloc(count * sizeof(int));
+    coo->x = (double *) malloc(count * sizeof(double));
+    coo->value_map = (int *) malloc(count * sizeof(int));
+
+    /* Pass 2: fill arrays */
+    int idx = 0;
+    for (int r = 0; r < A->m; r++)
+    {
+        for (int k = A->p[r]; k < A->p[r + 1]; k++)
+        {
+            if (A->i[k] <= r)
+            {
+                coo->rows[idx] = r;
+                coo->cols[idx] = A->i[k];
+                coo->x[idx] = A->x[k];
+                coo->value_map[idx] = k;
+                idx++;
+            }
+        }
+    }
+
+    return coo;
+}
+
+void refresh_lower_triangular_coo(COO_Matrix *coo, const double *vals)
+{
+    for (int i = 0; i < coo->nnz; i++)
+    {
+        coo->x[i] = vals[coo->value_map[i]];
+    }
+}
+
+void free_coo_matrix(COO_Matrix *matrix)
+{
+    if (matrix)
+    {
+        free(matrix->rows);
+        free(matrix->cols);
+        free(matrix->x);
+        free(matrix->value_map);
+        free(matrix);
+    }
+}

tests/all_tests.c

@@ -43,6 +43,7 @@
 #include "jacobian_tests/test_trace.h"
 #include "jacobian_tests/test_transpose.h"
 #include "problem/test_problem.h"
+#include "utils/test_coo_matrix.h"
 #include "utils/test_csc_matrix.h"
 #include "utils/test_csr_csc_conversion.h"
 #include "utils/test_csr_matrix.h"
@@ -271,6 +272,9 @@ int main(void)
     mu_run_test(test_ATA_alloc_random, tests_run);
     mu_run_test(test_ATA_alloc_random2, tests_run);
     mu_run_test(test_BTA_alloc_and_BTDA_fill, tests_run);
+    mu_run_test(test_csr_to_coo, tests_run);
+    mu_run_test(test_csr_to_coo_lower_triangular, tests_run);
+    mu_run_test(test_refresh_lower_triangular_coo, tests_run);
 
     printf("\n--- Problem Struct Tests ---\n");
     mu_run_test(test_problem_new_free, tests_run);

tests/utils/test_coo_matrix.h

@@ -0,0 +1,104 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "minunit.h"
+#include "test_helpers.h"
+#include "utils/COO_Matrix.h"
+
+const char *test_csr_to_coo()
+{
+    /* Create a 3x3 CSR matrix A:
+     * [1.0  2.0  0.0]
+     * [0.0  3.0  4.0]
+     * [5.0  0.0  6.0]
+     */
+    CSR_Matrix *A = new_csr_matrix(3, 3, 6);
+    double Ax[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+    int Ai[6] = {0, 1, 1, 2, 0, 2};
+    int Ap[4] = {0, 2, 4, 6};
+    memcpy(A->x, Ax, 6 * sizeof(double));
+    memcpy(A->i, Ai, 6 * sizeof(int));
+    memcpy(A->p, Ap, 4 * sizeof(int));
+
+    COO_Matrix *coo = new_coo_matrix(A);
+
+    mu_assert("m incorrect", coo->m == 3);
+    mu_assert("n incorrect", coo->n == 3);
+    mu_assert("nnz incorrect", coo->nnz == 6);
+
+    int expected_rows[6] = {0, 0, 1, 1, 2, 2};
+    int expected_cols[6] = {0, 1, 1, 2, 0, 2};
+    double expected_x[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+
+    mu_assert("rows incorrect", cmp_int_array(coo->rows, expected_rows, 6));
+    mu_assert("cols incorrect", cmp_int_array(coo->cols, expected_cols, 6));
+    mu_assert("vals incorrect", cmp_double_array(coo->x, expected_x, 6));
+
+    free_coo_matrix(coo);
+    free_csr_matrix(A);
+
+    return 0;
+}
+
+const char *test_csr_to_coo_lower_triangular()
+{
+    /* Symmetric 3x3 matrix:
+     * [1  2  3]
+     * [2  5  6]
+     * [3  6  9]
+     */
+    CSR_Matrix *A = new_csr_matrix(3, 3, 9);
+    int Ap[4] = {0, 3, 6, 9};
+    int Ai[9] = {0, 1, 2, 0, 1, 2, 0, 1, 2};
+    double Ax[9] = {1, 2, 3, 2, 5, 6, 3, 6, 9};
+    memcpy(A->p, Ap, 4 * sizeof(int));
+    memcpy(A->i, Ai, 9 * sizeof(int));
+    memcpy(A->x, Ax, 9 * sizeof(double));
+
+    COO_Matrix *coo = new_coo_matrix_lower_triangular(A);
+
+    mu_assert("ltri m incorrect", coo->m == 3);
+    mu_assert("ltri n incorrect", coo->n == 3);
+    mu_assert("ltri nnz incorrect", coo->nnz == 6);
+
+    int expected_rows[6] = {0, 1, 1, 2, 2, 2};
+    int expected_cols[6] = {0, 0, 1, 0, 1, 2};
+    double expected_x[6] = {1, 2, 5, 3, 6, 9};
+    int expected_map[6] = {0, 3, 4, 6, 7, 8};
+
+    mu_assert("ltri rows incorrect", cmp_int_array(coo->rows, expected_rows, 6));
+    mu_assert("ltri cols incorrect", cmp_int_array(coo->cols, expected_cols, 6));
+    mu_assert("ltri vals incorrect", cmp_double_array(coo->x, expected_x, 6));
+    mu_assert("ltri value_map incorrect",
+              cmp_int_array(coo->value_map, expected_map, 6));
+
+    free_coo_matrix(coo);
+    free_csr_matrix(A);
+
+    return 0;
+}
+
+const char *test_refresh_lower_triangular_coo()
+{
+    CSR_Matrix *A = new_csr_matrix(3, 3, 9);
+    int Ap[4] = {0, 3, 6, 9};
+    int Ai[9] = {0, 1, 2, 0, 1, 2, 0, 1, 2};
+    double Ax[9] = {1, 2, 3, 2, 5, 6, 3, 6, 9};
+    memcpy(A->p, Ap, 4 * sizeof(int));
+    memcpy(A->i, Ai, 9 * sizeof(int));
+    memcpy(A->x, Ax, 9 * sizeof(double));
+
+    COO_Matrix *coo = new_coo_matrix_lower_triangular(A);
+
+    double vals2[9] = {10, 20, 30, 20, 50, 60, 30, 60, 90};
+    refresh_lower_triangular_coo(coo, vals2);
+
+    double expected_x[6] = {10, 20, 50, 30, 60, 90};
+    mu_assert("refresh vals incorrect", cmp_double_array(coo->x, expected_x, 6));
+
+    free_coo_matrix(coo);
+    free_csr_matrix(A);
+
+    return 0;
+}