4 * Copyright: (c) Universitaet Karlsruhe
5 * Licence: This file protected by GPL - GNU GENERAL PUBLIC LICENSE.
7 * Sparse matrix storage with linked lists for rows and cols.
8 * I did not need floats, so this is all integer.
25 #include "sp_matrix.h"
29 #define MAX(a,b) ((a<b)?b:a)
31 typedef enum _iter_direction_t {
35 typedef struct _entry_t {
36 struct list_head col_chain; /**< points to next element in same column */
37 struct list_head row_chain; /**< points to next element in same row */
42 /* These specify the dimensions of the matrix.
43 * They equal the largest values ever used in matrix_set */
45 /* These are the dimensions of allocated arrays below.
46 * rowc >= maxrow and colc >= maxcol hold. */
48 /* number of entries */
50 /* arrays of list_head* as entry-points to rows and cols */
51 struct list_head **rows, **cols;
52 /* for iteration: first is to remember start-point;
53 * last was returned just before
54 * next is used in case last was removed from list */
56 struct list_head *first, *last, *next;
57 int iter_row; /* used for iteration over all elements */
60 #define _offsetof(type,member) ((char *) &(((type *) 0)->member) - (char *) 0)
61 #define _container_of(ptr,type,member) ((type *) ((char *) (ptr) - _offsetof(type, member)))
62 #define is_empty_row(row) (row>m->maxrow || list_empty(m->rows[row]))
63 #define is_empty_col(col) (col>m->maxcol || list_empty(m->cols[col]))
64 #define list_entry_by_col(h) (&list_entry(h, entry_t, col_chain)->e)
65 #define list_entry_by_row(h) (&list_entry(h, entry_t, row_chain)->e)
68 * Returns the new size for an array of size old_size,
69 * which must at least store an entry at position min.
71 static INLINE int _m_new_size(int old_size, int min) {
73 assert(min>=old_size);
78 assert(bits < sizeof(min)*8-1);
83 * Allocates space for @p count entries in the rows array and
84 * inititializes all entries from @p start to the end.
86 static INLINE void _m_alloc_row(sp_matrix_t *m, int start, int count) {
89 m->rows = realloc(m->rows, m->rowc * sizeof(*m->rows));
91 for (p=start; p<m->rowc; ++p) {
92 m->rows[p] = malloc(sizeof(*m->rows[p]));
94 INIT_LIST_HEAD(m->rows[p]);
99 * Allocates space for @p count entries in the cols array and
100 * inititializes all entries from @p start to the end.
102 static INLINE void _m_alloc_col(sp_matrix_t *m, int start, int count) {
105 m->cols = realloc(m->cols, m->colc * sizeof(*m->cols));
107 for (p=start; p<m->colc; ++p) {
108 m->cols[p] = malloc(sizeof(*m->cols[p]));
110 INIT_LIST_HEAD(m->cols[p]);
115 * Searches in row @p row for the matrix element m[row, col].
116 * @return If the element exists: Element m[row, col] and @p prev points to the list_head in the entry_t holding the element.
117 * Else: NULL and @p prev points to the list_head after which the element would be inserted.
119 static INLINE matrix_elem_t *_m_search_in_row(const sp_matrix_t *m, int row, int col, struct list_head **prev) {
120 struct list_head *start;
121 start = *prev = m->rows[row];
122 while ((*prev)->next != start && list_entry_by_row((*prev)->next)->col <= col)
123 *prev = (*prev)->next;
124 if (*prev != start) {
125 matrix_elem_t *me = list_entry_by_row(*prev);
126 if (me->row == row && me->col == col)
133 * Searches in col @p col for the matrix element m[row, col].
134 * @return If the element exists: Element m[row, col] and @p prev points to the list_head in the entry_t holding the element.
135 * Else: NULL and @p prev points to the list_head after which the element would be inserted.
137 static INLINE matrix_elem_t *_m_search_in_col(const sp_matrix_t *m, int row, int col, struct list_head **prev) {
138 struct list_head *start;
139 start = *prev = m->cols[col];
140 while ((*prev)->next != start && list_entry_by_col((*prev)->next)->row <= row)
141 *prev = (*prev)->next;
142 if (*prev != start) {
143 matrix_elem_t *me = list_entry_by_col(*prev);
144 if (me->row == row && me->col == col)
150 sp_matrix_t *new_matrix(int row_init, int col_init) {
151 sp_matrix_t *res = calloc(1, sizeof(*res));
154 _m_alloc_row(res, 0, MAX(0, row_init));
155 _m_alloc_col(res, 0, MAX(0, col_init));
159 void del_matrix(sp_matrix_t *m) {
163 for (i=0; i<m->rowc; ++i) {
164 list_for_each_entry_safe(entry_t, e, tmp, m->rows[i], row_chain)
168 for (i=0; i<m->colc; ++i)
175 void matrix_set(sp_matrix_t *m, int row, int col, int val) {
176 matrix_elem_t *me = NULL;
178 struct list_head *leftof, *above;
180 /* if necessary enlarge the matrix */
184 _m_alloc_row(m, m->rowc, _m_new_size(m->rowc, row));
189 _m_alloc_col(m, m->colc, _m_new_size(m->colc, col));
192 /* search for existing entry */
193 if (m->maxrow < m->maxcol)
194 me = _m_search_in_col(m, row, col, &above);
196 me = _m_search_in_row(m, row, col, &leftof);
198 /* if it exists, set the value and return */
203 entr = _container_of(me, entry_t, e);
204 list_del(&entr->row_chain);
205 list_del(&entr->col_chain);
212 /* if it does not exist and 0 should be set just quit */
216 /* if it does not exist and val != 0 search the other direction */
217 if (m->maxrow >= m->maxcol)
218 _m_search_in_col(m, row, col, &above);
220 _m_search_in_row(m, row, col, &leftof);
221 /* now leftof and above are the entry_t's prior the new one in each direction */
223 /* insert new entry */
224 entr = malloc(sizeof(*entr));
228 list_add(&entr->row_chain, leftof);
229 list_add(&entr->col_chain, above);
233 int matrix_get(const sp_matrix_t *m, int row, int col) {
234 struct list_head *dummy;
237 if (is_empty_row(row) || is_empty_col(col))
240 if (m->maxrow < m->maxcol)
241 me = _m_search_in_col(m, row, col, &dummy);
243 me = _m_search_in_row(m, row, col, &dummy);
246 assert(me->col == col && me->row == row);
251 int matrix_get_entries(const sp_matrix_t *m) {
255 int matrix_get_rowcount(const sp_matrix_t *m) {
259 int matrix_get_colcount(const sp_matrix_t *m) {
263 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row) {
264 if (is_empty_row(row))
267 m->first = m->rows[row];
268 m->last = m->first->next;
269 m->next = m->last->next;
270 assert (list_entry_by_row(m->last)->row == row);
271 return list_entry_by_row(m->last);
274 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col) {
275 if (is_empty_col(col))
278 m->first = m->cols[col];
279 m->last = m->first->next;
280 m->next = m->last->next;
281 assert (list_entry_by_col(m->last)->col == col);
282 return list_entry_by_col(m->last);
285 static INLINE const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow) {
286 const matrix_elem_t *res;
288 for (i=startrow; i<=m->maxrow; ++i)
289 if ((res = matrix_row_first(m, i))) {
297 const matrix_elem_t *matrix_first(sp_matrix_t *m) {
298 return matrix_first_from(m, 0);
301 const matrix_elem_t *matrix_next(sp_matrix_t *m) {
302 assert(m->last && "Start iteration with matrix_???_first, before calling me!");
305 if (m->next == m->first) {
307 return matrix_first_from(m, ++m->iter_row);
312 m->next = m->next->next;
314 return list_entry_by_col(m->last);
315 else /* right or all */
316 return list_entry_by_row(m->last);
319 static int cmp_count(const void *e1, const void *e2) {
320 return (int *)e2 - (int *)e1;
323 static INLINE void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow) {
324 const matrix_elem_t *e;
325 bitset_set(fullrow, row);
326 matrix_foreach_in_col(m, row, e)
327 if (!bitset_is_set(fullrow, e->row)) {
328 assert(0 == matrix_get(m, e->col, e->row));
329 matrix_set(m, e->col, e->row, e->val);
330 matrix_set(m, e->row, e->col, 0);
334 void matrix_optimize(sp_matrix_t *m) {
337 const matrix_elem_t *e;
340 size = MAX(m->maxcol, m->maxrow)+1;
342 /* kill all double-entries (Mij and Mji are set) */
343 matrix_foreach(m, e) {
346 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will alwazs have the same color as root");
347 t_val = matrix_get(m, e->col, e->row);
349 matrix_set(m, e->col, e->row, 0);
350 matrix_set(m, e->row, e->col, e->val + t_val);
354 c = alloca(size * sizeof(*c));
356 fullrow = bitset_alloca(size);
357 /* kill 'all' rows containing only 1 entry */
360 /* count elements in rows */
361 memset(c, 0, size * sizeof(*c));
364 for (i=0; i<size; ++i)
365 if (c[i] == 1 && !bitset_is_set(fullrow, i)) {
367 /* if the other row isn't empty move the e in there, else fill e's row */
368 if (e = matrix_row_first(m, i), e) {
370 matrix_fill_row(m, e->col, fullrow);
372 matrix_fill_row(m, e->row, fullrow);
378 memset(c, 0, size * sizeof(*c));
381 qsort(c, size, sizeof(*c), cmp_count);
384 for (i=0; i<size; ++i)
385 if (!bitset_is_set(fullrow, i))
386 matrix_fill_row(m, i, fullrow);
389 void matrix_dump(sp_matrix_t *m, FILE *out, int factor) {
391 const matrix_elem_t *e;
393 for (i=0; i <= m->maxrow; ++i) {
395 matrix_foreach_in_row(m, i, e) {
396 for (o=last_idx+1; o<e->col; ++o)
398 fprintf(out, "%d", factor*e->val);
401 for (o=last_idx+1; o<=m->maxcol; ++o)
407 void matrix_self_test(int d) {
409 const matrix_elem_t *e;
410 sp_matrix_t *m = new_matrix(10, 10);
414 matrix_set(m, i, o, i*o);
418 assert(matrix_get(m, i, o) == i*o);
421 matrix_foreach_in_row(m,1,e) {
425 assert(!matrix_next(m)); /*iter must finish */
428 matrix_foreach_in_col(m,d-1,e) {
432 assert(!matrix_next(m));
434 m = new_matrix(16,16);
435 matrix_set(m, 1,1,9);
436 matrix_set(m, 1,2,8);
437 matrix_set(m, 1,3,7);
439 matrix_set(m, 1,3,6);
440 matrix_set(m, 1,2,5);
441 matrix_set(m, 1,1,4);
443 matrix_foreach_in_row(m, 1, e) {
444 assert(e->row == 1 && e->col == i && e->val == i+3);
451 matrix_set(m, 1,1,1);
452 matrix_set(m, 2,2,2);
453 matrix_set(m, 3,3,3);
454 matrix_set(m, 3,5,4);
455 matrix_set(m, 4,4,5);
456 matrix_set(m, 5,5,6);
457 for (i=1, e = matrix_first(m); e; ++i, e=matrix_next(m))
460 matrix_set(m, 1,1,0);
461 assert(5 == matrix_get_entries(m));