copy sparsity function (#61)

* add work struct for expression

* run formatter

* copy sparsity function

Author: dance858

include/utils/CSR_Matrix.h

@@ -28,6 +28,7 @@ typedef struct CSR_Matrix
 /* constructors and destructors */
 CSR_Matrix *new_csr_matrix(int m, int n, int nnz);
 CSR_Matrix *new_csr(const CSR_Matrix *A);
+CSR_Matrix *new_csr_copy_sparsity(const CSR_Matrix *A);
 void free_csr_matrix(CSR_Matrix *matrix);
 void copy_csr_matrix(const CSR_Matrix *A, CSR_Matrix *C);
 

src/affine/broadcast.c

@@ -191,9 +191,7 @@ static void wsum_hess_init(expr *node)
     x->wsum_hess_init(x);
 
     /* Same sparsity as child - weights get summed */
-    node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, x->wsum_hess->nnz);
-    memcpy(node->wsum_hess->p, x->wsum_hess->p, (x->wsum_hess->m + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, x->wsum_hess->i, x->wsum_hess->nnz * sizeof(int));
+    node->wsum_hess = new_csr_copy_sparsity(x->wsum_hess);
 
     /* allocate space for weight vector */
     node->work->dwork = malloc(node->size * sizeof(double));

src/affine/diag_vec.c

@@ -105,9 +105,7 @@ static void wsum_hess_init(expr *node)
     /* Copy child's Hessian structure (diag_vec is linear, so its own Hessian is
      * zero) */
     CSR_Matrix *Hx = x->wsum_hess;
-    node->wsum_hess = new_csr_matrix(Hx->m, Hx->n, Hx->nnz);
-    memcpy(node->wsum_hess->p, Hx->p, (Hx->m + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, Hx->i, Hx->nnz * sizeof(int));
+    node->wsum_hess = new_csr_copy_sparsity(Hx);
 }
 
 static void eval_wsum_hess(expr *node, const double *w)

src/affine/index.c

@@ -113,9 +113,7 @@ static void wsum_hess_init(expr *node)
        structural zeros, but we do not try to exploit that sparsity
        right now. */
     CSR_Matrix *Hx = x->wsum_hess;
-    node->wsum_hess = new_csr_matrix(Hx->m, Hx->n, Hx->nnz);
-    memcpy(node->wsum_hess->p, Hx->p, (Hx->m + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, Hx->i, Hx->nnz * sizeof(int));
+    node->wsum_hess = new_csr_copy_sparsity(Hx);
 }
 
 static void eval_wsum_hess(expr *node, const double *w)

src/affine/neg.c

@@ -39,11 +39,7 @@ static void jacobian_init(expr *node)
     x->jacobian_init(x);
 
     /* same sparsity pattern as child */
-    node->jacobian = new_csr_matrix(node->size, node->n_vars, x->jacobian->nnz);
-
-    /* copy row pointers and column indices (sparsity pattern is constant) */
-    memcpy(node->jacobian->p, x->jacobian->p, (x->jacobian->m + 1) * sizeof(int));
-    memcpy(node->jacobian->i, x->jacobian->i, x->jacobian->nnz * sizeof(int));
+    node->jacobian = new_csr_copy_sparsity(x->jacobian);
 }
 
 static void eval_jacobian(expr *node)
@@ -68,11 +64,7 @@ static void wsum_hess_init(expr *node)
 
     /* same sparsity pattern as child */
     CSR_Matrix *child_hess = x->wsum_hess;
-    node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, child_hess->nnz);
-
-    /* copy row pointers and column indices (sparsity pattern is constant) */
-    memcpy(node->wsum_hess->p, child_hess->p, (child_hess->m + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, child_hess->i, child_hess->nnz * sizeof(int));
+    node->wsum_hess = new_csr_copy_sparsity(child_hess);
 }
 
 static void eval_wsum_hess(expr *node, const double *w)

src/affine/promote.c

@@ -78,12 +78,7 @@ static void wsum_hess_init(expr *node)
 
     /* same sparsity as child since we're summing weights */
     CSR_Matrix *child_hess = node->left->wsum_hess;
-    node->wsum_hess = new_csr_matrix(child_hess->m, child_hess->n, child_hess->nnz);
-
-    /* copy sparsity pattern */
-    memcpy(node->wsum_hess->p, child_hess->p, (child_hess->m + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, child_hess->i, child_hess->nnz * sizeof(int));
-    node->wsum_hess->nnz = child_hess->nnz;
+    node->wsum_hess = new_csr_copy_sparsity(child_hess);
 }
 
 static void eval_wsum_hess(expr *node, const double *w)

src/affine/reshape.c

@@ -35,10 +35,7 @@ static void jacobian_init(expr *node)
 {
     expr *x = node->left;
     x->jacobian_init(x);
-    node->jacobian = new_csr_matrix(node->size, node->n_vars, x->jacobian->nnz);
-    CSR_Matrix *jac = node->jacobian;
-    memcpy(jac->p, x->jacobian->p, (x->size + 1) * sizeof(int));
-    memcpy(jac->i, x->jacobian->i, x->jacobian->nnz * sizeof(int));
+    node->jacobian = new_csr_copy_sparsity(x->jacobian);
 }
 
 static void eval_jacobian(expr *node)
@@ -52,10 +49,7 @@ static void wsum_hess_init(expr *node)
 {
     expr *x = node->left;
     x->wsum_hess_init(x);
-    node->wsum_hess =
-        new_csr_matrix(x->wsum_hess->m, x->wsum_hess->n, x->wsum_hess->nnz);
-    memcpy(node->wsum_hess->p, x->wsum_hess->p, (x->wsum_hess->m + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, x->wsum_hess->i, x->wsum_hess->nnz * sizeof(int));
+    node->wsum_hess = new_csr_copy_sparsity(x->wsum_hess);
 }
 
 static void eval_wsum_hess(expr *node, const double *w)

src/affine/sum.c

@@ -136,12 +136,8 @@ static void wsum_hess_init(expr *node)
     x->wsum_hess_init(x);
 
     /* we never have to store more than the child's nnz */
-    node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, x->wsum_hess->nnz);
+    node->wsum_hess = new_csr_copy_sparsity(x->wsum_hess);
     node->work->dwork = malloc(x->size * sizeof(double));
-
-    /* copy sparsity pattern */
-    memcpy(node->wsum_hess->p, x->wsum_hess->p, (x->n_vars + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, x->wsum_hess->i, x->wsum_hess->nnz * sizeof(int));
 }
 
 static void eval_wsum_hess(expr *node, const double *w)

src/affine/trace.c

@@ -102,16 +102,14 @@ static void wsum_hess_init(expr *node)
 
     /* initialize child's hessian */
     x->wsum_hess_init(x);
-
-    node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, x->wsum_hess->nnz);
+  
     node->work->dwork = (double *) calloc(x->size, sizeof(double));
 
     /* We copy over the sparsity pattern from the child. This also includes the
        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. */
-    memcpy(node->wsum_hess->p, x->wsum_hess->p, (x->n_vars + 1) * sizeof(int));
-    memcpy(node->wsum_hess->i, x->wsum_hess->i, x->wsum_hess->nnz * sizeof(int));
+    node->wsum_hess = new_csr_copy_sparsity(x->wsum_hess);
 }
 
 static void eval_wsum_hess(expr *node, const double *w)

src/affine/transpose.c

@@ -92,11 +92,7 @@ static void wsum_hess_init(expr *node)
     x->wsum_hess_init(x);
 
     /* same sparsity pattern as child */
-    CSR_Matrix *H = node->wsum_hess;
-    H = new_csr_matrix(x->wsum_hess->m, node->n_vars, x->wsum_hess->nnz);
-    memcpy(H->p, x->wsum_hess->p, (H->m + 1) * sizeof(int));
-    memcpy(H->i, x->wsum_hess->i, H->nnz * sizeof(int));
-    node->wsum_hess = H;
+    node->wsum_hess = new_csr_copy_sparsity(x->wsum_hess);
 
     /* for computing Kw where K is the commutation matrix */
     node->work->dwork = (double *) malloc(node->size * sizeof(double));