move old code into old-code-folder

include/affine.h → include/atoms/affine.h

@@ -22,8 +22,6 @@
 #include "subexpr.h"
 #include "utils/CSR_Matrix.h"
 
-expr *new_linear(expr *u, const CSR_Matrix *A, const double *b);
-
 expr *new_add(expr *left, expr *right);
 expr *new_neg(expr *child);
 

include/other.h → include/atoms/non_elementwise_full_dom.h

@@ -15,8 +15,8 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-#ifndef OTHER_H
-#define OTHER_H
+#ifndef NON_ELEMENTWISE_FULL_DOM_H
+#define NON_ELEMENTWISE_FULL_DOM_H
 
 #include "expr.h"
 #include "subexpr.h"
@@ -33,4 +33,4 @@ expr *new_prod_axis_zero(expr *child);
 /* product of entries along axis=1 (rowwise products) */
 expr *new_prod_axis_one(expr *child);
 
-#endif /* OTHER_H */
+#endif /* NON_ELEMENTWISE_FULL_DOM_H */

include/old-code/old_CSR.h

@@ -0,0 +1,29 @@
+#ifndef OLD_CSR_H
+#define OLD_CSR_H
+
+#include "utils/CSR_Matrix.h"
+
+/* Build (I_p kron A) = blkdiag(A, A, ..., A) of size (p*A->m) x (p*A->n) */
+CSR_Matrix *block_diag_repeat_csr(const CSR_Matrix *A, int p);
+
+/* Build (A kron I_p) of size (A->m * p) x (A->n * p) with nnz = A->nnz * p. */
+CSR_Matrix *kron_identity_csr(const CSR_Matrix *A, int p);
+
+/* Computes values of the row matrix C = z^T A (column indices must have been
+   pre-computed) and transposed matrix AT must be provided) */
+void Ax_csr_fill_values(const CSR_Matrix *AT, const double *z, CSR_Matrix *C);
+
+/* Insert value into CSR matrix A with just one row at col_idx. Assumes that A
+has enough space and that A does not have an element at col_idx. It does update
+nnz. */
+void csr_insert_value(CSR_Matrix *A, int col_idx, double value);
+
+/* Compute C = diag(d) * A where d is an array and A, C are CSR matrices
+ * d must have length m
+ * C must be pre-allocated with same dimensions as A */
+void diag_csr_mult(const double *d, const CSR_Matrix *A, CSR_Matrix *C);
+
+/* y = Ax, where y is returned as dense (no column offset) */
+void Ax_csr_wo_offset(const CSR_Matrix *A, const double *x, double *y);
+
+#endif /* OLD_CSR_H */

include/old-code/old_CSR_sum.h

@@ -0,0 +1,44 @@
+#ifndef OLD_CSR_SUM_H
+#define OLD_CSR_SUM_H
+
+#include "utils/CSR_Matrix.h"
+
+/* Compute C = A + B where A, B, C are CSR matrices
+ * A and B must have same dimensions
+ * C must be pre-allocated with sufficient nnz capacity.
+ * C must be different from A and B */
+void sum_csr_matrices(const CSR_Matrix *A, const CSR_Matrix *B, CSR_Matrix *C);
+
+/* Compute C = diag(d1) * A + diag(d2) * B where A, B, C are CSR matrices */
+void sum_scaled_csr_matrices(const CSR_Matrix *A, const CSR_Matrix *B, CSR_Matrix *C,
+                             const double *d1, const double *d2);
+
+/* forward declaration */
+struct int_double_pair;
+
+/* Sum all rows of A into a single row matrix C */
+void sum_all_rows_csr(const CSR_Matrix *A, CSR_Matrix *C,
+                      struct int_double_pair *pairs);
+
+/* Sum blocks of rows of A into a matrix C */
+void sum_block_of_rows_csr(const CSR_Matrix *A, CSR_Matrix *C,
+                           struct int_double_pair *pairs, int row_block_size);
+
+/* Sum evenly spaced rows of A into a matrix C */
+void sum_evenly_spaced_rows_csr(const CSR_Matrix *A, CSR_Matrix *C,
+                                struct int_double_pair *pairs, int row_spacing);
+
+/* Sum evenly spaced rows of A starting at offset into a row matrix C */
+void sum_spaced_rows_into_row_csr(const CSR_Matrix *A, CSR_Matrix *C,
+                                  struct int_double_pair *pairs, int offset,
+                                  int spacing);
+
+/* Fill values of summed rows using precomputed idx_map and sparsity of C */
+void sum_all_rows_csr_fill_values(const CSR_Matrix *A, CSR_Matrix *C,
+                                  const int *idx_map);
+
+/* Fill values of summed block rows using precomputed idx_map */
+void sum_block_of_rows_csr_fill_values(const CSR_Matrix *A, CSR_Matrix *C,
+                                       const int *idx_map);
+
+#endif /* OLD_CSR_SUM_H */

include/old-code/old_affine.h

@@ -0,0 +1,9 @@
+#ifndef OLD_AFFINE_H
+#define OLD_AFFINE_H
+
+#include "expr.h"
+#include "utils/CSR_Matrix.h"
+
+expr *new_linear(expr *u, const CSR_Matrix *A, const double *b);
+
+#endif /* OLD_AFFINE_H */

include/utils/CSC_Matrix.h

@@ -23,10 +23,8 @@ typedef struct CSC_Matrix
     int nnz;
 } CSC_Matrix;
 
-/* Allocate a new CSC matrix with given dimensions and nnz */
+/* constructor and destructor */
 CSC_Matrix *new_csc_matrix(int m, int n, int nnz);
-
-/* Free a CSC matrix */
 void free_csc_matrix(CSC_Matrix *matrix);
 
 /* Fill sparsity of C = A^T D A for diagonal D */

include/utils/CSR_Matrix.h

@@ -2,9 +2,6 @@
 #define CSR_MATRIX_H
 #include <stdbool.h>
 
-/* forward declaration */
-struct int_double_pair;
-
 /* CSR (Compressed Sparse Row) Matrix Format
  *
  * For an m x n matrix with nnz nonzeros:
@@ -37,37 +34,19 @@ CSR_Matrix *transpose(const CSR_Matrix *A, int *iwork);
 CSR_Matrix *AT_alloc(const CSR_Matrix *A, int *iwork);
 void AT_fill_values(const CSR_Matrix *A, CSR_Matrix *AT, int *iwork);
 
-/* Build (I_p kron A) = blkdiag(A, A, ..., A) of size (p*A->m) x (p*A->n) */
-CSR_Matrix *block_diag_repeat_csr(const CSR_Matrix *A, int p);
-
-/* Build (A kron I_p) of size (A->m * p) x (A->n * p) with nnz = A->nnz * p. */
-CSR_Matrix *kron_identity_csr(const CSR_Matrix *A, int p);
-
-/* y = Ax, where y is returned as dense */
-void csr_matvec(const CSR_Matrix *A, const double *x, double *y, int col_offset);
-void csr_matvec_wo_offset(const CSR_Matrix *A, const double *x, double *y);
-
-/* Computes values of the row matrix C = z^T A (column indices must have been
-   pre-computed) and transposed matrix AT must be provided) */
-void csr_matvec_fill_values(const CSR_Matrix *AT, const double *z, CSR_Matrix *C);
-
-/* Insert value into CSR matrix A with just one row at col_idx. Assumes that A
-has enough space and that A does not have an element at col_idx. It does update
-nnz. */
-void csr_insert_value(CSR_Matrix *A, int col_idx, double value);
+/* computes dense y = Ax */
+void Ax_csr(const CSR_Matrix *A, const double *x, double *y, int col_offset);
 
-/* Compute C = diag(d) * A where d is an array and A, C are CSR matrices
- * d must have length m
- * C must be pre-allocated with same dimensions as A */
-void diag_csr_mult(const double *d, const CSR_Matrix *A, CSR_Matrix *C);
-void diag_csr_mult_fill_values(const double *d, const CSR_Matrix *A, CSR_Matrix *C);
+/* fills values of C = diag(d) @ A */
+void DA_fill_values(const double *d, const CSR_Matrix *A, CSR_Matrix *C);
 
-/* Count number of columns with nonzero entries */
+/* Count number of columns with nonzero entries in A and marks them in col_nz */
 int count_nonzero_cols(const CSR_Matrix *A, bool *col_nz);
 
 /* inserts 'idx' into array 'arr' in sorted order, and moves the other elements */
 void insert_idx(int idx, int *arr, int len);
 
+/* get value at position (row, col) in A */
 double csr_get_value(const CSR_Matrix *A, int row, int col);
 
 /* Expand symmetric CSR matrix A to full matrix C. A is assumed to store

include/utils/CSR_sum.h

@@ -6,93 +6,55 @@
 /* forward declaration */
 struct int_double_pair;
 
-/* Compute C = A + B where A, B, C are CSR matrices
- * A and B must have same dimensions
- * C must be pre-allocated with sufficient nnz capacity.
- * C must be different from A and B */
-void sum_csr_matrices(const CSR_Matrix *A, const CSR_Matrix *B, CSR_Matrix *C);
-
-/* Compute sparsity pattern of A + B where A, B, C are CSR matrices.
- * Fills C->p, C->i, and C->nnz; does not touch C->x. */
+/* Compute sparsity pattern of C = A + B (and sets C->nnz) */
 void sum_csr_alloc(const CSR_Matrix *A, const CSR_Matrix *B, CSR_Matrix *C);
 
-/* Fill only the values of C = A + B, assuming C's sparsity pattern (p and i)
- * is already filled and matches the union of A and B per row. Does not modify
- * C->p, C->i, or C->nnz. */
+/* Fills values of C = A + B (assuming C's sparsity pattern is set) */
 void sum_csr_fill_values(const CSR_Matrix *A, const CSR_Matrix *B, CSR_Matrix *C);
 
-/* Compute C = diag(d1) * A + diag(d2) * B where A, B, C are CSR matrices */
-void sum_scaled_csr_matrices(const CSR_Matrix *A, const CSR_Matrix *B, CSR_Matrix *C,
-                             const double *d1, const double *d2);
-
-/* Fill only the values of C = diag(d1) * A + diag(d2) * B, assuming C's sparsity
- * pattern (p and i) is already filled and matches the union of A and B per row.
- * Does not modify C->p, C->i, or C->nnz. */
+/* Fills values of C = diag(d1) * A + diag(d2) * B (assuming C's sparsity is set)*/
 void sum_scaled_csr_matrices_fill_values(const CSR_Matrix *A, const CSR_Matrix *B,
                                          CSR_Matrix *C, const double *d1,
                                          const double *d2);
 
-/* Sum all rows of A into a single row matrix C */
-void sum_all_rows_csr(const CSR_Matrix *A, CSR_Matrix *C,
-                      struct int_double_pair *pairs);
-
-/* iwork must have size max(C->n, A->nnz), and idx_map must have size A->nnz. */
-void sum_all_rows_csr_fill_sparsity_and_idx_map(const CSR_Matrix *A, CSR_Matrix *C,
-                                                int *iwork, int *idx_map);
-
-/* Fill values of summed rows using precomputed idx_map and sparsity of C */
-// void sum_all_rows_csr_fill_values(const CSR_Matrix *A, CSR_Matrix *C,
-//                                  const int *idx_map);
-
-/* Fill accumulator for summing rows using precomputed idx_map for each nnz of A.
-   Must memset accumulator to zero before calling. */
-void idx_map_accumulator(const CSR_Matrix *A, const int *idx_map,
-                         double *accumulator);
-void idx_map_accumulator_with_spacing(const CSR_Matrix *A, const int *idx_map,
-                                      double *accumulator, int spacing);
-
-/* Sum blocks of rows of A into a matrix C */
-void sum_block_of_rows_csr(const CSR_Matrix *A, CSR_Matrix *C,
-                           struct int_double_pair *pairs, int row_block_size);
-
-/* Build sparsity and index map for summing blocks of rows.
- * iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz. */
-void sum_block_of_rows_csr_fill_sparsity_and_idx_map(const CSR_Matrix *A,
-                                                     CSR_Matrix *C,
-                                                     int row_block_size, int *iwork,
-                                                     int *idx_map);
-
-/* Sum evenly spaced rows of A into a matrix C */
-void sum_evenly_spaced_rows_csr(const CSR_Matrix *A, CSR_Matrix *C,
-                                struct int_double_pair *pairs, int row_spacing);
-
-/* Build sparsity and index map for summing evenly spaced rows.
- * iwork must have size max(A->n, A->nnz), and idx_map must have size A->nnz. */
-void sum_evenly_spaced_rows_csr_fill_sparsity_and_idx_map(const CSR_Matrix *A,
-                                                          CSR_Matrix *C,
-                                                          int row_spacing,
-                                                          int *iwork, int *idx_map);
-
-/* Sum evenly spaced rows of A starting at offset into a row matrix C */
-void sum_spaced_rows_into_row_csr(const CSR_Matrix *A, CSR_Matrix *C,
-                                  struct int_double_pair *pairs, int offset,
-                                  int spacing);
-
-/* Fills the sparsity and index map for summing spaced rows into a row matrix */
-void sum_spaced_rows_into_row_csr_fill_sparsity_and_idx_map(const CSR_Matrix *A,
-                                                            CSR_Matrix *C,
-                                                            int spacing, int *iwork,
-                                                            int *idx_map);
-
-/* 4-way sorted merge of CSR matrices A, B, C, D (same dimensions).
- * Allocates and returns the output CSR with the union sparsity pattern.
- * Allocates and fills idx_maps[0..3] (one per input, size input->nnz
- * each) mapping each input entry to its position in the output.
- * Caller owns the returned CSR and all 4 idx_map arrays. */
-CSR_Matrix *sum_4_csr_fill_sparsity_and_idx_maps(const CSR_Matrix *A,
-                                                 const CSR_Matrix *B,
-                                                 const CSR_Matrix *C,
-                                                 const CSR_Matrix *D,
-                                                 int *idx_maps[4]);
+/* The following five functions are used for summing either more than two CSR
+   matrices or rows of CSR 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 output
+   matrix. This allows us to fill the values with a single pass of the output matrix
+   through the input matrices, without needing to search for the position of each
+   entry in the output matrix. So each idx_map should have size equal to the number
+   of nonzeros in the corresponding input matrix, and idx_map[j] should give the
+   index in the output matrix of the entry (in the value array of the output 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. */
+// ------------------------------------------------------------------------------------
+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);
+
+/* Compute sparsity pattern of out = A + B + C + D */
+CSR_Matrix *sum_4_csr_alloc(const CSR_Matrix *A, const CSR_Matrix *B,
+                            const CSR_Matrix *C, const CSR_Matrix *D,
+                            int *idx_maps[4]);
+// ------------------------------------------------------------------------------------
+
+/* Accumulates values from A according to map. Must memset to zero before calling. */
+void accumulator(const CSR_Matrix *A, const int *idx_map, double *out);
+
+/* Accumulates values from A according to map with spacing. Must memset to zero
+ * before calling. */
+void accumulator_with_spacing(const CSR_Matrix *A, const int *idx_map, double *out,
+                              int spacing);
 
 #endif /* CSR_SUM_H */

include/utils/mini_numpy.h

@@ -1,36 +1,27 @@
 #ifndef MINI_NUMPY_H
 #define MINI_NUMPY_H
 
-/* Repeat each element of array 'a' 'repeats' times
- * Example: a = [1, 2], len = 2, repeats = 3
- *          result = [1, 1, 1, 2, 2, 2]
- */
+/* Example: a = [1, 2], len = 2, repeats = 3, result = [1, 1, 1, 2, 2, 2] */
 void repeat(double *result, const double *a, int len, int repeats);
 
-/* Tile array 'a' 'tiles' times
- * Example: a = [1, 2], len = 2, tiles = 3
- *          result = [1, 2, 1, 2, 1, 2]
- */
+/* Example: a = [1, 2], len = 2, tiles = 3, result = [1, 2, 1, 2, 1, 2] */
 void tile_double(double *result, const double *a, int len, int tiles);
 void tile_int(int *result, const int *a, int len, int tiles);
 
-/* Fill array with 'size' copies of 'value'
- * Example: size = 5, value = 3.0
- *          result = [3.0, 3.0, 3.0, 3.0, 3.0]
- */
+/* Example: size = 5, value = 3.0, result = [3.0, 3.0, 3.0, 3.0, 3.0] */
 void scaled_ones(double *result, int size, double value);
 
 /* Naive implementation of Z = X @ Y, X is m x k, Y is k x n, Z is m x n */
 void mat_mat_mult(const double *X, const double *Y, double *Z, int m, int k, int n);
 
-/* Compute v = (Y kron I_m) @ w where Y is k x n (col-major), w has
-   length m*n, and v has length m*k.  Equivalently, reshape w as the
-   m x n matrix W (col-major) and compute v = vec(W @ Y^T). */
+/* Compute v = (Y kron I_m) @ w where Y is k x n (col-major), len(w) = m * n, and
+   len(v) = m * k.  Equivalently, reshape w as the m x n matrix W (col-major) and
+   compute v = vec(W @ Y^T). */
 void Y_kron_I_vec(int m, int k, int n, const double *Y, const double *w, double *v);
 
-/* Compute v = (I_n kron X^T) @ w where X is m x k (col-major), w has
-   length m*n, and v has length k*n.  Equivalently, reshape w as the
-   m x n matrix W (col-major) and compute v = vec(X^T @ W). */
+/* Compute v = (I_n kron X^T) @ w where X is m x k (col-major), len(w) = m * n, and
+   len(v) = k * n.  Equivalently, reshape w as the m x n matrix W (col-major) and
+   compute v = vec(X^T @ W). */
 void I_kron_XT_vec(int m, int k, int n, const double *X, const double *w, double *v);
 
 #endif /* MINI_NUMPY_H */

src/affine/add.c → src/atoms/affine/add.c

@@ -15,7 +15,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-#include "affine.h"
+#include "atoms/affine.h"
 #include "utils/CSR_sum.h"
 #include <assert.h>
 #include <stdio.h>