Skip to content

Commit eba6858

Browse files
dance858dance858
committed
A^T A alloc function, hessian infrastructure preparation
1 parent c89efa3 commit eba6858

8 files changed

Lines changed: 454 additions & 0 deletions

File tree

‎include/utils/CSC_Matrix.h‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
#ifndef CSC_MATRIX_H
2+
#define CSC_MATRIX_H
3+
4+
#include "CSR_Matrix.h"
5+
6+
/* CSC (Compressed Sparse Column) Matrix Format
7+
*
8+
* For an m x n matrix with nnz nonzeros:
9+
* - p: array of size (n + 1) indicating start of each column
10+
* - i: array of size nnz containing row indices
11+
* - x: array of size nnz containing values
12+
* - m: number of rows
13+
* - n: number of columns
14+
* - nnz: number of nonzero entries
15+
*/
16+
typedef struct CSC_Matrix
17+
{
18+
int *p;
19+
int *i;
20+
double *x;
21+
int m;
22+
int n;
23+
int nnz;
24+
} CSC_Matrix;
25+
26+
/* Allocate a new CSC matrix with given dimensions and nnz */
27+
CSC_Matrix *new_csc_matrix(int m, int n, int nnz);
28+
29+
/* Free a CSC matrix */
30+
void free_csc_matrix(CSC_Matrix *matrix);
31+
32+
/* Allocate sparsity pattern for C = A^T D A or C = A^T A
33+
*/
34+
CSR_Matrix *ATA_alloc(const CSC_Matrix *A);
35+
36+
#endif /* CSC_MATRIX_H */

‎include/utils/CSR_Matrix.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,4 +82,8 @@ void insert_idx(int idx, int *arr, int len);
8282

8383
CSR_Matrix *transpose(const CSR_Matrix *A, int *iwork);
8484

85+
/* Expand symmetric CSR matrix A to full matrix C. A is assumed to store
86+
only upper triangle. C must be pre-allocated with sufficient nnz */
87+
void symmetrize_csr(const int *Ap, const int *Ai, int m, CSR_Matrix *C);
88+
8589
#endif /* CSR_MATRIX_H */

‎include/utils/Vec_macros.h‎

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
/*
2+
* Copyright 2025 Daniel Cederberg
3+
*
4+
* This file is part of the PSLP project (LP Presolver).
5+
*
6+
* Licensed under the Apache License, Version 2.0 (the "License");
7+
* you may not use this file except in compliance with the License.
8+
* You may obtain a copy of the License at
9+
*
10+
* http://www.apache.org/licenses/LICENSE-2.0
11+
*
12+
* Unless required by applicable law or agreed to in writing, software
13+
* distributed under the License is distributed on an "AS IS" BASIS,
14+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15+
* See the License for the specific language governing permissions and
16+
* limitations under the License.
17+
*/
18+
19+
#ifndef VEC_MACROS_H
20+
#define VEC_MACROS_H
21+
22+
#include <assert.h>
23+
#include <stdio.h>
24+
#include <stdlib.h>
25+
#include <string.h>
26+
27+
// Macro to define a generic vector class
28+
#define DEFINE_VECTOR(TYPE, TYPE_NAME) \
29+
typedef struct TYPE_NAME##Vec \
30+
{ \
31+
TYPE *data; \
32+
size_t len; \
33+
size_t capacity; \
34+
} TYPE_NAME##Vec; \
35+
\
36+
__attribute__((unused)) static TYPE_NAME##Vec *TYPE_NAME##Vec_new( \
37+
size_t capacity) \
38+
{ \
39+
assert(capacity > 0); \
40+
TYPE_NAME##Vec *vec = (TYPE_NAME##Vec *) malloc(sizeof(TYPE_NAME##Vec)); \
41+
if (vec == NULL) return NULL; \
42+
vec->data = (TYPE *) malloc(capacity * sizeof(TYPE)); \
43+
if (vec->data == NULL) \
44+
{ \
45+
free(vec); \
46+
return NULL; \
47+
} \
48+
\
49+
vec->len = 0; \
50+
vec->capacity = capacity; \
51+
return vec; \
52+
} \
53+
\
54+
static inline void TYPE_NAME##Vec_free(TYPE_NAME##Vec *vec) \
55+
{ \
56+
free(vec->data); \
57+
free(vec); \
58+
} \
59+
\
60+
static inline void TYPE_NAME##Vec_clear_no_resize(TYPE_NAME##Vec *vec) \
61+
{ \
62+
vec->len = 0; \
63+
} \
64+
\
65+
static inline void TYPE_NAME##Vec_append(TYPE_NAME##Vec *vec, TYPE value) \
66+
{ \
67+
if (vec->len >= vec->capacity) \
68+
{ \
69+
vec->capacity *= 2; \
70+
assert(vec->capacity > 0); \
71+
TYPE *temp = (TYPE *) realloc(vec->data, vec->capacity * sizeof(TYPE)); \
72+
if (temp == NULL) \
73+
{ \
74+
TYPE_NAME##Vec_free(vec); \
75+
fprintf(stderr, "Error: realloc failed\n"); \
76+
exit(1); \
77+
} \
78+
\
79+
vec->data = temp; \
80+
} \
81+
\
82+
vec->data[vec->len++] = value; \
83+
} \
84+
\
85+
static inline void TYPE_NAME##Vec_append_array(TYPE_NAME##Vec *vec, \
86+
const TYPE *values, size_t n) \
87+
{ \
88+
if (vec->len + n > vec->capacity) \
89+
{ \
90+
size_t new_capacity = vec->capacity > 0 ? vec->capacity : 1; \
91+
while (vec->len + n > new_capacity) \
92+
{ \
93+
new_capacity *= 2; \
94+
} \
95+
\
96+
TYPE *temp = (TYPE *) realloc(vec->data, new_capacity * sizeof(TYPE)); \
97+
if (temp == NULL) \
98+
{ \
99+
TYPE_NAME##Vec_free(vec); \
100+
fprintf(stderr, "Error: realloc failed\n"); \
101+
exit(1); \
102+
} \
103+
\
104+
vec->data = temp; \
105+
vec->capacity = new_capacity; \
106+
} \
107+
\
108+
memcpy(vec->data + vec->len, values, n * sizeof(TYPE)); \
109+
vec->len += n; \
110+
} \
111+
__attribute__((unused)) static int TYPE_NAME##Vec_contains( \
112+
const TYPE_NAME##Vec *vec, TYPE value) \
113+
{ \
114+
for (size_t i = 0; i < vec->len; ++i) \
115+
{ \
116+
if (vec->data[i] == value) \
117+
{ \
118+
return 1; /* Element found */ \
119+
} \
120+
} \
121+
return 0; /* Element not found */ \
122+
}
123+
124+
#endif

‎include/utils/iVec.h‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
/*
2+
* Copyright 2025 Daniel Cederberg
3+
*
4+
* This file is part of the PSLP project (LP Presolver).
5+
*
6+
* Licensed under the Apache License, Version 2.0 (the "License");
7+
* you may not use this file except in compliance with the License.
8+
* You may obtain a copy of the License at
9+
*
10+
* http://www.apache.org/licenses/LICENSE-2.0
11+
*
12+
* Unless required by applicable law or agreed to in writing, software
13+
* distributed under the License is distributed on an "AS IS" BASIS,
14+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15+
* See the License for the specific language governing permissions and
16+
* limitations under the License.
17+
*/
18+
19+
#ifndef IVEC_H
20+
#define IVEC_H
21+
22+
#include <stdio.h>
23+
24+
#include "Vec_macros.h"
25+
26+
// This macro defines a vector of integers.
27+
DEFINE_VECTOR(int, i)
28+
29+
/*
30+
__attribute__((unused)) static void print_ivec(iVec *vec)
31+
{
32+
for (size_t i = 0; i < vec->len; ++i)
33+
{
34+
printf("%d ", vec->data[i]);
35+
}
36+
printf("\n");
37+
}
38+
*/
39+
#endif // IVEC_H

‎src/utils/CSC_Matrix.c‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
#include "utils/CSC_Matrix.h"
2+
#include "utils/iVec.h"
3+
#include <stdlib.h>
4+
#include <string.h>
5+
6+
CSC_Matrix *new_csc_matrix(int m, int n, int nnz)
7+
{
8+
CSC_Matrix *matrix = (CSC_Matrix *) malloc(sizeof(CSC_Matrix));
9+
if (!matrix) return NULL;
10+
11+
matrix->p = (int *) malloc((n + 1) * sizeof(int));
12+
matrix->i = (int *) malloc(nnz * sizeof(int));
13+
matrix->x = (double *) malloc(nnz * sizeof(double));
14+
15+
if (!matrix->p || !matrix->i || !matrix->x)
16+
{
17+
free(matrix->p);
18+
free(matrix->i);
19+
free(matrix->x);
20+
free(matrix);
21+
return NULL;
22+
}
23+
24+
matrix->m = m;
25+
matrix->n = n;
26+
matrix->nnz = nnz;
27+
28+
return matrix;
29+
}
30+
31+
void free_csc_matrix(CSC_Matrix *matrix)
32+
{
33+
if (matrix)
34+
{
35+
free(matrix->p);
36+
free(matrix->i);
37+
free(matrix->x);
38+
free(matrix);
39+
}
40+
}
41+
42+
CSR_Matrix *ATA_alloc(const CSC_Matrix *A)
43+
{
44+
/* A is m x n, A^T A is n x n */
45+
int n = A->n;
46+
int m = A->m;
47+
int nnz = 0;
48+
int i, j, ii, jj;
49+
50+
/* row ptr and column idxs for upper triangular part of C = A^T A */
51+
int *Cp = (int *) malloc((n + 1) * sizeof(int));
52+
iVec *Ci = iVec_new(m);
53+
Cp[0] = 0;
54+
55+
/* compute sparsity pattern, only storing upper triangular part */
56+
for (i = 0; i < n; i++)
57+
{
58+
/* check if Cij != 0 */
59+
for (j = i; j < n; j++)
60+
{
61+
ii = A->p[i];
62+
jj = A->p[j];
63+
64+
while (ii < A->p[i + 1] && jj < A->p[j + 1])
65+
{
66+
if (A->i[ii] == A->i[jj])
67+
{
68+
nnz += (j != i) ? 2 : 1;
69+
iVec_append(Ci, j);
70+
break;
71+
}
72+
else if (A->i[ii] < A->i[jj])
73+
{
74+
ii++;
75+
}
76+
else
77+
{
78+
jj++;
79+
}
80+
}
81+
}
82+
Cp[i + 1] = Ci->len;
83+
}
84+
85+
/* Allocate C and symmetrize it */
86+
CSR_Matrix *C = new_csr_matrix(n, n, nnz);
87+
88+
/* TODO: do we need to symmetrize here? If we are a bit careful with symmetry
89+
throughout the implementation we can skip this step. */
90+
symmetrize_csr(Cp, Ci->data, n, C);
91+
92+
/* free workspace */
93+
free(Cp);
94+
iVec_free(Ci);
95+
96+
return C;
97+
}

‎src/utils/CSR_Matrix.c‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -474,3 +474,49 @@ void csr_matvec_fill_values(const CSR_Matrix *AT, const double *z, CSR_Matrix *C
474474
// C->x[]
475475
// }
476476
// }
477+
478+
void symmetrize_csr(const int *Ap, const int *Ai, int m, CSR_Matrix *C)
479+
{
480+
int i, j, col;
481+
482+
/* Count entries per row */
483+
int *counts = (int *) calloc(m, sizeof(int));
484+
for (i = 0; i < m; i++)
485+
{
486+
for (j = Ap[i]; j < Ap[i + 1]; j++)
487+
{
488+
counts[i]++;
489+
if (Ai[j] != i) counts[Ai[j]]++;
490+
}
491+
}
492+
493+
/* Build row pointers */
494+
C->p[0] = 0;
495+
for (i = 0; i < m; i++)
496+
{
497+
C->p[i + 1] = C->p[i] + counts[i];
498+
}
499+
500+
/* Fill column indices */
501+
memset(counts, 0, m * sizeof(int));
502+
for (i = 0; i < m; i++)
503+
{
504+
for (j = Ap[i]; j < Ap[i + 1]; j++)
505+
{
506+
col = Ai[j];
507+
508+
/* Add to row i */
509+
C->i[C->p[i] + counts[i]] = col;
510+
counts[i]++;
511+
512+
/* Add symmetric entry to row col */
513+
if (col != i)
514+
{
515+
C->i[C->p[col] + counts[col]] = i;
516+
counts[col]++;
517+
}
518+
}
519+
}
520+
521+
free(counts);
522+
}

‎tests/all_tests.c‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
#include "jacobian_tests/test_quad_over_lin.h"
2020
#include "jacobian_tests/test_rel_entr.h"
2121
#include "jacobian_tests/test_sum.h"
22+
#include "utils/test_csc_matrix.h"
2223
#include "utils/test_csr_matrix.h"
2324

2425
int main(void)
@@ -75,6 +76,9 @@ int main(void)
7576
mu_run_test(test_sum_all_rows_csr, tests_run);
7677
mu_run_test(test_sum_block_of_rows_csr, tests_run);
7778
mu_run_test(test_sum_evenly_spaced_rows_csr, tests_run);
79+
mu_run_test(test_ATA_alloc_simple, tests_run);
80+
mu_run_test(test_ATA_alloc_diagonal_like, tests_run);
81+
mu_run_test(test_ATA_alloc_random, tests_run);
7882

7983
printf("\n=== All %d tests passed ===\n", tests_run);
8084

0 commit comments

Comments
 (0)