Bump sparsediff_linalg to v0.5.0

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));

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));
 }
 

src/atoms/affine/diag_vec.c

@@ -50,16 +50,25 @@ 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)
 {
     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)

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)

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));
 }

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;

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)

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)

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)

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)

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));
 }
 

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)

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));

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));

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)

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

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,