diff --git a/include/old-code/old_CSR_sum.h b/include/old-code/old_CSR_sum.h
index 3582f315..92592670 100644
--- a/include/old-code/old_CSR_sum.h
+++ b/include/old-code/old_CSR_sum.h
@@ -58,4 +58,23 @@ void sum_all_rows_csr_fill_values(const CSR_matrix *A, CSR_matrix *C,
 void sum_block_of_rows_csr_fill_values(const CSR_matrix *A, CSR_matrix *C,
                                        const int *idx_map);
 
+/* Row-sum kernels that also produce an idx_map (input nnz -> position in C->x),
+   so values can be filled by zeroing C->x and calling accumulator() from
+   utils/CSR_sum.h. C must be pre-allocated with capacity >= A->nnz; iwork must
+   have size max(A->n, A->nnz); idx_map must have size A->nnz. C->nnz is set. */
+
+/* All rows of A into the single row of C (C->m == 1). */
+void sum_all_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C, int *iwork,
+                            int *idx_map);
+
+/* Consecutive blocks of row_block_size rows of A into one row of C each
+   (C->m == A->m / row_block_size). */
+void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
+                                 int row_block_size, int *iwork, int *idx_map);
+
+/* Rows of A at stride row_spacing into one row of C each (C->m == row_spacing):
+   C[j, :] = sum_{i : i % row_spacing == j} A[i, :]. */
+void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
+                                      int row_spacing, int *iwork, int *idx_map);
+
 #endif /* OLD_CSR_SUM_H */
diff --git a/include/old-code/old_mini_numpy.h b/include/old-code/old_mini_numpy.h
new file mode 100644
index 00000000..f8d7bbbb
--- /dev/null
+++ b/include/old-code/old_mini_numpy.h
@@ -0,0 +1,26 @@
+/*
+ * Copyright 2026 Daniel Cederberg and William Zhang
+ *
+ * This file is part of the SparseDiffEngine 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.
+ */
+#ifndef OLD_MINI_NUMPY_H
+#define OLD_MINI_NUMPY_H
+
+/* Example: a = [1, 2], len = 2, tiles = 3, result = [1, 2, 1, 2, 1, 2].
+   Retired from utils/mini_numpy.h with the CSR broadcast kernel, its last
+   caller. */
+void tile_int(int *result, const int *a, int len, int tiles);
+
+#endif /* OLD_MINI_NUMPY_H */
diff --git a/include/subexpr.h b/include/subexpr.h
index b7ffd078..db639c97 100644
--- a/include/subexpr.h
+++ b/include/subexpr.h
@@ -76,7 +76,6 @@ typedef struct sum_expr
 {
     expr base;
     int axis;
-    int *idx_map; /* maps child nnz to summed-row positions */
 } sum_expr;
 
 /* trace */
diff --git a/include/utils/CSR_sum.h b/include/utils/CSR_sum.h
index 62f49d9e..d1708293 100644
--- a/include/utils/CSR_sum.h
+++ b/include/utils/CSR_sum.h
@@ -34,7 +34,7 @@ void sum_scaled_csr_matrices_fill_values(const CSR_matrix *A, const CSR_matrix *
                                          CSR_matrix *C, const double *d1,
                                          const double *d2);
 
-/* The following five functions are used for summing either more than two CSR_matrix
+/* The following two functions are used for summing either more than two CSR_matrix
    matrices or rows of CSR_matrix matrices. To implement the filling of values
    efficiently, we compute an idx_map when we fill the sparsity pattern of the output
    matrix, which maps each nonzero entry in the input matrix to its position in the
@@ -46,17 +46,8 @@ void sum_scaled_csr_matrices_fill_values(const CSR_matrix *A, const CSR_matrix *
    matrix) corresponding to the j-th nonzero in the input matrix.
 
    Output matrix C, input matrix A, iwork->size = max(A->n, A->nnz) for the first
-   four functions. The last function allocates the output matrix and returns it. */
+   function. The last function allocates the output matrix and returns it. */
 // ------------------------------------------------------------------------------------
-void sum_all_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C, int *iwork,
-                            int *idx_map);
-
-void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
-                                 int row_block_size, int *iwork, int *idx_map);
-
-void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
-                                      int row_spacing, int *iwork, int *idx_map);
-
 void sum_spaced_rows_into_row_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
                                         int spacing, int *iwork, int *idx_map);
 
diff --git a/include/utils/matrix.h b/include/utils/matrix.h
index 33d6bc20..2dfe06fc 100644
--- a/include/utils/matrix.h
+++ b/include/utils/matrix.h
@@ -87,18 +87,14 @@ typedef void (*matrix_row_gather_fill_values_fn)(const matrix *A, matrix *C);
 typedef matrix *(*matrix_diag_vec_alloc_fn)(matrix *A);
 typedef void (*matrix_diag_vec_fill_values_fn)(matrix *A, matrix *out);
 
-/* Allocate C as a row-wise reduction of A. The reduction pattern is chosen by
-   axis:
-     - axis = -1: sum all rows of A. C has shape (1, A->n).
-     - axis = 0:  block-sum rows in consecutive groups of d1. C has shape (A->m
-                  / d1, A->n). C[j, :] = sum_{i in [j*d1, (j+1)*d1)} A[i, :].
-     - axis = 1:  stride-sum rows at spacing d1. C has shape (d1, A->n). C[j, :]
-                  = sum_{i : i % d1 == j} A[i, :].
-
-   Caller pre-allocates idx_map of size A->nnz that can be used to compute the
-   numerical result of the operation using via accumulation. */
-typedef matrix *(*matrix_sum_row_partition_alloc_fn)(matrix *A, int axis, int d1,
-                                                     int *idx_map);
+/* Allocate C = row-reduce of A: C[j, :] = sum of rows i with group[i] == j,
+   group[i] in [0, m_out). C stores the reduction map internally, so the fill
+   takes none. */
+typedef matrix *(*matrix_row_reduce_alloc_fn)(const matrix *A, const int *group,
+                                              int m_out);
+
+/* Fill values of C = row-reduce of A. */
+typedef void (*matrix_row_reduce_fill_values_fn)(const matrix *A, matrix *C);
 
 typedef void (*matrix_free_fn)(matrix *self);
 
@@ -139,7 +135,8 @@ struct matrix
     matrix_row_gather_fill_values_fn row_gather_fill_values;
     matrix_diag_vec_alloc_fn diag_vec_alloc;
     matrix_diag_vec_fill_values_fn diag_vec_fill_values;
-    matrix_sum_row_partition_alloc_fn sum_row_partition_alloc;
+    matrix_row_reduce_alloc_fn row_reduce_alloc;
+    matrix_row_reduce_fill_values_fn row_reduce_fill_values;
 
     /* Lifecycle */
     matrix_free_fn free_fn;
diff --git a/include/utils/permuted_dense.h b/include/utils/permuted_dense.h
index a8f73215..13f10c54 100644
--- a/include/utils/permuted_dense.h
+++ b/include/utils/permuted_dense.h
@@ -68,8 +68,10 @@ typedef struct permuted_dense
 
     /* 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). */
+       dense row of the source it copies (length m0); row_reduce_alloc stores,
+       per dense row of the source, the dense row of this PD it adds into
+       (length source 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
@@ -101,6 +103,14 @@ matrix *row_gather_pd_alloc(const permuted_dense *A, const int *map, int m_out);
 /* 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 = row-reduce of A (C[j, :] = sum of rows i with group[i] == j,
+   group[i] in [0, m_out)), where A and C are permuted dense. C stores the
+   reduction map internally, so the fill takes none. */
+matrix *row_reduce_pd_alloc(const permuted_dense *A, const int *group, int m_out);
+
+/* Fill values of C = row-reduce of A, where A and C are permuted dense. */
+void row_reduce_pd_fill_values(const permuted_dense *A, permuted_dense *C);
+
 /* Allocate C = diag_vec(A), where A and C are permuted dense. */
 matrix *diag_vec_pd_alloc(const permuted_dense *A);
 
diff --git a/include/utils/sparse_matrix.h b/include/utils/sparse_matrix.h
index c281a4ce..039ca839 100644
--- a/include/utils/sparse_matrix.h
+++ b/include/utils/sparse_matrix.h
@@ -33,8 +33,9 @@ typedef struct sparse_matrix
     /* 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. */
+       for row_gather_fill_values; row_reduce_alloc stores the source-nnz-long
+       position map for row_reduce_fill_values. NULL otherwise; never touched
+       by any other kernel. */
     int *bound_iwork;
 } sparse_matrix;
 
diff --git a/include/utils/stacked_pd_linalg.h b/include/utils/stacked_pd_linalg.h
index 2d50b274..add55038 100644
--- a/include/utils/stacked_pd_linalg.h
+++ b/include/utils/stacked_pd_linalg.h
@@ -32,6 +32,14 @@ matrix *transpose_spd_alloc(const stacked_pd *A);
 /* Fill values of C = transpose(A). */
 void transpose_spd_fill_values(const stacked_pd *A, stacked_pd *C);
 
+/* Allocate C = row-reduce of A (C[j, :] = sum of rows i with group[i] == j,
+   group[i] in [0, m_out)); C is a stacked_pd. C stores the reduction map
+   internally, so the fill takes none. */
+matrix *row_reduce_spd_alloc(const stacked_pd *A, const int *group, int m_out);
+
+/* Fill values of C = row-reduce of A. */
+void row_reduce_spd_fill_values(const stacked_pd *A, stacked_pd *C);
+
 /* Fill values of C = diag(d) @ A, where 'd' is 'global' with length A->m. */
 void DA_spd_fill_values(const double *d, const stacked_pd *A, stacked_pd *C);
 
diff --git a/src/atoms/affine/sum.c b/src/atoms/affine/sum.c
index 42aa2ecd..420a34cd 100644
--- a/src/atoms/affine/sum.c
+++ b/src/atoms/affine/sum.c
@@ -17,15 +17,9 @@
  */
 #include "atoms/affine.h"
 #include "subexpr.h"
-#include "utils/CSR_sum.h"
-#include "utils/int_double_pair.h"
 #include "utils/mini_numpy.h"
-#include "utils/sparse_matrix.h"
-#include "utils/stacked_pd.h"
 #include "utils/tracked_alloc.h"
-#include "utils/utils.h"
 #include <assert.h>
-#include <stdlib.h>
 #include <string.h>
 
 static void forward(expr *node, const double *u)
@@ -86,26 +80,32 @@ static void jacobian_init_impl(expr *node)
     expr *x = node->left;
     sum_expr *snode = (sum_expr *) node;
     jacobian_init(x);
-
-    /* sum_row_partition_alloc fills idx_map so eval_jacobian can accumulate from
-       child->jacobian->x. */
-    snode->idx_map = sp_malloc(x->jacobian->nnz * sizeof(int));
-    node->jacobian = x->jacobian->sum_row_partition_alloc(x->jacobian, snode->axis,
-                                                          x->d1, snode->idx_map);
+    assert(x->jacobian->m == x->size);
+
+    /* Child rows are column-major: row i = r + c * x->d1. Row i of the child
+       Jacobian is summed into output row 0 (axis -1), c (axis 0) or r (axis 1).
+       The reduction map is bound to node->jacobian, so group is not kept. */
+    int d1 = x->d1;
+    int m_out = 1;
+    if (snode->axis == 0) m_out = x->d2;
+    if (snode->axis == 1) m_out = d1;
+    int *group = (int *) sp_malloc(x->size * sizeof(int));
+    for (int i = 0; i < x->size; i++)
+    {
+        int g = 0;
+        if (snode->axis == 0) g = i / d1;
+        if (snode->axis == 1) g = i % d1;
+        group[i] = g;
+    }
+    node->jacobian = x->jacobian->row_reduce_alloc(x->jacobian, group, m_out);
+    sp_free(group);
 }
 
 static void eval_jacobian_impl(expr *node)
 {
     expr *child = node->left;
-
-    /* evaluate child's jacobian */
     eval_jacobian(child);
-
-    /* we have precomputed an idx map between the nonzeros of the child's jacobian
-       and this node's jacobian, so we just accumulate accordingly */
-    memset(node->jacobian->x, 0, node->jacobian->nnz * sizeof(double));
-    accumulator(child->jacobian->x, child->jacobian->nnz,
-                ((sum_expr *) node)->idx_map, node->jacobian->x);
+    child->jacobian->row_reduce_fill_values(child->jacobian, node->jacobian);
 }
 
 static void wsum_hess_init_impl(expr *node)
@@ -149,12 +149,6 @@ static bool is_affine(const expr *node)
     return node->left->is_affine(node->left);
 }
 
-static void free_type_data(expr *node)
-{
-    sum_expr *snode = (sum_expr *) node;
-    sp_free(snode->idx_map);
-}
-
 expr *new_sum(expr *child, int axis)
 {
     int d2 = 0;
@@ -183,7 +177,7 @@ expr *new_sum(expr *child, int axis)
        sum with an axis argument as a row vector */
     init_expr(node, 1, d2, child->n_vars, forward, jacobian_init_impl,
               eval_jacobian_impl, is_affine, wsum_hess_init_impl,
-              eval_wsum_hess_impl, free_type_data);
+              eval_wsum_hess_impl, NULL);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/old-code/old_CSR_sum.c b/src/old-code/old_CSR_sum.c
index c6a659d4..b1a00cde 100644
--- a/src/old-code/old_CSR_sum.c
+++ b/src/old-code/old_CSR_sum.c
@@ -18,6 +18,7 @@
 #include "old-code/old_CSR_sum.h"
 #include "utils/CSR_matrix.h"
 #include "utils/int_double_pair.h"
+#include "utils/utils.h"
 #include <assert.h>
 #include <stdlib.h>
 #include <string.h>
@@ -330,3 +331,198 @@ void sum_spaced_rows_into_row_csr(const CSR_matrix *A, CSR_matrix *C,
 
     C->p[1] = C->nnz;
 }
+
+// ------------------------------------------------------------------------------------
+// Row-sum kernels with an idx_map (input nnz -> output position), retired from the
+// engine when the sum atom moved onto the generic row_reduce primitive. Kept as
+// standalone CSR operations. Fill values with accumulator() from utils/CSR_sum.h
+// after zeroing C->x.
+// ------------------------------------------------------------------------------------
+
+/* iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz */
+void sum_all_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C, int *iwork,
+                            int *idx_map)
+{
+    // -------------------------------------------------------------------
+    //           Build sparsity pattern of the summed row
+    // -------------------------------------------------------------------
+    int *cols = iwork;
+    memcpy(cols, A->i, A->nnz * sizeof(int));
+    sort_int_array(cols, A->nnz);
+
+    int unique_nnz = 0;
+    int prev_col = -1;
+    for (int j = 0; j < A->nnz; j++)
+    {
+        if (cols[j] != prev_col)
+        {
+            C->i[unique_nnz] = cols[j];
+            prev_col = cols[j];
+            unique_nnz++;
+        }
+    }
+
+    C->p[0] = 0;
+    C->p[1] = unique_nnz;
+    C->nnz = unique_nnz;
+
+    // -------------------------------------------------------------------
+    //  Map child values to summed-row positions. col_to_pos maps
+    //  column indices to positions in C's row.
+    // -------------------------------------------------------------------
+    int *col_to_pos = iwork;
+    for (int idx = 0; idx < unique_nnz; idx++)
+    {
+        col_to_pos[C->i[idx]] = idx;
+    }
+
+    for (int i = 0; i < A->m; i++)
+    {
+        for (int j = A->p[i]; j < A->p[i + 1]; j++)
+        {
+            idx_map[j] = col_to_pos[A->i[j]];
+        }
+    }
+}
+
+/* iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz */
+void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
+                                 int row_block_size, int *iwork, int *idx_map)
+{
+    assert(A->m % row_block_size == 0);
+    int n_blocks = A->m / row_block_size;
+    assert(C->m == n_blocks);
+
+    C->n = A->n;
+    C->p[0] = 0;
+    int cursor = 0;
+
+    int *cols = iwork;
+    int *col_to_pos = iwork;
+
+    for (int block = 0; block < n_blocks; block++)
+    {
+        int start_row = block * row_block_size;
+        int end_row = start_row + row_block_size;
+
+        // -----------------------------------------------------------------
+        // Build sparsity pattern of the row resulting from summing
+        // the block of rows from A
+        // -----------------------------------------------------------------
+        C->p[block] = cursor;
+        int count = 0;
+        for (int row = start_row; row < end_row; row++)
+        {
+            for (int j = A->p[row]; j < A->p[row + 1]; j++)
+            {
+                cols[count++] = A->i[j];
+            }
+        }
+
+        /* Sort columns and write unique pattern into C->i */
+        sort_int_array(cols, count);
+
+        int unique_nnz = 0;
+        int prev_col = -1;
+        for (int t = 0; t < count; t++)
+        {
+            int col = cols[t];
+            if (t == 0 || col != prev_col)
+            {
+                C->i[cursor + unique_nnz] = col;
+                prev_col = col;
+                unique_nnz++;
+            }
+        }
+
+        cursor += unique_nnz;
+        C->p[block + 1] = cursor;
+
+        // -----------------------------------------------------------------
+        //         Build idx_map for all entries in this block
+        // -----------------------------------------------------------------
+        int row_start = C->p[block];
+        for (int idx = 0; idx < unique_nnz; idx++)
+        {
+            col_to_pos[C->i[row_start + idx]] = row_start + idx;
+        }
+
+        for (int row = start_row; row < end_row; row++)
+        {
+            for (int j = A->p[row]; j < A->p[row + 1]; j++)
+            {
+                idx_map[j] = col_to_pos[A->i[j]];
+            }
+        }
+    }
+
+    C->nnz = cursor;
+}
+
+/* iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz */
+void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
+                                      int row_spacing, int *iwork, int *idx_map)
+{
+    assert(C->m == row_spacing);
+    C->n = A->n;
+    C->p[0] = 0;
+    int cursor = 0;
+
+    int *cols = iwork;
+    int *col_to_pos = iwork;
+
+    for (int C_row = 0; C_row < C->m; C_row++)
+    {
+        // -----------------------------------------------------------------
+        // Build sparsity pattern of the row resulting from summing
+        // evenly spaced rows from A
+        // -----------------------------------------------------------------
+        C->p[C_row] = cursor;
+        int count = 0;
+        for (int row = C_row; row < A->m; row += row_spacing)
+        {
+            for (int j = A->p[row]; j < A->p[row + 1]; j++)
+            {
+                cols[count++] = A->i[j];
+            }
+        }
+
+        /* Sort columns and write unique pattern into C->i */
+        sort_int_array(cols, count);
+
+        int unique_nnz = 0;
+        int prev_col = -1;
+        for (int t = 0; t < count; t++)
+        {
+            int col = cols[t];
+            if (t == 0 || col != prev_col)
+            {
+                C->i[cursor + unique_nnz] = col;
+                prev_col = col;
+                unique_nnz++;
+            }
+        }
+
+        cursor += unique_nnz;
+        C->p[C_row + 1] = cursor;
+
+        // -----------------------------------------------------------------
+        //         Build idx_map for all entries in evenly spaced rows
+        // -----------------------------------------------------------------
+        int row_start = C->p[C_row];
+        for (int idx = 0; idx < unique_nnz; idx++)
+        {
+            col_to_pos[C->i[row_start + idx]] = row_start + idx;
+        }
+
+        for (int row = C_row; row < A->m; row += row_spacing)
+        {
+            for (int j = A->p[row]; j < A->p[row + 1]; j++)
+            {
+                idx_map[j] = col_to_pos[A->i[j]];
+            }
+        }
+    }
+
+    C->nnz = cursor;
+}
diff --git a/src/old-code/old_mini_numpy.c b/src/old-code/old_mini_numpy.c
new file mode 100644
index 00000000..ec6a6ae5
--- /dev/null
+++ b/src/old-code/old_mini_numpy.c
@@ -0,0 +1,27 @@
+/*
+ * Copyright 2026 Daniel Cederberg and William Zhang
+ *
+ * This file is part of the SparseDiffEngine 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 "old-code/old_mini_numpy.h"
+#include <string.h>
+
+void tile_int(int *result, const int *a, int len, int tiles)
+{
+    for (int i = 0; i < tiles; i++)
+    {
+        memcpy(result + i * len, a, len * sizeof(int));
+    }
+}
diff --git a/src/utils/CSR_sum.c b/src/utils/CSR_sum.c
index e4330a98..dcdb5615 100644
--- a/src/utils/CSR_sum.c
+++ b/src/utils/CSR_sum.c
@@ -148,148 +148,6 @@ void sum_scaled_csr_matrices_fill_values(const CSR_matrix *A, const CSR_matrix *
     }
 }
 
-/* iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz */
-void sum_block_of_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
-                                 int row_block_size, int *iwork, int *idx_map)
-{
-    assert(A->m % row_block_size == 0);
-    int n_blocks = A->m / row_block_size;
-    assert(C->m == n_blocks);
-
-    C->n = A->n;
-    C->p[0] = 0;
-    int cursor = 0;
-
-    int *cols = iwork;
-    int *col_to_pos = iwork;
-
-    for (int block = 0; block < n_blocks; block++)
-    {
-        int start_row = block * row_block_size;
-        int end_row = start_row + row_block_size;
-
-        // -----------------------------------------------------------------
-        // Build sparsity pattern of the row resulting from summing
-        // the block of rows from A
-        // -----------------------------------------------------------------
-        C->p[block] = cursor;
-        int count = 0;
-        for (int row = start_row; row < end_row; row++)
-        {
-            for (int j = A->p[row]; j < A->p[row + 1]; j++)
-            {
-                cols[count++] = A->i[j];
-            }
-        }
-
-        /* Sort columns and write unique pattern into C->i */
-        sort_int_array(cols, count);
-
-        int unique_nnz = 0;
-        int prev_col = -1;
-        for (int t = 0; t < count; t++)
-        {
-            int col = cols[t];
-            if (t == 0 || col != prev_col)
-            {
-                C->i[cursor + unique_nnz] = col;
-                prev_col = col;
-                unique_nnz++;
-            }
-        }
-
-        cursor += unique_nnz;
-        C->p[block + 1] = cursor;
-
-        // -----------------------------------------------------------------
-        //         Build idx_map for all entries in this block
-        // -----------------------------------------------------------------
-        int row_start = C->p[block];
-        for (int idx = 0; idx < unique_nnz; idx++)
-        {
-            col_to_pos[C->i[row_start + idx]] = row_start + idx;
-        }
-
-        for (int row = start_row; row < end_row; row++)
-        {
-            for (int j = A->p[row]; j < A->p[row + 1]; j++)
-            {
-                idx_map[j] = col_to_pos[A->i[j]];
-            }
-        }
-    }
-
-    C->nnz = cursor;
-}
-
-/* iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz */
-void sum_evenly_spaced_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C,
-                                      int row_spacing, int *iwork, int *idx_map)
-{
-    assert(C->m == row_spacing);
-    C->n = A->n;
-    C->p[0] = 0;
-    int cursor = 0;
-
-    int *cols = iwork;
-    int *col_to_pos = iwork;
-
-    for (int C_row = 0; C_row < C->m; C_row++)
-    {
-        // -----------------------------------------------------------------
-        // Build sparsity pattern of the row resulting from summing
-        // evenly spaced rows from A
-        // -----------------------------------------------------------------
-        C->p[C_row] = cursor;
-        int count = 0;
-        for (int row = C_row; row < A->m; row += row_spacing)
-        {
-            for (int j = A->p[row]; j < A->p[row + 1]; j++)
-            {
-                cols[count++] = A->i[j];
-            }
-        }
-
-        /* Sort columns and write unique pattern into C->i */
-        sort_int_array(cols, count);
-
-        int unique_nnz = 0;
-        int prev_col = -1;
-        for (int t = 0; t < count; t++)
-        {
-            int col = cols[t];
-            if (t == 0 || col != prev_col)
-            {
-                C->i[cursor + unique_nnz] = col;
-                prev_col = col;
-                unique_nnz++;
-            }
-        }
-
-        cursor += unique_nnz;
-        C->p[C_row + 1] = cursor;
-
-        // -----------------------------------------------------------------
-        //         Build idx_map for all entries in evenly spaced rows
-        // -----------------------------------------------------------------
-        int row_start = C->p[C_row];
-        for (int idx = 0; idx < unique_nnz; idx++)
-        {
-            col_to_pos[C->i[row_start + idx]] = row_start + idx;
-        }
-
-        for (int row = C_row; row < A->m; row += row_spacing)
-        {
-            for (int j = A->p[row]; j < A->p[row + 1]; j++)
-            {
-                idx_map[j] = col_to_pos[A->i[j]];
-            }
-        }
-    }
-
-    C->nnz = cursor;
-}
-
 void accumulator(const double *vals, int nnz, const int *idx_map, double *out)
 {
     /* don't forget to initialize accumulator to 0 before calling this */
@@ -312,51 +170,6 @@ void accumulator_with_spacing(const CSR_matrix *A, const int *idx_map, double *o
     }
 }
 
-void sum_all_rows_csr_alloc(const CSR_matrix *A, CSR_matrix *C, int *iwork,
-                            int *idx_map)
-{
-    // -------------------------------------------------------------------
-    //           Build sparsity pattern of the summed row
-    // -------------------------------------------------------------------
-    int *cols = iwork;
-    memcpy(cols, A->i, A->nnz * sizeof(int));
-    sort_int_array(cols, A->nnz);
-
-    int unique_nnz = 0;
-    int prev_col = -1;
-    for (int j = 0; j < A->nnz; j++)
-    {
-        if (cols[j] != prev_col)
-        {
-            C->i[unique_nnz] = cols[j];
-            prev_col = cols[j];
-            unique_nnz++;
-        }
-    }
-
-    C->p[0] = 0;
-    C->p[1] = unique_nnz;
-    C->nnz = unique_nnz;
-
-    // -------------------------------------------------------------------
-    //  Map child values to summed-row positions. col_to_pos maps
-    //  column indices to positions in C's row.
-    // -------------------------------------------------------------------
-    int *col_to_pos = iwork;
-    for (int idx = 0; idx < unique_nnz; idx++)
-    {
-        col_to_pos[C->i[idx]] = idx;
-    }
-
-    for (int i = 0; i < A->m; i++)
-    {
-        for (int j = A->p[i]; j < A->p[i + 1]; j++)
-        {
-            idx_map[j] = col_to_pos[A->i[j]];
-        }
-    }
-}
-
 /*
  * Sums evenly spaced rows from A into a single row in C and fills an index map.
  * A: input CSR_matrix matrix
diff --git a/src/utils/permuted_dense.c b/src/utils/permuted_dense.c
index 637263e3..9268eddb 100644
--- a/src/utils/permuted_dense.c
+++ b/src/utils/permuted_dense.c
@@ -253,146 +253,73 @@ static void permuted_dense_vtable_block_left_mult_values(const matrix *A,
     I_kron_A_fill_values(A, J, C, pd->kernel_dwork);
 }
 
-/* C = sum-all-rows of A. */
-static matrix *sum_all_rows_pd_alloc(matrix *A, int *idx_map)
+matrix *row_reduce_pd_alloc(const permuted_dense *A, const int *group, int m_out)
 {
-    permuted_dense *pd = (permuted_dense *) A;
-
-    /* allocate C */
-    int zero = 0;
-    matrix *C = new_permuted_dense(1, A->n, 1, pd->n0, &zero, pd->col_perm, NULL);
-
-    /* fill idx_map */
-    for (int i = 0; i < pd->m0; i++)
+    /* C's dense rows are the groups hit by A's dense rows, in increasing order.
+       bound_iwork[ii] is the dense row of C that source dense row ii adds into. */
+    bool *seen = (bool *) sp_calloc(m_out, sizeof(bool));
+    for (int ii = 0; ii < A->m0; ii++)
     {
-        int *idx_base = idx_map + i * pd->n0;
-        for (int j = 0; j < pd->n0; j++)
-        {
-            idx_base[j] = j;
-        }
+        int g = group[A->row_perm[ii]];
+        assert(g >= 0 && g < m_out);
+        seen[g] = true;
     }
-    return C;
-}
-
-/* C = block-sum of A's rows in consecutive groups of d1. */
-static matrix *sum_block_of_rows_pd_alloc(matrix *A_matrix, int d1, int *idx_map)
-{
-    permuted_dense *A = (permuted_dense *) A_matrix;
-    int C_m0 = 0;
-    int last_bucket = -1;
-    int *C_row_perm = (int *) sp_malloc(A->m0 * sizeof(int));
-
-    /* per input dense row ii, the index of its bucket within C_row_perm */
-    int *row_to_out = (int *) sp_malloc(A->m0 * sizeof(int));
-
-    // ---------------------------------------------------------------------------
-    //                          determine C's row_perm
-    // ---------------------------------------------------------------------------
-
-    /* for every row in the dense block */
-    for (int ii = 0; ii < A->m0; ii++)
+    int *C_row_perm = (int *) sp_malloc(m_out * sizeof(int));
+    int *group_to_out = (int *) sp_malloc(m_out * sizeof(int));
+    int new_m0 = 0;
+    for (int g = 0; g < m_out; g++)
     {
-        /* find the bucket to which this row belongs */
-        int bucket = A->row_perm[ii] / d1;
-
-        /* add the bucket to C if it's new */
-        if (bucket != last_bucket)
+        if (seen[g])
         {
-            C_row_perm[C_m0++] = bucket;
-            last_bucket = bucket;
+            group_to_out[g] = new_m0;
+            C_row_perm[new_m0++] = g;
         }
-
-        /* map the input row of A to its row in C */
-        row_to_out[ii] = C_m0 - 1;
     }
 
-    matrix *C = new_permuted_dense(A_matrix->m / d1, A_matrix->n, C_m0, A->n0,
-                                   C_row_perm, A->col_perm, NULL);
-
-    // ---------------------------------------------------------------------------
-    //                          fill idx_map
-    // ---------------------------------------------------------------------------
-    for (int ii = 0; ii < A->m0; ii++)
+    matrix *out = new_permuted_dense(m_out, A->base.n, new_m0, A->n0, C_row_perm,
+                                     A->col_perm, NULL);
+    if (A->m0 > 0)
     {
-        int offset = row_to_out[ii] * A->n0;
-        int *idx_base = idx_map + ii * A->n0;
-        for (int jj = 0; jj < A->n0; jj++)
+        permuted_dense *C = (permuted_dense *) out;
+        C->bound_iwork = (int *) sp_malloc(A->m0 * sizeof(int));
+        for (int ii = 0; ii < A->m0; ii++)
         {
-            idx_base[jj] = offset + jj;
+            C->bound_iwork[ii] = group_to_out[group[A->row_perm[ii]]];
         }
     }
-
-    sp_free(row_to_out);
+    sp_free(group_to_out);
     sp_free(C_row_perm);
-    return C;
+    sp_free(seen);
+    return out;
 }
 
-/* C = stride-sum of A's rows at modular spacing d1. C has shape (d1, A->n);
-   C[j, :] = sum_{i : i % d1 == j} A[i, :]. */
-static matrix *sum_evenly_spaced_rows_pd_alloc(matrix *self, int d1, int *idx_map)
+void row_reduce_pd_fill_values(const permuted_dense *A, permuted_dense *C)
 {
-    permuted_dense *A = (permuted_dense *) self;
-
-    // ---------------------------------------------------------------------------
-    //          which buckets of [0, d1) are hit by A->row_perm?
-    // ---------------------------------------------------------------------------
-    bool *seen = (bool *) sp_calloc(d1, sizeof(bool));
-    for (int ii = 0; ii < A->m0; ii++)
-    {
-        seen[A->row_perm[ii] % d1] = true;
-    }
-
-    // ---------------------------------------------------------------------------
-    //              determine C's row_perm (guarantees sorted order)
-    // ---------------------------------------------------------------------------
-    int *C_row_perm = (int *) sp_malloc(A->m0 * sizeof(int));
-    int *bucket_to_out_idx = (int *) sp_malloc(d1 * sizeof(int));
-    int C_m0 = 0;
-    for (int ii = 0; ii < d1; ii++)
-    {
-        if (seen[ii])
-        {
-            bucket_to_out_idx[ii] = C_m0;
-            C_row_perm[C_m0++] = ii;
-        }
-    }
-    sp_free(seen);
-
-    matrix *C =
-        new_permuted_dense(d1, self->n, C_m0, A->n0, C_row_perm, A->col_perm, NULL);
-
-    // ---------------------------------------------------------------------------
-    //                          fill idx_map
-    // ---------------------------------------------------------------------------
+    if (C->base.nnz == 0) return;
+    assert(A->base.n == C->base.n && C->bound_iwork != NULL);
+    int n0 = A->n0;
+    memset(C->X, 0, (size_t) C->m0 * n0 * sizeof(double));
     for (int ii = 0; ii < A->m0; ii++)
     {
-        int base = bucket_to_out_idx[A->row_perm[ii] % d1] * A->n0;
-        int *idx_base = idx_map + ii * A->n0;
-        for (int jj = 0; jj < A->n0; jj++)
+        const double *src = A->X + ii * n0;
+        double *dst = C->X + C->bound_iwork[ii] * n0;
+        for (int jj = 0; jj < n0; jj++)
         {
-            idx_base[jj] = base + jj;
+            dst[jj] += src[jj];
         }
     }
-
-    sp_free(bucket_to_out_idx);
-    sp_free(C_row_perm);
-    return C;
 }
 
-static matrix *permuted_dense_vtable_sum_row_partition_alloc(matrix *self, int axis,
-                                                             int d1, int *idx_map)
+static matrix *permuted_dense_vtable_row_reduce_alloc(const matrix *self,
+                                                      const int *group, int m_out)
 {
-    if (axis == -1)
-    {
-        return sum_all_rows_pd_alloc(self, idx_map);
-    }
-
-    if (axis == 0)
-    {
-        return sum_block_of_rows_pd_alloc(self, d1, idx_map);
-    }
+    return row_reduce_pd_alloc((const permuted_dense *) self, group, m_out);
+}
 
-    return sum_evenly_spaced_rows_pd_alloc(self, d1, idx_map); /* axis == 1 */
+static void permuted_dense_vtable_row_reduce_fill_values(const matrix *self,
+                                                         matrix *out)
+{
+    row_reduce_pd_fill_values((const permuted_dense *) self, (permuted_dense *) out);
 }
 
 static void wire_vtable(permuted_dense *pd)
@@ -414,7 +341,8 @@ static void wire_vtable(permuted_dense *pd)
     pd->base.row_gather_fill_values = permuted_dense_vtable_row_gather_fill_values;
     pd->base.diag_vec_alloc = permuted_dense_vtable_diag_vec_alloc;
     pd->base.diag_vec_fill_values = permuted_dense_vtable_diag_vec_fill_values;
-    pd->base.sum_row_partition_alloc = permuted_dense_vtable_sum_row_partition_alloc;
+    pd->base.row_reduce_alloc = permuted_dense_vtable_row_reduce_alloc;
+    pd->base.row_reduce_fill_values = permuted_dense_vtable_row_reduce_fill_values;
     pd->base.refresh_csc_values = permuted_dense_refresh_csc_values;
 }
 
diff --git a/src/utils/sparse_matrix.c b/src/utils/sparse_matrix.c
index c1421beb..0bf2355f 100644
--- a/src/utils/sparse_matrix.c
+++ b/src/utils/sparse_matrix.c
@@ -235,43 +235,107 @@ static void sparse_refresh_csc_values(matrix *self)
     sm->csc_seen = sm->base.values_version;
 }
 
-static matrix *sparse_sum_row_partition_alloc(matrix *self, int axis, int d1,
-                                              int *idx_map)
+static matrix *sparse_row_reduce_alloc(const matrix *self, const int *group,
+                                       int m_out)
 {
-    CSR_matrix *A = ((sparse_matrix *) self)->csr;
-    int m;
-    if (axis == -1)
+    const CSR_matrix *A = ((const sparse_matrix *) self)->csr;
+    int m = A->m;
+    int n = A->n;
+
+    /* Inverse of group by counting sort: rows[start[j] .. start[j+1]) lists the
+       source rows of output row j. */
+    int *start = (int *) sp_calloc(m_out + 1, sizeof(int));
+    int *fill = (int *) sp_malloc(m_out * sizeof(int));
+    int *rows = (int *) sp_malloc(m * sizeof(int));
+    for (int i = 0; i < m; i++)
     {
-        m = 1;
+        assert(group[i] >= 0 && group[i] < m_out);
+        start[group[i] + 1]++;
     }
-    else if (axis == 0)
+    for (int j = 0; j < m_out; j++)
     {
-        m = A->m / d1;
+        start[j + 1] += start[j];
+        fill[j] = start[j];
     }
-    else
+    for (int i = 0; i < m; i++)
     {
-        m = d1;
+        rows[fill[group[i]]++] = i;
     }
-    int max_nnz = MIN(A->nnz, sat_mul_int(m, A->n));
-    CSR_matrix *out = new_CSR_matrix(m, A->n, max_nnz);
-    int *iwork = (int *) sp_malloc(MAX(A->n, A->nnz) * sizeof(int));
 
-    if (axis == -1)
+    /* Output row j is the sorted union of its source rows' columns. Summing can
+       only merge entries, so A->nnz bounds the output nnz. pos_of[c] is the
+       position in J of column c within the output row being built; anything
+       below row_start is left over from an earlier row and means "not yet in
+       this row". */
+    int cap = MIN(A->nnz, sat_mul_int(m_out, n));
+    CSR_matrix *J = new_CSR_matrix(m_out, n, cap);
+    int *map = (int *) sp_malloc(A->nnz * sizeof(int));
+    int *pos_of = (int *) sp_malloc(n * sizeof(int));
+    for (int c = 0; c < n; c++)
     {
-        sum_all_rows_csr_alloc(A, out, iwork, idx_map);
+        pos_of[c] = -1;
     }
-    else if (axis == 0)
-    {
-        sum_block_of_rows_csr_alloc(A, out, d1, iwork, idx_map);
-    }
-    else
+
+    int nnz = 0;
+    J->p[0] = 0;
+    for (int j = 0; j < m_out; j++)
     {
-        sum_evenly_spaced_rows_csr_alloc(A, out, m, iwork, idx_map);
+        int row_start = nnz;
+        for (int ii = start[j]; ii < start[j + 1]; ii++)
+        {
+            int i = rows[ii];
+            for (int jj = A->p[i]; jj < A->p[i + 1]; jj++)
+            {
+                int c = A->i[jj];
+                if (pos_of[c] < row_start)
+                {
+                    pos_of[c] = nnz;
+                    J->i[nnz++] = c;
+                }
+            }
+        }
+        if (nnz > row_start)
+        {
+            sort_int_array(J->i + row_start, nnz - row_start);
+        }
+        J->p[j + 1] = nnz;
+
+        /* sorting moved the columns: record final positions, then map every
+           source entry of this output row to its output position */
+        for (int pos = row_start; pos < nnz; pos++)
+        {
+            pos_of[J->i[pos]] = pos;
+        }
+        for (int ii = start[j]; ii < start[j + 1]; ii++)
+        {
+            int i = rows[ii];
+            for (int jj = A->p[i]; jj < A->p[i + 1]; jj++)
+            {
+                map[jj] = pos_of[A->i[jj]];
+            }
+        }
     }
+    J->nnz = nnz;
+    CSR_trim(J);
+
+    sp_free(pos_of);
+    sp_free(rows);
+    sp_free(fill);
+    sp_free(start);
 
-    sp_free(iwork);
-    CSR_trim(out);
-    return new_sparse_matrix(out);
+    sparse_matrix *out = (sparse_matrix *) new_sparse_matrix(J);
+    out->bound_iwork = map;
+    return &out->base;
+}
+
+static void sparse_row_reduce_fill_values(const matrix *self, matrix *out)
+{
+    if (out->nnz == 0) return;
+    const CSR_matrix *A = ((const sparse_matrix *) self)->csr;
+    sparse_matrix *sm_out = (sparse_matrix *) out;
+    assert(sm_out->bound_iwork != NULL && self->n == out->n);
+    memset(out->x, 0, out->nnz * sizeof(double));
+    accumulator(A->x, A->nnz, sm_out->bound_iwork, out->x);
 }
 
 static void wire_vtable(sparse_matrix *sm)
@@ -290,7 +354,8 @@ static void wire_vtable(sparse_matrix *sm)
     sm->base.row_gather_fill_values = sparse_row_gather_fill_values;
     sm->base.diag_vec_alloc = sparse_diag_vec_alloc;
     sm->base.diag_vec_fill_values = sparse_diag_vec_fill_values;
-    sm->base.sum_row_partition_alloc = sparse_sum_row_partition_alloc;
+    sm->base.row_reduce_alloc = sparse_row_reduce_alloc;
+    sm->base.row_reduce_fill_values = sparse_row_reduce_fill_values;
     sm->base.refresh_csc_values = sparse_refresh_csc_values;
     sm->base.free_fn = sparse_free;
 }
diff --git a/src/utils/stacked_pd.c b/src/utils/stacked_pd.c
index 4771d814..5f40ac67 100644
--- a/src/utils/stacked_pd.c
+++ b/src/utils/stacked_pd.c
@@ -17,8 +17,6 @@
  */
 #include "utils/stacked_pd.h"
 
-#include "utils/CSR_sum.h"
-#include "utils/iVec.h"
 #include "utils/matrix.h"
 #include "utils/permuted_dense.h"
 #include "utils/sparse_matrix.h"
@@ -331,92 +329,15 @@ static void assert_disjoint_row_perms(int n_blocks, permuted_dense *const *block
 }
 #endif
 
-/* C = sum(A, axis) for a stacked_pd A.
-
-   axis == -1: native union path. C is a 1-row permuted_dense whose col_perm
-   is the sorted union of all blocks' col_perms (every input cell collapses
-   to exactly one output column). idx_map is filled in block-major order
-   matching A's base.x layout.
-
-   axis == 0 or 1: internal CSR fallback. A row-reduction of a stacked_pd
-   doesn't in general fit a single PD output (different output rows may have
-   different col footprints), so we materialize via to_csr, dispatch to the
-   existing CSR helper, then re-index idx_map from CSR ordering into A's
-   block-major base.x ordering so a downstream values-fill pass can read
-   directly from A->base.x. */
-static matrix *stacked_pd_vtable_sum_row_partition_alloc(matrix *self, int axis,
-                                                         int d1, int *idx_map)
+static matrix *stacked_pd_vtable_row_reduce_alloc(const matrix *self,
+                                                  const int *group, int m_out)
 {
-    stacked_pd *spd = (stacked_pd *) self;
-
-    if (axis == -1)
-    {
-        /* multi-way sorted union of each block's col_perm */
-        iVec *col_union = iVec_new(8);
-        const int **col_arrs =
-            (const int **) sp_malloc(spd->n_blocks * sizeof(int *));
-        int *col_lens = (int *) sp_malloc(spd->n_blocks * sizeof(int));
-        for (int k = 0; k < spd->n_blocks; k++)
-        {
-            col_arrs[k] = spd->blocks[k]->col_perm;
-            col_lens[k] = spd->blocks[k]->n0;
-        }
-        sorted_union_int_arrays(col_arrs, col_lens, spd->n_blocks, col_union);
-        sp_free(col_arrs);
-        sp_free(col_lens);
-
-        /* inverse map: column-id (in [0, self->n)) → its position in col_union */
-        int *col_to_pos = (int *) sp_malloc(self->n * sizeof(int));
-        for (int p = 0; p < col_union->len; p++)
-        {
-            col_to_pos[col_union->data[p]] = p;
-        }
-
-        int row_zero = 0;
-        matrix *out = new_permuted_dense(1, self->n, 1, col_union->len, &row_zero,
-                                         col_union->data, NULL);
-        iVec_free(col_union);
-
-        /* fill idx_map in block-major order matching spd->base.x layout */
-        int native_pos = 0;
-        for (int k = 0; k < spd->n_blocks; k++)
-        {
-            permuted_dense *blk = spd->blocks[k];
-            for (int i = 0; i < blk->m0; i++)
-            {
-                for (int j = 0; j < blk->n0; j++)
-                {
-                    idx_map[native_pos++] = col_to_pos[blk->col_perm[j]];
-                }
-            }
-        }
-        sp_free(col_to_pos);
-        return out;
-    }
-
-    /* axis == 0 or 1: CSR fallback */
-    CSR_matrix *A = self->to_csr(self);
-    int m_out = (axis == 0) ? A->m / d1 : d1;
-    int max_out_nnz = MIN(A->nnz, sat_mul_int(m_out, A->n));
-    CSR_matrix *out = new_CSR_matrix(m_out, A->n, max_out_nnz);
-    int *iwork = (int *) sp_malloc(MAX(A->n, A->nnz) * sizeof(int));
-
-    if (axis == 0)
-    {
-        sum_block_of_rows_csr_alloc(A, out, d1, iwork, idx_map);
-    }
-    else
-    {
-        sum_evenly_spaced_rows_csr_alloc(A, out, m_out, iwork, idx_map);
-    }
-    sp_free(iwork);
-
-    /* idx_map is currently in CSR row order; re-index to block-major base.x
-       order so eval_jacobian reads child->jacobian->x directly. */
-    compose_csr_idx_map_for_spd(spd, A, idx_map);
+    return row_reduce_spd_alloc((const stacked_pd *) self, group, m_out);
+}
 
-    CSR_trim(out);
-    return new_sparse_matrix(out);
+static void stacked_pd_vtable_row_reduce_fill_values(const matrix *self, matrix *out)
+{
+    row_reduce_spd_fill_values((const stacked_pd *) self, (stacked_pd *) out);
 }
 
 static void wire_vtable(stacked_pd *spd)
@@ -435,7 +356,8 @@ static void wire_vtable(stacked_pd *spd)
     spd->base.row_gather_fill_values = stacked_pd_vtable_row_gather_fill_values;
     spd->base.diag_vec_alloc = stacked_pd_vtable_diag_vec_alloc;
     spd->base.diag_vec_fill_values = stacked_pd_vtable_diag_vec_fill_values;
-    spd->base.sum_row_partition_alloc = stacked_pd_vtable_sum_row_partition_alloc;
+    spd->base.row_reduce_alloc = stacked_pd_vtable_row_reduce_alloc;
+    spd->base.row_reduce_fill_values = stacked_pd_vtable_row_reduce_fill_values;
 }
 
 matrix *new_stacked_pd_unchecked(int m, int n, int n_blocks, permuted_dense **blocks,
diff --git a/src/utils/stacked_pd_linalg.c b/src/utils/stacked_pd_linalg.c
index abdc99be..da6ecada 100644
--- a/src/utils/stacked_pd_linalg.c
+++ b/src/utils/stacked_pd_linalg.c
@@ -154,6 +154,44 @@ static void spd_blockwise_fill_coalesce_accumulate(const stacked_pd *spd_iter,
     coalesce_spd_fill_values_accumulate(raw, C);
 }
 
+// ----------------------------------------------------------------------------------
+// C = row-reduce of stacked_pd A: C[j, :] = sum of rows i with group[i] == j.
+// Each block is reduced on its own into a PD (rows of one block that share a
+// group are summed there); the partials may overlap in rows and cells across
+// blocks, so they go through the coalesce-accumulate skeleton above. The result
+// is always a stacked_pd.
+// ----------------------------------------------------------------------------------
+typedef struct
+{
+    const int *group;
+    int m_out;
+} row_reduce_ctx;
+
+static matrix *row_reduce_partial_alloc(const permuted_dense *blk, const void *ctx)
+{
+    const row_reduce_ctx *c = (const row_reduce_ctx *) ctx;
+    return row_reduce_pd_alloc(blk, c->group, c->m_out);
+}
+
+matrix *row_reduce_spd_alloc(const stacked_pd *A, const int *group, int m_out)
+{
+    row_reduce_ctx ctx = {group, m_out};
+    return spd_blockwise_alloc_coalesce(A, m_out, A->base.n,
+                                        row_reduce_partial_alloc, &ctx);
+}
+
+void row_reduce_spd_fill_values(const stacked_pd *A, stacked_pd *C)
+{
+    if (C->base.nnz == 0) return;
+    stacked_pd *raw = C->pre_coalesce;
+    for (int k = 0; k < A->n_blocks; k++)
+    {
+        row_reduce_pd_fill_values(A->blocks[k], raw->blocks[k]);
+    }
+    memset(C->base.x, 0, C->base.nnz * sizeof(double));
+    coalesce_spd_fill_values_accumulate(raw, C);
+}
+
 // ------------------------------------------------------------------------------------
 // C = ATDA for stacked_pd A. Let A = [A1; A2; A3] where Ai has n columns (the same
 // number as A). Then ATDA = A1^T D1 A1 + A2^T D2 A2 + A3^T D3 A3. Term i and j
diff --git a/tests/all_tests.c b/tests/all_tests.c
index 98eb767b..dd7853e3 100644
--- a/tests/all_tests.c
+++ b/tests/all_tests.c
@@ -75,6 +75,7 @@
 #include "utils/test_matrix.h"
 #include "utils/test_permuted_dense.h"
 #include "utils/test_row_gather.h"
+#include "utils/test_row_reduce.h"
 #include "utils/test_stacked_pd.h"
 #include "wsum_hess/affine/test_broadcast.h"
 #include "wsum_hess/affine/test_convolve.h"
@@ -232,6 +233,9 @@ int main(void)
     mu_run_test(test_jacobian_sum_add_log_axis_0, tests_run);
     mu_run_test(test_jacobian_sum_log_axis_1, tests_run);
     mu_run_test(test_jacobian_sum_axis_minus_one_pd_child, tests_run);
+    mu_run_test(test_jacobian_sum_spd_child_axis_minus_one, tests_run);
+    mu_run_test(test_jacobian_sum_spd_child_axis_0, tests_run);
+    mu_run_test(test_jacobian_sum_spd_child_axis_1, tests_run);
     mu_run_test(test_jacobian_hstack_vectors, tests_run);
     mu_run_test(test_jacobian_hstack_matrix, tests_run);
     mu_run_test(test_jacobian_vstack_vectors, tests_run);
@@ -446,6 +450,14 @@ int main(void)
     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_row_reduce_sparse, tests_run);
+    mu_run_test(test_row_reduce_pd, tests_run);
+    mu_run_test(test_row_reduce_spd_cross_block_accumulate, tests_run);
+    mu_run_test(test_row_reduce_spd_within_block, tests_run);
+    mu_run_test(test_row_reduce_spd_all_to_one, tests_run);
+#ifdef SP_TRACK_MEMORY
+    mu_run_test(test_row_reduce_spd_fill_no_transient_alloc, tests_run);
 #endif
     mu_run_test(test_permuted_dense_diag_vec, tests_run);
     mu_run_test(test_permuted_dense_BTA_matching_row_perm, tests_run);
@@ -454,9 +466,6 @@ int main(void)
     mu_run_test(test_permuted_dense_BTDA_decomposition, tests_run);
     mu_run_test(test_permuted_dense_BTDA_matching_row_perm, tests_run);
     mu_run_test(test_permuted_dense_BTDA_partial_overlap, tests_run);
-    mu_run_test(test_permuted_dense_sum_all_rows, tests_run);
-    mu_run_test(test_permuted_dense_sum_block_of_rows, tests_run);
-    mu_run_test(test_permuted_dense_sum_evenly_spaced_rows, tests_run);
     mu_run_test(test_BTA_pd_csc_matches_csr, tests_run);
     mu_run_test(test_BA_pd_matrices_pd_pd_full_block_B, tests_run);
     mu_run_test(test_BA_pd_matrices_pd_pd_general_B, tests_run);
@@ -575,6 +584,7 @@ int main(void)
     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_jacobian_spd_sum_constraint, 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/jacobian_tests/affine/test_sum.h b/tests/jacobian_tests/affine/test_sum.h
index 12cd22b4..ef12c414 100644
--- a/tests/jacobian_tests/affine/test_sum.h
+++ b/tests/jacobian_tests/affine/test_sum.h
@@ -8,6 +8,7 @@
 #include "minunit.h"
 #include "numerical_diff.h"
 #include "test_helpers.h"
+#include "utils/stacked_pd.h"
 
 const char *test_jacobian_sum_log(void)
 {
@@ -230,3 +231,79 @@ const char *test_jacobian_sum_axis_minus_one_pd_child(void)
     free_expr(sum_node);
     return 0;
 }
+
+/* sum of a stacked_pd child stays stacked_pd for every axis. L = A @ X with A
+   2x3 and X a 3x2 variable has a 2-block spd Jacobian (one block per column of
+   X: rows {0, 1} and {2, 3}, columns 3j .. 3j+2). */
+static expr *sum_spd_child_fixture(int axis, expr **L_out)
+{
+    static const 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);
+    *L_out = L;
+    return new_sum(L, axis);
+}
+
+const char *test_jacobian_sum_spd_child_axis_minus_one(void)
+{
+    expr *L;
+    expr *sum_node = sum_spd_child_fixture(-1, &L);
+    double u[6] = {0.1, 0.2, 0.3, -0.1, -0.2, -0.3};
+    jacobian_init(sum_node);
+    mu_assert("child Jacobian should be spd", L->jacobian->is_stacked_pd);
+    mu_assert("sum Jacobian should be spd", sum_node->jacobian->is_stacked_pd);
+    mu_assert("check_jacobian failed",
+              check_jacobian_num(sum_node, u, NUMERICAL_DIFF_DEFAULT_H));
+    free_expr(sum_node);
+    return 0;
+}
+
+/* axis 0 groups rows {0, 1} and {2, 3}: each output row is fed by one block. */
+const char *test_jacobian_sum_spd_child_axis_0(void)
+{
+    expr *L;
+    expr *sum_node = sum_spd_child_fixture(0, &L);
+    double u[6] = {0.1, 0.2, 0.3, -0.1, -0.2, -0.3};
+    jacobian_init(sum_node);
+    mu_assert("sum Jacobian should be spd", sum_node->jacobian->is_stacked_pd);
+    mu_assert("check_jacobian failed",
+              check_jacobian_num(sum_node, u, NUMERICAL_DIFF_DEFAULT_H));
+    free_expr(sum_node);
+    return 0;
+}
+
+/* axis 1 groups rows {0, 2} and {1, 3}: each output row sums one row from each
+   block, so output row r is [A[r, :], A[r, :]]. Checked exactly through the
+   CSR view on top of the numerical check. */
+const char *test_jacobian_sum_spd_child_axis_1(void)
+{
+    expr *L;
+    expr *sum_node = sum_spd_child_fixture(1, &L);
+    double u[6] = {0.1, 0.2, 0.3, -0.1, -0.2, -0.3};
+    jacobian_init(sum_node);
+    mu_assert("sum Jacobian should be spd", sum_node->jacobian->is_stacked_pd);
+    mu_assert("check_jacobian failed",
+              check_jacobian_num(sum_node, u, NUMERICAL_DIFF_DEFAULT_H));
+
+    sum_node->forward(sum_node, u);
+    eval_jacobian(sum_node);
+    CSR_matrix *J = sum_node->jacobian->to_csr(sum_node->jacobian);
+    mu_assert("shape", J->m == 2 && J->n == 6 && J->nnz == 12);
+    double expected[2][6] = {{1.0, -0.5, 2.0, 1.0, -0.5, 2.0},
+                             {0.5, 1.5, -1.0, 0.5, 1.5, -1.0}};
+    double dense[2][6] = {{0}};
+    for (int r = 0; r < 2; 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 < 2; r++)
+    {
+        mu_assert("axis 1 rows must sum across blocks",
+                  cmp_double_array(dense[r], expected[r], 6));
+    }
+    free_expr(sum_node);
+    return 0;
+}
diff --git a/tests/problem/test_problem.h b/tests/problem/test_problem.h
index 8a385960..c401ae82 100644
--- a/tests/problem/test_problem.h
+++ b/tests/problem/test_problem.h
@@ -470,4 +470,50 @@ const char *test_problem_jacobian_spd_constraint_interleaved(void)
     return 0;
 }
 
+/* A sum over axis 1 of a 2-block spd child is itself a stacked_pd; used as a
+   constraint it must reach the aggregated Jacobian in CSR row order. Output
+   row r of sum(A @ X, axis=1) is [A[r, :], A[r, :]]. */
+const char *test_problem_jacobian_spd_sum_constraint(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 *S = new_sum(L, 1);
+    expr *constraints[1] = {S};
+
+    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", S->jacobian->is_stacked_pd);
+
+    double expected[2][6] = {{1.0, 2.0, 3.0, 1.0, 2.0, 3.0},
+                             {4.0, 5.0, 6.0, 4.0, 5.0, 6.0}};
+    CSR_matrix *J = prob->jacobian;
+    mu_assert("shape", J->m == 2 && J->n == 6 && J->nnz == 12);
+    double dense[2][6] = {{0}};
+    for (int r = 0; r < 2; 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 < 2; r++)
+    {
+        mu_assert("problem Jacobian row of sum constraint",
+                  cmp_double_array(dense[r], expected[r], 6));
+    }
+
+    free_problem(prob);
+    return 0;
+}
+
 #endif /* TEST_PROBLEM_H */
diff --git a/tests/profiling/profile_log_reg.h b/tests/profiling/profile_log_reg.h
index 7fcc69e6..d4dd4736 100644
--- a/tests/profiling/profile_log_reg.h
+++ b/tests/profiling/profile_log_reg.h
@@ -10,7 +10,6 @@
 #include "atoms/elementwise_full_dom.h"
 #include "expr.h"
 #include "minunit.h"
-#include "utils/CSR_sum.h"
 #include "utils/Timer.h"
 #include "utils/permuted_dense.h"
 #include "utils/permuted_dense_linalg.h"
@@ -80,13 +79,11 @@ const char *profile_log_reg(void)
     free(full_rows);
     free(full_cols);
 
-    /* CSR_matrix scaffolding for the row-sum step (PD owns the cached CSR_matrix
-     * view). */
-    CSR_matrix *Jlog_csr = Jlog_M->to_csr(Jlog_M);
-    CSR_matrix *Jobj_csr = new_CSR_matrix(1, n, n);
-    int *iwork = (int *) malloc((size_t) m * n * sizeof(int));
-    int *idx_map = (int *) malloc((size_t) m * n * sizeof(int));
-    sum_all_rows_csr_alloc(Jlog_csr, Jobj_csr, iwork, idx_map);
+    /* Row-sum step through the row_reduce primitive on the PD: every row into
+       output row 0, so Jobj is a 1-row PD. */
+    int *group = (int *) calloc((size_t) m, sizeof(int));
+    matrix *Jobj = Jlog_M->row_reduce_alloc(Jlog_M, group, 1);
+    free(group);
 
     double *d2 = (double *) malloc(m * sizeof(double));
     double *w_ones = (double *) malloc(m * sizeof(double));
@@ -100,8 +97,7 @@ const char *profile_log_reg(void)
     clock_gettime(CLOCK_MONOTONIC, &t_b_jac.start);
     log_obj->local_jacobian(log_obj, log_obj->work->dwork);
     DA_pd_fill_values(log_obj->work->dwork, A_pd, Jlog_pd);
-    memset(Jobj_csr->x, 0, Jobj_csr->nnz * sizeof(double));
-    accumulator(Jlog_csr->x, Jlog_csr->nnz, idx_map, Jobj_csr->x);
+    Jlog_M->row_reduce_fill_values(Jlog_M, Jobj);
     clock_gettime(CLOCK_MONOTONIC, &t_b_jac.end);
     clock_gettime(CLOCK_MONOTONIC, &t_b_hess.start);
     log_obj->local_wsum_hess(log_obj, d2, w_ones);
@@ -122,12 +118,13 @@ const char *profile_log_reg(void)
 
     /* ---- Compare Jacobian (1 x n, both have full sparsity) ---- */
     CSR_matrix *J_a = obj->jacobian->to_csr(obj->jacobian);
-    mu_assert("J n mismatch", J_a->n == Jobj_csr->n);
-    mu_assert("J nnz mismatch", J_a->nnz == Jobj_csr->nnz);
+    CSR_matrix *J_b = Jobj->to_csr(Jobj);
+    mu_assert("J n mismatch", J_a->n == J_b->n);
+    mu_assert("J nnz mismatch", J_a->nnz == J_b->nnz);
     double max_J_diff = 0.0;
     for (int j = 0; j < J_a->nnz; j++)
     {
-        double diff = fabs(J_a->x[j] - Jobj_csr->x[j]);
+        double diff = fabs(J_a->x[j] - J_b->x[j]);
         if (diff > max_J_diff) max_J_diff = diff;
     }
     printf("  Jacobian max abs diff:   %10.3e\n", max_J_diff);
@@ -161,10 +158,7 @@ const char *profile_log_reg(void)
     free(H_a_dense);
     free(d2);
     free(w_ones);
-    free(iwork);
-    free(idx_map);
-    free_CSR_matrix(Jobj_csr);
-    /* Jlog_csr is owned by Jlog_M's cache; released by free_matrix below. */
+    free_matrix(Jobj);
     free_matrix(H_pd_M);
     free_matrix(Jlog_M);
     free_matrix(A_pd_M);
diff --git a/tests/utils/test_permuted_dense.h b/tests/utils/test_permuted_dense.h
index 2dc6954f..53f9833c 100644
--- a/tests/utils/test_permuted_dense.h
+++ b/tests/utils/test_permuted_dense.h
@@ -995,100 +995,4 @@ const char *test_BA_pd_matrices_fast_path(void)
     return 0;
 }
 
-/* Direct vtable tests for sum_row_partition_alloc. The test PD represents a (6, 4)
-   matrix with a (3, 2) dense block at rows {0, 3, 4}, cols {1, 3}. We exercise all
-   three axes; for axis=0 (d1=2) the buckets {0/2, 3/2, 4/2} = {0, 1, 2}
-   are all distinct and non-decreasing (linear-scan dedupe); for axis=1
-   (d1=3) the buckets {0%3, 3%3, 4%3} = {0, 0, 1} collapse two input rows
-   onto the same output row (bitmap dedupe), so multiple idx_map entries
-   point to the same output position — a values-fill pass would scatter-add. */
-const char *test_permuted_dense_sum_all_rows(void)
-{
-    int row_perm[3] = {0, 3, 4};
-    int col_perm[2] = {1, 3};
-    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);
-
-    int idx_map[6];
-    matrix *out = M->sum_row_partition_alloc(M, -1, /*d1 unused*/ 0, idx_map);
-
-    mu_assert("output is PD", out->is_permuted_dense);
-    permuted_dense *opd = (permuted_dense *) out;
-    mu_assert("out m", out->m == 1);
-    mu_assert("out n", out->n == 4);
-    mu_assert("out m0", opd->m0 == 1);
-    mu_assert("out n0", opd->n0 == 2);
-    int expected_row_perm[1] = {0};
-    int expected_col_perm[2] = {1, 3};
-    mu_assert("out row_perm", cmp_int_array(opd->row_perm, expected_row_perm, 1));
-    mu_assert("out col_perm", cmp_int_array(opd->col_perm, expected_col_perm, 2));
-    /* every input row collapses to row 0, so idx_map[i*n0+j] = j */
-    int expected_idx_map[6] = {0, 1, 0, 1, 0, 1};
-    mu_assert("idx_map", cmp_int_array(idx_map, expected_idx_map, 6));
-
-    free_matrix(out);
-    free_matrix(M);
-    return 0;
-}
-
-const char *test_permuted_dense_sum_block_of_rows(void)
-{
-    int row_perm[3] = {0, 3, 4};
-    int col_perm[2] = {1, 3};
-    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);
-
-    int d1 = 2; /* child shape (d1, d2) = (2, 3); output rows = d2 = 3 */
-    int idx_map[6];
-    matrix *out = M->sum_row_partition_alloc(M, 0, d1, idx_map);
-
-    mu_assert("output is PD", out->is_permuted_dense);
-    permuted_dense *opd = (permuted_dense *) out;
-    mu_assert("out m", out->m == 3);
-    mu_assert("out n", out->n == 4);
-    mu_assert("out m0", opd->m0 == 3);
-    mu_assert("out n0", opd->n0 == 2);
-    int expected_row_perm[3] = {0, 1, 2}; /* {0/2, 3/2, 4/2} */
-    int expected_col_perm[2] = {1, 3};
-    mu_assert("out row_perm", cmp_int_array(opd->row_perm, expected_row_perm, 3));
-    mu_assert("out col_perm", cmp_int_array(opd->col_perm, expected_col_perm, 2));
-    int expected_idx_map[6] = {0, 1, 2, 3, 4, 5};
-    mu_assert("idx_map", cmp_int_array(idx_map, expected_idx_map, 6));
-
-    free_matrix(out);
-    free_matrix(M);
-    return 0;
-}
-
-const char *test_permuted_dense_sum_evenly_spaced_rows(void)
-{
-    int row_perm[3] = {0, 3, 4};
-    int col_perm[2] = {1, 3};
-    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);
-
-    int d1 = 3; /* output rows = d1 = 3, buckets = {0%3, 3%3, 4%3} = {0, 0, 1} */
-    int idx_map[6];
-    matrix *out = M->sum_row_partition_alloc(M, 1, d1, idx_map);
-
-    mu_assert("output is PD", out->is_permuted_dense);
-    permuted_dense *opd = (permuted_dense *) out;
-    mu_assert("out m", out->m == 3);
-    mu_assert("out n", out->n == 4);
-    mu_assert("out m0", opd->m0 == 2); /* two distinct buckets {0, 1} */
-    mu_assert("out n0", opd->n0 == 2);
-    int expected_row_perm[2] = {0, 1};
-    int expected_col_perm[2] = {1, 3};
-    mu_assert("out row_perm", cmp_int_array(opd->row_perm, expected_row_perm, 2));
-    mu_assert("out col_perm", cmp_int_array(opd->col_perm, expected_col_perm, 2));
-    /* input rows 0 and 3 both bucket to output row 0 (idx_map → 0, 1);
-       input row 4 buckets to output row 1 (idx_map → 2, 3) */
-    int expected_idx_map[6] = {0, 1, 0, 1, 2, 3};
-    mu_assert("idx_map", cmp_int_array(idx_map, expected_idx_map, 6));
-
-    free_matrix(out);
-    free_matrix(M);
-    return 0;
-}
-
 #endif /* TEST_PERMUTED_DENSE_H */
diff --git a/tests/utils/test_row_reduce.h b/tests/utils/test_row_reduce.h
new file mode 100644
index 00000000..7a468361
--- /dev/null
+++ b/tests/utils/test_row_reduce.h
@@ -0,0 +1,286 @@
+#ifndef TEST_ROW_REDUCE_H
+#define TEST_ROW_REDUCE_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 <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+/* row_reduce_alloc / row_reduce_fill_values across the three matrix kinds:
+   C[j, :] = sum of rows i of A with group[i] == j. The reduction map is bound to
+   the result at alloc time; the fill zeroes C and accumulates from A's current
+   values. */
+
+/* 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_reduce_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);
+}
+
+/* 4x4 stacked_pd whose blocks share column 2:
+     block 0: rows {0, 1}, cols {0, 2}, X = [[1, 2], [3, 4]]
+     block 1: rows {2, 3}, cols {1, 2}, X = [[5, 6], [7, 8]] */
+static matrix *row_reduce_spd_fixture(void)
+{
+    int rp0[2] = {0, 1}, 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};
+    permuted_dense *blocks[2] = {
+        (permuted_dense *) new_permuted_dense(4, 4, 2, 2, rp0, cp0, X0),
+        (permuted_dense *) new_permuted_dense(4, 4, 2, 2, rp1, cp1, X1)};
+    return new_stacked_pd(4, 4, 2, blocks, NULL, NULL);
+}
+
+static void row_reduce_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_reduce_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_reduce_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);
+}
+
+/* Reduce A and its sparse twin with the same group, compare through to_csr
+   after a poisoned fill and again after scaling both sources. Returns the
+   failing assertion message or NULL. */
+static const char *row_reduce_check_against_twin(matrix *A, const int *group,
+                                                 int m_out, matrix **C_out)
+{
+    matrix *A_tw = row_reduce_sparse_twin(A);
+    matrix *C = A->row_reduce_alloc(A, group, m_out);
+    matrix *C_tw = A_tw->row_reduce_alloc(A_tw, group, m_out);
+
+    row_reduce_poison(C);
+    A->row_reduce_fill_values(A, C);
+    A_tw->row_reduce_fill_values(A_tw, C_tw);
+    matrix_values_changed(C);
+    int ok = row_reduce_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);
+    A->row_reduce_fill_values(A, C);
+    A_tw->row_reduce_fill_values(A_tw, C_tw);
+    matrix_values_changed(C);
+    int ok_refill = row_reduce_same_csr(C, C_tw);
+
+    free_matrix(C_tw);
+    free_matrix(A_tw);
+    *C_out = C;
+    if (!ok) return "twin match";
+    if (!ok_refill) return "twin match after refill";
+    return NULL;
+}
+
+/* group = [1, 0, 1, 0], m_out = 3: output row 0 = rows 1 + 3 (columns merge
+   on 4), output row 1 = rows 0 + 2 (row 2 is empty), output row 2 has no
+   members. */
+const char *test_row_reduce_sparse(void)
+{
+    matrix *A = row_reduce_sparse_fixture();
+    int group[4] = {1, 0, 1, 0};
+    matrix *C = A->row_reduce_alloc(A, group, 3);
+
+    int exp_p[4] = {0, 4, 6, 6};
+    int exp_i[6] = {0, 1, 2, 4, 0, 3};
+    mu_assert("shape", C->m == 3 && C->n == 5);
+    mu_assert("sparsity", cmp_sparsity(C, exp_p, exp_i, 3, 6));
+
+    /* the group is bound at alloc time: clobbering the caller's copy is fine */
+    for (int k = 0; k < 4; k++) group[k] = -7;
+
+    /* poisoned output proves the fill zeroes before accumulating */
+    row_reduce_poison(C);
+    A->row_reduce_fill_values(A, C);
+    double exp_x[6] = {5.0, 3.0, 6.0, 11.0, 1.0, 2.0};
+    mu_assert("values", cmp_values(C, exp_x, 6));
+
+    for (int k = 0; k < A->nnz; k++) A->x[k] *= 10.0;
+    for (int k = 0; k < 6; k++) exp_x[k] *= 10.0;
+    A->row_reduce_fill_values(A, C);
+    mu_assert("values after refill", cmp_values(C, exp_x, 6));
+
+    free_matrix(C);
+    free_matrix(A);
+    return 0;
+}
+
+/* (6, 4) pd with dense rows {0, 3, 4} x cols {1, 3}; group = [0, 0, 1, 0, 1, 0]
+   sends dense rows 0 and 3 to output 0 and dense row 4 to output 1
+   (non-monotone in the row index, two rows collapse). */
+const char *test_row_reduce_pd(void)
+{
+    int row_perm[3] = {0, 3, 4};
+    int col_perm[2] = {1, 3};
+    double X[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+    matrix *A = new_permuted_dense(6, 4, 3, 2, row_perm, col_perm, X);
+
+    int group[6] = {0, 0, 1, 0, 1, 0};
+    matrix *C;
+    const char *msg = row_reduce_check_against_twin(A, group, 2, &C);
+    mu_assert(msg ? msg : "", msg == NULL);
+
+    mu_assert("kind preserved", C->is_permuted_dense);
+    permuted_dense *pd = (permuted_dense *) C;
+    mu_assert("shape", C->m == 2 && C->n == 4 && pd->m0 == 2 && pd->n0 == 2);
+    int exp_row_perm[2] = {0, 1};
+    int exp_col_perm[2] = {1, 3};
+    int exp_bound[3] = {0, 0, 1};
+    mu_assert("row_perm", cmp_int_array(pd->row_perm, exp_row_perm, 2));
+    mu_assert("col_perm", cmp_int_array(pd->col_perm, exp_col_perm, 2));
+    mu_assert("bound_iwork", cmp_int_array(pd->bound_iwork, exp_bound, 3));
+    /* sources were scaled by 2 inside the twin check */
+    double exp_X[4] = {8.0, 12.0, 10.0, 12.0};
+    mu_assert("X", cmp_double_array(pd->X, exp_X, 4));
+
+    free_matrix(C);
+    free_matrix(A);
+    return 0;
+}
+
+/* group = [0, 1, 0, 1]: every output row draws one row from each block, so
+   the blocks' shared column 2 must be summed across blocks. One output block
+   with the column union {0, 1, 2}. */
+const char *test_row_reduce_spd_cross_block_accumulate(void)
+{
+    matrix *A = row_reduce_spd_fixture();
+    int group[4] = {0, 1, 0, 1};
+    matrix *C;
+    const char *msg = row_reduce_check_against_twin(A, group, 2, &C);
+    mu_assert(msg ? msg : "", msg == NULL);
+
+    mu_assert("kind preserved", C->is_stacked_pd);
+    stacked_pd *spd = (stacked_pd *) C;
+    mu_assert("one output block", spd->n_blocks == 1);
+    int exp_row_perm[2] = {0, 1};
+    int exp_col_perm[3] = {0, 1, 2};
+    mu_assert("row_perm", cmp_int_array(spd->blocks[0]->row_perm, exp_row_perm, 2));
+    mu_assert("col_perm", cmp_int_array(spd->blocks[0]->col_perm, exp_col_perm, 3));
+    /* sources were scaled by 2: row 0 = [1, 5, 2 + 6], row 1 = [3, 7, 4 + 8] */
+    double exp_X[6] = {2.0, 10.0, 16.0, 6.0, 14.0, 24.0};
+    mu_assert("X", cmp_double_array(spd->blocks[0]->X, exp_X, 6));
+
+    free_matrix(C);
+    free_matrix(A);
+    return 0;
+}
+
+/* group = [0, 0, 1, 1]: each output row is fed by a single block, so the
+   result keeps two blocks with the source column footprints. */
+const char *test_row_reduce_spd_within_block(void)
+{
+    matrix *A = row_reduce_spd_fixture();
+    int group[4] = {0, 0, 1, 1};
+    matrix *C;
+    const char *msg = row_reduce_check_against_twin(A, group, 2, &C);
+    mu_assert(msg ? msg : "", msg == NULL);
+
+    mu_assert("kind preserved", C->is_stacked_pd);
+    stacked_pd *spd = (stacked_pd *) C;
+    mu_assert("two output blocks", spd->n_blocks == 2);
+    int exp_cp0[2] = {0, 2};
+    int exp_cp1[2] = {1, 2};
+    mu_assert("block 0", spd->blocks[0]->m0 == 1 &&
+                             spd->blocks[0]->row_perm[0] == 0 &&
+                             cmp_int_array(spd->blocks[0]->col_perm, exp_cp0, 2));
+    mu_assert("block 1", spd->blocks[1]->m0 == 1 &&
+                             spd->blocks[1]->row_perm[0] == 1 &&
+                             cmp_int_array(spd->blocks[1]->col_perm, exp_cp1, 2));
+    double exp_X0[2] = {8.0, 12.0};  /* 2 * ([1, 2] + [3, 4]) */
+    double exp_X1[2] = {24.0, 28.0}; /* 2 * ([5, 6] + [7, 8]) */
+    mu_assert("X0", cmp_double_array(spd->blocks[0]->X, exp_X0, 2));
+    mu_assert("X1", cmp_double_array(spd->blocks[1]->X, exp_X1, 2));
+
+    free_matrix(C);
+    free_matrix(A);
+    return 0;
+}
+
+/* All rows into one: a single 1-row block over the column union. */
+const char *test_row_reduce_spd_all_to_one(void)
+{
+    matrix *A = row_reduce_spd_fixture();
+    int group[4] = {0, 0, 0, 0};
+    matrix *C;
+    const char *msg = row_reduce_check_against_twin(A, group, 1, &C);
+    mu_assert(msg ? msg : "", msg == NULL);
+
+    mu_assert("kind preserved", C->is_stacked_pd);
+    stacked_pd *spd = (stacked_pd *) C;
+    mu_assert("one output block", spd->n_blocks == 1 && spd->blocks[0]->m0 == 1);
+    int exp_col_perm[3] = {0, 1, 2};
+    mu_assert("col_perm", cmp_int_array(spd->blocks[0]->col_perm, exp_col_perm, 3));
+    double exp_X[3] = {8.0, 24.0, 40.0}; /* 2 * [1 + 3, 5 + 7, 2 + 4 + 6 + 8] */
+    mu_assert("X", cmp_double_array(spd->blocks[0]->X, exp_X, 3));
+
+    free_matrix(C);
+    free_matrix(A);
+    return 0;
+}
+
+#ifdef SP_TRACK_MEMORY
+typedef struct
+{
+    matrix *A;
+    matrix *C;
+} row_reduce_fill_args;
+
+static void run_row_reduce_fill(const void *ctx)
+{
+    const row_reduce_fill_args *a = (const row_reduce_fill_args *) ctx;
+    a->A->row_reduce_fill_values(a->A, a->C);
+}
+
+/* The spd fill reduces per block into the raw spd, then coalesce-accumulates;
+   none of it may allocate once the alloc phase is done. */
+const char *test_row_reduce_spd_fill_no_transient_alloc(void)
+{
+    matrix *A = row_reduce_spd_fixture();
+    int group[4] = {0, 1, 0, 1};
+    matrix *C = A->row_reduce_alloc(A, group, 2);
+    row_reduce_fill_args args = {A, C};
+    run_row_reduce_fill(&args); /* warm-up */
+    mu_assert("spd row_reduce fill must not allocate",
+              fill_is_alloc_free(run_row_reduce_fill, &args));
+    free_matrix(C);
+    free_matrix(A);
+    return 0;
+}
+#endif
+
+#endif /* TEST_ROW_REDUCE_H */