2 * Sparse matrix storage with linked lists for rows and cols.
3 * I did not need floats, so this is all integer.
16 #include "sp_matrix.h"
20 #define MAX(a,b) ((a<b)?b:a)
22 typedef enum _iter_direction_t {
26 typedef struct _entry_t {
27 struct list_head col_chain; /**< points to next element in same column */
28 struct list_head row_chain; /**< points to next element in same row */
33 /* These specify the dimensions of the matrix.
34 * They equal the largest values ever used in matrix_set */
36 /* These are the dimensions of allocated arrays below.
37 * rowc >= maxrow and colc >= maxcol hold. */
39 /* number of entries */
41 /* arrays of list_head* as entry-points to rows and cols */
42 struct list_head **rows, **cols;
43 /* for iteration: first is to remember start-point;
44 * last was returned just before
45 * next is used in case last was removed from list */
47 struct list_head *first, *last, *next;
48 int iter_row; /* used for iteration over all elements */
51 #define _offsetof(type,member) ((char *) &(((type *) 0)->member) - (char *) 0)
52 #define _container_of(ptr,type,member) ((type *) ((char *) (ptr) - _offsetof(type, member)))
53 #define is_empty_row(row) (row>m->maxrow || list_empty(m->rows[row]))
54 #define is_empty_col(col) (col>m->maxcol || list_empty(m->cols[col]))
55 #define list_entry_by_col(h) (&list_entry(h, entry_t, col_chain)->e)
56 #define list_entry_by_row(h) (&list_entry(h, entry_t, row_chain)->e)
59 * Returns the new size for an array of size old_size,
60 * which must at least store an entry at position min.
62 static INLINE int _m_new_size(int old_size, int min) {
64 assert(min>=old_size);
69 assert(bits < sizeof(min)*8-1);
74 * Allocates space for @p count entries in the rows array and
75 * intitializes all entries from @p start to the end.
77 static INLINE void _m_alloc_row(sp_matrix_t *m, int start, int count) {
80 m->rows = realloc(m->rows, m->rowc * sizeof(*m->rows));
82 for (p=start; p<m->rowc; ++p) {
83 m->rows[p] = malloc(sizeof(*m->rows[p]));
85 INIT_LIST_HEAD(m->rows[p]);
90 * Allocates space for @p count entries in the cols array and
91 * intitializes all entries from @p start to the end.
93 static INLINE void _m_alloc_col(sp_matrix_t *m, int start, int count) {
96 m->cols = realloc(m->cols, m->colc * sizeof(*m->cols));
98 for (p=start; p<m->colc; ++p) {
99 m->cols[p] = malloc(sizeof(*m->cols[p]));
101 INIT_LIST_HEAD(m->cols[p]);
106 * Searches in row @p row for the matrix element m[row, col].
107 * @return If the element exists: Element m[row, col] and @p prev points to the list_head in the entry_t holding the element.
108 * Else: NULL and @p prev points to the list_head after which the element would be inserted.
110 static INLINE matrix_elem_t *_m_search_in_row(const sp_matrix_t *m, int row, int col, struct list_head **prev) {
111 struct list_head *start;
112 start = *prev = m->rows[row];
113 while ((*prev)->next != start && list_entry_by_row((*prev)->next)->col <= col)
114 *prev = (*prev)->next;
115 if (*prev != start) {
116 matrix_elem_t *me = list_entry_by_row(*prev);
117 if (me->row == row && me->col == col)
124 * Searches in col @p col for the matrix element m[row, col].
125 * @return If the element exists: Element m[row, col] and @p prev points to the list_head in the entry_t holding the element.
126 * Else: NULL and @p prev points to the list_head after which the element would be inserted.
128 static INLINE matrix_elem_t *_m_search_in_col(const sp_matrix_t *m, int row, int col, struct list_head **prev) {
129 struct list_head *start;
130 start = *prev = m->cols[col];
131 while ((*prev)->next != start && list_entry_by_col((*prev)->next)->row <= row)
132 *prev = (*prev)->next;
133 if (*prev != start) {
134 matrix_elem_t *me = list_entry_by_col(*prev);
135 if (me->row == row && me->col == col)
141 sp_matrix_t *new_matrix(int row_init, int col_init) {
142 sp_matrix_t *res = calloc(1, sizeof(*res));
145 _m_alloc_row(res, 0, MAX(0, row_init));
146 _m_alloc_col(res, 0, MAX(0, col_init));
150 void del_matrix(sp_matrix_t *m) {
154 for (i=0; i<m->rowc; ++i) {
155 list_for_each_entry_safe(entry_t, e, tmp, m->rows[i], row_chain)
159 for (i=0; i<m->colc; ++i)
166 void matrix_set(sp_matrix_t *m, int row, int col, int val) {
167 matrix_elem_t *me = NULL;
169 struct list_head *leftof, *above;
171 /* if necessary enlarge the matrix */
175 _m_alloc_row(m, m->rowc, _m_new_size(m->rowc, row));
180 _m_alloc_col(m, m->colc, _m_new_size(m->colc, col));
183 /* search for existing entry */
184 if (m->maxrow < m->maxcol)
185 me = _m_search_in_col(m, row, col, &above);
187 me = _m_search_in_row(m, row, col, &leftof);
189 /* if it exists, set the value and return */
194 entr = _container_of(me, entry_t, e);
195 list_del(&entr->row_chain);
196 list_del(&entr->col_chain);
203 /* if it does not exist and 0 should be set just quit */
207 /* if it does not exist and val != 0 search the other direction */
208 if (m->maxrow >= m->maxcol)
209 _m_search_in_col(m, row, col, &above);
211 _m_search_in_row(m, row, col, &leftof);
212 /* now leftof and above are the entry_t's prior the new one in each direction */
214 /* insert new entry */
215 entr = malloc(sizeof(*entr));
219 list_add(&entr->row_chain, leftof);
220 list_add(&entr->col_chain, above);
224 int matrix_get(const sp_matrix_t *m, int row, int col) {
225 struct list_head *dummy;
228 if (is_empty_row(row) || is_empty_col(col))
231 if (m->maxrow < m->maxcol)
232 me = _m_search_in_col(m, row, col, &dummy);
234 me = _m_search_in_row(m, row, col, &dummy);
237 assert(me->col == col && me->row == row);
242 int matrix_get_entries(const sp_matrix_t *m) {
246 int matrix_get_rowcount(const sp_matrix_t *m) {
250 int matrix_get_colcount(const sp_matrix_t *m) {
254 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row) {
255 if (is_empty_row(row))
258 m->first = m->rows[row];
259 m->last = m->first->next;
260 m->next = m->last->next;
261 assert (list_entry_by_row(m->last)->row == row);
262 return list_entry_by_row(m->last);
265 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col) {
266 if (is_empty_col(col))
269 m->first = m->cols[col];
270 m->last = m->first->next;
271 m->next = m->last->next;
272 assert (list_entry_by_col(m->last)->col == col);
273 return list_entry_by_col(m->last);
276 static INLINE const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow) {
277 const matrix_elem_t *res;
279 for (i=startrow; i<=m->maxrow; ++i)
280 if ((res = matrix_row_first(m, i))) {
288 const matrix_elem_t *matrix_first(sp_matrix_t *m) {
289 return matrix_first_from(m, 0);
292 const matrix_elem_t *matrix_next(sp_matrix_t *m) {
293 assert(m->last && "Start iteration with matrix_???_first, before calling me!");
296 if (m->next == m->first) {
298 return matrix_first_from(m, ++m->iter_row);
303 m->next = m->next->next;
305 return list_entry_by_col(m->last);
306 else /* right or all */
307 return list_entry_by_row(m->last);
310 static int cmp_count(const void *e1, const void *e2) {
311 return (int *)e2 - (int *)e1;
314 static INLINE void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow) {
315 const matrix_elem_t *e;
316 bitset_set(fullrow, row);
317 matrix_foreach_in_col(m, row, e)
318 if (!bitset_is_set(fullrow, e->row)) {
319 assert(0 == matrix_get(m, e->col, e->row));
320 matrix_set(m, e->col, e->row, e->val);
321 matrix_set(m, e->row, e->col, 0);
325 void matrix_optimize(sp_matrix_t *m) {
328 const matrix_elem_t *e;
331 size = MAX(m->maxcol, m->maxrow)+1;
333 /* kill all double-entries (Mij and Mji are set) */
334 matrix_foreach(m, e) {
335 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will alwazs have the same color as root");
336 int t_val = matrix_get(m, e->col, e->row);
338 matrix_set(m, e->col, e->row, 0);
339 matrix_set(m, e->row, e->col, e->val + t_val);
343 c = alloca(size * sizeof(*c));
345 fullrow = bitset_alloca(size);
346 /* kill 'all' rows containing only 1 entry */
349 /* count elements in rows */
350 memset(c, 0, size * sizeof(*c));
353 for (i=0; i<size; ++i)
354 if (c[i] == 1 && !bitset_is_set(fullrow, i)) {
356 /* if the other row isn't empty move the e in there, else fill e's row */
357 if (e = matrix_row_first(m, i), e) {
359 matrix_fill_row(m, e->col, fullrow);
361 matrix_fill_row(m, e->row, fullrow);
367 memset(c, 0, size * sizeof(*c));
370 qsort(c, size, sizeof(*c), cmp_count);
373 for (i=0; i<size; ++i)
374 if (!bitset_is_set(fullrow, i))
375 matrix_fill_row(m, i, fullrow);
378 void matrix_dump(sp_matrix_t *m, FILE *out, int factor) {
380 const matrix_elem_t *e;
382 for (i=0; i <= m->maxrow; ++i) {
384 matrix_foreach_in_row(m, i, e) {
385 for (o=last_idx+1; o<e->col; ++o)
387 fprintf(out, "%d", factor*e->val);
390 for (o=last_idx+1; o<=m->maxcol; ++o)
396 void matrix_self_test(int d) {
398 const matrix_elem_t *e;
399 sp_matrix_t *m = new_matrix(10, 10);
403 matrix_set(m, i, o, i*o);
407 assert(matrix_get(m, i, o) == i*o);
410 matrix_foreach_in_row(m,1,e) {
414 assert(!matrix_next(m)); /*iter must finish */
417 matrix_foreach_in_col(m,d-1,e) {
421 assert(!matrix_next(m));
423 m = new_matrix(16,16);
424 matrix_set(m, 1,1,9);
425 matrix_set(m, 1,2,8);
426 matrix_set(m, 1,3,7);
428 matrix_set(m, 1,3,6);
429 matrix_set(m, 1,2,5);
430 matrix_set(m, 1,1,4);
432 matrix_foreach_in_row(m, 1, e) {
433 assert(e->row == 1 && e->col == i && e->val == i+3);
440 matrix_set(m, 1,1,1);
441 matrix_set(m, 2,2,2);
442 matrix_set(m, 3,3,3);
443 matrix_set(m, 3,5,4);
444 matrix_set(m, 4,4,5);
445 matrix_set(m, 5,5,6);
446 for (i=1, e = matrix_first(m); e; ++i, e=matrix_next(m))
449 matrix_set(m, 1,1,0);
450 assert(5 == matrix_get_entries(m));