From 6064f74f6964e8da3ad8407c31039789f96320fb Mon Sep 17 00:00:00 2001 From: dance858 Date: Fri, 11 Sep 2026 06:37:49 -0700 Subject: [PATCH 1/4] Replace index_alloc/fill_values with a row_gather primitive (lesson 3, M1) row_gather_alloc(A, map, m_out) / row_gather_fill_values(A, C) replace the atom-named index slots. map[i] in [-1, A->m), -1 = empty output row, repeats allowed; the map is bound to the result at alloc time (bound_iwork on sparse_matrix, generalizing transpose_iwork, and new on permuted_dense). index and transpose atoms migrated; new_permuted_dense checks m0*n0 overflow unconditionally. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_012QaXUu4VhNvNyBc2dCzpkp --- include/utils/matrix.h | 21 +- include/utils/permuted_dense.h | 18 +- include/utils/sparse_matrix.h | 10 +- src/atoms/affine/index.c | 9 +- src/atoms/affine/transpose.c | 10 +- src/utils/permuted_dense.c | 61 +++-- src/utils/sparse_matrix.c | 69 +++-- src/utils/stacked_pd.c | 38 +-- tests/all_tests.c | 15 +- tests/jacobian_tests/affine/test_transpose.h | 53 ++++ tests/test_helpers.h | 18 ++ tests/utils/test_alloc_overflow.h | 15 +- tests/utils/test_matmul_dispatchers.h | 17 +- tests/utils/test_permuted_dense.h | 46 ++-- tests/utils/test_row_gather.h | 253 +++++++++++++++++++ tests/utils/test_stacked_pd.h | 42 +-- 16 files changed, 536 insertions(+), 159 deletions(-) create mode 100644 tests/utils/test_row_gather.h diff --git a/include/utils/matrix.h b/include/utils/matrix.h index 1aad9322..f64d270b 100644 --- a/include/utils/matrix.h +++ b/include/utils/matrix.h @@ -82,12 +82,17 @@ typedef CSR_matrix *(*matrix_to_csr_fn)(matrix *A); cache already matches values_version, so it is cheap to call when fresh. */ typedef void (*matrix_refresh_csc_values_fn)(matrix *A); -/* Allocate C = A[indices, :] */ -typedef matrix *(*matrix_index_alloc_fn)(matrix *A, const int *indices, int n_idxs); - -/* Fill values of C = A[indices, :] */ -typedef void (*matrix_index_fill_values_fn)(matrix *A, const int *indices, - int n_idxs, matrix *C); +/* Row gather: allocate C of shape (m_out, A->n) with C[i, :] = A[map[i], :] + for map[i] in [0, A->m), and C[i, :] structurally empty when map[i] == -1. + Repeated entries are allowed, so C->nnz may exceed A->nnz. Whatever the fill + needs from map is copied onto C here; the caller may free map afterwards. */ +typedef matrix *(*matrix_row_gather_alloc_fn)(const matrix *A, const int *map, + int m_out); + +/* Fill values of C = A[map, :]. C must be the matrix returned by + row_gather_alloc(A, map, m_out): the gather state is bound to it at alloc + time, and copy_sparsity copies carry none. */ +typedef void (*matrix_row_gather_fill_values_fn)(const matrix *A, matrix *C); /* Row-tiling for the promote atom: A must be a 1-row matrix; returns a new matrix of shape (size, A->n) where every row is a copy of A's @@ -156,8 +161,8 @@ struct matrix matrix_refresh_csc_values_fn refresh_csc_values; /* Atom-specific ops */ - matrix_index_alloc_fn index_alloc; - matrix_index_fill_values_fn index_fill_values; + matrix_row_gather_alloc_fn row_gather_alloc; + matrix_row_gather_fill_values_fn row_gather_fill_values; matrix_promote_alloc_fn promote_alloc; matrix_promote_fill_values_fn promote_fill_values; matrix_broadcast_alloc_fn broadcast_alloc; diff --git a/include/utils/permuted_dense.h b/include/utils/permuted_dense.h index 671a0fc1..eab01e5c 100644 --- a/include/utils/permuted_dense.h +++ b/include/utils/permuted_dense.h @@ -66,6 +66,12 @@ typedef struct permuted_dense int *kernel_iwork; size_t kernel_iwork_size; + /* Int state bound by the alloc that produced this PD and read by the + matching fill: row_gather_alloc stores, per dense row of this PD, the + dense row of the source it copies (length m0). NULL otherwise; never + touched by any other kernel (unlike kernel_iwork). */ + int *bound_iwork; + /* Cached transpose of this PD as another permuted_dense, allocated lazily on first call to permuted_dense_ensure_transpose_cache. On the cache PD itself, transpose_seen records the source's base.values_version whose @@ -96,12 +102,14 @@ matrix *broadcast_pd_alloc(const permuted_dense *A, broadcast_type type, int d1, void broadcast_pd_fill_values(const permuted_dense *A, broadcast_type type, int d1, int d2, permuted_dense *C); -/* Allocate C = A[indices, :], where A and C are permuted dense. */ -matrix *index_pd_alloc(const permuted_dense *A, const int *indices, int n_idxs); +/* Allocate C = A[map, :], where A and C are permuted dense. See + matrix_row_gather_alloc_fn for the map contract (-1 = empty row, repeats + allowed). C's row_perm is the set of output positions whose map entry hits + A->row_perm; C->bound_iwork holds the matching source dense rows. */ +matrix *row_gather_pd_alloc(const permuted_dense *A, const int *map, int m_out); -/* Fill values of C = A[indices, :]. */ -void index_pd_fill_values(const permuted_dense *A, const int *indices, int n_idxs, - permuted_dense *C); +/* Fill values of C = A[map, :]; C must come from row_gather_pd_alloc(A, ...). */ +void row_gather_pd_fill_values(const permuted_dense *A, permuted_dense *C); /* Allocate C = promote(A, size), where A and C are permuted dense. */ matrix *promote_pd_alloc(const permuted_dense *A, int size); diff --git a/include/utils/sparse_matrix.h b/include/utils/sparse_matrix.h index 85ee3f35..c281a4ce 100644 --- a/include/utils/sparse_matrix.h +++ b/include/utils/sparse_matrix.h @@ -30,10 +30,12 @@ typedef struct sparse_matrix CSC_matrix *csc_cache; uint64_t csc_seen; /* base.values_version the csc_cache values reflect */ int *csc_iwork; - int *transpose_iwork; /* sized csr->n; allocated by sparse_transpose_alloc - on the output sm and reused by - sparse_transpose_fill_values. NULL when this - sm wasn't produced by transpose_alloc. */ + /* Int state bound by the alloc that produced this matrix and read by the + matching fill: transpose_alloc stores csr->n scratch for + transpose_fill_values; row_gather_alloc stores the base.m-long row map + for row_gather_fill_values. NULL otherwise; never touched by any other + kernel. */ + int *bound_iwork; } sparse_matrix; /* Constructor. Takes ownership of A; the caller must not free A separately diff --git a/src/atoms/affine/index.c b/src/atoms/affine/index.c index 28baa1d0..776e02e2 100644 --- a/src/atoms/affine/index.c +++ b/src/atoms/affine/index.c @@ -65,20 +65,19 @@ static void jacobian_init_impl(expr *node) jacobian_init(x); /* allocate sparsity pattern for the matrix consisting of rows - 'idx->indices' of the child's Jacobian */ + 'idx->indices' of the child's Jacobian; the gather map is bound to + node->jacobian, so eval only needs the two matrices */ node->jacobian = - x->jacobian->index_alloc(x->jacobian, idx->indices, idx->n_idxs); + x->jacobian->row_gather_alloc(x->jacobian, idx->indices, idx->n_idxs); } static void eval_jacobian_impl(expr *node) { expr *x = node->left; - index_expr *idx = (index_expr *) node; eval_jacobian(x); /* copy values of the selected rows into the preallocated output */ - x->jacobian->index_fill_values(x->jacobian, idx->indices, idx->n_idxs, - node->jacobian); + x->jacobian->row_gather_fill_values(x->jacobian, node->jacobian); } static void wsum_hess_init_impl(expr *node) diff --git a/src/atoms/affine/transpose.c b/src/atoms/affine/transpose.c index 40bf61de..95695f11 100644 --- a/src/atoms/affine/transpose.c +++ b/src/atoms/affine/transpose.c @@ -50,18 +50,18 @@ static void jacobian_init_impl(expr *node) indices[r] = (r / d1) + (r % d1) * d2; } - node->jacobian = child->jacobian->index_alloc(child->jacobian, indices, n_out); + node->jacobian = + child->jacobian->row_gather_alloc(child->jacobian, indices, n_out); - /* save indices for eval_jacobian */ - node->work->iwork = indices; + /* the gather map is bound to node->jacobian; nothing to keep */ + sp_free(indices); } static void eval_jacobian_impl(expr *node) { expr *child = node->left; eval_jacobian(child); - child->jacobian->index_fill_values(child->jacobian, node->work->iwork, - node->size, node->jacobian); + child->jacobian->row_gather_fill_values(child->jacobian, node->jacobian); } static void wsum_hess_init_impl(expr *node) diff --git a/src/utils/permuted_dense.c b/src/utils/permuted_dense.c index ed7dca50..f2fbdd4e 100644 --- a/src/utils/permuted_dense.c +++ b/src/utils/permuted_dense.c @@ -23,6 +23,7 @@ #include "utils/tracked_alloc.h" #include "utils/utils.h" #include +#include #include #include #include @@ -47,6 +48,7 @@ static void permuted_dense_free(matrix *self) } sp_free(pd->kernel_dwork); sp_free(pd->kernel_iwork); + sp_free(pd->bound_iwork); free_matrix((matrix *) pd->transpose_cache); sp_free(pd); } @@ -108,49 +110,58 @@ static void permuted_dense_vtable_transpose_fill_values(const matrix *self, transpose_pd_fill_values((const permuted_dense *) self, (permuted_dense *) out); } -matrix *index_pd_alloc(const permuted_dense *A, const int *indices, int n_idxs) +matrix *row_gather_pd_alloc(const permuted_dense *A, const int *map, int m_out) { - /* Scan indices: which output positions i hit a row in A->row_perm? */ - int *new_row_perm = (int *) sp_malloc(n_idxs * sizeof(int)); + /* Output position i is dense iff map[i] hits a row in A->row_perm. The kept + positions form C's row_perm (strictly increasing by construction, repeats + in map included); src[k] is the dense row of A that C's row k copies. */ + int *new_row_perm = (int *) sp_malloc(m_out * sizeof(int)); + int *src = (int *) sp_malloc(m_out * sizeof(int)); int new_m0 = 0; - for (int i = 0; i < n_idxs; i++) + for (int i = 0; i < m_out; i++) { - if (A->row_inv[indices[i]] >= 0) + int ii = map[i] < 0 ? -1 : A->row_inv[map[i]]; + if (ii >= 0) { - new_row_perm[new_m0++] = i; + new_row_perm[new_m0] = i; + src[new_m0] = ii; + new_m0++; } } - matrix *out = new_permuted_dense(n_idxs, A->base.n, new_m0, A->n0, new_row_perm, + matrix *out = new_permuted_dense(m_out, A->base.n, new_m0, A->n0, new_row_perm, A->col_perm, NULL); + if (new_m0 > 0) + { + permuted_dense *C = (permuted_dense *) out; + C->bound_iwork = (int *) sp_malloc(new_m0 * sizeof(int)); + memcpy(C->bound_iwork, src, new_m0 * sizeof(int)); + } sp_free(new_row_perm); + sp_free(src); return out; } -void index_pd_fill_values(const permuted_dense *A, const int *indices, int n_idxs, - permuted_dense *C) +void row_gather_pd_fill_values(const permuted_dense *A, permuted_dense *C) { - (void) n_idxs; + assert(A->base.n == C->base.n && (C->m0 == 0 || C->bound_iwork != NULL)); int n0 = A->n0; for (int k = 0; k < C->m0; k++) { - int i = C->row_perm[k]; - int old_ii = A->row_inv[indices[i]]; - memcpy(C->X + k * n0, A->X + old_ii * n0, n0 * sizeof(double)); + memcpy(C->X + k * n0, A->X + C->bound_iwork[k] * n0, n0 * sizeof(double)); } } -static matrix *permuted_dense_vtable_index_alloc(matrix *self, const int *indices, - int n_idxs) +static matrix *permuted_dense_vtable_row_gather_alloc(const matrix *self, + const int *map, int m_out) { - return index_pd_alloc((const permuted_dense *) self, indices, n_idxs); + return row_gather_pd_alloc((const permuted_dense *) self, map, m_out); } -static void permuted_dense_vtable_index_fill_values(matrix *self, const int *indices, - int n_idxs, matrix *out) +static void permuted_dense_vtable_row_gather_fill_values(const matrix *self, + matrix *out) { - index_pd_fill_values((const permuted_dense *) self, indices, n_idxs, - (permuted_dense *) out); + row_gather_pd_fill_values((const permuted_dense *) self, (permuted_dense *) out); } matrix *promote_pd_alloc(const permuted_dense *A, int size) @@ -554,8 +565,8 @@ static void wire_vtable(permuted_dense *pd) pd->base.to_csr = permuted_dense_to_csr; pd->base.transpose_alloc = permuted_dense_vtable_transpose_alloc; pd->base.transpose_fill_values = permuted_dense_vtable_transpose_fill_values; - pd->base.index_alloc = permuted_dense_vtable_index_alloc; - pd->base.index_fill_values = permuted_dense_vtable_index_fill_values; + pd->base.row_gather_alloc = permuted_dense_vtable_row_gather_alloc; + pd->base.row_gather_fill_values = permuted_dense_vtable_row_gather_fill_values; pd->base.promote_alloc = permuted_dense_vtable_promote_alloc; pd->base.promote_fill_values = permuted_dense_vtable_promote_fill_values; pd->base.broadcast_alloc = permuted_dense_vtable_broadcast_alloc; @@ -590,6 +601,12 @@ matrix *new_permuted_dense(int m, int n, int m0, int n0, const int *row_perm, permuted_dense *pd = (permuted_dense *) sp_calloc(1, sizeof(permuted_dense)); pd->base.m = m; pd->base.n = n; + if (m0 > 0 && n0 > INT_MAX / m0) + { + fprintf(stderr, "Error in new_permuted_dense: dense block m0 * n0 " + "exceeds INT_MAX.\n"); + exit(1); + } pd->base.nnz = m0 * n0; wire_vtable(pd); diff --git a/src/utils/sparse_matrix.c b/src/utils/sparse_matrix.c index 8a7162f2..ff4f956f 100644 --- a/src/utils/sparse_matrix.c +++ b/src/utils/sparse_matrix.c @@ -24,6 +24,7 @@ #include "utils/mini_numpy.h" #include "utils/tracked_alloc.h" #include "utils/utils.h" +#include #include #include #include @@ -55,7 +56,7 @@ static void sparse_free(matrix *self) free_CSR_matrix(sm->csr); free_CSC_matrix(sm->csc_cache); sp_free(sm->csc_iwork); - sp_free(sm->transpose_iwork); + sp_free(sm->bound_iwork); sp_free(sm); } @@ -117,7 +118,7 @@ static matrix *sparse_transpose_alloc(const matrix *self) int *iwork = (int *) sp_malloc(sm->csr->n * sizeof(int)); CSR_matrix *AT = AT_alloc(sm->csr, iwork); sparse_matrix *out = (sparse_matrix *) new_sparse_matrix(AT); - out->transpose_iwork = iwork; + out->bound_iwork = iwork; return &out->base; } @@ -125,52 +126,68 @@ static void sparse_transpose_fill_values(const matrix *self, matrix *out) { const sparse_matrix *sm_in = (const sparse_matrix *) self; sparse_matrix *sm_out = (sparse_matrix *) out; - AT_fill_values(sm_in->csr, sm_out->csr, sm_out->transpose_iwork); + AT_fill_values(sm_in->csr, sm_out->csr, sm_out->bound_iwork); } -static matrix *sparse_index_alloc(matrix *self, const int *indices, int n_idxs) +static matrix *sparse_row_gather_alloc(const matrix *self, const int *map, int m_out) { - CSR_matrix *Jx = ((sparse_matrix *) self)->csr; + const CSR_matrix *Jx = ((const sparse_matrix *) self)->csr; /* Exact output nnz: sum the selected rows' nnz. Jx->nnz is NOT an upper - bound — duplicated indices select the same source row more than once - (cvxpy#3442). Duplicated gathers of dense rows can push the true count - past INT_MAX, which a CSR cannot represent, so fail before wrapping. */ + bound — repeated map entries select the same source row more than once + (cvxpy#3442). Repeated gathers of dense rows can push the true count + past INT_MAX, which a CSR cannot represent, so fail before wrapping. + map[i] == -1 selects nothing and must be branched before touching p. */ int nnz = 0; - for (int i = 0; i < n_idxs; i++) + for (int i = 0; i < m_out; i++) { - int len = Jx->p[indices[i] + 1] - Jx->p[indices[i]]; + int row = map[i]; + int len = row < 0 ? 0 : Jx->p[row + 1] - Jx->p[row]; if (len > INT_MAX - nnz) { - fprintf(stderr, "Error in sparse_index_alloc: gathered nnz " + fprintf(stderr, "Error in sparse_row_gather_alloc: gathered nnz " "exceeds INT_MAX.\n"); exit(1); } nnz += len; } - CSR_matrix *J = new_CSR_matrix(n_idxs, self->n, nnz); + CSR_matrix *J = new_CSR_matrix(m_out, self->n, nnz); J->p[0] = 0; - for (int i = 0; i < n_idxs; i++) + for (int i = 0; i < m_out; i++) { - int row = indices[i]; - int len = Jx->p[row + 1] - Jx->p[row]; - memcpy(J->i + J->p[i], Jx->i + Jx->p[row], len * sizeof(int)); + int row = map[i]; + int len = row < 0 ? 0 : Jx->p[row + 1] - Jx->p[row]; + if (len > 0) + { + memcpy(J->i + J->p[i], Jx->i + Jx->p[row], len * sizeof(int)); + } J->p[i + 1] = J->p[i] + len; } - J->nnz = J->p[n_idxs]; - return new_sparse_matrix(J); + J->nnz = J->p[m_out]; + + sparse_matrix *out = (sparse_matrix *) new_sparse_matrix(J); + out->bound_iwork = (int *) sp_malloc(m_out * sizeof(int)); + if (m_out > 0) + { + memcpy(out->bound_iwork, map, m_out * sizeof(int)); + } + return &out->base; } -static void sparse_index_fill_values(matrix *self, const int *indices, int n_idxs, - matrix *out) +static void sparse_row_gather_fill_values(const matrix *self, matrix *out) { - CSR_matrix *Jx = ((sparse_matrix *) self)->csr; - CSR_matrix *J = ((sparse_matrix *) out)->csr; - for (int i = 0; i < n_idxs; i++) + const CSR_matrix *Jx = ((const sparse_matrix *) self)->csr; + sparse_matrix *sm_out = (sparse_matrix *) out; + CSR_matrix *J = sm_out->csr; + const int *map = sm_out->bound_iwork; + assert(map != NULL && self->n == out->n); + for (int i = 0; i < J->m; i++) { + int row = map[i]; + if (row < 0) continue; int len = J->p[i + 1] - J->p[i]; - memcpy(J->x + J->p[i], Jx->x + Jx->p[indices[i]], len * sizeof(double)); + memcpy(J->x + J->p[i], Jx->x + Jx->p[row], len * sizeof(double)); } } @@ -389,8 +406,8 @@ static void wire_vtable(sparse_matrix *sm) sm->base.to_csr = sparse_to_csr; sm->base.transpose_alloc = sparse_transpose_alloc; sm->base.transpose_fill_values = sparse_transpose_fill_values; - sm->base.index_alloc = sparse_index_alloc; - sm->base.index_fill_values = sparse_index_fill_values; + sm->base.row_gather_alloc = sparse_row_gather_alloc; + sm->base.row_gather_fill_values = sparse_row_gather_fill_values; sm->base.promote_alloc = sparse_promote_alloc; sm->base.promote_fill_values = sparse_promote_fill_values; sm->base.broadcast_alloc = sparse_broadcast_alloc; diff --git a/src/utils/stacked_pd.c b/src/utils/stacked_pd.c index f29bcc2f..a178279a 100644 --- a/src/utils/stacked_pd.c +++ b/src/utils/stacked_pd.c @@ -231,37 +231,39 @@ void compose_csr_idx_map_for_spd(const stacked_pd *spd, const CSR_matrix *csr, } // ----------------------------------------------------------------------------- -// index of stacked_pd: C = A[indices, :] where A is stacked_pd +// row gather of stacked_pd: C = A[map, :] where A is stacked_pd // ----------------------------------------------------------------------------- typedef struct { - const int *indices; - int n_idxs; -} pd_index_ctx; + const int *map; + int m_out; +} pd_row_gather_ctx; -static matrix *wrapper_pd_index(permuted_dense *Bk, const void *ctx) +static matrix *wrapper_pd_row_gather(permuted_dense *Bk, const void *ctx) { - const pd_index_ctx *c = (const pd_index_ctx *) ctx; - return index_pd_alloc(Bk, c->indices, c->n_idxs); + const pd_row_gather_ctx *c = (const pd_row_gather_ctx *) ctx; + return row_gather_pd_alloc(Bk, c->map, c->m_out); } -static matrix *stacked_pd_vtable_index_alloc(matrix *self, const int *indices, - int n_idxs) +/* Every output row has exactly one source row, which lives in exactly one + block, so the per-block results stay row-disjoint even when map repeats. */ +static matrix *stacked_pd_vtable_row_gather_alloc(const matrix *self, const int *map, + int m_out) { - stacked_pd *src = (stacked_pd *) self; - pd_index_ctx ctx = {indices, n_idxs}; - return spd_map_filter_blocks(src, n_idxs, src->base.n, wrapper_pd_index, &ctx); + const stacked_pd *src = (const stacked_pd *) self; + pd_row_gather_ctx ctx = {map, m_out}; + return spd_map_filter_blocks(src, m_out, src->base.n, wrapper_pd_row_gather, + &ctx); } -static void stacked_pd_vtable_index_fill_values(matrix *self, const int *indices, - int n_idxs, matrix *out) +static void stacked_pd_vtable_row_gather_fill_values(const matrix *self, matrix *out) { - stacked_pd *src = (stacked_pd *) self; + const stacked_pd *src = (const stacked_pd *) self; stacked_pd *out_spd = (stacked_pd *) out; for (int k = 0; k < out_spd->n_blocks; k++) { int sk = out_spd->src_block_idx[k]; - index_pd_fill_values(src->blocks[sk], indices, n_idxs, out_spd->blocks[k]); + row_gather_pd_fill_values(src->blocks[sk], out_spd->blocks[k]); } } @@ -492,8 +494,8 @@ static void wire_vtable(stacked_pd *spd) spd->base.transpose_fill_values = stacked_pd_vtable_transpose_fill_values; spd->base.refresh_csc_values = stacked_pd_vtable_refresh_csc_values; spd->base.to_csr = stacked_pd_to_csr; - spd->base.index_alloc = stacked_pd_vtable_index_alloc; - spd->base.index_fill_values = stacked_pd_vtable_index_fill_values; + spd->base.row_gather_alloc = stacked_pd_vtable_row_gather_alloc; + spd->base.row_gather_fill_values = stacked_pd_vtable_row_gather_fill_values; spd->base.promote_alloc = stacked_pd_vtable_promote_alloc; spd->base.promote_fill_values = stacked_pd_vtable_promote_fill_values; spd->base.diag_vec_alloc = stacked_pd_vtable_diag_vec_alloc; diff --git a/tests/all_tests.c b/tests/all_tests.c index 5e0a9e40..1d504b97 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -74,6 +74,7 @@ #include "utils/test_matmul_dispatchers.h" #include "utils/test_matrix.h" #include "utils/test_permuted_dense.h" +#include "utils/test_row_gather.h" #include "utils/test_stacked_pd.h" #include "wsum_hess/affine/test_broadcast.h" #include "wsum_hess/affine/test_convolve.h" @@ -268,6 +269,7 @@ int main(void) mu_run_test(test_jacobian_kron_composite, tests_run); mu_run_test(test_jacobian_transpose, tests_run); mu_run_test(test_jacobian_transpose_pd_preserved, tests_run); + mu_run_test(test_jacobian_transpose_spd_preserved, tests_run); mu_run_test(test_diag_mat_jacobian_variable, tests_run); mu_run_test(test_diag_mat_jacobian_of_log, tests_run); mu_run_test(test_upper_tri_jacobian_variable, tests_run); @@ -380,7 +382,7 @@ int main(void) printf("\n--- Utility Tests ---\n"); mu_run_test(test_sat_mul_int_clamps_on_overflow, tests_run); - mu_run_test(test_sparse_index_alloc_no_int_overflow, tests_run); + mu_run_test(test_row_gather_alloc_no_int_overflow, tests_run); mu_run_test(test_cblas_ddot, tests_run); mu_run_test(test_diag_csr_mult, tests_run); mu_run_test(test_csr_sum, tests_run); @@ -435,7 +437,14 @@ int main(void) mu_run_test(test_permuted_dense_times_csc_no_active, tests_run); mu_run_test(test_permuted_dense_to_csr_lazy, tests_run); mu_run_test(test_permuted_dense_col_inv, tests_run); - mu_run_test(test_permuted_dense_index, tests_run); + mu_run_test(test_permuted_dense_row_gather, tests_run); + mu_run_test(test_row_gather_sparse, tests_run); + mu_run_test(test_row_gather_sparse_all_empty, tests_run); + mu_run_test(test_row_gather_pd_vs_sparse_twin, tests_run); + mu_run_test(test_row_gather_spd_vs_sparse_twin, tests_run); +#ifdef SP_TRACK_MEMORY + mu_run_test(test_row_gather_spd_fill_no_transient_alloc, tests_run); +#endif mu_run_test(test_permuted_dense_promote, tests_run); mu_run_test(test_permuted_dense_broadcast_scalar, tests_run); mu_run_test(test_permuted_dense_broadcast_row, tests_run); @@ -550,7 +559,7 @@ int main(void) mu_run_test(test_spd_vtable_ATDA_fill_values, tests_run); mu_run_test(test_spd_vtable_transpose, tests_run); mu_run_test(test_spd_vtable_refresh_csc_values_noop, tests_run); - mu_run_test(test_spd_vtable_index, tests_run); + mu_run_test(test_spd_vtable_row_gather, tests_run); mu_run_test(test_spd_vtable_promote, tests_run); mu_run_test(test_spd_vtable_diag_vec, tests_run); mu_run_test(test_spd_vtable_broadcast_row, tests_run); diff --git a/tests/jacobian_tests/affine/test_transpose.h b/tests/jacobian_tests/affine/test_transpose.h index 40b1b763..244b7065 100644 --- a/tests/jacobian_tests/affine/test_transpose.h +++ b/tests/jacobian_tests/affine/test_transpose.h @@ -90,4 +90,57 @@ const char *test_jacobian_transpose_pd_preserved(void) return 0; } +/* When the child Jacobian is a stacked_pd (left_matmul_dense with a multi-column + variable), the transpose must stay stacked_pd and equal the row-permuted child + Jacobian. X is 3x2, A is 2x3, L = A @ X is 2x2, so d1 = d2 = 2 and + k(r) = (r/2) + (r%2)*2 = [0, 2, 1, 3]. */ +const char *test_jacobian_transpose_spd_preserved(void) +{ + double A[6] = {1.0, -0.5, 2.0, 0.5, 1.5, -1.0}; + expr *X = new_variable(3, 2, 0, 6); + expr *L = new_left_matmul_dense(NULL, X, 2, 3, A); + expr *T = new_transpose(L); + + double u[6] = {0.1, 0.2, 0.3, -0.1, -0.2, -0.3}; + jacobian_init(T); + T->forward(T, u); + eval_jacobian(T); + mu_assert("child Jacobian should be spd", L->jacobian->is_stacked_pd); + mu_assert("transpose Jacobian should be spd", T->jacobian->is_stacked_pd); + + /* Dense child Jacobian: output entry (i, j) (col-major row i + 2j) carries + A[i, :] in the columns of X's j-th column (vars 3j .. 3j+2). */ + double JL[4][6] = {{0}}; + for (int j = 0; j < 2; j++) + { + for (int i = 0; i < 2; i++) + { + for (int c = 0; c < 3; c++) + { + JL[i + 2 * j][3 * j + c] = A[3 * i + c]; + } + } + } + int k[4] = {0, 2, 1, 3}; + + CSR_matrix *JT = T->jacobian->to_csr(T->jacobian); + mu_assert("nnz", JT->nnz == 12); + double dense[4][6] = {{0}}; + for (int r = 0; r < 4; r++) + { + for (int jj = JT->p[r]; jj < JT->p[r + 1]; jj++) + { + dense[r][JT->i[jj]] += JT->x[jj]; + } + } + for (int r = 0; r < 4; r++) + { + mu_assert("transpose rows must be permuted child rows", + cmp_double_array(dense[r], JL[k[r]], 6)); + } + + free_expr(T); + return 0; +} + #endif // TEST_TRANSPOSE_H diff --git a/tests/test_helpers.h b/tests/test_helpers.h index a16ec8d5..cead42cd 100644 --- a/tests/test_helpers.h +++ b/tests/test_helpers.h @@ -31,4 +31,22 @@ int csr_is_valid(const CSR_matrix *A); * in [0, 1]. Nonzero values are standard Gaussian (Box-Muller). */ CSR_matrix *new_csr_random(int m, int n, double density); +/* Only available with -DSP_TRACK_MEMORY=ON: reads the tracked allocator + * counters, which do not exist in a default build. */ +#ifdef SP_TRACK_MEMORY +#include "utils/tracked_alloc.h" + +/* No-alloc-in-fill contract: after alloc and one warm-up fill, a second fill + * must not touch the tracked allocator at all. Any transient sp_malloc inside + * the fill raises g_peak_bytes above the baseline even if freed before + * returning; a permanent one raises g_allocated_bytes. */ +static inline int fill_is_alloc_free(void (*fill)(const void *ctx), const void *ctx) +{ + size_t base = g_allocated_bytes; + g_peak_bytes = base; + fill(ctx); + return g_allocated_bytes == base && g_peak_bytes == base; +} +#endif + #endif /* TEST_HELPERS_H */ diff --git a/tests/utils/test_alloc_overflow.h b/tests/utils/test_alloc_overflow.h index 5e783ee1..b403512b 100644 --- a/tests/utils/test_alloc_overflow.h +++ b/tests/utils/test_alloc_overflow.h @@ -40,12 +40,13 @@ const char *test_sat_mul_int_clamps_on_overflow(void) return 0; } -/* sparse_index_alloc sized its allocation as MIN(Jx->nnz, n_idxs * self->n). For a - large jacobian the dense product n_idxs * self->n overflows int (e.g. a transpose - of a 250000-variable matrix), wrapping negative so MIN selected it -> calloc - overflow -> SIGSEGV. This builds an index op whose product overflows and checks - the result is a valid CSR with the true (subset) nnz. */ -const char *test_sparse_index_alloc_no_int_overflow(void) +/* The sparse row gather (then sparse_index_alloc) once sized its allocation as + MIN(Jx->nnz, n_idxs * self->n). For a large jacobian the dense product + n_idxs * self->n overflows int (e.g. a transpose of a 250000-variable matrix), + wrapping negative so MIN selected it -> calloc overflow -> SIGSEGV. This builds + a gather whose product overflows and checks the result is a valid CSR with the + true (subset) nnz. */ +const char *test_row_gather_alloc_no_int_overflow(void) { /* n_idxs * self->n = 25000 * 100000 = 2.5e9 overflows int32. */ const int n_idxs = 25000; @@ -67,7 +68,7 @@ const char *test_sparse_index_alloc_no_int_overflow(void) for (int i = 0; i < n_idxs; i++) indices[i] = i; /* pre-fix: crashes here in new_CSR_matrix(... negative nnz) */ - matrix *out = mat->index_alloc(mat, indices, n_idxs); + matrix *out = mat->row_gather_alloc(mat, indices, n_idxs); CSR_matrix *out_csr = out->to_csr(out); mu_assert("nnz must equal selected-row total", out_csr->nnz == n_idxs); diff --git a/tests/utils/test_matmul_dispatchers.h b/tests/utils/test_matmul_dispatchers.h index ebd5eb3c..9b09eab4 100644 --- a/tests/utils/test_matmul_dispatchers.h +++ b/tests/utils/test_matmul_dispatchers.h @@ -2441,20 +2441,9 @@ const char *test_BA_pd_kron_spd_no_cache_staleness(void) allocator counters, which do not exist in a default build. */ #ifdef SP_TRACK_MEMORY -/* No-alloc-in-fill contract for the BTDA kernels: after alloc and one warm-up - fill, a second fill must not touch the tracked allocator at all. Any - transient sp_malloc inside the fill raises g_peak_bytes above the baseline - even if freed before returning; a permanent one raises g_allocated_bytes. - Covers BTDA_pd_pd (both the matching-row_perm and gather paths), - BTDA_pd_spd, and the blockwise BTDA_spd_pd. */ -static int fill_is_alloc_free(void (*fill)(const void *ctx), const void *ctx) -{ - size_t base = g_allocated_bytes; - g_peak_bytes = base; - fill(ctx); - return g_allocated_bytes == base && g_peak_bytes == base; -} - +/* No-alloc-in-fill contract for the BTDA kernels via fill_is_alloc_free + (tests/test_helpers.h). Covers BTDA_pd_pd (both the matching-row_perm and + gather paths), BTDA_pd_spd, and the blockwise BTDA_spd_pd. */ typedef struct { const permuted_dense *B; diff --git a/tests/utils/test_permuted_dense.h b/tests/utils/test_permuted_dense.h index 761708b7..585808fb 100644 --- a/tests/utils/test_permuted_dense.h +++ b/tests/utils/test_permuted_dense.h @@ -341,10 +341,10 @@ const char *test_permuted_dense_col_inv(void) return 0; } -/* PD index_alloc / index_fill_values: select rows from a PD; output must be - another PD with row_perm equal to the output positions where indices[i] - hit the source row_perm. */ -const char *test_permuted_dense_index(void) +/* PD row_gather_alloc / row_gather_fill_values: output must be another PD whose + row_perm is the set of output positions where map[i] hits the source + row_perm, with repeats and -1 entries handled. */ +const char *test_permuted_dense_row_gather(void) { /* Source PD, shape (6, 4), dense block at rows {1, 3, 4} x cols {0, 2}. */ int row_perm[3] = {1, 3, 4}; @@ -352,34 +352,36 @@ const char *test_permuted_dense_index(void) double X[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; matrix *M = new_permuted_dense(6, 4, 3, 2, row_perm, col_perm, X); - /* Index by [0, 3, 1, 5, 4]: + /* map = [0, 3, 1, -1, 4, 3, 5]: - position 0 -> source row 0 (not in row_perm, zero) - - position 1 -> source row 3 (in row_perm at ii=1, dense) - - position 2 -> source row 1 (in row_perm at ii=0, dense) - - position 3 -> source row 5 (not in row_perm, zero) - - position 4 -> source row 4 (in row_perm at ii=2, dense) */ - int indices[5] = {0, 3, 1, 5, 4}; - matrix *out = M->index_alloc(M, indices, 5); + - position 1 -> source row 3 (ii=1, dense) + - position 2 -> source row 1 (ii=0, dense) + - position 3 -> -1 (structurally empty) + - position 4 -> source row 4 (ii=2, dense) + - position 5 -> source row 3 again (ii=1, dense) + - position 6 -> source row 5 (not in row_perm, zero) */ + int map[7] = {0, 3, 1, -1, 4, 3, 5}; + matrix *out = M->row_gather_alloc(M, map, 7); permuted_dense *out_pd = (permuted_dense *) out; - mu_assert("out m", out->m == 5); + mu_assert("out m", out->m == 7); mu_assert("out n", out->n == 4); - mu_assert("out nnz", out->nnz == 6); /* m0=3 * n0=2 */ - mu_assert("m0", out_pd->m0 == 3); + mu_assert("out nnz", out->nnz == 8); /* m0=4 * n0=2 */ + mu_assert("m0", out_pd->m0 == 4); mu_assert("n0", out_pd->n0 == 2); - int expected_row_perm[3] = {1, 2, 4}; - mu_assert("row_perm", cmp_int_array(out_pd->row_perm, expected_row_perm, 3)); + int expected_row_perm[4] = {1, 2, 4, 5}; + mu_assert("row_perm", cmp_int_array(out_pd->row_perm, expected_row_perm, 4)); int expected_col_perm[2] = {0, 2}; mu_assert("col_perm", cmp_int_array(out_pd->col_perm, expected_col_perm, 2)); + int expected_src[4] = {1, 0, 2, 1}; + mu_assert("bound_iwork", cmp_int_array(out_pd->bound_iwork, expected_src, 4)); - M->index_fill_values(M, indices, 5, out); + M->row_gather_fill_values(M, out); - /* Row 0 of out (i=1) = source row 3 = X[1, :] = {3, 4}. - Row 1 of out (i=2) = source row 1 = X[0, :] = {1, 2}. - Row 2 of out (i=4) = source row 4 = X[2, :] = {5, 6}. */ - double expected_X[6] = {3.0, 4.0, 1.0, 2.0, 5.0, 6.0}; - mu_assert("values", cmp_double_array(out_pd->X, expected_X, 6)); + /* dense rows of out = X[1, :], X[0, :], X[2, :], X[1, :] */ + double expected_X[8] = {3.0, 4.0, 1.0, 2.0, 5.0, 6.0, 3.0, 4.0}; + mu_assert("values", cmp_double_array(out_pd->X, expected_X, 8)); free_matrix(out); free_matrix(M); diff --git a/tests/utils/test_row_gather.h b/tests/utils/test_row_gather.h new file mode 100644 index 00000000..eceec09b --- /dev/null +++ b/tests/utils/test_row_gather.h @@ -0,0 +1,253 @@ +#ifndef TEST_ROW_GATHER_H +#define TEST_ROW_GATHER_H + +#include "minunit.h" +#include "test_helpers.h" +#include "utils/CSR_matrix.h" +#include "utils/permuted_dense.h" +#include "utils/sparse_matrix.h" +#include "utils/stacked_pd.h" +#include +#include +#include + +/* row_gather_alloc / row_gather_fill_values across the three matrix kinds. + The map contract under test: map[i] in [0, A->m) copies source row map[i], + map[i] == -1 leaves output row i structurally empty, repeats are allowed, + and the map is bound to the result at alloc time (the caller's copy is + dead afterwards). */ + +/* Shared 4x5 CSR source with an empty row 2: + row 0: (0: 1.0) (3: 2.0) + row 1: (1: 3.0) (4: 4.0) + row 2: empty + row 3: (0: 5.0) (2: 6.0) (4: 7.0) */ +static matrix *row_gather_sparse_fixture(void) +{ + CSR_matrix *A = new_CSR_matrix(4, 5, 7); + int p[5] = {0, 2, 4, 4, 7}; + int i[7] = {0, 3, 1, 4, 0, 2, 4}; + double x[7] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0}; + memcpy(A->p, p, 5 * sizeof(int)); + memcpy(A->i, i, 7 * sizeof(int)); + memcpy(A->x, x, 7 * sizeof(double)); + return new_sparse_matrix(A); +} + +/* 7x5 stacked_pd with three blocks; block 2 is never hit by the maps below so + the gather must drop it: + block 0: rows {0, 4}, cols {0, 2}, X = [[1, 2], [3, 4]] + block 1: rows {2, 3}, cols {1, 2}, X = [[5, 6], [7, 8]] + block 2: rows {6}, cols {4}, X = [[9]] + Row 1 and row 5 belong to no block. */ +static matrix *row_gather_spd_fixture(void) +{ + int rp0[2] = {0, 4}, cp0[2] = {0, 2}; + double X0[4] = {1.0, 2.0, 3.0, 4.0}; + int rp1[2] = {2, 3}, cp1[2] = {1, 2}; + double X1[4] = {5.0, 6.0, 7.0, 8.0}; + int rp2[1] = {6}, cp2[1] = {4}; + double X2[1] = {9.0}; + permuted_dense *blocks[3] = { + (permuted_dense *) new_permuted_dense(7, 5, 2, 2, rp0, cp0, X0), + (permuted_dense *) new_permuted_dense(7, 5, 2, 2, rp1, cp1, X1), + (permuted_dense *) new_permuted_dense(7, 5, 1, 1, rp2, cp2, X2)}; + return new_stacked_pd(7, 5, 3, blocks, NULL, NULL); +} + +static void row_gather_poison(matrix *M) +{ + for (int k = 0; k < M->nnz; k++) M->x[k] = (double) NAN; +} + +/* Deep-copy any matrix's CSR view into a fresh sparse_matrix. */ +static matrix *row_gather_sparse_twin(matrix *M) +{ + CSR_matrix *src = M->to_csr(M); + CSR_matrix *dst = new_CSR_matrix(src->m, src->n, src->nnz); + memcpy(dst->p, src->p, (size_t) (src->m + 1) * sizeof(int)); + memcpy(dst->i, src->i, (size_t) src->nnz * sizeof(int)); + memcpy(dst->x, src->x, (size_t) src->nnz * sizeof(double)); + dst->nnz = src->nnz; + return new_sparse_matrix(dst); +} + +/* Structure + values equality through the CSR views. */ +static int row_gather_same_csr(matrix *X, matrix *Y) +{ + CSR_matrix *a = X->to_csr(X); + CSR_matrix *b = Y->to_csr(Y); + if (a->m != b->m || a->n != b->n || a->nnz != b->nnz) return 0; + return cmp_int_array(a->p, b->p, a->m + 1) && + cmp_int_array(a->i, b->i, a->nnz) && cmp_double_array(a->x, b->x, a->nnz); +} + +/* map = [3, -1, 0, 3, 2, 1]: a permutation, a repeat (row 3 twice), a -1, and + the empty source row 2. Output is 6x5 with nnz 3+0+2+3+0+2 = 10. */ +const char *test_row_gather_sparse(void) +{ + matrix *A = row_gather_sparse_fixture(); + int map[6] = {3, -1, 0, 3, 2, 1}; + matrix *C = A->row_gather_alloc(A, map, 6); + + int exp_p[7] = {0, 3, 3, 5, 8, 8, 10}; + int exp_i[10] = {0, 2, 4, 0, 3, 0, 2, 4, 1, 4}; + mu_assert("shape", C->m == 6 && C->n == 5); + mu_assert("sparsity", cmp_sparsity(C, exp_p, exp_i, 6, 10)); + + /* the map is bound at alloc time: clobbering the caller's copy is fine */ + for (int k = 0; k < 6; k++) map[k] = -7; + + row_gather_poison(C); + A->row_gather_fill_values(A, C); + double exp_x[10] = {5.0, 6.0, 7.0, 1.0, 2.0, 5.0, 6.0, 7.0, 3.0, 4.0}; + mu_assert("values", cmp_values(C, exp_x, 10)); + + /* fill reads the source's current values, not a snapshot */ + for (int k = 0; k < A->nnz; k++) A->x[k] *= 10.0; + for (int k = 0; k < 10; k++) exp_x[k] *= 10.0; + A->row_gather_fill_values(A, C); + mu_assert("values after refill", cmp_values(C, exp_x, 10)); + + free_matrix(C); + free_matrix(A); + return 0; +} + +const char *test_row_gather_sparse_all_empty(void) +{ + matrix *A = row_gather_sparse_fixture(); + int map[3] = {-1, -1, -1}; + matrix *C = A->row_gather_alloc(A, map, 3); + int exp_p[4] = {0, 0, 0, 0}; + mu_assert("empty structure", cmp_sparsity(C, exp_p, NULL, 3, 0)); + A->row_gather_fill_values(A, C); /* nothing to write; must not touch A */ + free_matrix(C); + free_matrix(A); + return 0; +} + +/* pd 6x5 with rows {1, 3, 4} x cols {0, 2}. The map hits a row outside + row_perm (0 and 5), repeats row 3, and contains a -1. Result must stay a pd + and agree entrywise with the same gather on a sparse twin. */ +const char *test_row_gather_pd_vs_sparse_twin(void) +{ + int row_perm[3] = {1, 3, 4}; + int col_perm[2] = {0, 2}; + double X[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + matrix *A = new_permuted_dense(6, 5, 3, 2, row_perm, col_perm, X); + matrix *A_tw = row_gather_sparse_twin(A); + + int map[7] = {0, 3, -1, 1, 3, 5, 4}; + matrix *C = A->row_gather_alloc(A, map, 7); + matrix *C_tw = A_tw->row_gather_alloc(A_tw, map, 7); + mu_assert("kind preserved", C->is_permuted_dense); + mu_assert("m0", ((permuted_dense *) C)->m0 == 4); + mu_assert("nnz", C->nnz == 8); + + row_gather_poison(C); + A->row_gather_fill_values(A, C); + A_tw->row_gather_fill_values(A_tw, C_tw); + matrix_values_changed(C); + mu_assert("twin match", row_gather_same_csr(C, C_tw)); + + /* order-independent mutation of both sources, then refill */ + for (int k = 0; k < A->nnz; k++) A->x[k] *= 2.0; + for (int k = 0; k < A_tw->nnz; k++) A_tw->x[k] *= 2.0; + matrix_values_changed(A); + matrix_values_changed(A_tw); + A->row_gather_fill_values(A, C); + A_tw->row_gather_fill_values(A_tw, C_tw); + matrix_values_changed(C); + mu_assert("twin match after refill", row_gather_same_csr(C, C_tw)); + + free_matrix(C_tw); + free_matrix(C); + free_matrix(A_tw); + free_matrix(A); + return 0; +} + +/* spd fixture above; map = [4, 2, -1, 4, 0, 3, 1, 2]: + pos 0, 3 -> row 4 (block 0, twice) + pos 4 -> row 0 (block 0) + pos 1, 7 -> row 2 (block 1, twice) + pos 5 -> row 3 (block 1) + pos 2 -> -1 + pos 6 -> row 1 (no block) + Block 2 gets no hits and must be dropped; output blocks stay row-disjoint. */ +const char *test_row_gather_spd_vs_sparse_twin(void) +{ + matrix *A = row_gather_spd_fixture(); + matrix *A_tw = row_gather_sparse_twin(A); + + int map[8] = {4, 2, -1, 4, 0, 3, 1, 2}; + matrix *C = A->row_gather_alloc(A, map, 8); + matrix *C_tw = A_tw->row_gather_alloc(A_tw, map, 8); + mu_assert("kind preserved", C->is_stacked_pd); + stacked_pd *C_spd = (stacked_pd *) C; + mu_assert("empty block dropped", C_spd->n_blocks == 2); + mu_assert("src_block_idx", + C_spd->src_block_idx[0] == 0 && C_spd->src_block_idx[1] == 1); + int exp_rp0[3] = {0, 3, 4}; + int exp_rp1[3] = {1, 5, 7}; + mu_assert("block 0 row_perm", + C_spd->blocks[0]->m0 == 3 && + cmp_int_array(C_spd->blocks[0]->row_perm, exp_rp0, 3)); + mu_assert("block 1 row_perm", + C_spd->blocks[1]->m0 == 3 && + cmp_int_array(C_spd->blocks[1]->row_perm, exp_rp1, 3)); + + row_gather_poison(C); + A->row_gather_fill_values(A, C); + A_tw->row_gather_fill_values(A_tw, C_tw); + matrix_values_changed(C); + mu_assert("twin match", row_gather_same_csr(C, C_tw)); + + for (int k = 0; k < A->nnz; k++) A->x[k] *= 2.0; + for (int k = 0; k < A_tw->nnz; k++) A_tw->x[k] *= 2.0; + matrix_values_changed(A); + matrix_values_changed(A_tw); + A->row_gather_fill_values(A, C); + A_tw->row_gather_fill_values(A_tw, C_tw); + matrix_values_changed(C); + mu_assert("twin match after refill", row_gather_same_csr(C, C_tw)); + + free_matrix(C_tw); + free_matrix(C); + free_matrix(A_tw); + free_matrix(A); + return 0; +} + +#ifdef SP_TRACK_MEMORY +typedef struct +{ + matrix *A; + matrix *C; +} row_gather_fill_args; + +static void run_row_gather_fill(const void *ctx) +{ + const row_gather_fill_args *a = (const row_gather_fill_args *) ctx; + a->A->row_gather_fill_values(a->A, a->C); +} + +/* The spd fill dispatches through src_block_idx into per-block pd fills; none + of it may allocate once the alloc phase is done. */ +const char *test_row_gather_spd_fill_no_transient_alloc(void) +{ + matrix *A = row_gather_spd_fixture(); + int map[8] = {4, 2, -1, 4, 0, 3, 1, 2}; + matrix *C = A->row_gather_alloc(A, map, 8); + row_gather_fill_args args = {A, C}; + run_row_gather_fill(&args); /* warm-up */ + mu_assert("spd row_gather fill must not allocate", + fill_is_alloc_free(run_row_gather_fill, &args)); + free_matrix(C); + free_matrix(A); + return 0; +} +#endif + +#endif /* TEST_ROW_GATHER_H */ diff --git a/tests/utils/test_stacked_pd.h b/tests/utils/test_stacked_pd.h index 08e70f4e..5429e49d 100644 --- a/tests/utils/test_stacked_pd.h +++ b/tests/utils/test_stacked_pd.h @@ -1321,9 +1321,10 @@ const char *test_spd_vtable_refresh_csc_values_noop(void) return 0; } -/* index_* on spd: rows are routed to the source block (if any) that - carries them; per-block index_* is reused; empty blocks are dropped. */ -const char *test_spd_vtable_index(void) +/* row_gather on spd: each output row is routed to the source block (if any) + that carries its source row; per-block row_gather_pd_* is reused; empty + result blocks are dropped; repeats and -1 entries are allowed. */ +const char *test_spd_vtable_row_gather(void) { /* 6x4 spd: block 0: rows {0,1}, cols {0,2}, X = [[1,2],[3,4]] @@ -1341,23 +1342,24 @@ const char *test_spd_vtable_index(void) permuted_dense *blocks[2] = {(permuted_dense *) blk0, (permuted_dense *) blk1}; matrix *M = new_stacked_pd(6, 4, 2, blocks, NULL, NULL); - /* indices = [3, 0, 5, 1]: row 3 → blk1, row 0/1 → blk0, row 5 → none. - Expected output (4 rows, dropping output position 2): + /* map = [3, 0, 5, 1, -1, 3]: row 3 -> blk1 (twice, positions 0 and 5), + rows 0/1 -> blk0 (positions 1 and 3), row 5 -> none, position 4 = -1. + Expected output (6 rows, positions 2 and 4 empty): output row 0 = src row 3 = blk1 row [7, 8] in cols {1,3} output row 1 = src row 0 = blk0 row [1, 2] in cols {0,2} - output row 3 = src row 1 = blk0 row [3, 4] in cols {0,2} */ - int indices[4] = {3, 0, 5, 1}; - matrix *C_m = M->index_alloc(M, indices, 4); - M->index_fill_values(M, indices, 4, C_m); + output row 3 = src row 1 = blk0 row [3, 4] in cols {0,2} + output row 5 = src row 3 = blk1 row [7, 8] in cols {1,3} */ + int map[6] = {3, 0, 5, 1, -1, 3}; + matrix *C_m = M->row_gather_alloc(M, map, 6); + M->row_gather_fill_values(M, C_m); stacked_pd *C = (stacked_pd *) C_m; mu_assert("n_blocks", C->n_blocks == 2); - mu_assert("base.m", C_m->m == 4); + mu_assert("base.m", C_m->m == 6); mu_assert("base.n", C_m->n == 4); - /* Per the per-block alloc, blk0 yields output positions where indices - hits {0,1}: that's i=1 (indices[1]=0) and i=3 (indices[3]=1). So - output block 0 carries row_perm={1,3}, col_perm={0,2}, X=[1,2,3,4]. */ + /* blk0 is hit at positions 1 and 3: output block 0 carries + row_perm={1,3}, col_perm={0,2}, X=[1,2,3,4]. */ permuted_dense *out0 = C->blocks[0]; int expected_row_perm_0[2] = {1, 3}; int expected_col_perm_0[2] = {0, 2}; @@ -1370,19 +1372,19 @@ const char *test_spd_vtable_index(void) cmp_int_array(out0->col_perm, expected_col_perm_0, 2)); mu_assert("out0 X", cmp_double_array(out0->X, expected_X0, 4)); - /* blk1 hit only at i=0 (indices[0]=3 → blk1 row 1). Output block 1: - row_perm={0}, col_perm={1,3}, X=[7,8]. */ + /* blk1 is hit at positions 0 and 5, both from its row 1: output block 1 + carries row_perm={0,5}, col_perm={1,3}, X=[7,8,7,8]. */ permuted_dense *out1 = C->blocks[1]; - int expected_row_perm_1[1] = {0}; + int expected_row_perm_1[2] = {0, 5}; int expected_col_perm_1[2] = {1, 3}; - double expected_X1[2] = {7.0, 8.0}; - mu_assert("out1 m0", out1->m0 == 1); + double expected_X1[4] = {7.0, 8.0, 7.0, 8.0}; + mu_assert("out1 m0", out1->m0 == 2); mu_assert("out1 n0", out1->n0 == 2); mu_assert("out1 row_perm", - cmp_int_array(out1->row_perm, expected_row_perm_1, 1)); + cmp_int_array(out1->row_perm, expected_row_perm_1, 2)); mu_assert("out1 col_perm", cmp_int_array(out1->col_perm, expected_col_perm_1, 2)); - mu_assert("out1 X", cmp_double_array(out1->X, expected_X1, 2)); + mu_assert("out1 X", cmp_double_array(out1->X, expected_X1, 4)); /* src_block_idx: both source blocks survived, identity mapping. */ mu_assert("src_block_idx[0]", C->src_block_idx[0] == 0); From 26ce74b221b420eeb37e6164849dd33eaa7c070f Mon Sep 17 00:00:00 2001 From: dance858 Date: Sat, 12 Sep 2026 13:19:08 -0700 Subject: [PATCH 2/4] row_gather: drop -1 map entries to keep M1 minimal map[i] must be a valid source row in [0, A->m); the structurally-empty-row case (-1) has no caller until diag_vec migrates in M2 and is deferred there. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01Azb8o1AgiFRU4t9aAY7oDr --- include/utils/matrix.h | 6 +-- include/utils/permuted_dense.h | 6 +-- src/utils/permuted_dense.c | 2 +- src/utils/sparse_matrix.c | 16 +++---- tests/all_tests.c | 1 - tests/utils/test_permuted_dense.h | 19 ++++---- tests/utils/test_row_gather.h | 73 ++++++++++++------------------- tests/utils/test_stacked_pd.h | 22 +++++----- 8 files changed, 61 insertions(+), 84 deletions(-) diff --git a/include/utils/matrix.h b/include/utils/matrix.h index f64d270b..6f141b12 100644 --- a/include/utils/matrix.h +++ b/include/utils/matrix.h @@ -83,9 +83,9 @@ typedef CSR_matrix *(*matrix_to_csr_fn)(matrix *A); typedef void (*matrix_refresh_csc_values_fn)(matrix *A); /* Row gather: allocate C of shape (m_out, A->n) with C[i, :] = A[map[i], :] - for map[i] in [0, A->m), and C[i, :] structurally empty when map[i] == -1. - Repeated entries are allowed, so C->nnz may exceed A->nnz. Whatever the fill - needs from map is copied onto C here; the caller may free map afterwards. */ + for map[i] in [0, A->m). Repeated entries are allowed, so C->nnz may exceed + A->nnz. Whatever the fill needs from map is copied onto C here; the caller + may free map afterwards. */ typedef matrix *(*matrix_row_gather_alloc_fn)(const matrix *A, const int *map, int m_out); diff --git a/include/utils/permuted_dense.h b/include/utils/permuted_dense.h index eab01e5c..d94f8616 100644 --- a/include/utils/permuted_dense.h +++ b/include/utils/permuted_dense.h @@ -103,9 +103,9 @@ void broadcast_pd_fill_values(const permuted_dense *A, broadcast_type type, int int d2, permuted_dense *C); /* Allocate C = A[map, :], where A and C are permuted dense. See - matrix_row_gather_alloc_fn for the map contract (-1 = empty row, repeats - allowed). C's row_perm is the set of output positions whose map entry hits - A->row_perm; C->bound_iwork holds the matching source dense rows. */ + matrix_row_gather_alloc_fn for the map contract (repeats allowed). C's + row_perm is the set of output positions whose map entry hits A->row_perm; + C->bound_iwork holds the matching source dense rows. */ matrix *row_gather_pd_alloc(const permuted_dense *A, const int *map, int m_out); /* Fill values of C = A[map, :]; C must come from row_gather_pd_alloc(A, ...). */ diff --git a/src/utils/permuted_dense.c b/src/utils/permuted_dense.c index f2fbdd4e..f5ad9ead 100644 --- a/src/utils/permuted_dense.c +++ b/src/utils/permuted_dense.c @@ -120,7 +120,7 @@ matrix *row_gather_pd_alloc(const permuted_dense *A, const int *map, int m_out) int new_m0 = 0; for (int i = 0; i < m_out; i++) { - int ii = map[i] < 0 ? -1 : A->row_inv[map[i]]; + int ii = A->row_inv[map[i]]; if (ii >= 0) { new_row_perm[new_m0] = i; diff --git a/src/utils/sparse_matrix.c b/src/utils/sparse_matrix.c index ff4f956f..762107b5 100644 --- a/src/utils/sparse_matrix.c +++ b/src/utils/sparse_matrix.c @@ -136,13 +136,12 @@ static matrix *sparse_row_gather_alloc(const matrix *self, const int *map, int m /* Exact output nnz: sum the selected rows' nnz. Jx->nnz is NOT an upper bound — repeated map entries select the same source row more than once (cvxpy#3442). Repeated gathers of dense rows can push the true count - past INT_MAX, which a CSR cannot represent, so fail before wrapping. - map[i] == -1 selects nothing and must be branched before touching p. */ + past INT_MAX, which a CSR cannot represent, so fail before wrapping. */ int nnz = 0; for (int i = 0; i < m_out; i++) { int row = map[i]; - int len = row < 0 ? 0 : Jx->p[row + 1] - Jx->p[row]; + int len = Jx->p[row + 1] - Jx->p[row]; if (len > INT_MAX - nnz) { fprintf(stderr, "Error in sparse_row_gather_alloc: gathered nnz " @@ -157,11 +156,8 @@ static matrix *sparse_row_gather_alloc(const matrix *self, const int *map, int m for (int i = 0; i < m_out; i++) { int row = map[i]; - int len = row < 0 ? 0 : Jx->p[row + 1] - Jx->p[row]; - if (len > 0) - { - memcpy(J->i + J->p[i], Jx->i + Jx->p[row], len * sizeof(int)); - } + int len = Jx->p[row + 1] - Jx->p[row]; + memcpy(J->i + J->p[i], Jx->i + Jx->p[row], len * sizeof(int)); J->p[i + 1] = J->p[i] + len; } J->nnz = J->p[m_out]; @@ -184,10 +180,8 @@ static void sparse_row_gather_fill_values(const matrix *self, matrix *out) assert(map != NULL && self->n == out->n); for (int i = 0; i < J->m; i++) { - int row = map[i]; - if (row < 0) continue; int len = J->p[i + 1] - J->p[i]; - memcpy(J->x + J->p[i], Jx->x + Jx->p[row], len * sizeof(double)); + memcpy(J->x + J->p[i], Jx->x + Jx->p[map[i]], len * sizeof(double)); } } diff --git a/tests/all_tests.c b/tests/all_tests.c index 1d504b97..4be944e9 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -439,7 +439,6 @@ int main(void) mu_run_test(test_permuted_dense_col_inv, tests_run); mu_run_test(test_permuted_dense_row_gather, tests_run); mu_run_test(test_row_gather_sparse, tests_run); - mu_run_test(test_row_gather_sparse_all_empty, tests_run); mu_run_test(test_row_gather_pd_vs_sparse_twin, tests_run); mu_run_test(test_row_gather_spd_vs_sparse_twin, tests_run); #ifdef SP_TRACK_MEMORY diff --git a/tests/utils/test_permuted_dense.h b/tests/utils/test_permuted_dense.h index 585808fb..ffed9921 100644 --- a/tests/utils/test_permuted_dense.h +++ b/tests/utils/test_permuted_dense.h @@ -343,7 +343,7 @@ const char *test_permuted_dense_col_inv(void) /* PD row_gather_alloc / row_gather_fill_values: output must be another PD whose row_perm is the set of output positions where map[i] hits the source - row_perm, with repeats and -1 entries handled. */ + row_perm, with repeats handled. */ const char *test_permuted_dense_row_gather(void) { /* Source PD, shape (6, 4), dense block at rows {1, 3, 4} x cols {0, 2}. */ @@ -352,25 +352,24 @@ const char *test_permuted_dense_row_gather(void) double X[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; matrix *M = new_permuted_dense(6, 4, 3, 2, row_perm, col_perm, X); - /* map = [0, 3, 1, -1, 4, 3, 5]: + /* map = [0, 3, 1, 4, 3, 5]: - position 0 -> source row 0 (not in row_perm, zero) - position 1 -> source row 3 (ii=1, dense) - position 2 -> source row 1 (ii=0, dense) - - position 3 -> -1 (structurally empty) - - position 4 -> source row 4 (ii=2, dense) - - position 5 -> source row 3 again (ii=1, dense) - - position 6 -> source row 5 (not in row_perm, zero) */ - int map[7] = {0, 3, 1, -1, 4, 3, 5}; - matrix *out = M->row_gather_alloc(M, map, 7); + - position 3 -> source row 4 (ii=2, dense) + - position 4 -> source row 3 again (ii=1, dense) + - position 5 -> source row 5 (not in row_perm, zero) */ + int map[6] = {0, 3, 1, 4, 3, 5}; + matrix *out = M->row_gather_alloc(M, map, 6); permuted_dense *out_pd = (permuted_dense *) out; - mu_assert("out m", out->m == 7); + mu_assert("out m", out->m == 6); mu_assert("out n", out->n == 4); mu_assert("out nnz", out->nnz == 8); /* m0=4 * n0=2 */ mu_assert("m0", out_pd->m0 == 4); mu_assert("n0", out_pd->n0 == 2); - int expected_row_perm[4] = {1, 2, 4, 5}; + int expected_row_perm[4] = {1, 2, 3, 4}; mu_assert("row_perm", cmp_int_array(out_pd->row_perm, expected_row_perm, 4)); int expected_col_perm[2] = {0, 2}; mu_assert("col_perm", cmp_int_array(out_pd->col_perm, expected_col_perm, 2)); diff --git a/tests/utils/test_row_gather.h b/tests/utils/test_row_gather.h index eceec09b..9c74d73e 100644 --- a/tests/utils/test_row_gather.h +++ b/tests/utils/test_row_gather.h @@ -13,9 +13,8 @@ /* row_gather_alloc / row_gather_fill_values across the three matrix kinds. The map contract under test: map[i] in [0, A->m) copies source row map[i], - map[i] == -1 leaves output row i structurally empty, repeats are allowed, - and the map is bound to the result at alloc time (the caller's copy is - dead afterwards). */ + repeats are allowed, and the map is bound to the result at alloc time (the + caller's copy is dead afterwards). */ /* Shared 4x5 CSR source with an empty row 2: row 0: (0: 1.0) (3: 2.0) @@ -82,21 +81,21 @@ static int row_gather_same_csr(matrix *X, matrix *Y) cmp_int_array(a->i, b->i, a->nnz) && cmp_double_array(a->x, b->x, a->nnz); } -/* map = [3, -1, 0, 3, 2, 1]: a permutation, a repeat (row 3 twice), a -1, and - the empty source row 2. Output is 6x5 with nnz 3+0+2+3+0+2 = 10. */ +/* map = [3, 0, 3, 2, 1]: a permutation, a repeat (row 3 twice), and the empty + source row 2. Output is 5x5 with nnz 3+2+3+0+2 = 10. */ const char *test_row_gather_sparse(void) { matrix *A = row_gather_sparse_fixture(); - int map[6] = {3, -1, 0, 3, 2, 1}; - matrix *C = A->row_gather_alloc(A, map, 6); + int map[5] = {3, 0, 3, 2, 1}; + matrix *C = A->row_gather_alloc(A, map, 5); - int exp_p[7] = {0, 3, 3, 5, 8, 8, 10}; + int exp_p[6] = {0, 3, 5, 8, 8, 10}; int exp_i[10] = {0, 2, 4, 0, 3, 0, 2, 4, 1, 4}; - mu_assert("shape", C->m == 6 && C->n == 5); - mu_assert("sparsity", cmp_sparsity(C, exp_p, exp_i, 6, 10)); + mu_assert("shape", C->m == 5 && C->n == 5); + mu_assert("sparsity", cmp_sparsity(C, exp_p, exp_i, 5, 10)); /* the map is bound at alloc time: clobbering the caller's copy is fine */ - for (int k = 0; k < 6; k++) map[k] = -7; + for (int k = 0; k < 5; k++) map[k] = -7; row_gather_poison(C); A->row_gather_fill_values(A, C); @@ -114,22 +113,9 @@ const char *test_row_gather_sparse(void) return 0; } -const char *test_row_gather_sparse_all_empty(void) -{ - matrix *A = row_gather_sparse_fixture(); - int map[3] = {-1, -1, -1}; - matrix *C = A->row_gather_alloc(A, map, 3); - int exp_p[4] = {0, 0, 0, 0}; - mu_assert("empty structure", cmp_sparsity(C, exp_p, NULL, 3, 0)); - A->row_gather_fill_values(A, C); /* nothing to write; must not touch A */ - free_matrix(C); - free_matrix(A); - return 0; -} - -/* pd 6x5 with rows {1, 3, 4} x cols {0, 2}. The map hits a row outside - row_perm (0 and 5), repeats row 3, and contains a -1. Result must stay a pd - and agree entrywise with the same gather on a sparse twin. */ +/* pd 6x5 with rows {1, 3, 4} x cols {0, 2}. The map hits rows outside + row_perm (0 and 5) and repeats row 3. Result must stay a pd and agree + entrywise with the same gather on a sparse twin. */ const char *test_row_gather_pd_vs_sparse_twin(void) { int row_perm[3] = {1, 3, 4}; @@ -138,9 +124,9 @@ const char *test_row_gather_pd_vs_sparse_twin(void) matrix *A = new_permuted_dense(6, 5, 3, 2, row_perm, col_perm, X); matrix *A_tw = row_gather_sparse_twin(A); - int map[7] = {0, 3, -1, 1, 3, 5, 4}; - matrix *C = A->row_gather_alloc(A, map, 7); - matrix *C_tw = A_tw->row_gather_alloc(A_tw, map, 7); + int map[6] = {0, 3, 1, 3, 5, 4}; + matrix *C = A->row_gather_alloc(A, map, 6); + matrix *C_tw = A_tw->row_gather_alloc(A_tw, map, 6); mu_assert("kind preserved", C->is_permuted_dense); mu_assert("m0", ((permuted_dense *) C)->m0 == 4); mu_assert("nnz", C->nnz == 8); @@ -168,29 +154,28 @@ const char *test_row_gather_pd_vs_sparse_twin(void) return 0; } -/* spd fixture above; map = [4, 2, -1, 4, 0, 3, 1, 2]: - pos 0, 3 -> row 4 (block 0, twice) - pos 4 -> row 0 (block 0) - pos 1, 7 -> row 2 (block 1, twice) - pos 5 -> row 3 (block 1) - pos 2 -> -1 - pos 6 -> row 1 (no block) +/* spd fixture above; map = [4, 2, 4, 0, 3, 1, 2]: + pos 0, 2 -> row 4 (block 0, twice) + pos 3 -> row 0 (block 0) + pos 1, 6 -> row 2 (block 1, twice) + pos 4 -> row 3 (block 1) + pos 5 -> row 1 (no block) Block 2 gets no hits and must be dropped; output blocks stay row-disjoint. */ const char *test_row_gather_spd_vs_sparse_twin(void) { matrix *A = row_gather_spd_fixture(); matrix *A_tw = row_gather_sparse_twin(A); - int map[8] = {4, 2, -1, 4, 0, 3, 1, 2}; - matrix *C = A->row_gather_alloc(A, map, 8); - matrix *C_tw = A_tw->row_gather_alloc(A_tw, map, 8); + int map[7] = {4, 2, 4, 0, 3, 1, 2}; + matrix *C = A->row_gather_alloc(A, map, 7); + matrix *C_tw = A_tw->row_gather_alloc(A_tw, map, 7); mu_assert("kind preserved", C->is_stacked_pd); stacked_pd *C_spd = (stacked_pd *) C; mu_assert("empty block dropped", C_spd->n_blocks == 2); mu_assert("src_block_idx", C_spd->src_block_idx[0] == 0 && C_spd->src_block_idx[1] == 1); - int exp_rp0[3] = {0, 3, 4}; - int exp_rp1[3] = {1, 5, 7}; + int exp_rp0[3] = {0, 2, 3}; + int exp_rp1[3] = {1, 4, 6}; mu_assert("block 0 row_perm", C_spd->blocks[0]->m0 == 3 && cmp_int_array(C_spd->blocks[0]->row_perm, exp_rp0, 3)); @@ -238,8 +223,8 @@ static void run_row_gather_fill(const void *ctx) const char *test_row_gather_spd_fill_no_transient_alloc(void) { matrix *A = row_gather_spd_fixture(); - int map[8] = {4, 2, -1, 4, 0, 3, 1, 2}; - matrix *C = A->row_gather_alloc(A, map, 8); + int map[7] = {4, 2, 4, 0, 3, 1, 2}; + matrix *C = A->row_gather_alloc(A, map, 7); row_gather_fill_args args = {A, C}; run_row_gather_fill(&args); /* warm-up */ mu_assert("spd row_gather fill must not allocate", diff --git a/tests/utils/test_stacked_pd.h b/tests/utils/test_stacked_pd.h index 5429e49d..495b62f6 100644 --- a/tests/utils/test_stacked_pd.h +++ b/tests/utils/test_stacked_pd.h @@ -1323,7 +1323,7 @@ const char *test_spd_vtable_refresh_csc_values_noop(void) /* row_gather on spd: each output row is routed to the source block (if any) that carries its source row; per-block row_gather_pd_* is reused; empty - result blocks are dropped; repeats and -1 entries are allowed. */ + result blocks are dropped; repeats are allowed. */ const char *test_spd_vtable_row_gather(void) { /* 6x4 spd: @@ -1342,20 +1342,20 @@ const char *test_spd_vtable_row_gather(void) permuted_dense *blocks[2] = {(permuted_dense *) blk0, (permuted_dense *) blk1}; matrix *M = new_stacked_pd(6, 4, 2, blocks, NULL, NULL); - /* map = [3, 0, 5, 1, -1, 3]: row 3 -> blk1 (twice, positions 0 and 5), - rows 0/1 -> blk0 (positions 1 and 3), row 5 -> none, position 4 = -1. - Expected output (6 rows, positions 2 and 4 empty): + /* map = [3, 0, 5, 1, 3]: row 3 -> blk1 (twice, positions 0 and 4), + rows 0/1 -> blk0 (positions 1 and 3), row 5 -> none. + Expected output (5 rows, position 2 empty): output row 0 = src row 3 = blk1 row [7, 8] in cols {1,3} output row 1 = src row 0 = blk0 row [1, 2] in cols {0,2} output row 3 = src row 1 = blk0 row [3, 4] in cols {0,2} - output row 5 = src row 3 = blk1 row [7, 8] in cols {1,3} */ - int map[6] = {3, 0, 5, 1, -1, 3}; - matrix *C_m = M->row_gather_alloc(M, map, 6); + output row 4 = src row 3 = blk1 row [7, 8] in cols {1,3} */ + int map[5] = {3, 0, 5, 1, 3}; + matrix *C_m = M->row_gather_alloc(M, map, 5); M->row_gather_fill_values(M, C_m); stacked_pd *C = (stacked_pd *) C_m; mu_assert("n_blocks", C->n_blocks == 2); - mu_assert("base.m", C_m->m == 6); + mu_assert("base.m", C_m->m == 5); mu_assert("base.n", C_m->n == 4); /* blk0 is hit at positions 1 and 3: output block 0 carries @@ -1372,10 +1372,10 @@ const char *test_spd_vtable_row_gather(void) cmp_int_array(out0->col_perm, expected_col_perm_0, 2)); mu_assert("out0 X", cmp_double_array(out0->X, expected_X0, 4)); - /* blk1 is hit at positions 0 and 5, both from its row 1: output block 1 - carries row_perm={0,5}, col_perm={1,3}, X=[7,8,7,8]. */ + /* blk1 is hit at positions 0 and 4, both from its row 1: output block 1 + carries row_perm={0,4}, col_perm={1,3}, X=[7,8,7,8]. */ permuted_dense *out1 = C->blocks[1]; - int expected_row_perm_1[2] = {0, 5}; + int expected_row_perm_1[2] = {0, 4}; int expected_col_perm_1[2] = {1, 3}; double expected_X1[4] = {7.0, 8.0, 7.0, 8.0}; mu_assert("out1 m0", out1->m0 == 2); From ffec387540ae5e06ab647996b94d74b8b7d0cc37 Mon Sep 17 00:00:00 2001 From: dance858 Date: Sat, 12 Sep 2026 13:22:13 -0700 Subject: [PATCH 3/4] row_gather: terse slot comments, matching the rest of matrix.h Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01Azb8o1AgiFRU4t9aAY7oDr --- include/utils/matrix.h | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/include/utils/matrix.h b/include/utils/matrix.h index 6f141b12..5b3c01fe 100644 --- a/include/utils/matrix.h +++ b/include/utils/matrix.h @@ -82,16 +82,11 @@ typedef CSR_matrix *(*matrix_to_csr_fn)(matrix *A); cache already matches values_version, so it is cheap to call when fresh. */ typedef void (*matrix_refresh_csc_values_fn)(matrix *A); -/* Row gather: allocate C of shape (m_out, A->n) with C[i, :] = A[map[i], :] - for map[i] in [0, A->m). Repeated entries are allowed, so C->nnz may exceed - A->nnz. Whatever the fill needs from map is copied onto C here; the caller - may free map afterwards. */ +/* Allocate C = A[map, :] */ typedef matrix *(*matrix_row_gather_alloc_fn)(const matrix *A, const int *map, int m_out); -/* Fill values of C = A[map, :]. C must be the matrix returned by - row_gather_alloc(A, map, m_out): the gather state is bound to it at alloc - time, and copy_sparsity copies carry none. */ +/* Fill values of C = A[map, :] */ typedef void (*matrix_row_gather_fill_values_fn)(const matrix *A, matrix *C); /* Row-tiling for the promote atom: A must be a 1-row matrix; returns From ebe1456f6dc77c3b292f618f4bfc867f775c8686 Mon Sep 17 00:00:00 2001 From: dance858 Date: Sat, 12 Sep 2026 13:33:13 -0700 Subject: [PATCH 4/4] row_gather: trim header comments; note that C stores the map Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01Azb8o1AgiFRU4t9aAY7oDr --- include/utils/matrix.h | 2 +- include/utils/permuted_dense.h | 8 +++----- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/include/utils/matrix.h b/include/utils/matrix.h index 5b3c01fe..0fc72870 100644 --- a/include/utils/matrix.h +++ b/include/utils/matrix.h @@ -82,7 +82,7 @@ typedef CSR_matrix *(*matrix_to_csr_fn)(matrix *A); cache already matches values_version, so it is cheap to call when fresh. */ typedef void (*matrix_refresh_csc_values_fn)(matrix *A); -/* Allocate C = A[map, :] */ +/* Allocate C = A[map, :]. C stores map internally, so the fill takes none. */ typedef matrix *(*matrix_row_gather_alloc_fn)(const matrix *A, const int *map, int m_out); diff --git a/include/utils/permuted_dense.h b/include/utils/permuted_dense.h index d94f8616..bb01b880 100644 --- a/include/utils/permuted_dense.h +++ b/include/utils/permuted_dense.h @@ -102,13 +102,11 @@ matrix *broadcast_pd_alloc(const permuted_dense *A, broadcast_type type, int d1, void broadcast_pd_fill_values(const permuted_dense *A, broadcast_type type, int d1, int d2, permuted_dense *C); -/* Allocate C = A[map, :], where A and C are permuted dense. See - matrix_row_gather_alloc_fn for the map contract (repeats allowed). C's - row_perm is the set of output positions whose map entry hits A->row_perm; - C->bound_iwork holds the matching source dense rows. */ +/* Allocate C = A[map, :], where A and C are permuted dense. C stores map + internally, so the fill takes none. */ matrix *row_gather_pd_alloc(const permuted_dense *A, const int *map, int m_out); -/* Fill values of C = A[map, :]; C must come from row_gather_pd_alloc(A, ...). */ +/* Fill values of C = A[map, :], where A and C are permuted dense. */ void row_gather_pd_fill_values(const permuted_dense *A, permuted_dense *C); /* Allocate C = promote(A, size), where A and C are permuted dense. */