clean up CSC_Matrix class

Author: dance858

include/utils/CSC_Matrix.h

@@ -29,46 +29,39 @@ CSC_Matrix *new_csc_matrix(int m, int n, int nnz);
 /* Free a CSC matrix */
 void free_csc_matrix(CSC_Matrix *matrix);
 
-CSC_Matrix *csr_to_csc(const CSR_Matrix *A);
-
-/* Allocate sparsity pattern for C = A^T D A for diagonal D */
+/* Fill sparsity of C = A^T D A for diagonal D */
 CSR_Matrix *ATA_alloc(const CSC_Matrix *A);
 
-/* Allocate sparsity pattern for C = B^T D A for diagonal D */
+/* Fill sparsity of C = B^T D A for diagonal D */
 CSR_Matrix *BTA_alloc(const CSC_Matrix *A, const CSC_Matrix *B);
 
-/* Compute values for C = A^T D A. C must have precomputed sparsity pattern.
- * If d is NULL, D is treated as the identity (computes A^T A). */
+/* Fill sparsity of C = BA, where B is symmetric. */
+CSC_Matrix *symBA_alloc(const CSR_Matrix *B, const CSC_Matrix *A);
+
+/* Compute values for C = A^T D A (null d corresponds to D as identity) */
 void ATDA_fill_values(const CSC_Matrix *A, const double *d, CSR_Matrix *C);
 
-/* Compute values for C = B^T D A. C must have precomputed sparsity pattern.
- * If d is NULL, D is treated as the identity (computes B^T A). */
+/* Compute values for C = B^T D A (null d corresonds to D as identity) */
 void BTDA_fill_values(const CSC_Matrix *A, const CSC_Matrix *B, const double *d,
                       CSR_Matrix *C);
 
-/* C = z^T A where A is in CSC format and C is assumed to have one row.
- * C must have column indices pre-computed. Fills in values of C only.
- */
-void csc_matvec_fill_values(const CSC_Matrix *A, const double *z, CSR_Matrix *C);
-
-/* Allocate B = Q * A (sparsity only). Q must be symmetric.
- * Q is CSR, A is CSC, B is CSC. */
-CSC_Matrix *sym_csr_csc_multiply_fill_sparsity(const CSR_Matrix *Q,
-                                               const CSC_Matrix *A);
+/* Fill values of C = BA. The matrix B does not have to be symmetric */
+void BA_fill_values(const CSR_Matrix *B, const CSC_Matrix *A, CSC_Matrix *C);
 
-/* Fill values of B = Q * A. Q must be symmetric.
- * B must have sparsity from sym_csr_csc_multiply_fill_sparsity. */
-void sym_csr_csc_multiply_fill_values(const CSR_Matrix *Q,
-                                      const CSC_Matrix *A,
-                                      CSC_Matrix *B);
+/* Fill values of C = x^T A. The matrix C must have filled sparsity. */
+void yTA_fill_values(const CSC_Matrix *A, const double *x, CSR_Matrix *C);
 
 /* Count nonzero columns of a CSC matrix */
 int count_nonzero_cols_csc(const CSC_Matrix *A);
 
-CSC_Matrix *csr_to_csc_fill_sparsity(const CSR_Matrix *A, int *iwork);
+/* convert from CSR to CSC format */
+CSC_Matrix *csr_to_csc_alloc(const CSR_Matrix *A, int *iwork);
 void csr_to_csc_fill_values(const CSR_Matrix *A, CSC_Matrix *C, int *iwork);
 
-CSR_Matrix *csc_to_csr_fill_sparsity(const CSC_Matrix *A, int *iwork);
+
+/* convert from CSC to CSR format */
+CSR_Matrix *csc_to_csr_alloc(const CSC_Matrix *A, int *iwork);
 void csc_to_csr_fill_values(const CSC_Matrix *A, CSR_Matrix *C, int *iwork);
 
+
 #endif /* CSC_MATRIX_H */

src/affine/left_matmul.c

@@ -88,12 +88,12 @@ static void jacobian_init_impl(expr *node)
 
     /* initialize child's jacobian and precompute sparsity of its CSC */
     jacobian_init(x);
-    lnode->Jchild_CSC = csr_to_csc_fill_sparsity(x->jacobian, node->work->iwork);
+    lnode->Jchild_CSC = csr_to_csc_alloc(x->jacobian, node->work->iwork);
 
     /* precompute sparsity of this node's jacobian in CSC and CSR */
     lnode->J_CSC = lnode->A->block_left_mult_sparsity(lnode->A, lnode->Jchild_CSC,
                                                       lnode->n_blocks);
-    node->jacobian = csc_to_csr_fill_sparsity(lnode->J_CSC, lnode->csc_to_csr_work);
+    node->jacobian = csc_to_csr_alloc(lnode->J_CSC, lnode->csc_to_csr_work);
 }
 
 static void eval_jacobian(expr *node)

src/bivariate_restricted_dom/quad_over_lin.c

@@ -164,8 +164,7 @@ static void eval_jacobian(expr *node)
         }
 
         /* chain rule (no derivative wrt y) using CSC format */
-        csc_matvec_fill_values(x->work->jacobian_csc, node->work->dwork,
-                               node->jacobian);
+        yTA_fill_values(x->work->jacobian_csc, node->work->dwork, node->jacobian);
 
         /* insert derivative wrt y at right place (for correctness this assumes
            that y does not appear in the numerator, but this will always be

src/expr.c

@@ -51,7 +51,7 @@ void jacobian_csc_init(expr *node)
     }
     node->work->csc_work = (int *) malloc(node->n_vars * sizeof(int));
     node->work->jacobian_csc =
-        csr_to_csc_fill_sparsity(node->jacobian, node->work->csc_work);
+        csr_to_csc_alloc(node->jacobian, node->work->csc_work);
 }
 
 void free_expr(expr *node)

src/other/quad_form.c

@@ -99,8 +99,7 @@ static void eval_jacobian(expr *node)
 
         /* The jacobian has same values as the gradient, which is
            J_f^T (Q @ f(x)). Here, dwork stores Q @ f(x) from forward */
-        csc_matvec_fill_values(x->work->jacobian_csc, node->work->dwork,
-                               node->jacobian);
+        yTA_fill_values(x->work->jacobian_csc, node->work->dwork, node->jacobian);
 
         cblas_dscal(node->jacobian->nnz, 2.0, node->jacobian->x, 1);
     }
@@ -146,7 +145,7 @@ static void wsum_hess_init_impl(expr *node)
         CSC_Matrix *Jf = x->work->jacobian_csc;
 
         /* term1 = Jf^T W Jf = Jf^T B*/
-        CSC_Matrix *B = sym_csr_csc_multiply_fill_sparsity(Q, Jf);
+        CSC_Matrix *B = symBA_alloc(Q, Jf);
         qnode->QJf = B;
         node->work->hess_term1 = BTA_alloc(Jf, B);
 
@@ -194,7 +193,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         CSR_Matrix *term2 = node->work->hess_term2;
 
         /* term1 = J_f^T Q J_f = J_f^T B  */
-        sym_csr_csc_multiply_fill_values(Q, Jf, QJf);
+        BA_fill_values(Q, Jf, QJf);
         BTDA_fill_values(Jf, QJf, NULL, term1);
 
         /* term2 */

src/utils/CSC_Matrix.c

@@ -202,54 +202,7 @@ void ATDA_fill_values(const CSC_Matrix *A, const double *d, CSR_Matrix *C)
     }
 }
 
-CSC_Matrix *csr_to_csc(const CSR_Matrix *A)
-{
-    CSC_Matrix *C = new_csc_matrix(A->m, A->n, A->nnz);
-
-    int i, j;
-    int *count = malloc(A->n * sizeof(int));
-
-    memset(count, 0, A->n * sizeof(int));
-
-    // -------------------------------------------------------------------
-    //              compute nnz in each column of A
-    // -------------------------------------------------------------------
-    for (i = 0; i < A->m; ++i)
-    {
-        for (j = A->p[i]; j < A->p[i + 1]; ++j)
-        {
-            count[A->i[j]]++;
-        }
-    }
-
-    // ------------------------------------------------------------------
-    //                      compute column pointers
-    // ------------------------------------------------------------------
-    C->p[0] = 0;
-    for (i = 0; i < A->n; ++i)
-    {
-        C->p[i + 1] = C->p[i] + count[i];
-        count[i] = C->p[i];
-    }
-
-    // ------------------------------------------------------------------
-    //                         fill matrix
-    // ------------------------------------------------------------------
-    for (i = 0; i < A->m; ++i)
-    {
-        for (j = A->p[i]; j < A->p[i + 1]; ++j)
-        {
-            C->x[count[A->i[j]]] = A->x[j];
-            C->i[count[A->i[j]]] = i;
-            count[A->i[j]]++;
-        }
-    }
-
-    free(count);
-    return C;
-}
-
-CSC_Matrix *csr_to_csc_fill_sparsity(const CSR_Matrix *A, int *iwork)
+CSC_Matrix *csr_to_csc_alloc(const CSR_Matrix *A, int *iwork)
 {
     CSC_Matrix *C = new_csc_matrix(A->m, A->n, A->nnz);
 
@@ -312,7 +265,7 @@ void csr_to_csc_fill_values(const CSR_Matrix *A, CSC_Matrix *C, int *iwork)
     }
 }
 
-CSR_Matrix *csc_to_csr_fill_sparsity(const CSC_Matrix *A, int *iwork)
+CSR_Matrix *csc_to_csr_alloc(const CSC_Matrix *A, int *iwork)
 {
     CSR_Matrix *C = new_csr_matrix(A->m, A->n, A->nnz);
 
@@ -435,19 +388,14 @@ CSR_Matrix *BTA_alloc(const CSC_Matrix *A, const CSC_Matrix *B)
     return C;
 }
 
-void csc_matvec_fill_values(const CSC_Matrix *A, const double *z, CSR_Matrix *C)
+void yTA_fill_values(const CSC_Matrix *A, const double *y, CSR_Matrix *C)
 {
-    /* Compute C = z^T * A where A is in CSC format
-     * C is a single-row CSR matrix with column indices pre-computed
-     * This fills in the values of C only
-     */
-
     for (int col = 0; col < A->n; col++)
     {
         double val = 0;
         for (int j = A->p[col]; j < A->p[col + 1]; j++)
         {
-            val += z[A->i[j]] * A->x[j];
+            val += y[A->i[j]] * A->x[j];
         }
 
         if (A->p[col + 1] - A->p[col] == 0) continue;
@@ -464,6 +412,7 @@ void csc_matvec_fill_values(const CSC_Matrix *A, const double *z, CSR_Matrix *C)
     }
 }
 
+/* computes C = B^T * D * A in CSR */
 void BTDA_fill_values(const CSC_Matrix *A, const CSC_Matrix *B, const double *d,
                       CSR_Matrix *C)
 {
@@ -477,40 +426,62 @@ void BTDA_fill_values(const CSC_Matrix *A, const CSC_Matrix *B, const double *d,
             int nnz_bi = B->p[i + 1] - B->p[i];
             int nnz_aj = A->p[j + 1] - A->p[j];
 
-            double sum;
             if (d != NULL)
             {
-                sum = sparse_wdot(B->x + B->p[i], B->i + B->p[i], nnz_bi,
-                                  A->x + A->p[j], A->i + A->p[j], nnz_aj, d);
+                C->x[jj] = sparse_wdot(B->x + B->p[i], B->i + B->p[i], nnz_bi,
+                                       A->x + A->p[j], A->i + A->p[j], nnz_aj, d);
             }
             else
             {
-                sum = sparse_dot(B->x + B->p[i], B->i + B->p[i], nnz_bi,
-                                 A->x + A->p[j], A->i + A->p[j], nnz_aj);
+                C->x[jj] = sparse_dot(B->x + B->p[i], B->i + B->p[i], nnz_bi,
+                                      A->x + A->p[j], A->i + A->p[j], nnz_aj);
             }
+        }
+    }
+}
+
+/* NOTE: an alternative marker-based approach (scatter A_{k,j} * Q[k,:]
+ * into column j of C using a marker array for position lookup) may be
+ * faster when Q is dense, since it touches each Q entry exactly once.
+ * The sparse_dot approach below is simpler but redundantly scans
+ * column j of A for each nonzero row of C. */
+void BA_fill_values(const CSR_Matrix *Q, const CSC_Matrix *A, CSC_Matrix *C)
+{
+    /* fill values of C = Q * A, given the sparsity pattern of C. */
+    int i, j, ii;
 
-            C->x[jj] = sum;
+    /* for each column j of C */
+    for (j = 0; j < C->n; j++)
+    {
+        for (ii = C->p[j]; ii < C->p[j + 1]; ii++)
+        {
+            i = C->i[ii];
+            int nnz_q = Q->p[i + 1] - Q->p[i];
+            int nnz_a = A->p[j + 1] - A->p[j];
+
+            /* inner product between row i of Q and column j of A */
+            C->x[ii] = sparse_dot(Q->x + Q->p[i], Q->i + Q->p[i], nnz_q,
+                                  A->x + A->p[j], A->i + A->p[j], nnz_a);
         }
     }
 }
 
-CSC_Matrix *sym_csr_csc_multiply_fill_sparsity(const CSR_Matrix *Q,
-                                               const CSC_Matrix *A)
+CSC_Matrix *symBA_alloc(const CSR_Matrix *B, const CSC_Matrix *A)
 {
-    /* Allocate B = Q * A (sparsity only). Q must be symmetric.
-     * Q is CSR (m x m), A is CSC (m x n), B is CSC (m x n).
+    /* Allocate C = B * A (sparsity only). B must be symmetric.
+     * B is CSR (m x m), A is CSC (m x n), C is CSC (m x n).
      *
-     * Column j of B is Q * a_j = sum_k A_{k,j} Q[:, k], so the nonzero
-     * rows of column j of B are the union of the nonzero rows of Q[:, k].
+     * Column j of C is B * a_j = sum_k A_{k,j} B[:, k], so the nonzero
+     * rows of column j of C are the union of the nonzero rows of B[:, k].
      *
-     * Since Q is symmetric, we can find the nonzero rows of Q[:, k] by
-     * finding the nonzero columns of Q in row k.
+     * Since B is symmetric, we can find the nonzero rows of B[:, k] by
+     * finding the nonzero columns of B in row k.
      *
      * We use a marker array to avoid duplicates: marker[l] stores the
      * last column j that registered l as nonzero, so checking
      * marker[l] != j avoids duplicates. */
 
-    int m = Q->m;
+    int m = B->m;
     int n = A->n;
     int i, j, k, jj, ii, ell;
 
@@ -521,11 +492,11 @@ CSC_Matrix *sym_csr_csc_multiply_fill_sparsity(const CSR_Matrix *Q,
         marker[i] = -1;
     }
 
-    int *Bp = (int *) malloc((n + 1) * sizeof(int));
-    iVec *Bi = iVec_new(A->nnz);
-    Bp[0] = 0;
+    int *Cp = (int *) malloc((n + 1) * sizeof(int));
+    iVec *Ci = iVec_new(A->nnz);
+    Cp[0] = 0;
 
-    /* for each column j of B */
+    /* for each column j of C */
     for (j = 0; j < n; j++)
     {
         int col_nnz = 0;
@@ -535,66 +506,33 @@ CSC_Matrix *sym_csr_csc_multiply_fill_sparsity(const CSR_Matrix *Q,
         {
             k = A->i[ii];
 
-            /* find nonzero rows ell of column k of Q */
-            for (jj = Q->p[k]; jj < Q->p[k + 1]; jj++)
+            /* find nonzero rows ell of column k of B */
+            for (jj = B->p[k]; jj < B->p[k + 1]; jj++)
             {
-                ell = Q->i[jj];
+                ell = B->i[jj];
                 if (marker[ell] != j)
                 {
                     marker[ell] = j;
-                    iVec_append(Bi, ell);
+                    iVec_append(Ci, ell);
                     col_nnz++;
                 }
             }
         }
 
-        Bp[j + 1] = Bp[j] + col_nnz;
+        Cp[j + 1] = Cp[j] + col_nnz;
     }
 
-    /* allocate B and copy the computed structure */
-    int total_nnz = Bp[n];
-    CSC_Matrix *B = new_csc_matrix(m, n, total_nnz);
-    memcpy(B->p, Bp, (n + 1) * sizeof(int));
-    memcpy(B->i, Bi->data, total_nnz * sizeof(int));
+    /* allocate C and copy the computed structure */
+    int total_nnz = Cp[n];
+    CSC_Matrix *C = new_csc_matrix(m, n, total_nnz);
+    memcpy(C->p, Cp, (n + 1) * sizeof(int));
+    memcpy(C->i, Ci->data, total_nnz * sizeof(int));
 
     free(marker);
-    free(Bp);
-    iVec_free(Bi);
-
-    return B;
-}
-
-/* NOTE: an alternative marker-based approach (scatter A_{k,j} * Q[k,:]
- * into column j of B using a marker array for position lookup) may be
- * faster when Q is dense, since it touches each Q entry exactly once.
- * The sparse_dot approach below is simpler but redundantly scans
- * column j of A for each nonzero row of B. */
-void sym_csr_csc_multiply_fill_values(const CSR_Matrix *Q, const CSC_Matrix *A,
-                                      CSC_Matrix *B)
-{
-    /* Fill values of B = Q * A. Q must be symmetric.
-     * B must have sparsity from sym_csr_csc_multiply_fill_sparsity.
-     *
-     * B_{l,j} = sum_k Q_{l,k} * A_{k,j} = dot(Q[l,:], A[:,j]).
-     * Since Q is symmetric, row l of Q has the same entries as
-     * column l, so we iterate over row l of Q in CSR format. */
-
-    int i, j, ii;
-
-    /* for each column j of B */
-    for (j = 0; j < B->n; j++)
-    {
-        for (ii = B->p[j]; ii < B->p[j + 1]; ii++)
-        {
-            i = B->i[ii];
-            int nnz_q = Q->p[i + 1] - Q->p[i];
-            int nnz_a = A->p[j + 1] - A->p[j];
+    free(Cp);
+    iVec_free(Ci);
 
-            /* inner product between row i of Q and column j of A */
-            B->x[ii] = sparse_dot(Q->x + Q->p[i], Q->i + Q->p[i], nnz_q,
-                                  A->x + A->p[j], A->i + A->p[j], nnz_a);
-        }
-    }
+    return C;
 }
 
 int count_nonzero_cols_csc(const CSC_Matrix *A)