2 * Sparse matrix storage with linked lists for rows and cols.
3 * I did not need floats, so this is all integer.
22 #include "sp_matrix.h"
26 #define MAX(a,b) ((a<b)?b:a)
28 typedef enum _iter_direction_t {
32 typedef struct _entry_t {
33 struct list_head col_chain; /**< points to next element in same column */
34 struct list_head row_chain; /**< points to next element in same row */
39 /* These specify the dimensions of the matrix.
40 * They equal the largest values ever used in matrix_set */
42 /* These are the dimensions of allocated arrays below.
43 * rowc >= maxrow and colc >= maxcol hold. */
45 /* number of entries */
47 /* arrays of list_head* as entry-points to rows and cols */
48 struct list_head **rows, **cols;
49 /* for iteration: first is to remember start-point;
50 * last was returned just before
51 * next is used in case last was removed from list */
53 struct list_head *first, *last, *next;
54 int iter_row; /* used for iteration over all elements */
57 #define _offsetof(type,member) ((char *) &(((type *) 0)->member) - (char *) 0)
58 #define _container_of(ptr,type,member) ((type *) ((char *) (ptr) - _offsetof(type, member)))
59 #define is_empty_row(row) (row>m->maxrow || list_empty(m->rows[row]))
60 #define is_empty_col(col) (col>m->maxcol || list_empty(m->cols[col]))
61 #define list_entry_by_col(h) (&list_entry(h, entry_t, col_chain)->e)
62 #define list_entry_by_row(h) (&list_entry(h, entry_t, row_chain)->e)
65 * Returns the new size for an array of size old_size,
66 * which must at least store an entry at position min.
68 static INLINE int _m_new_size(int old_size, int min) {
70 assert(min>=old_size);
75 assert(bits < sizeof(min)*8-1);
80 * Allocates space for @p count entries in the rows array and
81 * inititializes all entries from @p start to the end.
83 static INLINE void _m_alloc_row(sp_matrix_t *m, int start, int count) {
86 m->rows = realloc(m->rows, m->rowc * sizeof(*m->rows));
88 for (p=start; p<m->rowc; ++p) {
89 m->rows[p] = malloc(sizeof(*m->rows[p]));
91 INIT_LIST_HEAD(m->rows[p]);
96 * Allocates space for @p count entries in the cols array and
97 * inititializes all entries from @p start to the end.
99 static INLINE void _m_alloc_col(sp_matrix_t *m, int start, int count) {
102 m->cols = realloc(m->cols, m->colc * sizeof(*m->cols));
104 for (p=start; p<m->colc; ++p) {
105 m->cols[p] = malloc(sizeof(*m->cols[p]));
107 INIT_LIST_HEAD(m->cols[p]);
112 * Searches in row @p row for the matrix element m[row, col].
113 * @return If the element exists: Element m[row, col] and @p prev points to the list_head in the entry_t holding the element.
114 * Else: NULL and @p prev points to the list_head after which the element would be inserted.
116 static INLINE matrix_elem_t *_m_search_in_row(const sp_matrix_t *m, int row, int col, struct list_head **prev) {
117 struct list_head *start;
118 start = *prev = m->rows[row];
119 while ((*prev)->next != start && list_entry_by_row((*prev)->next)->col <= col)
120 *prev = (*prev)->next;
121 if (*prev != start) {
122 matrix_elem_t *me = list_entry_by_row(*prev);
123 if (me->row == row && me->col == col)
130 * Searches in col @p col for the matrix element m[row, col].
131 * @return If the element exists: Element m[row, col] and @p prev points to the list_head in the entry_t holding the element.
132 * Else: NULL and @p prev points to the list_head after which the element would be inserted.
134 static INLINE matrix_elem_t *_m_search_in_col(const sp_matrix_t *m, int row, int col, struct list_head **prev) {
135 struct list_head *start;
136 start = *prev = m->cols[col];
137 while ((*prev)->next != start && list_entry_by_col((*prev)->next)->row <= row)
138 *prev = (*prev)->next;
139 if (*prev != start) {
140 matrix_elem_t *me = list_entry_by_col(*prev);
141 if (me->row == row && me->col == col)
147 sp_matrix_t *new_matrix(int row_init, int col_init) {
148 sp_matrix_t *res = calloc(1, sizeof(*res));
151 _m_alloc_row(res, 0, MAX(0, row_init));
152 _m_alloc_col(res, 0, MAX(0, col_init));
156 void del_matrix(sp_matrix_t *m) {
160 for (i=0; i<m->rowc; ++i) {
161 list_for_each_entry_safe(entry_t, e, tmp, m->rows[i], row_chain)
165 for (i=0; i<m->colc; ++i)
172 void matrix_set(sp_matrix_t *m, int row, int col, int val) {
173 matrix_elem_t *me = NULL;
175 struct list_head *leftof, *above;
177 /* if necessary enlarge the matrix */
181 _m_alloc_row(m, m->rowc, _m_new_size(m->rowc, row));
186 _m_alloc_col(m, m->colc, _m_new_size(m->colc, col));
189 /* search for existing entry */
190 if (m->maxrow < m->maxcol)
191 me = _m_search_in_col(m, row, col, &above);
193 me = _m_search_in_row(m, row, col, &leftof);
195 /* if it exists, set the value and return */
200 entr = _container_of(me, entry_t, e);
201 list_del(&entr->row_chain);
202 list_del(&entr->col_chain);
209 /* if it does not exist and 0 should be set just quit */
213 /* if it does not exist and val != 0 search the other direction */
214 if (m->maxrow >= m->maxcol)
215 _m_search_in_col(m, row, col, &above);
217 _m_search_in_row(m, row, col, &leftof);
218 /* now leftof and above are the entry_t's prior the new one in each direction */
220 /* insert new entry */
221 entr = malloc(sizeof(*entr));
225 list_add(&entr->row_chain, leftof);
226 list_add(&entr->col_chain, above);
230 int matrix_get(const sp_matrix_t *m, int row, int col) {
231 struct list_head *dummy;
234 if (is_empty_row(row) || is_empty_col(col))
237 if (m->maxrow < m->maxcol)
238 me = _m_search_in_col(m, row, col, &dummy);
240 me = _m_search_in_row(m, row, col, &dummy);
243 assert(me->col == col && me->row == row);
248 int matrix_get_entries(const sp_matrix_t *m) {
252 int matrix_get_rowcount(const sp_matrix_t *m) {
256 int matrix_get_colcount(const sp_matrix_t *m) {
260 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row) {
261 if (is_empty_row(row))
264 m->first = m->rows[row];
265 m->last = m->first->next;
266 m->next = m->last->next;
267 assert (list_entry_by_row(m->last)->row == row);
268 return list_entry_by_row(m->last);
271 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col) {
272 if (is_empty_col(col))
275 m->first = m->cols[col];
276 m->last = m->first->next;
277 m->next = m->last->next;
278 assert (list_entry_by_col(m->last)->col == col);
279 return list_entry_by_col(m->last);
282 static INLINE const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow) {
283 const matrix_elem_t *res;
285 for (i=startrow; i<=m->maxrow; ++i)
286 if ((res = matrix_row_first(m, i))) {
294 const matrix_elem_t *matrix_first(sp_matrix_t *m) {
295 return matrix_first_from(m, 0);
298 const matrix_elem_t *matrix_next(sp_matrix_t *m) {
299 assert(m->last && "Start iteration with matrix_???_first, before calling me!");
302 if (m->next == m->first) {
304 return matrix_first_from(m, ++m->iter_row);
309 m->next = m->next->next;
311 return list_entry_by_col(m->last);
312 else /* right or all */
313 return list_entry_by_row(m->last);
316 static int cmp_count(const void *e1, const void *e2) {
317 return (int *)e2 - (int *)e1;
320 static INLINE void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow) {
321 const matrix_elem_t *e;
322 bitset_set(fullrow, row);
323 matrix_foreach_in_col(m, row, e)
324 if (!bitset_is_set(fullrow, e->row)) {
325 assert(0 == matrix_get(m, e->col, e->row));
326 matrix_set(m, e->col, e->row, e->val);
327 matrix_set(m, e->row, e->col, 0);
331 void matrix_optimize(sp_matrix_t *m) {
334 const matrix_elem_t *e;
337 size = MAX(m->maxcol, m->maxrow)+1;
339 /* kill all double-entries (Mij and Mji are set) */
340 matrix_foreach(m, e) {
343 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will alwazs have the same color as root");
344 t_val = matrix_get(m, e->col, e->row);
346 matrix_set(m, e->col, e->row, 0);
347 matrix_set(m, e->row, e->col, e->val + t_val);
351 c = alloca(size * sizeof(*c));
353 fullrow = bitset_alloca(size);
354 /* kill 'all' rows containing only 1 entry */
357 /* count elements in rows */
358 memset(c, 0, size * sizeof(*c));
361 for (i=0; i<size; ++i)
362 if (c[i] == 1 && !bitset_is_set(fullrow, i)) {
364 /* if the other row isn't empty move the e in there, else fill e's row */
365 if (e = matrix_row_first(m, i), e) {
367 matrix_fill_row(m, e->col, fullrow);
369 matrix_fill_row(m, e->row, fullrow);
375 memset(c, 0, size * sizeof(*c));
378 qsort(c, size, sizeof(*c), cmp_count);
381 for (i=0; i<size; ++i)
382 if (!bitset_is_set(fullrow, i))
383 matrix_fill_row(m, i, fullrow);
386 void matrix_dump(sp_matrix_t *m, FILE *out, int factor) {
388 const matrix_elem_t *e;
390 for (i=0; i <= m->maxrow; ++i) {
392 matrix_foreach_in_row(m, i, e) {
393 for (o=last_idx+1; o<e->col; ++o)
395 fprintf(out, "%d", factor*e->val);
398 for (o=last_idx+1; o<=m->maxcol; ++o)
404 void matrix_self_test(int d) {
406 const matrix_elem_t *e;
407 sp_matrix_t *m = new_matrix(10, 10);
411 matrix_set(m, i, o, i*o);
415 assert(matrix_get(m, i, o) == i*o);
418 matrix_foreach_in_row(m,1,e) {
422 assert(!matrix_next(m)); /*iter must finish */
425 matrix_foreach_in_col(m,d-1,e) {
429 assert(!matrix_next(m));
431 m = new_matrix(16,16);
432 matrix_set(m, 1,1,9);
433 matrix_set(m, 1,2,8);
434 matrix_set(m, 1,3,7);
436 matrix_set(m, 1,3,6);
437 matrix_set(m, 1,2,5);
438 matrix_set(m, 1,1,4);
440 matrix_foreach_in_row(m, 1, e) {
441 assert(e->row == 1 && e->col == i && e->val == i+3);
448 matrix_set(m, 1,1,1);
449 matrix_set(m, 2,2,2);
450 matrix_set(m, 3,3,3);
451 matrix_set(m, 3,5,4);
452 matrix_set(m, 4,4,5);
453 matrix_set(m, 5,5,6);
454 for (i=1, e = matrix_first(m); e; ++i, e=matrix_next(m))
457 matrix_set(m, 1,1,0);
458 assert(5 == matrix_get_entries(m));