diff --git a/sparsediff_linalg b/sparsediff_linalg
index cba21acb..4c121859 160000
--- a/sparsediff_linalg
+++ b/sparsediff_linalg
@@ -1 +1 @@
-Subproject commit cba21acb5ecc18c497f70d7bb98522355c25d53f
+Subproject commit 4c1218595c0e4b5e5dc3af57dd1da3550f998e58
diff --git a/src/atoms/affine/broadcast.c b/src/atoms/affine/broadcast.c
index ced71017..cb6bcb1d 100644
--- a/src/atoms/affine/broadcast.c
+++ b/src/atoms/affine/broadcast.c
@@ -107,7 +107,7 @@ static void wsum_hess_init_impl(expr *node)
     wsum_hess_init(x);
 
     /* Same sparsity as child - weights get summed */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 
     /* allocate space for weight vector */
     node->work->dwork = sp_malloc(node->size * sizeof(double));
diff --git a/src/atoms/affine/convolve.c b/src/atoms/affine/convolve.c
index 61c91c78..6bbfe034 100644
--- a/src/atoms/affine/convolve.c
+++ b/src/atoms/affine/convolve.c
@@ -138,7 +138,7 @@ static void wsum_hess_init_impl(expr *node)
     convolve_expr *cnode = (convolve_expr *) node;
 
     wsum_hess_init(child);
-    node->wsum_hess = child->wsum_hess->copy_sparsity(child->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(child->wsum_hess);
     node->work->dwork = (double *) sp_malloc(cnode->n * sizeof(double));
 }
 
diff --git a/src/atoms/affine/diag_vec.c b/src/atoms/affine/diag_vec.c
index 4b48b244..80a687c3 100644
--- a/src/atoms/affine/diag_vec.c
+++ b/src/atoms/affine/diag_vec.c
@@ -50,8 +50,16 @@ static void jacobian_init_impl(expr *node)
     expr *x = node->left;
     jacobian_init(x);
 
-    /* output type matches child's; rows i*(n+1) hold child row i, others zero. */
-    node->jacobian = x->jacobian->diag_vec_alloc(x->jacobian);
+    /* row i of the child's jacobian lands at row i*(n+1) (the diagonal
+       positions), every other row is zero: a row_reduce with singleton groups */
+    int n = x->size;
+    int *group = (int *) sp_malloc(n * sizeof(int));
+    for (int i = 0; i < n; i++)
+    {
+        group[i] = i * (n + 1);
+    }
+    node->jacobian = x->jacobian->row_reduce_alloc(x->jacobian, group, n * n);
+    sp_free(group);
 }
 
 static void eval_jacobian_impl(expr *node)
@@ -59,7 +67,8 @@ static void eval_jacobian_impl(expr *node)
     eval_jacobian(node->left);
 
     /* fill the diagonal rows of the preallocated output. */
-    node->left->jacobian->diag_vec_fill_values(node->left->jacobian, node->jacobian);
+    node->left->jacobian->row_reduce_fill_values(node->left->jacobian,
+                                                 node->jacobian);
 }
 
 static void wsum_hess_init_impl(expr *node)
@@ -74,7 +83,7 @@ static void wsum_hess_init_impl(expr *node)
 
     /* Copy child's Hessian structure (diag_vec is linear, so its own Hessian is
      * zero) */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 }
 
 static void eval_wsum_hess_impl(expr *node, const double *w)
diff --git a/src/atoms/affine/index.c b/src/atoms/affine/index.c
index 0aac5882..b04c6234 100644
--- a/src/atoms/affine/index.c
+++ b/src/atoms/affine/index.c
@@ -96,7 +96,7 @@ static void wsum_hess_init_impl(expr *node)
        many numerical zeros in child->wsum_hess that are actually
        structural zeros, but we do not try to exploit that sparsity
        right now. */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 }
 
 static void eval_wsum_hess_impl(expr *node, const double *w)
diff --git a/src/atoms/affine/kron.c b/src/atoms/affine/kron.c
index c74523cf..180f6b7c 100644
--- a/src/atoms/affine/kron.c
+++ b/src/atoms/affine/kron.c
@@ -138,7 +138,7 @@ static void wsum_hess_init_impl(expr *node)
     expr *child = node->left;
 
     wsum_hess_init(child);
-    node->wsum_hess = child->wsum_hess->copy_sparsity(child->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(child->wsum_hess);
     /* backprop workspace: one weight per child entry */
     node->work->dwork = (double *) sp_malloc(child->size * sizeof(double));
 }
diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c
index 2088b617..8640e153 100644
--- a/src/atoms/affine/left_matmul.c
+++ b/src/atoms/affine/left_matmul.c
@@ -18,6 +18,7 @@
 #include "atoms/affine.h"
 #include "sparsediff_linalg/dense.h"
 #include "sparsediff_linalg/matmul.h"
+#include "sparsediff_linalg/matmul_kron.h"
 #include "sparsediff_linalg/matrix.h"
 #include "sparsediff_linalg/tracked_alloc.h"
 #include "subexpr.h"
@@ -81,7 +82,7 @@ static void forward(expr *node, const double *u)
     /* y = A_kron @ vec(f(x)) */
     matrix *A = lnode->A;
     int n_blocks = lnode->n_blocks;
-    A->block_left_mult_vec(A, x->value, node->value, n_blocks);
+    BA_kron_vec_matrices(A, x->value, node->value, n_blocks);
 }
 
 static bool is_affine(const expr *node)
@@ -168,13 +169,9 @@ static void jacobian_init_sparse(expr *node)
         return;
     }
 
-    /* initialize this node's jacobian */
-    lnode->Jchild_CSC =
-        csr_to_csc_alloc(x->jacobian->to_csr(x->jacobian), node->work->iwork);
-    lnode->J_CSC = lnode->A->block_left_mult_sparsity(lnode->A, lnode->Jchild_CSC,
-                                                      lnode->n_blocks);
+    /* J_y = kron(I_p, A) J_f */
     node->jacobian =
-        new_sparse_matrix(csc_to_csr_alloc(lnode->J_CSC, lnode->csc_to_csr_work));
+        BA_sparse_kron_matrices_alloc(lnode->A, lnode->n_blocks, x->jacobian);
 }
 
 /* eval_jacobian for sparse left matmul */
@@ -191,14 +188,8 @@ static void eval_jacobian_sparse(expr *node)
         return;
     }
 
-    /* evaluate this node's jacobian */
-    CSC_matrix *Jchild_CSC = lnode->Jchild_CSC;
-    CSC_matrix *J_CSC = lnode->J_CSC;
-    csr_to_csc_fill_values(x->jacobian->to_csr(x->jacobian), Jchild_CSC,
-                           node->work->iwork);
-    lnode->A->block_left_mult_values(lnode->A, Jchild_CSC, J_CSC);
-    csc_to_csr_fill_values(J_CSC, node->jacobian->to_csr(node->jacobian),
-                           lnode->csc_to_csr_work);
+    BA_sparse_kron_matrices_fill_values(lnode->A, lnode->n_blocks, x->jacobian,
+                                        node->jacobian);
 }
 
 static void wsum_hess_init_impl(expr *node)
@@ -208,7 +199,7 @@ static void wsum_hess_init_impl(expr *node)
     wsum_hess_init(x);
 
     /* allocate this node's hessian with the same sparsity as child's */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 
     /* work for computing A^T w*/
     int n_blocks = ((left_matmul_expr *) node)->n_blocks;
@@ -223,7 +214,7 @@ static void eval_wsum_hess_impl(expr *node, const double *w)
     /* compute A^T w*/
     matrix *AT = lnode->AT;
     int n_blocks = lnode->n_blocks;
-    AT->block_left_mult_vec(AT, w, node->work->dwork, n_blocks);
+    BA_kron_vec_matrices(AT, w, node->work->dwork, n_blocks);
 
     eval_wsum_hess(node->left, node->work->dwork);
     memcpy(node->wsum_hess->x, node->left->wsum_hess->x,
@@ -299,8 +290,8 @@ expr *new_left_matmul(expr *param_node, expr *u, const CSR_matrix *A)
 
     /* store A and AT. new_sparse_matrix takes ownership, so clone first. */
     lnode->A = new_sparse_matrix(new_csr(A));
-    lnode->AT = lnode->A->transpose_alloc(lnode->A);
-    lnode->A->transpose_fill_values(lnode->A, lnode->AT);
+    lnode->AT = transpose_matrices_alloc(lnode->A);
+    transpose_matrices_fill_values(lnode->A, lnode->AT);
 
     /* parameter support */
     lnode->param_source = param_node;
@@ -360,8 +351,8 @@ expr *new_left_matmul_dense(expr *param_node, expr *u, int m, int n,
         }
 
         lnode->A = new_permuted_dense_full(m, n, data);
-        lnode->AT = lnode->A->transpose_alloc(lnode->A);
-        lnode->A->transpose_fill_values(lnode->A, lnode->AT);
+        lnode->AT = transpose_matrices_alloc(lnode->A);
+        transpose_matrices_fill_values(lnode->A, lnode->AT);
     }
 
     return node;
diff --git a/src/atoms/affine/neg.c b/src/atoms/affine/neg.c
index 6367566e..a3e2c75a 100644
--- a/src/atoms/affine/neg.c
+++ b/src/atoms/affine/neg.c
@@ -40,7 +40,7 @@ static void jacobian_init_impl(expr *node)
     jacobian_init(x);
 
     /* same sparsity pattern as child */
-    node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
+    node->jacobian = copy_sparsity_matrices_alloc(x->jacobian);
 }
 
 static void eval_jacobian_impl(expr *node)
@@ -63,7 +63,7 @@ static void wsum_hess_init_impl(expr *node)
     wsum_hess_init(x);
 
     /* same sparsity pattern as child (polymorphic copy) */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 }
 
 static void eval_wsum_hess_impl(expr *node, const double *w)
diff --git a/src/atoms/affine/promote.c b/src/atoms/affine/promote.c
index 103c1bd6..1d70a445 100644
--- a/src/atoms/affine/promote.c
+++ b/src/atoms/affine/promote.c
@@ -59,7 +59,7 @@ static void eval_jacobian_impl(expr *node)
 static void wsum_hess_init_impl(expr *node)
 {
     wsum_hess_init(node->left);
-    node->wsum_hess = node->left->wsum_hess->copy_sparsity(node->left->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(node->left->wsum_hess);
 }
 
 static void eval_wsum_hess_impl(expr *node, const double *w)
diff --git a/src/atoms/affine/reshape.c b/src/atoms/affine/reshape.c
index 3e5719a3..c89f2957 100644
--- a/src/atoms/affine/reshape.c
+++ b/src/atoms/affine/reshape.c
@@ -36,7 +36,7 @@ static void jacobian_init_impl(expr *node)
 {
     expr *x = node->left;
     jacobian_init(x);
-    node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
+    node->jacobian = copy_sparsity_matrices_alloc(x->jacobian);
 }
 
 static void eval_jacobian_impl(expr *node)
@@ -50,7 +50,7 @@ static void wsum_hess_init_impl(expr *node)
 {
     expr *x = node->left;
     wsum_hess_init(x);
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 }
 
 static void eval_wsum_hess_impl(expr *node, const double *w)
diff --git a/src/atoms/affine/scalar_mult.c b/src/atoms/affine/scalar_mult.c
index bb48eb1b..617c7d54 100644
--- a/src/atoms/affine/scalar_mult.c
+++ b/src/atoms/affine/scalar_mult.c
@@ -59,7 +59,7 @@ static void jacobian_init_impl(expr *node)
     jacobian_init(x);
 
     /* same sparsity as child */
-    node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
+    node->jacobian = copy_sparsity_matrices_alloc(x->jacobian);
 }
 
 static void eval_jacobian_impl(expr *node)
@@ -85,7 +85,7 @@ static void wsum_hess_init_impl(expr *node)
     wsum_hess_init(x);
 
     /* same sparsity as child */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 }
 
 static void eval_wsum_hess_impl(expr *node, const double *w)
diff --git a/src/atoms/affine/sum.c b/src/atoms/affine/sum.c
index 5304ff1c..df4a8adc 100644
--- a/src/atoms/affine/sum.c
+++ b/src/atoms/affine/sum.c
@@ -129,7 +129,7 @@ static void wsum_hess_init_impl(expr *node)
     wsum_hess_init(child);
 
     /* we never have to store more than the child's nnz */
-    node->wsum_hess = child->wsum_hess->copy_sparsity(child->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(child->wsum_hess);
     node->work->dwork = sp_malloc(child->size * sizeof(double));
 }
 
diff --git a/src/atoms/affine/trace.c b/src/atoms/affine/trace.c
index 6c7c4ce4..175e1569 100644
--- a/src/atoms/affine/trace.c
+++ b/src/atoms/affine/trace.c
@@ -87,7 +87,7 @@ static void wsum_hess_init_impl(expr *node)
        contribution to wsum_hess of entries of the child that will always have
        zero weight in eval_wsum_hess. We do this for simplicity. But the Hessian
        can for sure be made more sophisticated. */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 }
 
 static void eval_wsum_hess_impl(expr *node, const double *w)
diff --git a/src/atoms/affine/transpose.c b/src/atoms/affine/transpose.c
index 6dd4ab28..9d6b2698 100644
--- a/src/atoms/affine/transpose.c
+++ b/src/atoms/affine/transpose.c
@@ -71,7 +71,7 @@ static void wsum_hess_init_impl(expr *node)
     wsum_hess_init(x);
 
     /* same sparsity pattern as child */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 
     /* for computing Kw where K is the commutation matrix */
     node->work->dwork = (double *) sp_malloc(node->size * sizeof(double));
diff --git a/src/atoms/affine/vector_mult.c b/src/atoms/affine/vector_mult.c
index c56cc426..3b598975 100644
--- a/src/atoms/affine/vector_mult.c
+++ b/src/atoms/affine/vector_mult.c
@@ -59,7 +59,7 @@ static void jacobian_init_impl(expr *node)
     jacobian_init(x);
 
     /* same sparsity as child */
-    node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
+    node->jacobian = copy_sparsity_matrices_alloc(x->jacobian);
 }
 
 static void eval_jacobian_impl(expr *node)
@@ -71,7 +71,7 @@ static void eval_jacobian_impl(expr *node)
     eval_jacobian(x);
 
     /* row-wise scale child's jacobian: diag(a) @ Jx */
-    x->jacobian->DA_fill_values(a, x->jacobian, node->jacobian);
+    DA_matrices_fill_values(a, x->jacobian, node->jacobian);
 }
 
 static void wsum_hess_init_impl(expr *node)
@@ -82,7 +82,7 @@ static void wsum_hess_init_impl(expr *node)
     wsum_hess_init(x);
 
     /* same sparsity as child */
-    node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->wsum_hess = copy_sparsity_matrices_alloc(x->wsum_hess);
 
     /* workspace for storing scaled weights */
     node->work->dwork = (double *) sp_malloc(node->size * sizeof(double));
diff --git a/src/atoms/bivariate_full_dom/matmul.c b/src/atoms/bivariate_full_dom/matmul.c
index 87b8001e..43493f6b 100644
--- a/src/atoms/bivariate_full_dom/matmul.c
+++ b/src/atoms/bivariate_full_dom/matmul.c
@@ -21,6 +21,7 @@
 #include "sparsediff_linalg/CSR_matrix.h"
 #include "sparsediff_linalg/cblas_wrapper.h"
 #include "sparsediff_linalg/dense.h"
+#include "sparsediff_linalg/matmul_kron.h"
 #include "sparsediff_linalg/matrix.h"
 #include "sparsediff_linalg/matrix_sum.h"
 #include "sparsediff_linalg/tracked_alloc.h"
@@ -428,7 +429,7 @@ static void wsum_hess_init_chain_rule(expr *node)
     mnode->B = new_sparse_matrix(build_cross_hessian_sparsity(m, k, n));
     mnode->BJg = BA_matrices_alloc(mnode->B, g->jacobian);
     mnode->C = BTA_matrices_alloc(mnode->BJg, f->jacobian);
-    mnode->CT = mnode->C->transpose_alloc(mnode->C);
+    mnode->CT = transpose_matrices_alloc(mnode->C);
 
     /* initialize Hessians of children */
     wsum_hess_init(f);
@@ -462,7 +463,7 @@ static void eval_wsum_hess_chain_rule(expr *node, const double *w)
     matrix_values_changed(mnode->B);
     BA_matrices_fill_values(mnode->B, g->jacobian, mnode->BJg);
     BTA_matrices_fill_values(mnode->BJg, f->jacobian, mnode->C);
-    mnode->C->transpose_fill_values(mnode->C, mnode->CT);
+    transpose_matrices_fill_values(mnode->C, mnode->CT);
 
     /* compute Hessian of f */
     if (!is_f_affine)
diff --git a/src/atoms/bivariate_full_dom/multiply.c b/src/atoms/bivariate_full_dom/multiply.c
index 7a44939d..2661e6de 100644
--- a/src/atoms/bivariate_full_dom/multiply.c
+++ b/src/atoms/bivariate_full_dom/multiply.c
@@ -152,7 +152,7 @@ static void wsum_hess_init_impl(expr *node)
 
         /* compute sparsity of C and prepare CT */
         matrix *C = BTA_matrices_alloc(x->jacobian, y->jacobian);
-        matrix *CT = C->transpose_alloc(C);
+        matrix *CT = transpose_matrices_alloc(C);
 
         /* initialize wsum_hessians of children */
         wsum_hess_init(x);
@@ -189,7 +189,7 @@ static void eval_wsum_hess_impl(expr *node, const double *w)
         // ---------------------------------------------------------------
         elementwise_mult_expr *mul_node = (elementwise_mult_expr *) node;
         BTDA_matrices_fill_values(x->jacobian, w, y->jacobian, mul_node->C);
-        mul_node->C->transpose_fill_values(mul_node->C, mul_node->CT);
+        transpose_matrices_fill_values(mul_node->C, mul_node->CT);
 
         // ---------------------------------------------------------------
         //              compute term2 and term 3
diff --git a/src/atoms/elementwise_full_dom/common.c b/src/atoms/elementwise_full_dom/common.c
index 46083486..5aec6f8c 100644
--- a/src/atoms/elementwise_full_dom/common.c
+++ b/src/atoms/elementwise_full_dom/common.c
@@ -18,6 +18,7 @@
 #include "atoms/elementwise_full_dom.h"
 #include "sparsediff_linalg/CSC_matrix.h"
 #include "sparsediff_linalg/CSR_matrix.h"
+#include "sparsediff_linalg/matmul.h"
 #include "sparsediff_linalg/matrix.h"
 #include "sparsediff_linalg/matrix_sum.h"
 #include "sparsediff_linalg/tracked_alloc.h"
@@ -46,7 +47,7 @@ void jacobian_init_elementwise(expr *node)
     {
         /* jacobian of h(x) = f(g(x)) is Jf @ Jg, and here Jf is diagonal */
         jacobian_init(child);
-        node->jacobian = child->jacobian->copy_sparsity(child->jacobian);
+        node->jacobian = copy_sparsity_matrices_alloc(child->jacobian);
         node->work->dwork = (double *) sp_malloc(node->size * sizeof(double));
         node->work->local_jac_diag =
             (double *) sp_malloc(node->size * sizeof(double));
@@ -68,8 +69,7 @@ void eval_jacobian_elementwise(expr *node)
         node->local_jacobian(node, node->work->local_jac_diag);
         memcpy(node->work->dwork, node->work->local_jac_diag,
                node->size * sizeof(double));
-        child->jacobian->DA_fill_values(node->work->dwork, child->jacobian,
-                                        node->jacobian);
+        DA_matrices_fill_values(node->work->dwork, child->jacobian, node->jacobian);
     }
 }
 
@@ -105,17 +105,16 @@ void wsum_hess_init_elementwise(expr *node)
             For elementwise functions, D is diagonal. */
         if (child->is_affine(child))
         {
-            node->wsum_hess = child->jacobian->ATA_alloc(child->jacobian);
+            node->wsum_hess = ATA_matrices_alloc(child->jacobian);
         }
         else
         {
             /* term1: Jg^T @ D @ Jg */
-            node->work->hess_term1 = child->jacobian->ATA_alloc(child->jacobian);
+            node->work->hess_term1 = ATA_matrices_alloc(child->jacobian);
 
             /* term2: child's Hessian (mirror its sparsity polymorphically) */
             wsum_hess_init(child);
-            node->work->hess_term2 =
-                child->wsum_hess->copy_sparsity(child->wsum_hess);
+            node->work->hess_term2 = copy_sparsity_matrices_alloc(child->wsum_hess);
 
             /* wsum_hess = term1 + term2 */
             int max_nnz = node->work->hess_term1->nnz + node->work->hess_term2->nnz;
@@ -140,15 +139,15 @@ void eval_wsum_hess_elementwise(expr *node, const double *w)
         if (child->is_affine(child))
         {
             node->local_wsum_hess(node, node->work->dwork, w);
-            child->jacobian->ATDA_fill_values(child->jacobian, node->work->dwork,
-                                              node->wsum_hess);
+            ATDA_matrices_fill_values(child->jacobian, node->work->dwork,
+                                      node->wsum_hess);
         }
         else
         {
             /* term1: Jg^T @ D @ Jg */
             node->local_wsum_hess(node, node->work->dwork, w);
-            child->jacobian->ATDA_fill_values(child->jacobian, node->work->dwork,
-                                              node->work->hess_term1);
+            ATDA_matrices_fill_values(child->jacobian, node->work->dwork,
+                                      node->work->hess_term1);
 
             /* term2: child Hessian with weight Jf^T w */
             memcpy(node->work->dwork, node->work->local_jac_diag,
diff --git a/src/atoms/other/quad_form.c b/src/atoms/other/quad_form.c
index dc3266ba..1d283a0d 100644
--- a/src/atoms/other/quad_form.c
+++ b/src/atoms/other/quad_form.c
@@ -18,6 +18,7 @@
 #include "atoms/non_elementwise_full_dom.h"
 #include "sparsediff_linalg/cblas_wrapper.h"
 #include "sparsediff_linalg/matmul.h"
+#include "sparsediff_linalg/matmul_kron.h"
 #include "sparsediff_linalg/matrix.h"
 #include "sparsediff_linalg/matrix_sum.h"
 #include "sparsediff_linalg/tracked_alloc.h"
@@ -60,7 +61,7 @@ static void forward(expr *node, const double *u)
 
     /* dwork = Q @ x; value = x' (Q x) */
     matrix *Q = qnode->Q;
-    Q->block_left_mult_vec(Q, x->value, node->work->dwork, 1);
+    BA_kron_vec_matrices(Q, x->value, node->work->dwork, 1);
     node->value[0] = cblas_ddot(qnode->n, x->value, 1, node->work->dwork, 1);
 }
 
@@ -168,7 +169,7 @@ static void wsum_hess_init_impl(expr *node)
     node->work->hess_term1 = BTA_matrices_alloc(x->jacobian, qnode->QJf);
 
     wsum_hess_init(x);
-    node->work->hess_term2 = x->wsum_hess->copy_sparsity(x->wsum_hess);
+    node->work->hess_term2 = copy_sparsity_matrices_alloc(x->wsum_hess);
 
     /* hess = term1 + term2 */
     int max_nnz = node->work->hess_term1->nnz + node->work->hess_term2->nnz;
diff --git a/tests/profiling/profile_hessian_exp_AX.h b/tests/profiling/profile_hessian_exp_AX.h
index 3ff48924..5e7359ff 100644
--- a/tests/profiling/profile_hessian_exp_AX.h
+++ b/tests/profiling/profile_hessian_exp_AX.h
@@ -12,6 +12,7 @@
 #include "expr.h"
 #include "minunit.h"
 #include "sparsediff_linalg/CSR_matrix.h"
+#include "sparsediff_linalg/matmul.h"
 #include "sparsediff_linalg/matrix.h"
 #include "sparsediff_linalg/tracked_alloc.h"
 
@@ -25,9 +26,9 @@
    - Approach 2: manually construct the Jacobian
      J = ∂vec(A @ X)/∂vec(X) = I_n ⊗ A as a stacked_pd with n blocks
      (each block carries A's values at the corresponding row/col
-     window). Call J->ATA_alloc (alloc) and J->ATDA_fill_values(J, d, H)
-     (fill), where d_k = w_k * exp((A@X)_k). The AX compute is
-     done OUTSIDE the timer per user request.
+     window). Call ATA_matrices_alloc(J) (alloc) and ATDA_matrices_fill_values(J, d,
+   H) (fill), where d_k = w_k * exp((A@X)_k). The AX compute is done OUTSIDE the
+   timer per user request.
 
    Both Hessians are densified to a 2500x2500 row-major buffer and
    compared via max-abs diff (asserted < 1e-9). */
@@ -126,12 +127,12 @@ const char *profile_hessian_exp_AX(void)
     free(row_perm);
     matrix *J = new_stacked_pd(n_vars, n_vars, n, j_blocks, NULL, NULL);
     free(j_blocks);
-    matrix *H2 = J->ATA_alloc(J);
+    matrix *H2 = ATA_matrices_alloc(J);
     clock_gettime(CLOCK_MONOTONIC, &t2a.end);
 
     Timer t2f;
     clock_gettime(CLOCK_MONOTONIC, &t2f.start);
-    J->ATDA_fill_values(J, d, H2);
+    ATDA_matrices_fill_values(J, d, H2);
     clock_gettime(CLOCK_MONOTONIC, &t2f.end);
 
     /* ------------------------------------------------------------ */
diff --git a/tests/profiling/profile_log_reg.h b/tests/profiling/profile_log_reg.h
index 7bfd45a3..f16e6db6 100644
--- a/tests/profiling/profile_log_reg.h
+++ b/tests/profiling/profile_log_reg.h
@@ -11,6 +11,7 @@
 #include "atoms/elementwise_full_dom.h"
 #include "expr.h"
 #include "minunit.h"
+#include "sparsediff_linalg/matmul.h"
 #include "sparsediff_linalg/matrix.h"
 
 /* Profile and validate Jacobian + Hessian of obj = sum(logistic(A x)).
@@ -70,7 +71,7 @@ const char *profile_log_reg(void)
 
     matrix *A_pd_M = new_permuted_dense(m, n, m, n, full_rows, full_cols, A_data);
     matrix *Jlog_M = new_permuted_dense(m, n, m, n, full_rows, full_cols, NULL);
-    matrix *H_pd_M = A_pd_M->ATA_alloc(A_pd_M);
+    matrix *H_pd_M = ATA_matrices_alloc(A_pd_M);
 
     free(full_rows);
     free(full_cols);
@@ -92,12 +93,12 @@ const char *profile_log_reg(void)
        dwork) as sigmas read by local_wsum_hess. */
     clock_gettime(CLOCK_MONOTONIC, &t_b_jac.start);
     log_obj->local_jacobian(log_obj, log_obj->work->dwork);
-    A_pd_M->DA_fill_values(log_obj->work->dwork, A_pd_M, Jlog_M);
+    DA_matrices_fill_values(log_obj->work->dwork, A_pd_M, Jlog_M);
     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);
-    A_pd_M->ATDA_fill_values(A_pd_M, d2, H_pd_M);
+    ATDA_matrices_fill_values(A_pd_M, d2, H_pd_M);
     clock_gettime(CLOCK_MONOTONIC, &t_b_hess.end);
     double sec_b_jac = GET_ELAPSED_SECONDS(t_b_jac);
     double sec_b_hess = GET_ELAPSED_SECONDS(t_b_hess);
diff --git a/tests/wsum_hess/affine/test_left_matmul.h b/tests/wsum_hess/affine/test_left_matmul.h
index 7833b971..2ff4c1ce 100644
--- a/tests/wsum_hess/affine/test_left_matmul.h
+++ b/tests/wsum_hess/affine/test_left_matmul.h
@@ -200,7 +200,7 @@ const char *test_wsum_hess_left_matmul_matrix(void)
 /* Regression test for the Phase 3 transpose_fill_values omission in
    new_left_matmul_dense. Mirrors the Python failure
        cp.sum(A @ cp.exp(X))   with X a (2,2) Variable and A a 2x2 numpy array.
-   eval_wsum_hess reads lnode->AT->X via AT->block_left_mult_vec; before the
+   eval_wsum_hess reads lnode->AT->X via BA_kron_vec_matrices(AT, ...); before the
    fix, AT->X was uninitialized memory (transpose_alloc allocates without
    filling) and the analytic Hessian disagreed with finite differences. */
 const char *test_wsum_hess_left_matmul_dense_matrix_exp(void)