From 96d483a94223a155c512bb591bb0d898aeacd930 Mon Sep 17 00:00:00 2001 From: dance858 Date: Sat, 12 Sep 2026 15:30:36 -0700 Subject: [PATCH] problem_jacobian: copy constraint values through the CSR view problem_init_jacobian lays out the aggregated Jacobian from each constraint's to_csr view, but problem_jacobian copied values from the constraint matrix's own buffer. For a stacked_pd that buffer is block-major, which equals CSR row order only when the blocks cover contiguous row ranges; a constraint like transpose(A @ X) with a multi-column X has interleaved blocks and its rows were copied in the wrong order. Copy from the to_csr view instead. Sparse and permuted_dense views alias their buffer, so nothing changes for them; stacked_pd refreshes its cache once per changed evaluation under the existing version guard. Regression test builds exactly that constraint and checks the assembled Jacobian row by row; it fails on main with rows 1 and 2 swapped. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01B6Bs765i1HTbqx3LUu6LMP --- src/problem.c | 10 +++++-- tests/all_tests.c | 1 + tests/problem/test_problem.h | 54 ++++++++++++++++++++++++++++++++++++ 3 files changed, 62 insertions(+), 3 deletions(-) diff --git a/src/problem.c b/src/problem.c index 71fa2ed1..dc97e3de 100644 --- a/src/problem.c +++ b/src/problem.c @@ -535,11 +535,15 @@ void problem_jacobian(problem *prob) /* copy only when the constraint's jacobian values actually changed (an affine constraint's eval is a no-op after its first call and - leaves the version untouched until the next parameter update) */ + leaves the version untouched until the next parameter update). + Copy through the CSR view: J's rows were laid out from to_csr in + problem_init_jacobian, and a stacked_pd's own value buffer is + block-major, which is not row order when its blocks interleave + rows. For sparse and permuted_dense the view aliases the buffer. */ if (prob->constraint_jac_seen[i] != c->jacobian->values_version) { - memcpy(J->x + nnz_offset, c->jacobian->x, - c->jacobian->nnz * sizeof(double)); + const CSR_matrix *Jc = c->jacobian->to_csr(c->jacobian); + memcpy(J->x + nnz_offset, Jc->x, Jc->nnz * sizeof(double)); prob->constraint_jac_seen[i] = c->jacobian->values_version; } nnz_offset += c->jacobian->nnz; diff --git a/tests/all_tests.c b/tests/all_tests.c index ab6235e9..98eb767b 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -574,6 +574,7 @@ int main(void) mu_run_test(test_problem_gradient, tests_run); mu_run_test(test_problem_jacobian, tests_run); mu_run_test(test_problem_jacobian_multi, tests_run); + mu_run_test(test_problem_jacobian_spd_constraint_interleaved, tests_run); mu_run_test(test_problem_constraint_forward, tests_run); mu_run_test(test_problem_hessian, tests_run); mu_run_test(test_problem_hessian_sum_exp_left_matmul_dense_transpose, tests_run); diff --git a/tests/problem/test_problem.h b/tests/problem/test_problem.h index 0944a628..8a385960 100644 --- a/tests/problem/test_problem.h +++ b/tests/problem/test_problem.h @@ -416,4 +416,58 @@ const char *test_problem_hessian_sum_exp_left_matmul_dense_transpose(void) return 0; } +/* Regression: problem_jacobian must copy constraint values in CSR row order, + not in the constraint Jacobian's native buffer order. A stacked_pd stores + values block-major, which only coincides with row order when its blocks + cover contiguous row ranges. Constraint T = transpose(A @ X) with A 2x3 and + X a 3x2 variable: A @ X has a 2-block spd Jacobian (one block per column of + X, rows {0,1} and {2,3}); the transpose gathers rows {0,2} into block 0 and + {1,3} into block 1, so block-major order swaps rows 1 and 2. */ +const char *test_problem_jacobian_spd_constraint_interleaved(void) +{ + double A[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + + expr *X_obj = new_variable(3, 2, 0, 6); + expr *objective = new_sum(X_obj, -1); + + expr *X = new_variable(3, 2, 0, 6); + expr *L = new_left_matmul_dense(NULL, X, 2, 3, A); + expr *T = new_transpose(L); + expr *constraints[1] = {T}; + + problem *prob = new_problem(objective, constraints, 1, false); + problem_init_jacobian(prob); + + double u[6] = {0.1, 0.2, 0.3, -0.1, -0.2, -0.3}; + problem_constraint_forward(prob, u); + problem_jacobian(prob); + + mu_assert("constraint Jacobian should be spd", T->jacobian->is_stacked_pd); + + /* T[r] = L[k(r)] with k = {0, 2, 1, 3}; L's Jacobian row (i, j) carries + A[i, :] in the columns of X's j-th column (vars 3j .. 3j+2). */ + double expected[4][6] = {{1.0, 2.0, 3.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 1.0, 2.0, 3.0}, + {4.0, 5.0, 6.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 4.0, 5.0, 6.0}}; + CSR_matrix *J = prob->jacobian; + mu_assert("shape", J->m == 4 && J->n == 6 && J->nnz == 12); + double dense[4][6] = {{0}}; + for (int r = 0; r < 4; r++) + { + for (int jj = J->p[r]; jj < J->p[r + 1]; jj++) + { + dense[r][J->i[jj]] += J->x[jj]; + } + } + for (int r = 0; r < 4; r++) + { + mu_assert("problem Jacobian row must match CSR row order", + cmp_double_array(dense[r], expected[r], 6)); + } + + free_problem(prob); + return 0; +} + #endif /* TEST_PROBLEM_H */