diff --git a/.github/workflows/formatting.yml b/.github/workflows/formatting.yml index 2dbf8e2..5db2ae9 100644 --- a/.github/workflows/formatting.yml +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 55181f5..3611755 100644 --- a/CMakeLists.txt +++ b/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}") diff --git a/include/problem.h b/include/problem.h index 2462ffd..9d7a446 100644 --- a/include/problem.h +++ b/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 @@ -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); diff --git a/include/utils/COO_Matrix.h b/include/utils/COO_Matrix.h new file mode 100644 index 0000000..16a9aad --- /dev/null +++ b/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 */ diff --git a/src/problem.c b/src/problem.c index 3413858..6dbd1a6 100644 --- a/src/problem.c +++ b/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) */ diff --git a/src/utils/COO_Matrix.c b/src/utils/COO_Matrix.c new file mode 100644 index 0000000..1a39df4 --- /dev/null +++ b/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 +#include + +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); + } +} diff --git a/tests/all_tests.c b/tests/all_tests.c index 6aa699a..fece8de 100644 --- a/tests/all_tests.c +++ b/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); diff --git a/tests/utils/test_coo_matrix.h b/tests/utils/test_coo_matrix.h new file mode 100644 index 0000000..8c6f14e --- /dev/null +++ b/tests/utils/test_coo_matrix.h @@ -0,0 +1,104 @@ +#include +#include +#include + +#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; +}